Message ID | 20210226174419.340501-1-Paul.Zimmermann@inria.fr |
---|---|
State | New |
Headers | show |
Series | Improve the accuracy of tgamma for negative inputs (BZ 26983) | expand |
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.
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 --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); }