From patchwork Fri Feb 26 17:44:19 2021 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Subject: Improve the accuracy of tgamma for negative inputs (BZ 26983) X-Patchwork-Submitter: Paul Zimmermann X-Patchwork-Id: 49822 Message-Id: <20210226174419.340501-1-Paul.Zimmermann@inria.fr> To: libc-alpha@sourceware.org Date: Fri, 26 Feb 2021 18:44:19 +0100 From: Paul Zimmermann List-Id: Libc-alpha mailing list --- sysdeps/ieee754/dbl-64/e_gamma_r.c | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) -- 2.30.0 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 #include #include +#include /* 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); }