diff mbox series

Improve the accuracy of tgamma for negative inputs (BZ 26983)

Message ID 20210226174419.340501-1-Paul.Zimmermann@inria.fr
State Superseded
Headers show
Series Improve the accuracy of tgamma for negative inputs (BZ 26983) | expand

Commit Message

Paul Zimmermann Feb. 26, 2021, 5:44 p.m. UTC
---
 sysdeps/ieee754/dbl-64/e_gamma_r.c | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

Comments

Joseph Myers Feb. 26, 2021, 6:22 p.m. UTC | #1
On Fri, 26 Feb 2021, Paul Zimmermann wrote:

> ---
>  sysdeps/ieee754/dbl-64/e_gamma_r.c | 12 ++++++++++--
>  1 file changed, 10 insertions(+), 2 deletions(-)

Please include a longer commit message explaining the change made and why 
it's an appropriate fix for the issue seen.

> @@ -188,8 +189,15 @@ __ieee754_gamma_r (double x, int *signgamp)
>  			       ? __sin (M_PI * frac)
>  			       : __cos (M_PI * (0.5 - frac)));
>  	      int exp2_adj;
> -	      double tret = M_PI / (-x * sinpix
> -				    * gamma_positive (-x, &exp2_adj));
> +              double h1, l1, h2, l2;
> +              mul_split (&h1, &l1, sinpix, gamma_positive (-x, &exp2_adj));
> +              /* sinpix*gamma_positive(.) = h1 + l1 */
> +              mul_split (&h2, &l2, h1, x);
> +              /* h1*x = h2 + l2 */
> +              /* (h1 + l1) * x = h1*x + l1*x = h2 + l2 + l1*x */
> +              l2 += l1 * x;
> +              h2 += l2;
> +	      double tret = M_PI / -h2;

Note that glibc uses tabs for the initial multiple-of-8-columns of 
indentation; this is inserting code using spaces instead.

I'd expect similar changes for ldbl-96 and ldbl-128 at least, in the 
absence of some reason it makes sense for the implementations to diverge.

The test input that justifies this change as being needed should also be 
added to auto-libm-test-in by the patch (or XFAILs for that test input 
removed if already present), with any corresponding libm-test-ulps updates 
also being included for at least one architecture on which the patch was 
tested.  Likewise the test input you found giving too-large ulps for 
binary128, unless the binary128 problem is still present after making a 
corresponding fix to the implementation for ldbl-128.
Andreas Schwab Feb. 26, 2021, 6:31 p.m. UTC | #2
On Feb 26 2021, Paul Zimmermann wrote:

> @@ -188,8 +189,15 @@ __ieee754_gamma_r (double x, int *signgamp)
>  			       ? __sin (M_PI * frac)
>  			       : __cos (M_PI * (0.5 - frac)));
>  	      int exp2_adj;
> -	      double tret = M_PI / (-x * sinpix
> -				    * gamma_positive (-x, &exp2_adj));
> +              double h1, l1, h2, l2;
> +              mul_split (&h1, &l1, sinpix, gamma_positive (-x, &exp2_adj));
> +              /* sinpix*gamma_positive(.) = h1 + l1 */
> +              mul_split (&h2, &l2, h1, x);
> +              /* h1*x = h2 + l2 */
> +              /* (h1 + l1) * x = h1*x + l1*x = h2 + l2 + l1*x */
> +              l2 += l1 * x;
> +              h2 += l2;
> +	      double tret = M_PI / -h2;

There is a mix of tab and space indentation. Please always use tab
indentation.

Andreas.
diff mbox series

Patch

diff --git a/sysdeps/ieee754/dbl-64/e_gamma_r.c b/sysdeps/ieee754/dbl-64/e_gamma_r.c
index fe69920c76..37e15e1428 100644
--- a/sysdeps/ieee754/dbl-64/e_gamma_r.c
+++ b/sysdeps/ieee754/dbl-64/e_gamma_r.c
@@ -24,6 +24,7 @@ 
 #include <math-underflow.h>
 #include <float.h>
 #include <libm-alias-finite.h>
+#include <mul_split.h>
 
 /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
    approximation to gamma function.  */
@@ -188,8 +189,15 @@  __ieee754_gamma_r (double x, int *signgamp)
 			       ? __sin (M_PI * frac)
 			       : __cos (M_PI * (0.5 - frac)));
 	      int exp2_adj;
-	      double tret = M_PI / (-x * sinpix
-				    * gamma_positive (-x, &exp2_adj));
+              double h1, l1, h2, l2;
+              mul_split (&h1, &l1, sinpix, gamma_positive (-x, &exp2_adj));
+              /* sinpix*gamma_positive(.) = h1 + l1 */
+              mul_split (&h2, &l2, h1, x);
+              /* h1*x = h2 + l2 */
+              /* (h1 + l1) * x = h1*x + l1*x = h2 + l2 + l1*x */
+              l2 += l1 * x;
+              h2 += l2;
+	      double tret = M_PI / -h2;
 	      ret = __scalbn (tret, -exp2_adj);
 	      math_check_force_underflow_nonneg (ret);
 	    }