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

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

[Paul Zimmermann]

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

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);
 	    }
-- 
2.30.0

Re: [PATCH] Improve the accuracy of tgamma for negative inputs (BZ 26983)

[Joseph Myers]

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.

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: [PATCH] Improve the accuracy of tgamma for negative inputs (BZ 26983)

[Andreas Schwab]

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.

-- 
Andreas Schwab, schwab@linux-m68k.org
GPG Key fingerprint = 7578 EB47 D4E5 4D69 2510  2552 DF73 E780 A9DA AEC1
"And now for something completely different."