diff options
Diffstat (limited to 'libc/sysdeps/ieee754/ldbl-128')
| -rw-r--r-- | libc/sysdeps/ieee754/ldbl-128/e_gammal_r.c | 145 | ||||
| -rw-r--r-- | libc/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c | 10 | ||||
| -rw-r--r-- | libc/sysdeps/ieee754/ldbl-128/gamma_productl.c | 75 |
3 files changed, 220 insertions, 10 deletions
diff --git a/libc/sysdeps/ieee754/ldbl-128/e_gammal_r.c b/libc/sysdeps/ieee754/ldbl-128/e_gammal_r.c index b6da31c13..e8d49e987 100644 --- a/libc/sysdeps/ieee754/ldbl-128/e_gammal_r.c +++ b/libc/sysdeps/ieee754/ldbl-128/e_gammal_r.c @@ -20,14 +20,108 @@ #include <math.h> #include <math_private.h> +#include <float.h> +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's + approximation to gamma function. */ + +static const long double gamma_coeff[] = + { + 0x1.5555555555555555555555555555p-4L, + -0xb.60b60b60b60b60b60b60b60b60b8p-12L, + 0x3.4034034034034034034034034034p-12L, + -0x2.7027027027027027027027027028p-12L, + 0x3.72a3c5631fe46ae1d4e700dca8f2p-12L, + -0x7.daac36664f1f207daac36664f1f4p-12L, + 0x1.a41a41a41a41a41a41a41a41a41ap-8L, + -0x7.90a1b2c3d4e5f708192a3b4c5d7p-8L, + 0x2.dfd2c703c0cfff430edfd2c703cp-4L, + -0x1.6476701181f39edbdb9ce625987dp+0L, + 0xd.672219167002d3a7a9c886459cp+0L, + -0x9.cd9292e6660d55b3f712eb9e07c8p+4L, + 0x8.911a740da740da740da740da741p+8L, + -0x8.d0cc570e255bf59ff6eec24b49p+12L, + }; + +#define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0])) + +/* Return gamma (X), for positive X less than 1775, in the form R * + 2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to + avoid overflow or underflow in intermediate calculations. */ + +static long double +gammal_positive (long double x, int *exp2_adj) +{ + int local_signgam; + if (x < 0.5L) + { + *exp2_adj = 0; + return __ieee754_expl (__ieee754_lgammal_r (x + 1, &local_signgam)) / x; + } + else if (x <= 1.5L) + { + *exp2_adj = 0; + return __ieee754_expl (__ieee754_lgammal_r (x, &local_signgam)); + } + else if (x < 12.5L) + { + /* Adjust into the range for using exp (lgamma). */ + *exp2_adj = 0; + long double n = __ceill (x - 1.5L); + long double x_adj = x - n; + long double eps; + long double prod = __gamma_productl (x_adj, 0, n, &eps); + return (__ieee754_expl (__ieee754_lgammal_r (x_adj, &local_signgam)) + * prod * (1.0L + eps)); + } + else + { + long double eps = 0; + long double x_eps = 0; + long double x_adj = x; + long double prod = 1; + if (x < 24.0L) + { + /* Adjust into the range for applying Stirling's + approximation. */ + long double n = __ceill (24.0L - x); + x_adj = x + n; + x_eps = (x - (x_adj - n)); + prod = __gamma_productl (x_adj - n, x_eps, n, &eps); + } + /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)). + 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. */ + long double exp_adj = -eps; + long double x_adj_int = __roundl (x_adj); + long double x_adj_frac = x_adj - x_adj_int; + int x_adj_log2; + long double x_adj_mant = __frexpl (x_adj, &x_adj_log2); + if (x_adj_mant < M_SQRT1_2l) + { + x_adj_log2--; + x_adj_mant *= 2.0L; + } + *exp2_adj = x_adj_log2 * (int) x_adj_int; + long double ret = (__ieee754_powl (x_adj_mant, x_adj) + * __ieee754_exp2l (x_adj_log2 * x_adj_frac) + * __ieee754_expl (-x_adj) + * __ieee754_sqrtl (2 * M_PIl / x_adj) + / prod); + exp_adj += x_eps * __ieee754_logl (x); + long double bsum = gamma_coeff[NCOEFF - 1]; + long 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 * __expm1l (exp_adj); + } +} long double __ieee754_gammal_r (long double x, int *signgamp) { - /* We don't have a real gamma implementation now. We'll use lgamma - and the exp function. But due to the required boundary - conditions we must check some values separately. */ int64_t hx; u_int64_t lx; @@ -51,8 +145,49 @@ __ieee754_gammal_r (long double x, int *signgamp) *signgamp = 0; return x - x; } + if ((hx & 0x7fff000000000000ULL) == 0x7fff000000000000ULL) + { + /* Positive infinity (return positive infinity) or NaN (return + NaN). */ + *signgamp = 0; + return x + x; + } - /* XXX FIXME. */ - return __ieee754_expl (__ieee754_lgammal_r (x, signgamp)); + if (x >= 1756.0L) + { + /* Overflow. */ + *signgamp = 0; + return LDBL_MAX * LDBL_MAX; + } + else if (x > 0.0L) + { + *signgamp = 0; + int exp2_adj; + long double ret = gammal_positive (x, &exp2_adj); + return __scalbnl (ret, exp2_adj); + } + else if (x >= -LDBL_EPSILON / 4.0L) + { + *signgamp = 0; + return 1.0f / x; + } + else + { + long double tx = __truncl (x); + *signgamp = (tx == 2.0L * __truncl (tx / 2.0L)) ? -1 : 1; + if (x <= -1775.0L) + /* Underflow. */ + return LDBL_MIN * LDBL_MIN; + long double frac = tx - x; + if (frac > 0.5L) + frac = 1.0L - frac; + long double sinpix = (frac <= 0.25L + ? __sinl (M_PIl * frac) + : __cosl (M_PIl * (0.5L - frac))); + int exp2_adj; + long double ret = M_PIl / (-x * sinpix + * gammal_positive (-x, &exp2_adj)); + return __scalbnl (ret, -exp2_adj); + } } strong_alias (__ieee754_gammal_r, __gammal_r_finite) diff --git a/libc/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c b/libc/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c index 84846fdc3..ee856acd5 100644 --- a/libc/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c +++ b/libc/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c @@ -184,13 +184,13 @@ static const int32_t two_over_pi[] = { }; static const long double c[] = { -/* 93 bits of pi/2 */ +/* 113 bits of pi/2 */ #define PI_2_1 c[0] - 1.57079632679489661923132169155131424e+00L, /* 3fff921fb54442d18469898cc5100000 */ + 0x1.921fb54442d18469898cc51701b8p+0L, /* pi/2 - PI_2_1 */ #define PI_2_1t c[1] - 8.84372056613570112025531863263659260e-29L, /* 3fa1c06e0e68948127044533e63a0106 */ + 0x3.9a252049c1114cf98e804177d4c8p-116L, }; int32_t __ieee754_rem_pio2l(long double x, long double *y) @@ -213,7 +213,7 @@ int32_t __ieee754_rem_pio2l(long double x, long double *y) { if (hx > 0) { - /* 113 + 93 bit PI is ok */ + /* 113 + 113 bit PI is ok */ z = x - PI_2_1; y[0] = z - PI_2_1t; y[1] = (z - y[0]) - PI_2_1t; @@ -221,7 +221,7 @@ int32_t __ieee754_rem_pio2l(long double x, long double *y) } else { - /* 113 + 93 bit PI is ok */ + /* 113 + 113 bit PI is ok */ z = x + PI_2_1; y[0] = z + PI_2_1t; y[1] = (z - y[0]) + PI_2_1t; diff --git a/libc/sysdeps/ieee754/ldbl-128/gamma_productl.c b/libc/sysdeps/ieee754/ldbl-128/gamma_productl.c new file mode 100644 index 000000000..157dbab9f --- /dev/null +++ b/libc/sysdeps/ieee754/ldbl-128/gamma_productl.c @@ -0,0 +1,75 @@ +/* Compute a product of X, X+1, ..., with an error estimate. + Copyright (C) 2013 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <float.h> + +/* Calculate X * Y exactly and store the result in *HI + *LO. It is + given that the values are small enough that no overflow occurs and + large enough (or zero) that no underflow occurs. */ + +static inline void +mul_split (long double *hi, long double *lo, long double x, long double y) +{ +#ifdef __FP_FAST_FMAL + /* Fast built-in fused multiply-add. */ + *hi = x * y; + *lo = __builtin_fmal (x, y, -*hi); +#elif defined FP_FAST_FMAL + /* Fast library fused multiply-add, compiler before GCC 4.6. */ + *hi = x * y; + *lo = __fmal (x, y, -*hi); +#else + /* Apply Dekker's algorithm. */ + *hi = x * y; +# define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1) + long double x1 = x * C; + long double y1 = y * C; +# undef C + x1 = (x - x1) + x1; + y1 = (y - y1) + y1; + long double x2 = x - x1; + long double y2 = y - y1; + *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2; +#endif +} + +/* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N + - 1, in the form R * (1 + *EPS) where the return value R is an + approximation to the product and *EPS is set to indicate the + approximate error in the return value. X is such that all the + values X + 1, ..., X + N - 1 are exactly representable, and X_EPS / + X is small enough that factors quadratic in it can be + neglected. */ + +long double +__gamma_productl (long double x, long double x_eps, int n, long double *eps) +{ + SET_RESTORE_ROUNDL (FE_TONEAREST); + long double ret = x; + *eps = x_eps / x; + for (int i = 1; i < n; i++) + { + *eps += x_eps / (x + i); + long double lo; + mul_split (&ret, &lo, ret, x + i); + *eps += lo / ret; + } + return ret; +} |