diff options
author | Joseph Myers <joseph@codesourcery.com> | 2015-09-10 22:27:58 +0000 |
---|---|---|
committer | Joseph Myers <joseph@codesourcery.com> | 2015-09-10 22:27:58 +0000 |
commit | 050f29c18873ec05ba04a4034bed8cb3f6ae4463 (patch) | |
tree | cc9721f5bb410543f63ad95a5913dc54eaba521d /sysdeps/ieee754/flt-32/lgamma_negf.c | |
parent | d18c36e6007b03533a38c890c68544daa78d301a (diff) | |
download | glibc-050f29c18873ec05ba04a4034bed8cb3f6ae4463.tar.gz glibc-050f29c18873ec05ba04a4034bed8cb3f6ae4463.tar.xz glibc-050f29c18873ec05ba04a4034bed8cb3f6ae4463.zip |
Fix lgamma (negative) inaccuracy (bug 2542, bug 2543, bug 2558).
The existing implementations of lgamma functions (except for the ia64 versions) use the reflection formula for negative arguments. This suffers large inaccuracy from cancellation near zeros of lgamma (near where the gamma function is +/- 1). This patch fixes this inaccuracy. For arguments above -2, there are no zeros and no large cancellation, while for sufficiently large negative arguments the zeros are so close to integers that even for integers +/- 1ulp the log(gamma(1-x)) term dominates and cancellation is not significant. Thus, it is only necessary to take special care about cancellation for arguments around a limited number of zeros. Accordingly, this patch uses precomputed tables of relevant zeros, expressed as the sum of two floating-point values. The log of the ratio of two sines can be computed accurately using log1p in cases where log would lose accuracy. The log of the ratio of two gamma(1-x) values can be computed using Stirling's approximation (the difference between two values of that approximation to lgamma being computable without computing the two values and then subtracting), with appropriate adjustments (which don't reduce accuracy too much) in cases where 1-x is too small to use Stirling's approximation directly. In the interval from -3 to -2, using the ratios of sines and of gamma(1-x) can still produce too much cancellation between those two parts of the computation (and that interval is also the worst interval for computing the ratio between gamma(1-x) values, which computation becomes more accurate, while being less critical for the final result, for larger 1-x). Because this can result in errors slightly above those accepted in glibc, this interval is instead dealt with by polynomial approximations. Separate polynomial approximations to (|gamma(x)|-1)(x-n)/(x-x0) are used for each interval of length 1/8 from -3 to -2, where n (-3 or -2) is the nearest integer to the 1/8-interval and x0 is the zero of lgamma in the relevant half-integer interval (-3 to -2.5 or -2.5 to -2). Together, the two approaches are intended to give sufficient accuracy for all negative arguments in the problem range. Outside that range, the previous implementation continues to be used. Tested for x86_64, x86, mips64 and powerpc. The mips64 and powerpc testing shows up pre-existing problems for ldbl-128 and ldbl-128ibm with large negative arguments giving spurious "invalid" exceptions (exposed by newly added tests for cases this patch doesn't affect the logic for); I'll address those problems separately. [BZ #2542] [BZ #2543] [BZ #2558] * sysdeps/ieee754/dbl-64/e_lgamma_r.c (__ieee754_lgamma_r): Call __lgamma_neg for arguments from -28.0 to -2.0. * sysdeps/ieee754/flt-32/e_lgammaf_r.c (__ieee754_lgammaf_r): Call __lgamma_negf for arguments from -15.0 to -2.0. * sysdeps/ieee754/ldbl-128/e_lgammal_r.c (__ieee754_lgammal_r): Call __lgamma_negl for arguments from -48.0 or -50.0 to -2.0. * sysdeps/ieee754/ldbl-96/e_lgammal_r.c (__ieee754_lgammal_r): Call __lgamma_negl for arguments from -33.0 to -2.0. * sysdeps/ieee754/dbl-64/lgamma_neg.c: New file. * sysdeps/ieee754/dbl-64/lgamma_product.c: Likewise. * sysdeps/ieee754/flt-32/lgamma_negf.c: Likewise. * sysdeps/ieee754/flt-32/lgamma_productf.c: Likewise. * sysdeps/ieee754/ldbl-128/lgamma_negl.c: Likewise. * sysdeps/ieee754/ldbl-128/lgamma_productl.c: Likewise. * sysdeps/ieee754/ldbl-128ibm/lgamma_negl.c: Likewise. * sysdeps/ieee754/ldbl-128ibm/lgamma_productl.c: Likewise. * sysdeps/ieee754/ldbl-96/lgamma_negl.c: Likewise. * sysdeps/ieee754/ldbl-96/lgamma_product.c: Likewise. * sysdeps/ieee754/ldbl-96/lgamma_productl.c: Likewise. * sysdeps/generic/math_private.h (__lgamma_negf): New prototype. (__lgamma_neg): Likewise. (__lgamma_negl): Likewise. (__lgamma_product): Likewise. (__lgamma_productl): Likewise. * math/Makefile (libm-calls): Add lgamma_neg and lgamma_product. * math/auto-libm-test-in: Add more tests of lgamma. * math/auto-libm-test-out: Regenerated. * sysdeps/i386/fpu/libm-test-ulps: Update. * sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
Diffstat (limited to 'sysdeps/ieee754/flt-32/lgamma_negf.c')
-rw-r--r-- | sysdeps/ieee754/flt-32/lgamma_negf.c | 288 |
1 files changed, 288 insertions, 0 deletions
diff --git a/sysdeps/ieee754/flt-32/lgamma_negf.c b/sysdeps/ieee754/flt-32/lgamma_negf.c new file mode 100644 index 0000000000..ed9659801e --- /dev/null +++ b/sysdeps/ieee754/flt-32/lgamma_negf.c @@ -0,0 +1,288 @@ +/* lgammaf expanding around zeros. + Copyright (C) 2015 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 <float.h> +#include <math.h> +#include <math_private.h> + +static const float lgamma_zeros[][2] = + { + { -0x2.74ff94p+0f, 0x1.3fe0f2p-24f }, + { -0x2.bf682p+0f, -0x1.437b2p-24f }, + { -0x3.24c1b8p+0f, 0x6.c34cap-28f }, + { -0x3.f48e2cp+0f, 0x1.707a04p-24f }, + { -0x4.0a13ap+0f, 0x1.e99aap-24f }, + { -0x4.fdd5ep+0f, 0x1.64454p-24f }, + { -0x5.021a98p+0f, 0x2.03d248p-24f }, + { -0x5.ffa4cp+0f, 0x2.9b82fcp-24f }, + { -0x6.005ac8p+0f, -0x1.625f24p-24f }, + { -0x6.fff3p+0f, 0x2.251e44p-24f }, + { -0x7.000dp+0f, 0x8.48078p-28f }, + { -0x7.fffe6p+0f, 0x1.fa98c4p-28f }, + { -0x8.0001ap+0f, -0x1.459fcap-28f }, + { -0x8.ffffdp+0f, -0x1.c425e8p-24f }, + { -0x9.00003p+0f, 0x1.c44b82p-24f }, + { -0xap+0f, 0x4.9f942p-24f }, + { -0xap+0f, -0x4.9f93b8p-24f }, + { -0xbp+0f, 0x6.b9916p-28f }, + { -0xbp+0f, -0x6.b9915p-28f }, + { -0xcp+0f, 0x8.f76c8p-32f }, + { -0xcp+0f, -0x8.f76c7p-32f }, + { -0xdp+0f, 0xb.09231p-36f }, + { -0xdp+0f, -0xb.09231p-36f }, + { -0xep+0f, 0xc.9cba5p-40f }, + { -0xep+0f, -0xc.9cba5p-40f }, + { -0xfp+0f, 0xd.73f9fp-44f }, + }; + +static const float e_hi = 0x2.b7e15p+0f, e_lo = 0x1.628aeep-24f; + +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's + approximation to lgamma function. */ + +static const float lgamma_coeff[] = + { + 0x1.555556p-4f, + -0xb.60b61p-12f, + 0x3.403404p-12f, + }; + +#define NCOEFF (sizeof (lgamma_coeff) / sizeof (lgamma_coeff[0])) + +/* Polynomial approximations to (|gamma(x)|-1)(x-n)/(x-x0), where n is + the integer end-point of the half-integer interval containing x and + x0 is the zero of lgamma in that half-integer interval. Each + polynomial is expressed in terms of x-xm, where xm is the midpoint + of the interval for which the polynomial applies. */ + +static const float poly_coeff[] = + { + /* Interval [-2.125, -2] (polynomial degree 5). */ + -0x1.0b71c6p+0f, + -0xc.73a1ep-4f, + -0x1.ec8462p-4f, + -0xe.37b93p-4f, + -0x1.02ed36p-4f, + -0xe.cbe26p-4f, + /* Interval [-2.25, -2.125] (polynomial degree 5). */ + -0xf.29309p-4f, + -0xc.a5cfep-4f, + 0x3.9c93fcp-4f, + -0x1.02a2fp+0f, + 0x9.896bep-4f, + -0x1.519704p+0f, + /* Interval [-2.375, -2.25] (polynomial degree 5). */ + -0xd.7d28dp-4f, + -0xe.6964cp-4f, + 0xb.0d4f1p-4f, + -0x1.9240aep+0f, + 0x1.dadabap+0f, + -0x3.1778c4p+0f, + /* Interval [-2.5, -2.375] (polynomial degree 6). */ + -0xb.74ea2p-4f, + -0x1.2a82cp+0f, + 0x1.880234p+0f, + -0x3.320c4p+0f, + 0x5.572a38p+0f, + -0x9.f92bap+0f, + 0x1.1c347ep+4f, + /* Interval [-2.625, -2.5] (polynomial degree 6). */ + -0x3.d10108p-4f, + 0x1.cd5584p+0f, + 0x3.819c24p+0f, + 0x6.84cbb8p+0f, + 0xb.bf269p+0f, + 0x1.57fb12p+4f, + 0x2.7b9854p+4f, + /* Interval [-2.75, -2.625] (polynomial degree 6). */ + -0x6.b5d25p-4f, + 0x1.28d604p+0f, + 0x1.db6526p+0f, + 0x2.e20b38p+0f, + 0x4.44c378p+0f, + 0x6.62a08p+0f, + 0x9.6db3ap+0f, + /* Interval [-2.875, -2.75] (polynomial degree 5). */ + -0x8.a41b2p-4f, + 0xc.da87fp-4f, + 0x1.147312p+0f, + 0x1.7617dap+0f, + 0x1.d6c13p+0f, + 0x2.57a358p+0f, + /* Interval [-3, -2.875] (polynomial degree 5). */ + -0xa.046d6p-4f, + 0x9.70b89p-4f, + 0xa.a89a6p-4f, + 0xd.2f2d8p-4f, + 0xd.e32b4p-4f, + 0xf.fb741p-4f, + }; + +static const size_t poly_deg[] = + { + 5, + 5, + 5, + 6, + 6, + 6, + 5, + 5, + }; + +static const size_t poly_end[] = + { + 5, + 11, + 17, + 24, + 31, + 38, + 44, + 50, + }; + +/* Compute sin (pi * X) for -0.25 <= X <= 0.5. */ + +static float +lg_sinpi (float x) +{ + if (x <= 0.25f) + return __sinf ((float) M_PI * x); + else + return __cosf ((float) M_PI * (0.5f - x)); +} + +/* Compute cos (pi * X) for -0.25 <= X <= 0.5. */ + +static float +lg_cospi (float x) +{ + if (x <= 0.25f) + return __cosf ((float) M_PI * x); + else + return __sinf ((float) M_PI * (0.5f - x)); +} + +/* Compute cot (pi * X) for -0.25 <= X <= 0.5. */ + +static float +lg_cotpi (float x) +{ + return lg_cospi (x) / lg_sinpi (x); +} + +/* Compute lgamma of a negative argument -15 < X < -2, setting + *SIGNGAMP accordingly. */ + +float +__lgamma_negf (float x, int *signgamp) +{ + /* Determine the half-integer region X lies in, handle exact + integers and determine the sign of the result. */ + int i = __floorf (-2 * x); + if ((i & 1) == 0 && i == -2 * x) + return 1.0f / 0.0f; + float xn = ((i & 1) == 0 ? -i / 2 : (-i - 1) / 2); + i -= 4; + *signgamp = ((i & 2) == 0 ? -1 : 1); + + SET_RESTORE_ROUNDF (FE_TONEAREST); + + /* Expand around the zero X0 = X0_HI + X0_LO. */ + float x0_hi = lgamma_zeros[i][0], x0_lo = lgamma_zeros[i][1]; + float xdiff = x - x0_hi - x0_lo; + + /* For arguments in the range -3 to -2, use polynomial + approximations to an adjusted version of the gamma function. */ + if (i < 2) + { + int j = __floorf (-8 * x) - 16; + float xm = (-33 - 2 * j) * 0.0625f; + float x_adj = x - xm; + size_t deg = poly_deg[j]; + size_t end = poly_end[j]; + float g = poly_coeff[end]; + for (size_t j = 1; j <= deg; j++) + g = g * x_adj + poly_coeff[end - j]; + return __log1pf (g * xdiff / (x - xn)); + } + + /* The result we want is log (sinpi (X0) / sinpi (X)) + + log (gamma (1 - X0) / gamma (1 - X)). */ + float x_idiff = fabsf (xn - x), x0_idiff = fabsf (xn - x0_hi - x0_lo); + float log_sinpi_ratio; + if (x0_idiff < x_idiff * 0.5f) + /* Use log not log1p to avoid inaccuracy from log1p of arguments + close to -1. */ + log_sinpi_ratio = __ieee754_logf (lg_sinpi (x0_idiff) + / lg_sinpi (x_idiff)); + else + { + /* Use log1p not log to avoid inaccuracy from log of arguments + close to 1. X0DIFF2 has positive sign if X0 is further from + XN than X is from XN, negative sign otherwise. */ + float x0diff2 = ((i & 1) == 0 ? xdiff : -xdiff) * 0.5f; + float sx0d2 = lg_sinpi (x0diff2); + float cx0d2 = lg_cospi (x0diff2); + log_sinpi_ratio = __log1pf (2 * sx0d2 + * (-sx0d2 + cx0d2 * lg_cotpi (x_idiff))); + } + + float log_gamma_ratio; +#if FLT_EVAL_METHOD != 0 + volatile +#endif + float y0_tmp = 1 - x0_hi; + float y0 = y0_tmp; + float y0_eps = -x0_hi + (1 - y0) - x0_lo; +#if FLT_EVAL_METHOD != 0 + volatile +#endif + float y_tmp = 1 - x; + float y = y_tmp; + float y_eps = -x + (1 - y); + /* We now wish to compute LOG_GAMMA_RATIO + = log (gamma (Y0 + Y0_EPS) / gamma (Y + Y_EPS)). XDIFF + accurately approximates the difference Y0 + Y0_EPS - Y - + Y_EPS. Use Stirling's approximation. */ + float log_gamma_high + = (xdiff * __log1pf ((y0 - e_hi - e_lo + y0_eps) / e_hi) + + (y - 0.5f + y_eps) * __log1pf (xdiff / y)); + /* Compute the sum of (B_2k / 2k(2k-1))(Y0^-(2k-1) - Y^-(2k-1)). */ + float y0r = 1 / y0, yr = 1 / y; + float y0r2 = y0r * y0r, yr2 = yr * yr; + float rdiff = -xdiff / (y * y0); + float bterm[NCOEFF]; + float dlast = rdiff, elast = rdiff * yr * (yr + y0r); + bterm[0] = dlast * lgamma_coeff[0]; + for (size_t j = 1; j < NCOEFF; j++) + { + float dnext = dlast * y0r2 + elast; + float enext = elast * yr2; + bterm[j] = dnext * lgamma_coeff[j]; + dlast = dnext; + elast = enext; + } + float log_gamma_low = 0; + for (size_t j = 0; j < NCOEFF; j++) + log_gamma_low += bterm[NCOEFF - 1 - j]; + log_gamma_ratio = log_gamma_high + log_gamma_low; + + return log_sinpi_ratio + log_gamma_ratio; +} |