diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/e_gamma_r.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_gamma_r.c | 37 |
1 files changed, 26 insertions, 11 deletions
diff --git a/sysdeps/ieee754/dbl-64/e_gamma_r.c b/sysdeps/ieee754/dbl-64/e_gamma_r.c index fe69920c76..eb36a9f464 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. */ @@ -88,7 +89,6 @@ gamma_positive (double x, int *exp2_adj) Compute gamma (X_ADJ + X_EPS) using Stirling's approximation, starting by computing pow (X_ADJ, X_ADJ) with a power of 2 factored out. */ - double exp_adj = -eps; double x_adj_int = round (x_adj); double x_adj_frac = x_adj - x_adj_int; int x_adj_log2; @@ -99,18 +99,22 @@ gamma_positive (double x, int *exp2_adj) x_adj_mant *= 2.0; } *exp2_adj = x_adj_log2 * (int) x_adj_int; - double ret = (__ieee754_pow (x_adj_mant, x_adj) - * __ieee754_exp2 (x_adj_log2 * x_adj_frac) - * __ieee754_exp (-x_adj) - * sqrt (2 * M_PI / x_adj) - / prod); - exp_adj += x_eps * __ieee754_log (x_adj); + double h1, l1, h2, l2; + mul_split (&h1, &l1, __ieee754_pow (x_adj_mant, x_adj), + __ieee754_exp2 (x_adj_log2 * x_adj_frac)); + mul_split (&h2, &l2, __ieee754_exp (-x_adj), sqrt (2 * M_PI / x_adj)); + mul_expansion (&h1, &l1, h1, l1, h2, l2); + /* Divide by prod * (1 + eps). */ + div_expansion (&h1, &l1, h1, l1, prod, prod * eps); + double exp_adj = x_eps * __ieee754_log (x_adj); double bsum = gamma_coeff[NCOEFF - 1]; double x_adj2 = x_adj * x_adj; for (size_t i = 1; i <= NCOEFF - 1; i++) bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i]; exp_adj += bsum / x_adj; - return ret + ret * __expm1 (exp_adj); + /* Now return (h1+l1) * exp(exp_adj), where exp_adj is small. */ + l1 += h1 * __expm1 (exp_adj); + return h1 + l1; } } @@ -188,9 +192,20 @@ __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)); - ret = __scalbn (tret, -exp2_adj); + double h1, l1, h2, l2; + h2 = gamma_positive (-x, &exp2_adj); + mul_split (&h1, &l1, sinpix, h2); + /* 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; + /* x*sinpix*gamma_positive(.) ~ h2 + l2 */ + h1 = 0x3.243f6a8885a3p+0; /* binary64 approximation of Pi */ + l1 = 0x8.d313198a2e038p-56; /* |h1+l1-Pi| < 3e-33 */ + /* Now we divide h1 + l1 by h2 + l2. */ + div_expansion (&h1, &l1, h1, l1, h2, l2); + ret = __scalbn (-h1, -exp2_adj); math_check_force_underflow_nonneg (ret); } } |