From: Adhemerval Zanella Date: Fri, 10 Oct 2025 18:15:27 +0000 (-0300) Subject: math: Use lgamma from CORE-MATH X-Git-Url: http://git.ipfire.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d67d2f468872c3fe9d3ba2b60eab0e421f906ff2;p=thirdparty%2Fglibc.git math: Use lgamma from CORE-MATH The current implementation precision shows the following accuracy, on one range ([-1,1]) with 10e9 uniform randomly generated numbers for each range (first column is the accuracy in ULP, with '0' being correctly rounded, second is the number of samples with the corresponding precision): * Range [-20, 20] * FE_TONEAREST 0: 6701254075 67.01% 1: 3230897408 32.31% 2: 63986940 0.64% 3: 3605417 0.04% 4: 233189 0.00% 5: 20973 0.00% 6: 1869 0.00% 7: 125 0.00% 8: 4 0.00% * FE_UPWARDA 0: 4207428861 42.07% 1: 5001137116 50.01% 2: 740542213 7.41% 3: 49116304 0.49% 4: 1715617 0.02% 5: 54464 0.00% 6: 4956 0.00% 7: 451 0.00% 8: 16 0.00% 9: 2 0.00% * FE_DOWNWARD 0: 4155925193 41.56% 1: 4989821364 49.90% 2: 770312796 7.70% 3: 72014726 0.72% 4: 11040522 0.11% 5: 872811 0.01% 6: 12480 0.00% 7: 106 0.00% 8: 2 0.00% * FE_TOWARDZERO 0: 4225861532 42.26% 1: 5027051105 50.27% 2: 706443411 7.06% 3: 39877908 0.40% 4: 713109 0.01% 5: 47513 0.00% 6: 4961 0.00% 7: 438 0.00% 8: 23 0.00% * Range [20, 0x5.d53649e2d4674p+1012] * FE_TONEAREST 0: 7262241995 72.62% 1: 2737758005 27.38% * FE_UPWARD 0: 4690392401 46.90% 1: 5143728216 51.44% 2: 165879383 1.66% * FE_DOWNWARD 0: 4690333331 46.90% 1: 5143794937 51.44% 2: 165871732 1.66% * FE_TOWARDZERO 0: 4690343071 46.90% 1: 5143786761 51.44% 2: 165870168 1.66% The CORE-MATH implementation is correctly rounded for any rounding mode. The code was adapted to glibc style and to use the definition of math_config.h (to handle errno, overflow, and underflow). Benchtest on x64_64 (Ryzen 9 5900X, gcc 14.2.1), aarch64 (Neoverse-N1, gcc 13.3.1), and powerpc (POWER10, gcc 13.2.1) shows: reciprocal-throughput master patched improvement x86_64 112.9740 135.8640 -20.26% x86_64v2 111.8910 131.7590 -17.76% x86_64v3 108.2800 68.0935 37.11% aarch64 61.3759 49.2403 19.77% power10 42.4483 24.1943 43.00% Latency master patched improvement x86_64 144.0090 167.9750 -16.64% x86_64v2 139.2690 167.1900 -20.05% x86_64v3 130.1320 96.9347 25.51% aarch64 66.8538 53.2747 20.31% power10 49.5076 29.6917 40.03% For x86_64/x86_64-v2, most performance hit came from the fma call through the ifunc mechanism. Checked on x86_64-linux-gnu, aarch64-linux-gnu, and powerpc64le-linux-gnu. Reviewed-by: DJ Delorie --- diff --git a/SHARED-FILES b/SHARED-FILES index 5de5b4d744..00dead58a6 100644 --- a/SHARED-FILES +++ b/SHARED-FILES @@ -243,6 +243,8 @@ core-math: sysdeps/ieee754/dbl-64/e_acosh.c # src/binary64/atanh/atanh.c, revision 4da7f241 sysdeps/ieee754/dbl-64/e_atanh.c + # src/binary64/lgamma/lgamma.c, revision 0413bb7e + sysdeps/ieee754/dbl-64/e_lgamma_r.c # src/binary64/asinh/asinh.c, revision fde815f8 sysdeps/ieee754/dbl-64/s_asinh.c # src/binary32/acos/acosf.c, revision 56dd347 diff --git a/math/Makefile b/math/Makefile index b8d5e1b514..8f4340addb 100644 --- a/math/Makefile +++ b/math/Makefile @@ -206,8 +206,6 @@ libm-calls = \ e_sqrtF \ gamma_productF \ k_tanF \ - lgamma_negF \ - lgamma_productF \ s_asinhF \ s_atanF \ s_cbrtF \ @@ -346,6 +344,8 @@ type-ldouble-routines := \ k_cosl \ k_sincosl \ k_sinl \ + lgamma_negl \ + lgamma_productl \ s_iscanonicall \ t_sincosl \ # type-ldouble-routines @@ -389,6 +389,8 @@ type-float128-routines := \ k_cosf128 \ k_sincosf128 \ k_sinf128 \ + lgamma_negf128 \ + lgamma_productf128 \ t_sincosf128 \ # type-float128-routines type-float128-yes := float128 diff --git a/sysdeps/i386/Makefile b/sysdeps/i386/Makefile index 4f96f8e191..ec1dfd98ee 100644 --- a/sysdeps/i386/Makefile +++ b/sysdeps/i386/Makefile @@ -11,6 +11,10 @@ ifeq ($(subdir),math) # the handling of MATH_SET_BOTH_ROUNDING_MODES in # sysdeps/i386/fpu/fenv_private.h. CFLAGS-e_gamma_r.c += -DMATH_SET_BOTH_ROUNDING_MODES + +# The CORE-MATH implementation assumes FLT_EVAL_METHOD == 0 to provide +# correctly rounded results. +CFLAGS-e_lgamma_r.c += -fexcess-precision=standard endif ifeq ($(subdir),gmon) diff --git a/sysdeps/ieee754/dbl-64/e_lgamma_r.c b/sysdeps/ieee754/dbl-64/e_lgamma_r.c index 0b0c4e6f4c..c630e0913f 100644 --- a/sysdeps/ieee754/dbl-64/e_lgamma_r.c +++ b/sysdeps/ieee754/dbl-64/e_lgamma_r.c @@ -1,312 +1,2059 @@ -/* @(#)er_lgamma.c 5.1 93/09/24 */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -/* __ieee754_lgamma_r(x, signgamp) - * Reentrant version of the logarithm of the Gamma function - * with user provide pointer for the sign of Gamma(x). - * - * Method: - * 1. Argument Reduction for 0 < x <= 8 - * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may - * reduce x to a number in [1.5,2.5] by - * lgamma(1+s) = log(s) + lgamma(s) - * for example, - * lgamma(7.3) = log(6.3) + lgamma(6.3) - * = log(6.3*5.3) + lgamma(5.3) - * = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3) - * 2. Polynomial approximation of lgamma around its - * minimum ymin=1.461632144968362245 to maintain monotonicity. - * On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use - * Let z = x-ymin; - * lgamma(x) = -1.214862905358496078218 + z^2*poly(z) - * where - * poly(z) is a 14 degree polynomial. - * 2. Rational approximation in the primary interval [2,3] - * We use the following approximation: - * s = x-2.0; - * lgamma(x) = 0.5*s + s*P(s)/Q(s) - * with accuracy - * |P/Q - (lgamma(x)-0.5s)| < 2**-61.71 - * Our algorithms are based on the following observation - * - * zeta(2)-1 2 zeta(3)-1 3 - * lgamma(2+s) = s*(1-Euler) + --------- * s - --------- * s + ... - * 2 3 - * - * where Euler = 0.5771... is the Euler constant, which is very - * close to 0.5. - * - * 3. For x>=8, we have - * lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+.... - * (better formula: - * lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...) - * Let z = 1/x, then we approximation - * f(z) = lgamma(x) - (x-0.5)(log(x)-1) - * by - * 3 5 11 - * w = w0 + w1*z + w2*z + w3*z + ... + w6*z - * where - * |w - f(z)| < 2**-58.74 - * - * 4. For negative x, since (G is gamma function) - * -x*G(-x)*G(x) = pi/sin(pi*x), - * we have - * G(x) = pi/(sin(pi*x)*(-x)*G(-x)) - * since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0 - * Hence, for x<0, signgam = sign(sin(pi*x)) and - * lgamma(x) = log(|Gamma(x)|) - * = log(pi/(|x*sin(pi*x)|)) - lgamma(-x); - * Note: one should avoid compute pi*(-x) directly in the - * computation of sin(pi*(-x)). - * - * 5. Special Cases - * lgamma(2+s) ~ s*(1-Euler) for tiny s - * lgamma(1)=lgamma(2)=0 - * lgamma(x) ~ -log(x) for tiny x - * lgamma(0) = lgamma(inf) = inf - * lgamma(-integer) = +-inf - * - */ +/* Correctly-rounded natural logarithm of the gamma function for binary64 +value. +Copyright (c) 2025 Alexei Sibidanov + +The original version of this file was copied from the CORE-MATH +project (file src/binary64/lgamma/lgamma.c, revision 0413bb7e). + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ + +#include +#include +#include #include -#include -#include -#include #include +#include "math_config.h" + +static inline double +fasttwosum (double x, double y, double *e) +{ + double s = x + y, z = s - x; + *e = y - z; + return s; +} + +static inline double +twosum (double x, double y, double *e) +{ + if (__glibc_likely (fabs (x) > fabs (y))) + return fasttwosum (x, y, e); + else + return fasttwosum (y, x, e); +} + +static inline double +fastsum (double xh, double xl, double yh, double yl, double *e) +{ + double sl, sh = fasttwosum (xh, yh, &sl); + *e = (xl + yl) + sl; + return sh; +} + +static inline double +sumdd (double xh, double xl, double yh, double yl, double *e) +{ + double sl, sh; + char o = fabs (xh) > fabs (yh); + if (__glibc_likely (o)) + sh = fasttwosum (xh, yh, &sl); + else + sh = fasttwosum (yh, xh, &sl); + sl += xl + yl; + *e = sl; + return sh; +} + +static inline double +muldd (double xh, double xl, double ch, double cl, double *l) +{ + double ahhh = ch * xh; + *l = (ch * xl + cl * xh) + fma (ch, xh, -ahhh); + return ahhh; +} + +static inline double +mulddd (double x, double ch, double cl, double *l) +{ + double ahhh = ch * x; + *l = cl * x + fma (ch, x, -ahhh); + return ahhh; +} + +static inline double +polydd (double xh, double xl, int n, const double c[][2], double *l) +{ + int i = n - 1; + double cl, ch = fasttwosum (c[i][0], *l, &cl); + cl += c[i][1]; + while (--i >= 0) + { + ch = muldd (xh, xl, ch, cl, &cl); + ch = fastsum (c[i][0], c[i][1], ch, cl, &cl); + } + *l = cl; + return ch; +} + +static inline double +polydddfst (double x, int n, const double c[][2], double *l) +{ + int i = n - 1; + double cl, ch = fasttwosum (c[i][0], *l, &cl); + cl += c[i][1]; + while (--i >= 0) + { + ch = mulddd (x, ch, cl, &cl); + ch = fastsum (c[i][0], c[i][1], ch, cl, &cl); + } + *l = cl; + return ch; +} + +static inline double +polyd (double x, int n, const double c[][2]) +{ + int i = n - 1; + double ch = c[i][0]; + while (--i >= 0) + ch = c[i][0] + x * ch; + return ch; +} + +static double __attribute__ ((noinline)) as_logd (double, double *); +static double __attribute__ ((noinline)) as_logd_accurate (double, double *, + double *); +static double __attribute__ ((noinline)) as_sinpipid (double, double *); +static double __attribute__ ((noinline)) as_sinpipid_accurate (double, + double *); +static double __attribute__ ((noinline)) +as_lgamma_asym_accurate (double, double *, double *); -static const double -two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ -half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ -one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ -pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ -a0 = 7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */ -a1 = 3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */ -a2 = 6.73523010531292681824e-02, /* 0x3FB13E00, 0x1A5562A7 */ -a3 = 2.05808084325167332806e-02, /* 0x3F951322, 0xAC92547B */ -a4 = 7.38555086081402883957e-03, /* 0x3F7E404F, 0xB68FEFE8 */ -a5 = 2.89051383673415629091e-03, /* 0x3F67ADD8, 0xCCB7926B */ -a6 = 1.19270763183362067845e-03, /* 0x3F538A94, 0x116F3F5D */ -a7 = 5.10069792153511336608e-04, /* 0x3F40B6C6, 0x89B99C00 */ -a8 = 2.20862790713908385557e-04, /* 0x3F2CF2EC, 0xED10E54D */ -a9 = 1.08011567247583939954e-04, /* 0x3F1C5088, 0x987DFB07 */ -a10 = 2.52144565451257326939e-05, /* 0x3EFA7074, 0x428CFA52 */ -a11 = 4.48640949618915160150e-05, /* 0x3F07858E, 0x90A45837 */ -tc = 1.46163214496836224576e+00, /* 0x3FF762D8, 0x6356BE3F */ -tf = -1.21486290535849611461e-01, /* 0xBFBF19B9, 0xBCC38A42 */ -/* tt = -(tail of tf) */ -tt = -3.63867699703950536541e-18, /* 0xBC50C7CA, 0xA48A971F */ -t0 = 4.83836122723810047042e-01, /* 0x3FDEF72B, 0xC8EE38A2 */ -t1 = -1.47587722994593911752e-01, /* 0xBFC2E427, 0x8DC6C509 */ -t2 = 6.46249402391333854778e-02, /* 0x3FB08B42, 0x94D5419B */ -t3 = -3.27885410759859649565e-02, /* 0xBFA0C9A8, 0xDF35B713 */ -t4 = 1.79706750811820387126e-02, /* 0x3F9266E7, 0x970AF9EC */ -t5 = -1.03142241298341437450e-02, /* 0xBF851F9F, 0xBA91EC6A */ -t6 = 6.10053870246291332635e-03, /* 0x3F78FCE0, 0xE370E344 */ -t7 = -3.68452016781138256760e-03, /* 0xBF6E2EFF, 0xB3E914D7 */ -t8 = 2.25964780900612472250e-03, /* 0x3F6282D3, 0x2E15C915 */ -t9 = -1.40346469989232843813e-03, /* 0xBF56FE8E, 0xBF2D1AF1 */ -t10 = 8.81081882437654011382e-04, /* 0x3F4CDF0C, 0xEF61A8E9 */ -t11 = -5.38595305356740546715e-04, /* 0xBF41A610, 0x9C73E0EC */ -t12 = 3.15632070903625950361e-04, /* 0x3F34AF6D, 0x6C0EBBF7 */ -t13 = -3.12754168375120860518e-04, /* 0xBF347F24, 0xECC38C38 */ -t14 = 3.35529192635519073543e-04, /* 0x3F35FD3E, 0xE8C2D3F4 */ -u0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */ -u1 = 6.32827064025093366517e-01, /* 0x3FE4401E, 0x8B005DFF */ -u2 = 1.45492250137234768737e+00, /* 0x3FF7475C, 0xD119BD6F */ -u3 = 9.77717527963372745603e-01, /* 0x3FEF4976, 0x44EA8450 */ -u4 = 2.28963728064692451092e-01, /* 0x3FCD4EAE, 0xF6010924 */ -u5 = 1.33810918536787660377e-02, /* 0x3F8B678B, 0xBF2BAB09 */ -v1 = 2.45597793713041134822e+00, /* 0x4003A5D7, 0xC2BD619C */ -v2 = 2.12848976379893395361e+00, /* 0x40010725, 0xA42B18F5 */ -v3 = 7.69285150456672783825e-01, /* 0x3FE89DFB, 0xE45050AF */ -v4 = 1.04222645593369134254e-01, /* 0x3FBAAE55, 0xD6537C88 */ -v5 = 3.21709242282423911810e-03, /* 0x3F6A5ABB, 0x57D0CF61 */ -s0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */ -s1 = 2.14982415960608852501e-01, /* 0x3FCB848B, 0x36E20878 */ -s2 = 3.25778796408930981787e-01, /* 0x3FD4D98F, 0x4F139F59 */ -s3 = 1.46350472652464452805e-01, /* 0x3FC2BB9C, 0xBEE5F2F7 */ -s4 = 2.66422703033638609560e-02, /* 0x3F9B481C, 0x7E939961 */ -s5 = 1.84028451407337715652e-03, /* 0x3F5E26B6, 0x7368F239 */ -s6 = 3.19475326584100867617e-05, /* 0x3F00BFEC, 0xDD17E945 */ -r1 = 1.39200533467621045958e+00, /* 0x3FF645A7, 0x62C4AB74 */ -r2 = 7.21935547567138069525e-01, /* 0x3FE71A18, 0x93D3DCDC */ -r3 = 1.71933865632803078993e-01, /* 0x3FC601ED, 0xCCFBDF27 */ -r4 = 1.86459191715652901344e-02, /* 0x3F9317EA, 0x742ED475 */ -r5 = 7.77942496381893596434e-04, /* 0x3F497DDA, 0xCA41A95B */ -r6 = 7.32668430744625636189e-06, /* 0x3EDEBAF7, 0xA5B38140 */ -w0 = 4.18938533204672725052e-01, /* 0x3FDACFE3, 0x90C97D69 */ -w1 = 8.33333333333329678849e-02, /* 0x3FB55555, 0x5555553B */ -w2 = -2.77777777728775536470e-03, /* 0xBF66C16C, 0x16B02E5C */ -w3 = 7.93650558643019558500e-04, /* 0x3F4A019F, 0x98CF38B6 */ -w4 = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */ -w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */ -w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */ - -static const double zero= 0.00000000000000000000e+00; - -static double -sin_pi(double x) +static __attribute__ ((noinline)) double +as_lgamma_database (double x, double f) { - double y,z; - int n,ix; - - GET_HIGH_WORD(ix,x); - ix &= 0x7fffffff; - - if(ix<0x3fd00000) return __sin(pi*x); - y = -x; /* x is assume negative */ - - /* - * argument reduction, make sure inexact flag not raised if input - * is an integer - */ - z = floor(y); - if(z!=y) { /* inexact anyway */ - y *= 0.5; - y = 2.0*(y - floor(y)); /* y = |x| mod 2.0 */ - n = (int) (y*4.0); - } else { - if(ix>=0x43400000) { - y = zero; n = 0; /* y must be even */ - } else { - if(ix<0x43300000) z = y+two52; /* exact */ - GET_LOW_WORD(n,z); - n &= 1; - y = n; - n<<= 2; + static const double db[][3] = { + { -0x1.1b649eb4316fbp+3, -0x1.50332a035af1fp+3, -0x1p-51 }, + { -0x1.808d3e2f56b4fp+2, -0x1.d779a9ab6cbffp+0, 0x1p-54 }, + { -0x1.d02b2008d5bf5p+1, -0x1.67cf93db4863ep+0, -0x1p-54 }, + { -0x1.a123d403647d7p+1, -0x1.530824ff0740ap-1, 0x1p-55 }, + { -0x1.90404f46978b6p+1, 0x1.189e211cebce7p-3, -0x1p-57 }, + { -0x1.75692939b3defp+1, 0x1.a11a8b4dd58bbp-1, 0x1p-55 }, + { -0x1.64cc652d46a2dp+1, 0x1.7c2fd07e26d78p-4, -0x1p-58 }, + { -0x1.593ae533139c1p+1, -0x1.316f7a9bb3e4ep-4, 0x1p-58 }, + { -0x1.549cde4f1fd16p+1, -0x1.ab1280c638adap-4, -0x1p-58 }, + { -0x1.53712a3a51156p+1, -0x1.bedf4564c976dp-4, 0x1p-58 }, + { -0x1.51333e5a494f7p+1, -0x1.d93a67f2ad756p-4, 0x1p-58 }, + { -0x1.504f0b5b9d1f9p+1, -0x1.dfa224bc3435ep-4, -0x1p-58 }, + { -0x1.4c213b7fa02dcp+1, -0x1.e051ae181baabp-4, 0x1p-58 }, + { -0x1.75824f4f0c7c8p+0, 0x1.cb23f1fc0296ep-1, 0x1p-55 }, + { -0x1.667cb87bcfa0fp+0, 0x1.f478a478204c5p-1, 0x1p-55 }, + { -0x1.51668b6af122dp+0, 0x1.2728422577ep+0, 0x1p-54 }, + { -0x1.50b41410187c1p+0, 0x1.290131f9d5feap+0, -0x1p-54 }, + { -0x1.307eb80d114afp-2, 0x1.7870d113b3febp+0, -0x1p-54 }, + { -0x1.26923ac1f7c1p-2, 0x1.7df2c32cf08a3p+0, 0x1p-54 }, + }; + int a = 0, b = array_length (db) - 1, m = (a + b) / 2; + while (a <= b) + { // binary search + if (db[m][0] < x) + a = m + 1; + else if (db[m][0] == x) + { + f = db[m][1] + db[m][2]; + break; + } + else + b = m - 1; + m = (a + b) / 2; + } + return f; +} + +static __attribute__ ((noinline)) double +as_lgamma_accurate (double x) +{ + // polynomial coefficients for (lgamma(x)+ln(|x|))/x where x in [-0.25,0.25] + static const double c0[][2] + = { { -0x1.2788cfc6fb619p-1, 0x1.6cb90701fbf92p-58 }, + { 0x1.a51a6625307d3p-1, 0x1.1873d891220dfp-56 }, + { -0x1.9a4d55beab2d7p-2, 0x1.4c26d1b4d5994p-59 }, + { 0x1.151322ac7d848p-2, 0x1.b5f9120fdaca6p-57 }, + { -0x1.a8b9c17aa6149p-3, -0x1.2e826b9f2e3f2p-58 }, + { 0x1.5b40cb100c306p-3, 0x1.4a79a5f96e16fp-59 }, + { -0x1.2703a1dcea3aep-3, -0x1.63066fbe425d6p-57 }, + { 0x1.010b36af86397p-3, -0x1.7438c20a423a8p-59 }, + { -0x1.c806706d57db4p-4, -0x1.5a891c586f017p-58 }, + { 0x1.9a01e385d5f8fp-4, 0x1.9eb2a3ca8e809p-59 }, + { -0x1.748c33114c6d5p-4, -0x1.505ae0736c854p-58 }, + { 0x1.556ad63243bc2p-4, -0x1.0cfd5c4f86f62p-58 }, + { -0x1.3b1d971fc59e2p-4, -0x1.9a19a5ba46076p-61 }, + { 0x1.2496df8320d56p-4, 0x1.51588f44e1564p-58 }, + { -0x1.11133476e5f97p-4, -0x1.fb27c7d952db7p-59 }, + { 0x1.00010064c948ap-4, 0x1.c2a9309d327a6p-60 }, + { -0x1.e1e2d312e9abp-5, 0x1.beee6edc135cfp-59 }, + { 0x1.c71ce3a414da2p-5, 0x1.7d6e71837bc1dp-62 }, + { -0x1.af28a18673c39p-5, 0x1.df2c138661f9bp-60 }, + { 0x1.9999b2dfca9bbp-5, -0x1.afb6a810fa7c2p-59 }, + { -0x1.861874180c0f2p-5, -0x1.c7900b8893e2bp-59 }, + { 0x1.745d2798dac6ap-5, 0x1.e962ae5d1c788p-61 }, + { -0x1.642be2fd1947dp-5, 0x1.94490a4e3d6b4p-61 }, + { 0x1.55545da95a50bp-5, 0x1.b56c4cc4dd96bp-63 }, + { -0x1.47ba80f965799p-5, 0x1.3bcac9fdf0d1ap-60 }, + { 0x1.3b24eb343ccefp-5, -0x1.83fb9b692b013p-61 }, + { -0x1.2eba18fc9675p-5, -0x1.3069a0cc7e12dp-60 }, + { 0x1.23b27295dc8b4p-5, 0x1.a6d42b2188101p-60 }, + { -0x1.2136ba2a595c4p-5, -0x1.c2acae5604f9dp-59 }, + { 0x1.191f0bce33ba9p-5, -0x1.8470cbefa841ap-61 }, + { -0x1.b86662fc0fb81p-6, -0x1.236cd67b781c2p-60 }, + { 0x1.9d6d718c336abp-6, 0x1.c632d9ae925d4p-60 }, + { -0x1.9ea8d7c263a8p-5, -0x1.bdcf945212655p-59 }, + { 0x1.9f461fd74cc3p-5, 0x1.5e5f99eb89ba4p-60 } }; + + // polynomial coefficients for lgamma(1.5 + x) where x in [-0.25,0.25] + static const double b[][2] + = { { -0x1.eeb95b094c191p-4, -0x1.346863f58b074p-58 }, + { 0x1.2aed059bd608ap-5, 0x1.cd3d2ca77b58cp-63 }, + { 0x1.de9e64df22ef3p-2, -0x1.6d48ec99346d4p-57 }, + { -0x1.1ae55b180726cp-3, -0x1.959aeebbdcd5bp-59 }, + { 0x1.e0f840dad61dap-5, -0x1.599fc3f5e3ccp-59 }, + { -0x1.da59d5374a543p-6, -0x1.0628c37280b13p-63 }, + { 0x1.f9ca39daa929cp-7, -0x1.69f468b73daddp-67 }, + { -0x1.1a8ba4f0ea597p-7, -0x1.7d1edde63f06cp-61 }, + { 0x1.456f1ad666a3bp-8, -0x1.31d73594fa521p-62 }, + { -0x1.7edb812f6426ep-9, -0x1.5f45d321e0ae2p-64 }, + { 0x1.c9735ae9db2c1p-10, -0x1.827b962d41f07p-64 }, + { -0x1.148a319eec639p-10, 0x1.a64298294f3a5p-64 }, + { 0x1.517c5a1579f4p-11, -0x1.055841dae38a6p-67 }, + { -0x1.9eff1d1c8be2dp-12, -0x1.fba56c797e258p-67 }, + { 0x1.00c41c13e3352p-12, -0x1.9c7d468e4cbd4p-73 }, + { -0x1.3f6dff22a8ffep-13, 0x1.79e9e76cea415p-68 }, + { 0x1.8f36195536e0cp-14, -0x1.9de6ad54f6997p-70 }, + { -0x1.f4ea079eaf734p-15, -0x1.fed2cdd5ce819p-69 }, + { 0x1.3b5e73a6bf0cep-15, -0x1.b8bf8cedf2f5dp-70 }, + { -0x1.8e583400720b6p-16, -0x1.58007f65918dcp-70 }, + { 0x1.f88ed0371a528p-17, 0x1.ca095e5525ba6p-75 }, + { -0x1.40597e1beb94dp-17, 0x1.99bf6a1da80a1p-71 }, + { 0x1.97b2f3649dc6ap-18, 0x1.39d919c548f05p-72 }, + { -0x1.03fa905165ffep-18, 0x1.0e61588355356p-72 }, + { 0x1.4c8df9a3d368p-19, -0x1.7a6689d0b7dd2p-73 }, + { -0x1.a9bafae87dbd3p-20, 0x1.5777bb7ae251ap-74 }, + { 0x1.0b2f2a9e11fe8p-20, -0x1.cf209a5b02c67p-75 }, + { -0x1.567358802c384p-21, 0x1.c44014666f1c9p-75 }, + { 0x1.1219755496b17p-21, -0x1.b113ac02e4428p-75 }, + { -0x1.63804ffdcd648p-22, 0x1.ba1993b4ec9b7p-76 } }; + + double sx = x, fh, fl = 0, fll = 0; + x = fabs (x); + if (x < 0x1p-100) + { + double lll, ll, lh = as_logd_accurate (x, &ll, &lll); + fh = -lh; + fl = -ll; + fll = -lll; + fh = fasttwosum (fh, fl, &fl); + fl = fasttwosum (fl, fll, &fll); + double e; + fasttwosum (fh, 2 * fl, &e); + if (e == 0 && 1 + 0x1p-54 == 1 - 0x1p-54) + fl *= 1 + copysign (0x1p-52, fl) * copysign (1, fll); + } + else if (x < 0x1p-2) + { + fh = polydddfst (sx, array_length (c0), c0, &fl); + fh = mulddd (sx, fh, fl, &fl); + double lll, ll, lh = as_logd_accurate (x, &ll, &lll); + fh = sumdd (fh, fl, -ll, -lll, &fl); + fh = twosum (-lh, fh, &fll); + fl = twosum (fll, fl, &fll); + double e; + fasttwosum (fh, 2 * fl, &e); + if (e == 0 && 1 + 0x1p-54 == 1 - 0x1p-54) + fl *= 1 + copysign (0x1p-52, fl) * copysign (1, fll); + } + else + { + if (fabs (x - 0.5) < 0x1p-2) + { + fh = polydddfst (x - 0.5, array_length(b), b, &fl); + if (sx > 0) + { + double lll, ll, lh = as_logd_accurate (x, &ll, &lll); + fl = sumdd (fl, 0, -ll, -lll, &fll); + fh = twosum (fh, -lh, &lh); + fl = sumdd (fl, fll, lh, 0, &fll); } } - switch (n) { - case 0: y = __sin(pi*y); break; - case 1: - case 2: y = __cos(pi*(0.5-y)); break; - case 3: - case 4: y = __sin(pi*(one-y)); break; - case 5: - case 6: y = -__cos(pi*(y-1.5)); break; - default: y = __sin(pi*(y-2.0)); break; + else if (fabs (x - 2.5) < 0x1p-2) + { + fh = polydddfst (x - 2.5, array_length (b), b, &fl); + double lll, ll, lh = as_logd_accurate (x - 1, &ll, &lll); + fl = sumdd (fl, 0, ll, lll, &fll); + fh = twosum (fh, lh, &lh); + fl = sumdd (fl, fll, lh, 0, &fll); + if (sx < 0) + { + lh = as_logd_accurate (x, &ll, &lll); + fl = sumdd (fl, fll, ll, lll, &fll); + fh = twosum (fh, lh, &lh); + fl = sumdd (fl, fll, lh, 0, &fll); } - return -y; -} + } + else if (fabs (x - 3.5) < 0x1p-2) + { + double l2ll, l2l, l2h = as_logd_accurate (x - 2, &l2l, &l2ll); + double l1ll, l1l, l1h = as_logd_accurate (x - 1, &l1l, &l1ll); + l1l = sumdd (l1l, l1ll, l2l, l2ll, &l1ll); + l1h = fasttwosum (l1h, l2h, &l2h); + l1l = sumdd (l1l, l1ll, l2h, 0, &l1ll); + fh = polydddfst (x - 3.5, array_length (b), b, &fl); + fl = sumdd (fl, 0, l1l, l1ll, &fll); + fh = twosum (fh, l1h, &l1h); + fl = sumdd (fl, fll, l1h, 0, &fll); + if (sx < 0) + { + double lll, ll, lh = as_logd_accurate (x, &ll, &lll); + fl = sumdd (fl, fll, ll, lll, &fll); + fh = twosum (fh, lh, &lh); + fl = sumdd (fl, fll, lh, 0, &fll); + } + } + else if (fabs (x - 1) < 0x1p-2) + { + fh = polydddfst (x - 1, array_length (c0), c0, &fl); + fh = mulddd (x - 1, fh, fl, &fl); + if (sx < 0) + { + double lll, ll, lh = as_logd_accurate (x, &ll, &lll); + fl = sumdd (fl, 0, ll, lll, &fll); + fh = twosum (fh, lh, &lh); + fl = sumdd (fl, fll, lh, 0, &fll); + } + } + else if (fabs (x - 1.5) < 0x1p-2) + { + fh = polydddfst (x - 1.5, array_length (b), b, &fl); + if (sx < 0) + { + double lll, ll, lh = as_logd_accurate (x, &ll, &lll); + fl = sumdd (fl, 0, ll, lll, &fll); + fh = twosum (fh, lh, &lh); + fl = sumdd (fl, fll, lh, 0, &fll); + } + } + else if (fabs (x - 2) < 0x1p-2) + { + double lll, ll, lh = as_logd_accurate (x - 1, &ll, &lll); + fh = polydddfst (x - 2, array_length (c0), c0, &fl); + fh = mulddd (x - 2, fh, fl, &fl); + fl = sumdd (fl, 0, ll, lll, &fll); + fh = twosum (fh, lh, &lh); + fl = sumdd (fl, fll, lh, 0, &fll); + if (sx < 0) + { + lh = as_logd_accurate (x, &ll, &lll); + fl = sumdd (fl, fll, ll, lll, &fll); + fh = twosum (fh, lh, &lh); + fl = sumdd (fl, fll, lh, 0, &fll); + } + } + else if (fabs (x - 3) < 0x1p-2) + { + double l2ll, l2l, l2h = as_logd_accurate (x - 2, &l2l, &l2ll); + double l1ll, l1l, l1h = as_logd_accurate (x - 1, &l1l, &l1ll); + l1l = sumdd (l1l, l1ll, l2l, l2ll, &l1ll); + l1h = fasttwosum (l1h, l2h, &l2h); + l1l = sumdd (l1l, l1ll, l2h, 0, &l1ll); + fh = polydddfst (x - 3, array_length (c0), c0, &fl); + fh = mulddd (x - 3, fh, fl, &fl); + fl = sumdd (fl, 0, l1l, l1ll, &fll); + fh = twosum (fh, l1h, &l1h); + fl = sumdd (fl, fll, l1h, 0, &fll); + if (sx < 0) + { + double lll, ll, lh = as_logd_accurate (x, &ll, &lll); + fl = sumdd (fl, fll, ll, lll, &fll); + fh = twosum (fh, lh, &lh); + fl = sumdd (fl, fll, lh, 0, &fll); + } + } + else + { // x>3.75 + if (sx < 0) + x = fasttwosum (x, 1, &fl); + fh = as_lgamma_asym_accurate (x, &fl, &fll); + } + if (sx < 0) + { + double phi = (sx < -0.5) ? sx - floor (sx) : -sx; + double sl, sh = as_sinpipid_accurate (phi, &sl); + double lll, ll, lh = as_logd_accurate (sh, &ll, &lll); + ll += sl / sh + lll; + fl = sumdd (fl, fll, ll, 0, &fll); + fh = twosum (fh, lh, &lh); + fl = sumdd (fl, fll, lh, 0, &fll); + fh = -fh; + fl = -fl; + fll = -fll; + } + fh = fasttwosum (fh, fl, &fl); + fl = fasttwosum (fl, fll, &fll); + fh = fasttwosum (fh, fl, &fl); + fl = fasttwosum (fl, fll, &fll); + double e; + fasttwosum (fh, 2 * fl, &e); + if (e == 0) + { + double dfl = 1 + copysign (0x1p-26, fl) * copysign (0x1p-26, fll); + fl *= dfl; + } + } +#define SX_BND -0x1.4e147ae147ae1p+1 // approximates -2.61 to nearest + + if (fabs (fh) < 0x1.8p-2) + { + if (fabs (fh) < 0x1.74p-4 && sx > SX_BND && sx < -2) + { + static const double x0[] + = { 0x1.3a7fc9600f86cp+1, 0x1.55f64f98af8dp-55, + 0x1.c4b0cd201366ap-110 }; + static const double c[][2] + = { { 0x1.83fe966af535fp+0, -0x1.775909a36a68cp-55 }, + { 0x1.36eebb002f55dp-1, -0x1.8d4b2124a39f8p-55 }, + { 0x1.694a6058a7858p-6, -0x1.1d8c8b9b4d80ep-61 }, + { 0x1.1718d7ca09e5bp-6, 0x1.83195b07fd25p-60 }, + { 0x1.7339fe04b2764p-10, -0x1.48648f4a4bf9ep-64 }, + { 0x1.8d32f682aa0bdp-11, -0x1.90953dadfba01p-66 }, + { 0x1.809f04ee6e0fap-14, -0x1.6100a3d177e7cp-68 }, + { 0x1.48eaa81657361p-15, 0x1.42b86f83f623dp-70 }, + { 0x1.9297adb2def5ap-18, -0x1.0c295492288fdp-73 }, + { 0x1.286fb8cbaebb5p-19, -0x1.ebf1f6a584b62p-77 }, + { 0x1.a92e0a5de4bc9p-22, -0x1.688db240a9a82p-76 }, + { 0x1.1a9d4d8c6284ep-23, 0x1.08624c7186639p-78 }, + { 0x1.c4cd2594e91c9p-26, 0x1.a225e59cfef5dp-81 }, + { 0x1.18737ec8e68aap-27, 0x1.6421ab9e4c553p-82 }, + { 0x1.e6028795be5c1p-30, -0x1.b4bd2f36bb0ddp-84 }, + { 0x1.1eacaeb800afdp-31, -0x1.cedc75cbd8cc4p-87 }, + { 0x1.06bce9e1f6b9p-33, -0x1.0fb524fda64ebp-87 }, + { 0x1.2bb110a516b79p-35, 0x1.2d7a224941f92p-90 }, + { 0x1.1e0024589b848p-37, 0x1.1a84fa416703dp-92 }, + { 0x1.3ec6ebb4ba73p-39, -0x1.eb9cd73c631d9p-94 }, + { 0x1.3936af34a015p-41, -0x1.83b442df73674p-95 }, + { 0x1.57d4ece262198p-43, 0x1.20b160178fca3p-97 }, + { 0x1.5b530d802ffffp-45, 0x1.8bb2ee582dbd1p-104 }, + { 0x1.7c323d20053dp-47, -0x1.32d5dd28cac13p-103 }, + { 0x1.6131339e2b76ap-49, 0x1.4422b26549563p-104 }, + { 0x1.b4f6aaa0d8886p-52, 0x1.796e66ee63a34p-106 } }; + const double sc = 0x1p+3; + double zl, zh = fasttwosum (x0[0] + sx, x0[1], &zl); + zl += x0[2]; + double sh = zh * sc, sl = zl * sc; + int n = array_length (c), k = 6; + fl = sh * polyd (sh, k, c + n - k); + fh = polydd (sh, sl, n - k, c, &fl); + fh = muldd (zh, zl, fh, fl, &fl); + } + else if (fabs (fh) < 0x1.168p-4 && sx > -3 && sx < SX_BND) + { + static const double x0[] + = { 0x1.5fb410a1bd901p+1, -0x1.a19a96d2e6f85p-54, + -0x1.140b4ff4b7d6p-108 }; + static const double c[][2] + = { { -0x1.ea12da904b18cp+0, -0x1.220130f99b2cbp-54 }, + { 0x1.3267f3c265a52p-1, -0x1.1c630ff19db83p-55 }, + { -0x1.4185ac30c8bf2p-4, 0x1.f161263693e12p-59 }, + { 0x1.f504accc9f19bp-7, -0x1.eacc0226c7208p-62 }, + { -0x1.8588458207eacp-9, 0x1.4b51668f3dff4p-63 }, + { 0x1.4373f7cc709b3p-11, -0x1.2474e20e777aep-66 }, + { -0x1.12239bdd6c013p-13, 0x1.46df69aa3032cp-69 }, + { 0x1.dba65e27421c4p-16, 0x1.36bfe12004625p-71 }, + { -0x1.a2d2504d7e987p-18, 0x1.a6b6dfe9fa6f4p-75 }, + { 0x1.7581739ee6087p-20, -0x1.23100aff1ab78p-77 }, + { -0x1.506c65fad6188p-22, -0x1.17d1749eb738fp-77 }, + { 0x1.318ef724f7814p-24, -0x1.6517edf27abc2p-82 }, + { -0x1.17767260d9825p-26, -0x1.03e683dfe1d87p-82 }, + { 0x1.011e34454ade3p-28, -0x1.2a844fd6e7622p-82 }, + { -0x1.db8b9e6b1e0e1p-31, -0x1.f73094d7d4e4p-85 }, + { 0x1.b9bab1fca321p-33, -0x1.89b7151db3448p-88 }, + { -0x1.9bed46f4d0aa2p-35, -0x1.1af984f4b8e4p-90 }, + { 0x1.81780d4242119p-37, -0x1.7576147eed5dcp-91 }, + { -0x1.69d3ce15a3f87p-39, 0x1.9823fc1b6be9ep-93 }, + { 0x1.5494a34ce847p-41, -0x1.0ec96907a37b4p-95 }, + { -0x1.4160eccd5982ep-43, 0x1.2b4a92d7b0c76p-98 }, + { 0x1.2fe0a640fb1ap-45, -0x1.b1f41cf9689edp-106 }, + { -0x1.1ff8659402df5p-47, 0x1.511b4097b0fccp-104 }, + { 0x1.1357a3e9e1cf1p-49, -0x1.bda0e37acdd97p-106 }, + { -0x1.0b145c5c5b9abp-51, -0x1.b52b4002dbc0cp-105 }, + { 0x1.dc1888c7036cap-54, 0x1.2a69b093c9a8bp-109 }, + { -0x1.060b52b5d6f68p-56, 0x1.d17de397a2b1p-110 } }; + const double sc = 0x1p+4; + double zl, zh = fasttwosum (x0[0] + sx, x0[1], &zl); + zl += x0[2]; + double sh = zh * sc, sl = zl * sc; + int n = array_length (c), k = 6; + fl = sh * polyd (sh, k, c + n - k); + fh = polydd (sh, sl, n - k, c, &fl); + fh = muldd (zh, zl, fh, fl, &fl); + } + else if (fabs (fh) < 0x1.2d4p-4 && sx > -3.5 && sx < -3) + { + static const double x0[] + = { 0x1.9260dbc9e59afp+1, 0x1.f717cd335a7b3p-53, + 0x1.d32a2a65bfd63p-107 }; + static const double c[][2] + = { { 0x1.f20a65f2fac55p+2, -0x1.1d258e4b0beb5p-53 }, + { 0x1.9d4d2977150efp-2, 0x1.a04089578ae88p-56 }, + { 0x1.c1137124d5c5bp-6, 0x1.d6c922d2512d6p-61 }, + { 0x1.267203d776b0ep-9, -0x1.aa6081d790e05p-63 }, + { 0x1.99a6337da39ddp-13, 0x1.49aeda9400147p-68 }, + { 0x1.293c3f78d3bdbp-16, 0x1.ee59e48e8b181p-73 }, + { 0x1.bb97aa0b71e45p-20, -0x1.592602ea73795p-78 }, + { 0x1.51ea3345f5349p-23, 0x1.ef9552f50f84ep-79 }, + { 0x1.057f65c64b21p-26, -0x1.36ca6c4396475p-80 }, + { 0x1.99c8650e3a4f9p-30, 0x1.7d3b9d2266a1p-84 }, + { 0x1.44520c3a4f841p-33, 0x1.8c284070d32dfp-89 }, + { 0x1.02d221961d8b5p-36, -0x1.0d56a1938ace2p-90 }, + { 0x1.9ffcd97899d22p-40, -0x1.d6040c8c03da8p-95 }, + { 0x1.50494b6c20faap-43, -0x1.644c84757ab17p-97 }, + { 0x1.113fe94711cb9p-46, -0x1.f7b0588af7045p-100 }, + { 0x1.be09d37bab523p-50, 0x1.0f8dbcb173a93p-106 }, + { 0x1.6d6b80656d502p-53, -0x1.9e5c21b32ebp-110 }, + { 0x1.2ca70c7ef211ep-56, 0x1.6511b75c60833p-111 }, + { 0x1.f9264b4df26e4p-60, -0x1.aefcabcc01fd2p-121 }, + { 0x1.9262efb57925p-63, -0x1.1196c645b5736p-117 } }; + const double sc = 0x1p+6; + double zl, zh = fasttwosum (x0[0] + sx, x0[1], &zl); + zl += x0[2]; + double sh = zh * sc, sl = zl * sc; + int n = array_length (c), k = 7; + fl = sh * polyd (sh, k, c + n - k); + fh = polydd (sh, sl, n - k, c, &fl); + fh = muldd (zh, zl, fh, fl, &fl); + } + else if (fabs (fh) < 0x1.efp-5 && sx > -4 && sx < -3.5) + { + static const double x0[] + = { 0x1.fa471547c2fe5p+1, 0x1.70d4561291237p-56, + -0x1.9e6fadbbc171ap-111 }; + static const double c[][2] + = { { -0x1.4b99d966c5647p+4, 0x1.9cba2450b0003p-50 }, + { 0x1.f76deae0436bep-1, -0x1.5af99a1af2d4bp-55 }, + { -0x1.d25359d4b2f38p-5, 0x1.10c02bb27b3e8p-60 }, + { 0x1.e8f829f141aa5p-9, 0x1.4b3ff24054d7dp-65 }, + { -0x1.116f7806d26d3p-12, -0x1.a2efbebc9bc45p-68 }, + { 0x1.3e8f3ab9fc1f4p-16, 0x1.e3a4940cdb9f6p-70 }, + { -0x1.7dbbe062ffd9ep-20, -0x1.2986222077472p-74 }, + { 0x1.d2f76de7bd027p-24, -0x1.6c0ccf4d8669cp-78 }, + { -0x1.2225fe4f84932p-27, 0x1.a38a5afc8eebap-82 }, + { 0x1.6d12ae1936ba4p-31, -0x1.de21fabe6423p-86 }, + { -0x1.cffc2a8f6492dp-35, 0x1.64881462031d6p-91 }, + { 0x1.294e1bdd838bcp-38, 0x1.dcbb67f373f7dp-92 }, + { -0x1.7fab625b2f2ffp-42, -0x1.e8527e1dc6098p-100 }, + { 0x1.f211abe57d76p-46, 0x1.d65f8610deeeep-102 }, + { -0x1.44f2a05f33c9dp-49, 0x1.3fb27bfc54f94p-106 }, + { 0x1.a9e4622460be6p-53, 0x1.2f5912b9785cp-107 }, + { -0x1.182776b50f5cfp-56, -0x1.d30d2e8d300a5p-110 }, + { 0x1.722b65d64ddcp-60, -0x1.0ee84a196d89cp-114 }, + { -0x1.f31305b4f263ep-64, 0x1.40fbdbf9d128ap-118 }, + { 0x1.3db9842b31607p-67, 0x1.ee1ac155bab3ap-122 } }; + const double sc = 0x1p+8; + double zl, zh = fasttwosum (x0[0] + sx, x0[1], &zl); + zl += x0[2]; + double sh = zh * sc, sl = zl * sc; + int n = array_length (c), k = 7; + fl = sh * polyd (sh, k, c + n - k); + fh = polydd (sh, sl, n - k, c, &fl); + fh = muldd (zh, zl, fh, fl, &fl); + } + else if (fabs (fh) < 0x1p-3 && sx > -4.5 && sx < -4) + { + static const double x0[] + = { 0x1.0284e78599581p+2, -0x1.e78c1e9e43cfep-53, + 0x1.2ac17bfd6be92p-108 }; + static const double c[][2] + = { { 0x1.aca5cf4921642p+4, 0x1.a46a2e0d8fdfdp-51 }, + { 0x1.44415cd813f8ep+1, 0x1.afdc2672876f4p-56 }, + { 0x1.559b11b2a9c7cp-2, 0x1.17b8ada982c18p-57 }, + { 0x1.96d18e21aebdbp-5, -0x1.c2f2da6bea629p-62 }, + { 0x1.0261eb5732e4p-7, 0x1.3910ed6591782p-61 }, + { 0x1.55e3dbf99eb3dp-10, -0x1.e2d04c05513dcp-64 }, + { 0x1.d14fe49c4e437p-13, -0x1.d3009aa4e3bb9p-67 }, + { 0x1.433dce282da6ep-15, -0x1.6a33c9bb3cb85p-70 }, + { 0x1.c8399c7588ccfp-18, 0x1.6b47c3c6f3b7ap-73 }, + { 0x1.45fbe666d9415p-20, 0x1.9aaef873df6f7p-75 }, + { 0x1.d68d794caf0cfp-23, -0x1.62f893b4b1c92p-80 }, + { 0x1.56729dc75d1c6p-25, -0x1.9e28ba9f5e686p-79 }, + { 0x1.f5ec3352af509p-28, 0x1.63d11f1c09f18p-82 }, + { 0x1.7205756353fb8p-30, -0x1.e4315fc28841bp-85 }, + { 0x1.122e7755f935bp-32, 0x1.e94449e4ce532p-90 }, + { 0x1.982503f30184cp-35, 0x1.db16eb0f7f6e4p-89 }, + { 0x1.30f8c448c65bcp-37, -0x1.ba5e1bf8c5346p-95 }, + { 0x1.c957f8be0461cp-40, 0x1.78aa87d33801bp-95 }, + { 0x1.57fe1f7e3454fp-42, -0x1.61a8ccec45dcap-96 }, + { 0x1.035fcc6efbe53p-44, 0x1.454e6ac634253p-98 }, + { 0x1.87856bbe691edp-47, -0x1.92c14fb55274p-105 }, + { 0x1.2b877971688d2p-49, 0x1.c72aa6bdd8e43p-104 }, + { 0x1.e30abef69866bp-52, 0x1.4e95354428576p-107 }, + { 0x1.3f5a80dfa357ep-54, 0x1.dc7c678b87e2ep-108 } }; + const double sc = 0x1p+7; + double zl, zh = fasttwosum (x0[0] + sx, x0[1], &zl); + zl += x0[2]; + double sh = zh * sc, sl = zl * sc; + int n = array_length (c), k = 7; + fl = sh * polyd (sh, k, c + n - k); + fh = polydd (sh, sl, n - k, c, &fl); + fh = muldd (zh, zl, fh, fl, &fl); + } + else if (fabs (fh) < 0x1.08p-4 && sx > -5 && sx < -4.5) + { + static const double x0[] + = { 0x1.3f7577a6eeafdp+2, -0x1.5de5eab7f12cfp-53, + 0x1.4075f5e0494a2p-110 }; + static const double c[][2] + = { { -0x1.d224a3ef9e41fp+6, -0x1.9be272a13bcb6p-48 }, + { 0x1.b533c678a3956p+2, -0x1.37da6a2b62e3cp-53 }, + { -0x1.0d3f7fee65d34p-1, 0x1.e68bf720db9bep-55 }, + { 0x1.752a6f5ac2726p-5, -0x1.16f2fc46d2a34p-62 }, + { -0x1.13d5d163bd3f7p-8, -0x1.813df737e6e21p-62 }, + { 0x1.a8c5c53458ca5p-12, 0x1.02a507fa1273dp-66 }, + { -0x1.5068b3ed69409p-15, -0x1.8fc69d833eb31p-72 }, + { 0x1.0ffa575ea7fe3p-18, 0x1.ab6c92ac78cc4p-75 }, + { -0x1.bec12dd78a23p-22, 0x1.81392cb5ee3ap-76 }, + { 0x1.7382570f0b3c5p-25, 0x1.6c165015a71c3p-79 }, + { -0x1.380ebf6161ccep-28, -0x1.efc3933cb1d2cp-83 }, + { 0x1.084de43d1a947p-31, -0x1.1fca46f4984d6p-85 }, + { -0x1.c2d90dead7a34p-35, 0x1.b9c200207ddp-89 }, + { 0x1.82d0afe1a76ap-38, -0x1.958789fbe4df9p-93 }, + { -0x1.4d93cbb2d04c5p-41, -0x1.411bc9bf3a7ccp-95 }, + { 0x1.20ea504f0c629p-44, -0x1.1c1a63ab9daddp-100 }, + { -0x1.f6cb7ebb79bc2p-48, 0x1.f2f3b682d06fep-106 }, + { 0x1.be565af661ea4p-51, 0x1.b54b6ecaea8d5p-108 }, + { -0x1.78510a36bad6cp-54, -0x1.6935e8745afe8p-109 } }; + const double sc = 0x1p+10; + double zl, zh = fasttwosum (x0[0] + sx, x0[1], &zl); + zl += x0[2]; + double sh = zh * sc, sl = zl * sc; + int n = array_length (c), k = 7; + fl = sh * polyd (sh, k, c + n - k); + fh = polydd (sh, sl, n - k, c, &fl); + fh = muldd (zh, zl, fh, fl, &fl); + } + else if (fabs (fh) < 0x1.1p-4 && sx > -5.5 && sx < -5) + { + static const double x0[] + = { 0x1.4086a57f0b6d9p+2, 0x1.95262b72ca9cap-55, + 0x1.bd98d5e0861aap-109 }; + static const double c[][2] + = { { 0x1.ed72e0829ae02p+6, -0x1.fdc1859ae9c53p-50 }, + { 0x1.cecc32ec22f9bp+2, 0x1.b6ecc779b86fcp-53 }, + { 0x1.253d8563f7264p-1, -0x1.5cd273fbd9f05p-55 }, + { 0x1.a225df2da6e63p-5, -0x1.fe9d0b1b4c66p-59 }, + { 0x1.3e01773762671p-8, -0x1.f0e1184c71c78p-62 }, + { 0x1.f7d8d5bdcb186p-12, -0x1.d3eb6b34618fbp-66 }, + { 0x1.9a8d00c77a92cp-15, -0x1.107ec223e7938p-70 }, + { 0x1.557fd8c490b43p-18, 0x1.d8ff4f1780ac1p-72 }, + { 0x1.209221a624132p-21, -0x1.990b9a98394afp-75 }, + { 0x1.edc98d3bbeedfp-25, 0x1.58dea3f0bc5b9p-79 }, + { 0x1.aabd28e6c9a81p-28, -0x1.da694af0e0ffp-82 }, + { 0x1.73de2dd1cfd61p-31, 0x1.4f223cad85a45p-85 }, + { 0x1.4651832f82e2ep-34, 0x1.98ad9bf65ef61p-92 }, + { 0x1.200d84be7e1fcp-37, -0x1.5a13d3eb64f54p-93 }, + { 0x1.ff278dbfe563p-41, 0x1.7119b25262c2dp-100 }, + { 0x1.c77e7644c5411p-44, 0x1.c21cc9e02459cp-99 }, + { 0x1.97c73d128923ap-47, -0x1.155bb925610d9p-102 }, + { 0x1.74798cd59653bp-50, -0x1.90790da400208p-104 }, + { 0x1.436d37cbd1ac2p-53, 0x1.34f1a6c134053p-107 } }; + const double sc = 0x1p+10; + double zl, zh = fasttwosum (x0[0] + sx, x0[1], &zl); + zl += x0[2]; + double sh = zh * sc, sl = zl * sc; + int n = array_length (c), k = 6; + fl = sh * polyd (sh, k, c + n - k); + fh = polydd (sh, sl, n - k, c, &fl); + fh = muldd (zh, zl, fh, fl, &fl); + } + else if (fabs (fh) < 0x1.13p-3 && sx > -6 && sx < -5.5) + { + static const double x0[] + = { 0x1.7fe92f591f40dp+2, 0x1.7dd4ed62cbd32p-52, + -0x1.2071c071a2146p-108 }; + static const double c[][2] + = { { -0x1.661f6a43a5e12p+9, -0x1.0c437b83bc0ddp-45 }, + { 0x1.f79dcb794f26fp+5, -0x1.ada8018ade6efp-52 }, + { -0x1.d6e8088a19ffep+2, -0x1.2c08708631269p-53 }, + { 0x1.ef5d308dbfc97p-1, 0x1.87cd5c0d349fap-58 }, + { -0x1.15ea6b0ab529ep-3, 0x1.4ac7004538519p-66 }, + { 0x1.44d54e9fe2397p-6, 0x1.f1ca84e1a96cp-62 }, + { -0x1.8684e40cebb3dp-9, -0x1.19803482caec7p-64 }, + { 0x1.df44c1d81c723p-12, -0x1.0ed2e1ad825e8p-67 }, + { -0x1.2ac3053f4ee18p-14, -0x1.fe50dd710dd6bp-68 }, + { 0x1.79226ae04a847p-17, -0x1.379839e0000d6p-75 }, + { -0x1.e0dffb5f77ccbp-20, 0x1.a4015df3482b4p-74 }, + { 0x1.35217890708e7p-22, 0x1.2875bd051d55ep-78 }, + { -0x1.903aa9af5cf03p-25, -0x1.a619f5e70dc34p-83 }, + { 0x1.04a103323cd6ap-27, 0x1.a260aa452be23p-83 }, + { -0x1.552efc1ecd295p-30, -0x1.88a83f178c25p-84 }, + { 0x1.c0a12b600b141p-33, -0x1.c3732154b52acp-89 }, + { -0x1.281ce75ae916fp-35, -0x1.a96d8e0798081p-89 }, + { 0x1.88411a65d3af4p-38, -0x1.0e8c83eb281a9p-92 }, + { -0x1.049e79a4b5baap-40, -0x1.64caa9cfd44a2p-94 }, + { 0x1.5b0f16f419f7bp-43, 0x1.74e160750d876p-102 }, + { -0x1.ce70ccad93bc1p-46, 0x1.1b340ff57e754p-100 }, + { 0x1.3a7c53f75547fp-48, -0x1.88a8d5ad25f3fp-102 }, + { -0x1.c41c8d3e4c2c9p-51, 0x1.e4dd9e959c88bp-105 }, + { 0x1.f82da393ec32ap-54, -0x1.65c64b7f7e325p-115 } }; + const double sc = 0x1p+12; + double zl, zh = fasttwosum (x0[0] + sx, x0[1], &zl); + zl += x0[2]; + double sh = zh * sc, sl = zl * sc; + int n = array_length (c), k = 6; + fl = sh * polyd (sh, k, c + n - k); + fh = polydd (sh, sl, n - k, c, &fl); + fh = muldd (zh, zl, fh, fl, &fl); + } + else if (fabs (fh) < 0x1.15p-3 && sx > -6.5 && sx < -6) + { + static const double x0[] + = { 0x1.8016b25897c8dp+2, -0x1.27e0f49a4ba72p-54, + 0x1.72e1ab15a4d03p-110 }; + static const double c[][2] + = { { 0x1.69de49e3af2aap+9, 0x1.954b690943afcp-47 }, + { 0x1.fce23484cfd1p+5, 0x1.8266e7580f08cp-49 }, + { 0x1.de503a3c37c4p+2, 0x1.9fa7459d62427p-53 }, + { 0x1.f9c7b52558abbp-1, 0x1.b6896749474c1p-55 }, + { 0x1.1d3d50714416ap-3, 0x1.56020c7693337p-58 }, + { 0x1.4f21e2fb9e06p-6, 0x1.9e1cfc6a7c13p-60 }, + { 0x1.9500994cd8a9ep-9, -0x1.ce1777b097bd6p-63 }, + { 0x1.f3a2c23c19d79p-12, -0x1.972494ba4cd4dp-70 }, + { 0x1.39152652eb3aap-14, 0x1.8f85c4250ca8p-69 }, + { 0x1.8d45f8be891dep-17, 0x1.41e12761b265ap-71 }, + { 0x1.fd3214a702b29p-20, -0x1.140076f8887f9p-74 }, + { 0x1.490b476815d46p-22, -0x1.3c92eb9376d39p-81 }, + { 0x1.ac3b965280088p-25, -0x1.468539ffe88b5p-82 }, + { 0x1.1851c43c40d43p-27, -0x1.6918e4fb64554p-82 }, + { 0x1.70dfb529f2c5ap-30, -0x1.b5dfdc2a8890dp-84 }, + { 0x1.e791ef4ef8602p-33, 0x1.b20188eabd9e6p-90 }, + { 0x1.437e60f850abep-35, 0x1.b5786abafa12ep-89 }, + { 0x1.aec29bcf6c8c1p-38, 0x1.774aed905f5cfp-92 }, + { 0x1.1fb22834000f1p-40, 0x1.705349812bcbap-95 }, + { 0x1.811ce4ec7b869p-43, 0x1.8be5e35bce704p-99 }, + { 0x1.01e72b00fd763p-45, 0x1.c8c867f1090e3p-99 }, + { 0x1.60a0a9eee14d8p-48, 0x1.ed307bf965e72p-106 }, + { 0x1.fdce60428ba59p-51, 0x1.268eaf9e15016p-106 }, + { 0x1.1dcf4e3ac0033p-53, 0x1.bdbdd6256584p-108 } }; + const double sc = 0x1p+12; + double zl, zh = fasttwosum (x0[0] + sx, x0[1], &zl); + zl += x0[2]; + double sh = zh * sc, sl = zl * sc; + int n = array_length (c), k = 6; + fl = sh * polyd (sh, k, c + n - k); + fh = polydd (sh, sl, n - k, c, &fl); + fh = muldd (zh, zl, fh, fl, &fl); + } + else if (fabs (fh) < 0x1.34p-3 && sx > -7.0 && sx < -6.5) + { + static const double x0[] + = { 0x1.bffcbf76b86fp+2, -0x1.853b29347b806p-57, + 0x1.0fa018051dd41p-111 }; + static const double c[][2] + = { { -0x1.3abf7a5cea91bp+12, -0x1.8257b8abd02cep-42 }, + { 0x1.8349a2550422dp+9, -0x1.c6f2ef4105b99p-45 }, + { -0x1.3d91dadc98428p+7, 0x1.46605ff1fc41cp-48 }, + { 0x1.24f3d636f3339p+5, 0x1.5966a8447ea56p-49 }, + { -0x1.20427df1b3492p+3, -0x1.e9801a53bb06ep-52 }, + { 0x1.2775e857fb69cp+1, 0x1.88bd494f9ea51p-53 }, + { -0x1.377e70b463c13p-1, -0x1.2524c9cdd93dbp-55 }, + { 0x1.4f3d28edba5cdp-3, 0x1.92d2ac40dac77p-60 }, + { -0x1.6e8557168cf8ap-5, 0x1.0b5cf2a57fc71p-59 }, + { 0x1.95bb17ce427a4p-7, -0x1.0f3a87df6aaa1p-66 }, + { -0x1.c5ac12d48f7e4p-9, 0x1.1d042aabfd8fbp-63 }, + { 0x1.ff816dad74e75p-11, 0x1.24594bf64f76dp-67 }, + { -0x1.225f4a6a0c5b4p-12, -0x1.e1dd93bdafddap-66 }, + { 0x1.4ba3e5c007209p-14, 0x1.a47689fcd5098p-70 }, + { -0x1.7cb738040762cp-16, -0x1.f100cc43c4b52p-70 }, + { 0x1.b701200fca8b8p-18, 0x1.d825b96ec4cebp-72 }, + { -0x1.fc33bee21552cp-20, 0x1.096b388c1821dp-75 }, + { 0x1.272cc14fa7a53p-21, 0x1.d22fcb28979f7p-76 }, + { -0x1.57f61caeb3e48p-23, -0x1.03971054f5156p-79 }, + { 0x1.91f102889ab68p-25, -0x1.b9949bdf154e2p-80 }, + { -0x1.d641c300c655ep-27, -0x1.47bf835c77fc5p-82 }, + { 0x1.1319b4831f0ep-28, -0x1.2688b8ef3d8f7p-84 }, + { -0x1.4c01654897d29p-30, -0x1.5c8d35e08cb51p-86 }, + { 0x1.a7ea9eb248f1ep-32, -0x1.346b798e3ef55p-86 }, + { -0x1.8afe498da415bp-34, 0x1.3f02ddea04afbp-88 } }; + const double sc = 0x1p+14; + double zl, zh = fasttwosum (x0[0] + sx, x0[1], &zl); + zl += x0[2]; + double sh = zh * sc, sl = zl * sc; + int n = array_length (c), k = 6; + fl = sh * polyd (sh, k, c + n - k); + fh = polydd (sh, sl, n - k, c, &fl); + fh = muldd (zh, zl, fh, fl, &fl); + } + else if (fabs (fh) < 0x1.145p-3 && sx > -7.5 && sx < -7) + { + static const double x0[] + = { 0x1.c0033fdedfe1fp+2, -0x1.20bb7d2324678p-52, + -0x1.f5536678d69d3p-106 }; + static const double c[][2] + = { { 0x1.3b407aa387bd1p+12, 0x1.da1e57343b1dbp-43 }, + { 0x1.83e85daafbad6p+8, -0x1.f37538d99bd29p-46 }, + { 0x1.3e552b5e3c226p+5, -0x1.07b1550d11cf4p-49 }, + { 0x1.25e42a45e905bp+2, 0x1.61a63054c47c5p-54 }, + { 0x1.216a3560743eep-1, 0x1.f5024887fdb9dp-57 }, + { 0x1.28e1c70ef5313p-4, 0x1.3234af5080105p-59 }, + { 0x1.393e2bc330081p-7, -0x1.45ef9ada326bap-62 }, + { 0x1.5164141f5ae6ap-10, -0x1.802404a8fab89p-65 }, + { 0x1.712b3a86e1bdfp-13, -0x1.b2f4283e8796ep-69 }, + { 0x1.98fd36b906e08p-16, -0x1.c84ad8d8481b7p-70 }, + { 0x1.c9ae6ef6262f6p-19, 0x1.ff0ee333033bfp-73 }, + { 0x1.02382a958b63dp-21, -0x1.99f49c2b45652p-75 }, + { 0x1.256845ec9b6b1p-24, 0x1.822bd9a29ab32p-80 }, + { 0x1.4f5ff36087486p-27, -0x1.dc262e92f86f3p-82 }, + { 0x1.814f9ce48442fp-30, 0x1.544e6fd9ac4bep-84 }, + { 0x1.bca89f844ac84p-33, -0x1.63ceb4f54c669p-87 }, + { 0x1.01946c9b7b26ap-35, -0x1.9b4cb6a88b6fcp-90 }, + { 0x1.2b75936f37bc3p-38, -0x1.d8809500cd622p-94 }, + { 0x1.5d3d0a248647p-41, -0x1.061c56c047c1dp-95 }, + { 0x1.98297f709854fp-44, 0x1.4cc4272ef2e04p-103 }, + { 0x1.dd4cc8e9d77b5p-47, 0x1.0103a0db6c049p-102 }, + { 0x1.1ce46c495fa62p-49, 0x1.622087cd95f2ap-103 }, + { 0x1.678ec65743609p-52, -0x1.850d52cb53a9fp-107 }, + { 0x1.5fe290153bcaap-55, 0x1.a8b6414621298p-109 } }; + const double sc = 0x1p+15; + double zl, zh = fasttwosum (x0[0] + sx, x0[1], &zl); + zl += x0[2]; + double sh = zh * sc, sl = zl * sc; + int n = array_length (c), k = 6; + fl = sh * polyd (sh, k, c + n - k); + fh = polydd (sh, sl, n - k, c, &fl); + fh = muldd (zh, zl, fh, fl, &fl); + } + else if (fabs (fh) < 0x1.76bp-3 && sx > -8 && sx < -7.5) + { + static const double x0[] + = { 0x1.ffff97f8159cfp+2, 0x1.e54f415a91586p-55, + 0x1.53a5d106f9a3ep-109 }; + static const double c[][2] + = { { -0x1.3af76fe4c2fabp+15, -0x1.7cc92f0b999ep-40 }, + { 0x1.838e76caaf123p+12, 0x1.292e15f502d1fp-42 }, + { -0x1.3de68b3256526p+10, 0x1.5456a490ea6e3p-44 }, + { 0x1.255c052530c71p+8, -0x1.6700425996a73p-48 }, + { -0x1.20c2a8418126ap+6, 0x1.1d5e3d82023c1p-50 }, + { 0x1.28139342cefp+4, 0x1.0238cfdd7475dp-51 }, + { -0x1.384066c322246p+2, 0x1.ce82f225b48d1p-53 }, + { 0x1.502bc4dad47d3p+0, 0x1.f45416fab08dap-58 }, + { -0x1.6faadfece0e3p-2, -0x1.6c5874c40cd48p-56 }, + { 0x1.9724323c89909p-4, -0x1.27460acd35293p-58 }, + { -0x1.c7684c96f243p-6, 0x1.5d749ebbd28fep-62 }, + { 0x1.00d1f48743fb8p-7, 0x1.ad60125333897p-61 }, + { -0x1.23af6dd470e5cp-9, -0x1.1d084246f0edap-64 }, + { 0x1.4d416958eec41p-11, -0x1.dfeda2a7c037ep-67 }, + { -0x1.7eb3eb405269fp-13, -0x1.95169511f695dp-69 }, + { 0x1.b972ec5a016fp-15, 0x1.0a9a114630056p-69 }, + { -0x1.ff35af33eaacp-17, -0x1.73922358a5b41p-71 }, + { 0x1.290662f9abea1p-18, -0x1.8a5267d544bb7p-72 }, + { -0x1.5a395633f9d69p-20, -0x1.377120c488d2ep-76 }, + { 0x1.94b2a8f6b566ep-22, 0x1.685bd8f7318c7p-79 }, + { -0x1.da517a6c7ab89p-24, -0x1.940f32a9a6413p-78 }, + { 0x1.168992efa74a9p-25, -0x1.bb9b7bf492795p-81 }, + { -0x1.4665c665b014ep-27, 0x1.91ef4857aca34p-81 }, + { 0x1.7f49e136d0e91p-29, 0x1.9d7da1584fe9p-85 }, + { -0x1.df48e4db635efp-31, -0x1.ac6dd06d9a112p-85 }, + { 0x1.3665697708ddp-32, -0x1.e15f660c03e5bp-90 }, + { -0x1.015f65348c21ap-34, -0x1.e45b91e0267cp-90 } }; + const double sc = 0x1p+17; + double zl, zh = fasttwosum (x0[0] + sx, x0[1], &zl); + zl += x0[2]; + double sh = zh * sc, sl = zl * sc; + int n = array_length (c), k = 6; + fl = sh * polyd (sh, k, c + n - k); + fh = polydd (sh, sl, n - k, c, &fl); + fh = muldd (zh, zl, fh, fl, &fl); + } + else if (fabs (fh) < 0x1.76cp-3 && sx > -0x1.10f5c28f5c28fp+3 + && sx < -8) + { + static const double x0[] + = { 0x1.000034028b3f9p+3, 0x1.f60cb3cec1cedp-52, + -0x1.ea26620d6b1cap-106 }; + static const double c[][2] + = { { 0x1.3b088fed67718p+15, -0x1.505613ba29893p-39 }, + { 0x1.83a3893550edcp+12, 0x1.f52e3b240db3p-42 }, + { 0x1.3e0078db8ada4p+10, 0x1.506573eecc6e3p-44 }, + { 0x1.257bec9464251p+8, 0x1.8c4ea8394aa1ep-49 }, + { 0x1.20e9ea0755a47p+6, -0x1.978e59e1a3d9fp-48 }, + { 0x1.2843e1313c83bp+4, -0x1.50610d6737717p-55 }, + { 0x1.387bd6a785478p+2, -0x1.1f1fff3fa1b4bp-52 }, + { 0x1.5074e788de77p+0, 0x1.aeeb83580ea25p-55 }, + { 0x1.7004dd990d7dap-2, 0x1.3f7f554ffc67p-59 }, + { 0x1.9792ed5f6dfb4p-4, -0x1.1f23ec3a02887p-58 }, + { 0x1.c7f08cdaef32fp-6, -0x1.1e18e4c767dep-62 }, + { 0x1.0125c81121f59p-7, -0x1.a6731270f4d3fp-61 }, + { 0x1.2416931f22d22p-9, -0x1.2bc2867243e45p-64 }, + { 0x1.4dc0543bea659p-11, 0x1.635ca3a9b1a8dp-66 }, + { 0x1.7f501645b2fdep-13, -0x1.e45b9db4c768p-67 }, + { 0x1.ba331549d971dp-15, 0x1.97379411be4a6p-69 }, + { 0x1.001110caa210fp-16, 0x1.68ce3f9f57bc4p-70 }, + { 0x1.2997db55cba74p-18, -0x1.0f0fa5ac41b0ap-74 }, + { 0x1.5aec550050044p-20, 0x1.7fdca933bee56p-74 }, + { 0x1.958ee8d69504p-22, 0x1.c28e192bc3ba4p-76 }, + { 0x1.db608c6ccccc8p-24, -0x1.71a440ef759bep-80 }, + { 0x1.173059c7e38cfp-25, 0x1.40d18009e343fp-79 }, + { 0x1.4731fe682342cp-27, -0x1.535cdcf9e608bp-83 }, + { 0x1.8043ec54aba15p-29, -0x1.62d6d36f1d19ep-84 }, + { 0x1.e08fc2d5d4ap-31, -0x1.6e4b1587c874ep-85 }, + { 0x1.3742f9e468b63p-32, 0x1.4c5c5de4449bfp-88 }, + { 0x1.021df5d9d312dp-34, 0x1.b2da076b6955fp-90 } }; + const double sc = 0x1p+17; + double zl, zh = fasttwosum (x0[0] + sx, x0[1], &zl); + zl += x0[2]; + double sh = zh * sc, sl = zl * sc; + int n = array_length (c), k = 6; + fl = sh * polyd (sh, k, c + n - k); + fh = polydd (sh, sl, n - k, c, &fl); + fh = muldd (zh, zl, fh, fl, &fl); + } + else if (fabs (fh) < 0x1.99p-3 && sx > -9 && sx < -8.5) + { + static const double x0[] + = { 0x1.1ffffa3884bdp+3, 0x1.ff90c9d2ae925p-53, + -0x1.30c0efef78c04p-107 }; + static const double c[][2] + = { { -0x1.625edfc63db2fp+18, 0x1.da7fc3ed68f99p-37 }, + { 0x1.ea8c150480a7ap+15, 0x1.344e4cbf6d2a1p-39 }, + { -0x1.c4b30e4bc55c1p+13, -0x1.9ec40fdd810d7p-41 }, + { 0x1.d5fe468dbbf03p+11, -0x1.80705f1856688p-43 }, + { -0x1.043d21bc24decp+10, -0x1.b0db95c85f28bp-45 }, + { 0x1.2c334ae535e1dp+8, 0x1.530b0b83f8e2p-46 }, + { -0x1.64314b431cd64p+6, -0x1.5d7e4374b5c05p-49 }, + { 0x1.af6ed589b3a86p+4, -0x1.622e836ea2287p-51 }, + { -0x1.096e446edcfb4p+3, 0x1.f849c504338f9p-51 }, + { 0x1.4aaf49e713c0cp+1, 0x1.bbab832e792a1p-53 }, + { -0x1.a0246d9c1b57fp-1, -0x1.18f4c5d12b8b5p-56 }, + { 0x1.0806315c1a365p-2, 0x1.074a0bafc599fp-58 }, + { -0x1.515dd6b891094p-4, -0x1.bcca3edb6dc58p-59 }, + { 0x1.b1a5fe76de3d4p-6, -0x1.f8d7b8e44c315p-63 }, + { -0x1.182229539b98cp-7, 0x1.5c2ae71e11bfep-62 }, + { 0x1.6b8b31630a98cp-9, -0x1.c5c500af742dap-65 }, + { -0x1.d9a3c6789b8e8p-11, 0x1.7160288934ba5p-65 }, + { 0x1.359c2ea889b3ep-12, 0x1.e44b62289b60bp-69 }, + { -0x1.9606a9f683511p-14, -0x1.3a3c7889c7af5p-68 }, + { 0x1.0af84fb66238ep-15, -0x1.3d487a70215dap-69 }, + { -0x1.5ff88107d908ap-17, 0x1.0c4172ceada86p-73 }, + { 0x1.d14070a537453p-19, 0x1.964ddb3dba2bap-73 }, + { -0x1.33eaa64e4654ep-20, 0x1.e6deb19e51b4ep-74 }, + { 0x1.957ef375ea2abp-22, 0x1.065699585b36cp-77 }, + { -0x1.0cd5c9d34a242p-23, -0x1.6f2ed041b79a1p-77 }, + { 0x1.836d31552e5d6p-25, -0x1.789b28b32b093p-79 }, + { -0x1.19aeee4b34ef6p-26, -0x1.c45d882aaaep-81 }, + { 0x1.eb1912eaced35p-29, 0x1.e9842f288e554p-84 } }; + const double sc = 0x1p+20; + double zl, zh = fasttwosum (x0[0] + sx, x0[1], &zl); + zl += x0[2]; + double sh = zh * sc, sl = zl * sc; + int n = array_length (c), k = 6; + fl = sh * polyd (sh, k, c + n - k); + fh = polydd (sh, sl, n - k, c, &fl); + fh = muldd (zh, zl, fh, fl, &fl); + } + else if (fabs (fh) < 0x1.99p-3 && sx > -9.5 && sx < -9) + { + static const double x0[] + = { 0x1.200005c7768fbp+3, 0x1.b5b610ffb70d4p-54, + 0x1.deb7ad09ec5eap-108 }; + static const double c[][2] + = { { 0x1.626120391944p+18, 0x1.7d5e8272ce41dp-38 }, + { 0x1.ea8f32fb7f586p+15, -0x1.345b1cc20d4a5p-39 }, + { 0x1.c4b75ee68e2bap+13, -0x1.812d7bcedccb5p-42 }, + { 0x1.d6043fa1ffaa5p+11, -0x1.5a4eaff5f008ap-43 }, + { 0x1.04414411db7f4p+10, 0x1.d742c54011402p-44 }, + { 0x1.2c3903ec9c90cp+8, 0x1.73da27a4b1b2p-46 }, + { 0x1.64393744bb9bdp+6, -0x1.482bcb016f26fp-48 }, + { 0x1.af79ccdc71d33p+4, 0x1.0a7a439456745p-50 }, + { 0x1.0975db7d71fc6p+3, 0x1.bd477ffdd39f6p-51 }, + { 0x1.4ab9cba1e3478p+1, 0x1.1c929f5a50618p-54 }, + { 0x1.a032f8f114635p-1, 0x1.92279b6d57eb2p-55 }, + { 0x1.0810426bfa5a1p-2, 0x1.0cd3cd3381e06p-56 }, + { 0x1.516bc616eaf53p-4, -0x1.0ebc405b4871ap-58 }, + { 0x1.b1b948b11a02dp-6, 0x1.8b0ffc5af8ce1p-60 }, + { 0x1.182f8343c9e61p-7, -0x1.de95c6fadf8bbp-61 }, + { 0x1.6b9dacc2e40dp-9, 0x1.7d2afb82825d1p-63 }, + { 0x1.d9bd5c016a05ap-11, 0x1.afdd41ec7d841p-65 }, + { 0x1.35ade3d85cd0ep-12, -0x1.5080702a55ab7p-66 }, + { 0x1.961f2d2659b02p-14, 0x1.f377d9fd7df3ap-68 }, + { 0x1.0b0946f97b203p-15, 0x1.9ecb24ffa7674p-69 }, + { 0x1.600ffd56f5733p-17, 0x1.2ffa0ff103bfdp-71 }, + { 0x1.d160f5b944255p-19, -0x1.6f04901227877p-75 }, + { 0x1.3401289d10771p-20, -0x1.b2fca86443a23p-76 }, + { 0x1.959dec69f88edp-22, 0x1.1945258ad2d8p-76 }, + { 0x1.0ceb1d5fdcfdbp-23, -0x1.e1bf43eb8d61fp-77 }, + { 0x1.838cdbbcfd467p-25, -0x1.948ada050fa52p-80 }, + { 0x1.19c71e9accac3p-26, 0x1.74cb4459f6b9p-81 }, + { 0x1.eb4678e3fcf16p-29, -0x1.66454030f804fp-83 } }; + const double sc = 0x1p+20; + double zl, zh = fasttwosum (x0[0] + sx, x0[1], &zl); + zl += x0[2]; + double sh = zh * sc, sl = zl * sc; + int n = array_length (c), k = 6; + fl = sh * polyd (sh, k, c + n - k); + fh = polydd (sh, sl, n - k, c, &fl); + fh = muldd (zh, zl, fh, fl, &fl); + } + else if (fabs (fh) < 0x1.76bp-3 && sx > -10 && sx < -9.5) + { + static const double x0[] + = { 0x1.3fffff6c0d7cp+3, -0x1.197cea8c42d7dp-51, + -0x1.7072c5a292198p-105 }; + static const double c[][2] + = { { -0x1.baf7da5f3795dp+21, -0x1.16a79518c8367p-33 }, + { 0x1.7f3e8791fa0d2p+18, -0x1.2aec811c9609ep-36 }, + { -0x1.ba18befcaaa63p+15, -0x1.d18c4e260465p-39 }, + { 0x1.1ede14765dc0cp+13, 0x1.13bc952384f0fp-41 }, + { -0x1.8d1a9ab5a505p+10, -0x1.904ee4f738385p-46 }, + { 0x1.1e4d8c35d22ccp+8, -0x1.3b56f0fecec54p-48 }, + { -0x1.a8a191db109p+5, -0x1.11eec467fa814p-51 }, + { 0x1.4174f65ff868p+3, 0x1.c939d4bfd1404p-51 }, + { -0x1.ee6d90f2332c7p+0, 0x1.3b18f44783de9p-55 }, + { 0x1.80fd3420fba09p-2, 0x1.55812b31052fbp-57 }, + { -0x1.2ecd481762eaep-4, -0x1.d974aa28e0c7p-62 }, + { 0x1.e04a0b28db4c5p-7, -0x1.f0d1f410ed5c5p-61 }, + { -0x1.7f91af3f1200fp-9, 0x1.da2dc770582ddp-63 }, + { 0x1.342652fcfe1b6p-11, -0x1.7744fdcaebcfap-66 }, + { -0x1.f1a88f796a348p-14, 0x1.a8bb92c660ef4p-74 }, + { 0x1.93a68769f8a0fp-16, -0x1.f24bb59f0e5dap-70 }, + { -0x1.48af467f6d621p-18, 0x1.7948f005903cap-74 }, + { 0x1.0c92181e40e32p-20, 0x1.65951755e7adfp-74 }, + { -0x1.b842459ce40bfp-23, 0x1.4503ca5bb9422p-77 }, + { 0x1.69db729e3b047p-25, -0x1.c0423483dea32p-79 }, + { -0x1.2a3762889a1ebp-27, -0x1.c935e18265a7p-81 }, + { 0x1.ec8fbfde50a33p-30, 0x1.98aa5e946ac0fp-89 }, + { -0x1.95dcfae4547bcp-32, -0x1.c04aac4c4ba01p-87 }, + { 0x1.4f21457a90c36p-34, -0x1.e01d7e3277c24p-88 }, + { -0x1.26aacbe779418p-36, 0x1.e6b067bac93acp-92 }, + { 0x1.0c603e5994262p-38, 0x1.d6144d3f2ef17p-93 }, + { -0x1.38f713b1343d3p-41, -0x1.0c3eef336946bp-95 } }; + const double sc = 0x1p+24; + double zl, zh = fasttwosum (x0[0] + sx, x0[1], &zl); + zl += x0[2]; + double sh = zh * sc, sl = zl * sc; + int n = array_length (c), k = 6; + fl = sh * polyd (sh, k, c + n - k); + fh = polydd (sh, sl, n - k, c, &fl); + fh = muldd (zh, zl, fh, fl, &fl); + } + else if (fabs (fh) < 0x1.76bp-3 && sx > -10.5 && sx < -10) + { + static const double x0[] + = { 0x1.40000093f2777p+3, 0x1.927b45d95e154p-52, + 0x1.0780c21b6e452p-106 }; + static const double c[][2] + = { { 0x1.baf825a0c63b2p+21, -0x1.20323f1015cdfp-35 }, + { 0x1.7f3ec8ae05f2ep+18, 0x1.2aec80d23ccb8p-36 }, + { 0x1.ba192fa62a5c8p+15, -0x1.25660ae99a81cp-39 }, + { 0x1.1ede75ef431bp+13, -0x1.a691cbd9c9bebp-41 }, + { 0x1.8d1b435ece20fp+10, 0x1.5b052a017c6d3p-47 }, + { 0x1.1e4e1e218c99cp+8, 0x1.775895e089534p-46 }, + { 0x1.a8a28e596cccep+5, -0x1.aa094ee6e63a1p-49 }, + { 0x1.4175d0d35b3d4p+3, -0x1.0cf2eb638feb3p-51 }, + { 0x1.ee6f0af10b985p+0, -0x1.7d3173e84c277p-54 }, + { 0x1.80fe7b2913e68p-2, 0x1.58aa6aa5a2a1dp-56 }, + { 0x1.2ece6307c7cbp-4, 0x1.812bbc72b9bf4p-58 }, + { 0x1.e04bf4be0259p-7, 0x1.b2fdc34bbffbep-62 }, + { 0x1.7f9356d1f8f5dp-9, 0x1.200577f2d649bp-63 }, + { 0x1.3427c173faa4ep-11, 0x1.963e715185901p-65 }, + { 0x1.f1ab0995dd7aap-14, -0x1.982b6a256e922p-70 }, + { 0x1.93a8ac07adde7p-16, 0x1.1938f0da933adp-70 }, + { 0x1.48b1212550a17p-18, -0x1.73810e7394decp-72 }, + { 0x1.0c93b2c55f9b3p-20, 0x1.b11c8bd5d36cap-75 }, + { 0x1.b8450c2eb24f2p-23, 0x1.b77d60943d51ap-83 }, + { 0x1.69ddd961f0eb9p-25, 0x1.26555b649e2b9p-80 }, + { 0x1.2a39767b036bdp-27, -0x1.ff547cf250c2cp-82 }, + { 0x1.ec935873525f3p-30, -0x1.56d39c71ef428p-84 }, + { 0x1.95e014a90402ep-32, 0x1.ff041edd80219p-88 }, + { 0x1.4f23f078ebd25p-34, 0x1.22e3191f74383p-89 }, + { 0x1.26ad387c7d522p-36, -0x1.0581274dbf44p-91 }, + { 0x1.0c628e2dc8b58p-38, 0x1.4473a911adb21p-97 }, + { 0x1.38f9f8f1d99a8p-41, -0x1.6893c733787dp-95 } }; + const double sc = 0x1p+24; + double zl, zh = fasttwosum (x0[0] + sx, x0[1], &zl); + zl += x0[2]; + double sh = zh * sc, sl = zl * sc; + int n = array_length (c), k = 6; + fl = sh * polyd (sh, k, c + n - k); + fh = polydd (sh, sl, n - k, c, &fl); + fh = muldd (zh, zl, fh, fl, &fl); + } + } + + unsigned ft = (asuint64 (fl) + 2) & (~UINT64_C(0) >> 12); + if (ft <= 2u) + return as_lgamma_database (sx, fh + fl); + return fh + fl; +} double -__ieee754_lgamma_r(double x, int *signgamp) +__ieee754_lgamma_r (double x, int *signgamp) { - double t,y,z,nadj,p,p1,p2,p3,q,r,w; - int i,hx,lx,ix; - - EXTRACT_WORDS(hx,lx,x); - - /* purge off +-inf, NaN, +-0, and negative arguments */ - *signgamp = 1; - ix = hx&0x7fffffff; - if(__builtin_expect(ix>=0x7ff00000, 0)) return x*x; - if(__builtin_expect((ix|lx)==0, 0)) - { - if (hx < 0) - *signgamp = -1; - return one/fabs(x); - } - if(__builtin_expect(ix<0x3b900000, 0)) { - /* |x|<2**-70, return -log(|x|) */ - if(hx<0) { - *signgamp = -1; - return -__ieee754_log(-x); - } else return -__ieee754_log(x); - } - if(hx<0) { - if(__builtin_expect(ix>=0x43300000, 0)) - /* |x|>=2**52, must be -integer */ - return fabs (x)/zero; - if (x < -2.0 && x > -28.0) - return __lgamma_neg (x, signgamp); - t = sin_pi(x); - if(t==zero) return one/fabs(t); /* -integer */ - nadj = __ieee754_log(pi/fabs(t*x)); - if(t=0x3FE76944) {y = one-x; i= 0;} - else if(ix>=0x3FCDA661) {y= x-(tc-one); i=1;} - else {y = x; i=2;} - } else { - r = zero; - if(ix>=0x3FFBB4C3) {y=2.0-x;i=0;} /* [1.7316,2] */ - else if(ix>=0x3FF3B4C4) {y=x-tc;i=1;} /* [1.23,1.73] */ - else {y=x-one;i=2;} + // piece-wise polynomial approximation in [0.5, 8.29541] range + // range borders + static const unsigned ubrd[20] + = { 0x1ff0000, 0x1ff146c, 0x1ff2b7b, 0x1ff4532, 0x1ff614c, + 0x1ff8310, 0x1ff93f7, 0x1ffa880, 0x1ffc05e, 0x1ffdb73, + 0x1fff8a5, 0x2001147, 0x2002703, 0x20041ac, 0x200622a, + 0x20084d9, 0x2009ce7, 0x200ba2c, 0x200ddd7, 0x20104ba }; + // the region offset i.e. z = x-off for P(z) + static const double offs[19] + = { 0x1.146cd8p-1, 0x1.3fe898p-1, 0x1.70aea8p-1, 0x1.a67fcp-1, + 0x1.e76db8p-1, 0x1.170838p+0, 0x1.3c78a8p+0, 0x1.68df2p+0, + 0x1.9bd14p+0, 0x1.d41868p+0, 0x1.0d9a64p+1, 0x1.384b8p+1, + 0x1.68b06p+1, 0x1.a3d6dp+1, 0x1.ebdd9p+1, 0x1.21c1p+2, + 0x1.571368p+2, 0x1.9803e8p+2, 0x1.e74cc8p+2 }; + // polynomial coefficients low part + static const double cl[19][8] = { + { -0x1.18ad63ca097e9p+2, 0x1.af8e15b715c51p+2, -0x1.56213b7191ba4p+3, + 0x1.151f165a9425fp+4, -0x1.c826426e4b7cdp+4, 0x1.7c313095e4b75p+5, + -0x1.44f3d7d848e78p+6, 0x1.13384c97ea99dp+7 }, + { -0x1.0f58e76c8d235p+1, 0x1.67c3f6b7124f6p+1, -0x1.ec78d7d8185a3p+1, + 0x1.588d6487de574p+2, -0x1.e9fbe8564220dp+2, 0x1.60dd913b80b5ep+3, + -0x1.0465db7c895a6p+4, 0x1.7ca34d903fc3p+4 }, + { -0x1.0c505555a86b2p+0, 0x1.33d3a22d1bb51p+0, -0x1.6d2a2457f05d4p+0, + 0x1.bb1c77fad8b03p+0, -0x1.115210a553746p+1, 0x1.558e305fd694p+1, + -0x1.b4f9a0654679fp+1, 0x1.1489e269cbf39p+2 }, + { -0x1.1170ead9585bap-1, 0x1.10c67b04495d7p-1, -0x1.19de3af9dd349p-1, + 0x1.2a34cd6e66472p-1, -0x1.40e2f93066eb8p-1, 0x1.5ddc21a559735p-1, + -0x1.860a742d837aap-1, 0x1.ace7238771e7fp-1 }, + { -0x1.be31df8f7d605p-3, 0x1.8e33a32c94cf7p-3, -0x1.6c62efd534bf3p-3, + 0x1.53719e404a7d6p-3, -0x1.4074d1d083331p-3, 0x1.31c6c5226f2b3p-3, + -0x1.2b062b8eedd9fp-3, 0x1.219431f82fbfcp-3 }, + { -0x1.168e45409b785p-3, 0x1.a04e5759477fp-4, -0x1.43c620bb1d77fp-4, + 0x1.027d79414ff7dp-4, -0x1.a46afc0776356p-5, 0x1.5a92c3ddb75f5p-5, + -0x1.23d37b3b3e3b6p-5, 0x1.f66a6169fe8efp-6 }, + { -0x1.2db051283fb7ap-4, 0x1.8afce072c9222p-5, -0x1.0dcc84e49a658p-5, + 0x1.7af09a263459bp-6, -0x1.0f51524188551p-6, 0x1.8a1f8bbd73d04p-7, + -0x1.24e389a08ab23p-7, 0x1.b5c2228d32783p-8 }, + { -0x1.3fbb4a9e75e6dp-5, 0x1.6c40332da72cp-6, -0x1.b235e2a6ed724p-7, + 0x1.0a9023dcf81a5p-7, -0x1.4e128ec28e27ap-8, 0x1.a90e4421d59e7p-9, + -0x1.14e6becdb3889p-9, 0x1.68e6a1727f763p-10 }, + { -0x1.5365e61675f08p-6, 0x1.4fd143859dc2cp-7, -0x1.5cb1d911bcabbp-8, + 0x1.75a869d793508p-9, -0x1.9941834996ea5p-10, 0x1.c782075b40e2cp-11, + -0x1.03ab2e70c9df2p-11, 0x1.26c9636a06719p-12 }, + { -0x1.7145b3bd2da75p-7, 0x1.3e6749a0fe63p-8, -0x1.20ea7ae3d0208p-9, + 0x1.0f17a25bb48f3p-10, -0x1.045c4de3c2101p-11, 0x1.fcc9430345441p-13, + -0x1.fccc917990854p-14, 0x1.f41fbb5026b5p-15 }, + { 0x1.253f3fc844189p-9, -0x1.cadf5cc04da1bp-11, 0x1.73e9dbf6ed988p-12, + -0x1.34f75abb9acfdp-13, 0x1.05502104fc072p-14, -0x1.c015daee9145bp-16, + 0x1.8af9ccda4578cp-17, -0x1.5b973bad98b6bp-18 }, + { -0x1.7f80bfa6d705ep-9, 0x1.e416c7e5d3bb3p-11, -0x1.4361a69711e0fp-12, + 0x1.c0beed7451d56p-14, -0x1.3fc0c220552b8p-15, 0x1.d09a0850c9ad6p-17, + -0x1.5b403dca4645p-18, 0x1.07c1349c4989ap-19 }, + { -0x1.8bd8d36b6b68p-10, 0x1.ab590101636b6p-12, -0x1.e980083cba776p-14, + 0x1.23c1bb53d55c4p-15, -0x1.65be06922ac5p-17, 0x1.bfd7d06279b09p-19, + -0x1.2127a54c7c981p-20, 0x1.79ae1f9c24de8p-22 }, + { -0x1.8db1211cc179cp-11, 0x1.6c24488bdc8e9p-13, -0x1.62882b2ca3c96p-15, + 0x1.67e36a9a5f89cp-17, -0x1.785b8294e4cf2p-19, 0x1.925c40ccc4611p-21, + -0x1.bccb0b78110bcp-23, 0x1.f05d365676624p-25 }, + { -0x1.879379df4c28cp-12, 0x1.2e193030ccfd1p-14, -0x1.f0900c0b3c1fcp-17, + 0x1.aa304a80f3ce4p-19, -0x1.795cdc082b2dfp-21, 0x1.56025546a45a8p-23, + -0x1.4137abecddaa9p-25, 0x1.303d852294977p-27 }, + { -0x1.7b7a8dbd38635p-13, 0x1.eab932219e072p-16, -0x1.528291d0efb42p-18, + 0x1.e86163ded7066p-21, -0x1.6be3b6446506ap-23, 0x1.15d4d64fd204ap-25, + -0x1.b8820763832efp-28, 0x1.5ff781e6e4e19p-30 }, + { -0x1.6b1f37f261621p-14, 0x1.87d88455d6443p-17, -0x1.c3a69901a7d7cp-20, + 0x1.107e1d8456499p-22, -0x1.53f4b8e0d8be8p-25, 0x1.b3016ba9fadffp-28, + -0x1.2175c3eb445d8p-30, 0x1.841e8df7f84e7p-33 }, + { -0x1.57cce7fdc9fe7p-15, 0x1.347c6b65ace16p-18, -0x1.27eb9df26d911p-21, + 0x1.296c0dcf0b476p-24, -0x1.354fd38659786p-27, 0x1.4a2dbe1c4af19p-30, + -0x1.6f1651636db1p-33, 0x1.9b1457d56445ap-36 }, + { -0x1.423d487c99e54p-16, 0x1.df4d1f34022a3p-20, -0x1.7d54e9cd7e7eap-23, + 0x1.3e131c44c6382p-26, -0x1.12afb6bfa8c14p-29, 0x1.e7412bd9ebd87p-33, + -0x1.c2ab005ebc13bp-36, 0x1.a3bbacb6ee6b7p-39 }, + }; + // polynomial coefficients high part + static const double ch[19][13][2] = { + { { 0x1.fdbd7c56b02b5p-2, -0x1.9f8c66985b6f3p-56 }, + { -0x1.c771ed8981f3ep+0, 0x1.8d8b72ce9b19dp-54 }, + { 0x1.1558ba7c0144dp+1, 0x1.4fc1fa0f0451cp-53 }, + { -0x1.1fa938f4d4b53p+1, -0x1.f29beb3ca3738p-53 }, + { 0x1.7f7469f6781efp+1, -0x1.b59ce1aa03545p-53 } }, + { { 0x1.71c14e711391ep-2, 0x1.2ad5eb4fb4f59p-60 }, + { -0x1.740c890bd54d3p+0, -0x1.6978dab8a116p-55 }, + { 0x1.b38de2e957c18p+0, -0x1.aba2b91749902p-55 }, + { -0x1.7ab358c51c087p+0, 0x1.46a8f1bc5883bp-55 }, + { 0x1.af1b63b322b6dp+0, 0x1.2d98d261df8f3p-55 } }, + { { 0x1.e53b12b3407e2p-3, 0x1.97cb2965d31b5p-57 }, + { -0x1.2a144e9a8b92ep+0, -0x1.bbf90d2717ba5p-54 }, + { 0x1.5adc4ef58621ep+0, 0x1.d41b3282f1d5bp-54 }, + { -0x1.fb259e2817239p-1, 0x1.a19b744867ccbp-55 }, + { 0x1.ee43256a6bfd3p-1, 0x1.880c7ca4d6687p-55 } }, + { { 0x1.0719312af823cp-3, -0x1.77ca1d8b99601p-57 }, + { -0x1.d11f75dc5be7dp-1, 0x1.997295e7f58d5p-57 }, + { 0x1.18a58180335ddp+0, -0x1.9e4f675e9e244p-58 }, + { -0x1.5aea0e9166a08p-1, 0x1.4799eb996a78bp-55 }, + { 0x1.22b448094c052p-1, -0x1.221db12561423p-56 } }, + { { 0x1.3c3b637596f8dp-1, -0x1.051b18f5744bap-56 }, + { -0x1.b9ccef0d71197p-1, -0x1.cf98e73bfb3d7p-55 }, + { 0x1.c55517304ef35p-2, 0x1.dfe2299217a1ap-57 }, + { -0x1.4230fb2a20b13p-2, -0x1.8eb1c5690348fp-57 }, + { 0x1.03aa1691c1841p-2, 0x1.a0e14e4b5a96cp-57 } }, + { { -0x1.752403c835a4dp-5, 0x1.a3a43faf6ecccp-59 }, + { -0x1.c0be76051e3a5p-2, -0x1.c737cd3ea73d9p-57 }, + { 0x1.73c36ef7bf402p-1, -0x1.40c4dff8e4c1ep-56 }, + { -0x1.458cec1d1393dp-2, -0x1.c7f148cf356efp-56 }, + { 0x1.8ec2d305516c4p-3, 0x1.9566535c9eabp-57 } }, + { { -0x1.85361b993719fp-4, -0x1.dc41ac35a716fp-58 }, + { -0x1.f3e2bae2cdf7dp-3, 0x1.6d5cae27956a4p-57 }, + { 0x1.3745220b46975p-1, 0x1.56d68f9018bb8p-60 }, + { -0x1.d29172b1a4407p-3, -0x1.93fc4238117bdp-58 }, + { 0x1.ef0f914e4a75bp-4, -0x1.6f0339a5cbb3ap-58 } }, + { { -0x1.ec2ab5aa5843ap-4, -0x1.adc658df2c1c1p-62 }, + { -0x1.a6243a7f3534cp-5, 0x1.0dc0b707b85abp-59 }, + { 0x1.04116f85f23a3p-1, 0x1.517c0b25b9233p-57 }, + { -0x1.4bb33f1abe408p-3, 0x1.cc0c1f637cea4p-58 }, + { 0x1.2ecfafae59f8fp-4, 0x1.3c57c7651ae8ap-58 } }, + { { -0x1.c8928613eb4f5p-4, 0x1.55f36a43c02bcp-62 }, + { 0x1.1151b40dad4e9p-3, 0x1.67907a753aa66p-57 }, + { 0x1.b46b0b78660acp-2, -0x1.9bcdfa3bbcd41p-56 }, + { -0x1.d9cd6009ac89dp-4, -0x1.69c4d18a5c993p-59 }, + { 0x1.73b079d35c37cp-5, -0x1.4d3891ecef09ep-59 } }, + { { -0x1.00ad2093da6e4p-4, -0x1.cbf7cf885033p-58 }, + { 0x1.391f431d39831p-2, 0x1.8fb94bb0e7df5p-56 }, + { 0x1.71d5a6e677f1cp-2, 0x1.d1dc12aaa3806p-59 }, + { -0x1.57f6fbf9108c1p-4, 0x1.4e341fb4cef78p-61 }, + { 0x1.d1e33efae7a1dp-6, 0x1.c4938a6deffbep-60 } }, + { { 0x1.d344dabcc201ep-2, 0x1.574f453e55614p-56 }, + { 0x1.3c3a02b015763p-2, -0x1.342e3d6a27dfap-56 }, + { -0x1.f5d49f62ecfd6p-5, -0x1.07444b43ab601p-60 }, + { 0x1.22abe7bbdf628p-6, 0x1.2cb184651725ap-63 }, + { -0x1.8b52066552f48p-8, 0x1.bc2dbb1b8365dp-62 } }, + { { 0x1.f22e8b160e053p-3, 0x1.89c03c62a66d7p-57 }, + { 0x1.58ae0ae32162p-1, -0x1.594df075ee813p-56 }, + { 0x1.028e87f2859fdp-2, 0x1.bf1ead4dde3d4p-58 }, + { -0x1.55b4949f3971ap-5, -0x1.2cfd594571487p-59 }, + { 0x1.4cfe08a2baa09p-7, 0x1.495ab3aeecafp-62 } }, + { { 0x1.104861734d948p-1, 0x1.32e74856dbad8p-56 }, + { 0x1.b2445e9d82006p-1, 0x1.6e48e474ddfbfp-55 }, + { 0x1.b352d20042182p-3, 0x1.a8ac4f9b7c938p-60 }, + { -0x1.e6c5b3585790ep-6, -0x1.31a8ef26cbf2ep-60 }, + { 0x1.93111b206dab4p-8, -0x1.aa3ae79b1707p-63 } }, + { { 0x1.eed49cf014c0bp-1, 0x1.bca14c01f79aep-55 }, + { 0x1.0718fe597659bp+0, 0x1.7d14012138c17p-55 }, + { 0x1.6c89e19ff8e58p-3, 0x1.12dfe29d6e296p-59 }, + { -0x1.56a9890298c3ap-6, -0x1.2181516eb15d6p-61 }, + { 0x1.deaa0ec93f6d9p-9, 0x1.e4d7a3e816168p-63 } }, + { { 0x1.990530fe5fa37p+0, -0x1.cf639a3a54f76p-56 }, + { 0x1.35e029ece68dp+0, -0x1.e3db2cbb514ebp-60 }, + { 0x1.301f23426a05fp-3, -0x1.b5ec346a456bcp-57 }, + { -0x1.de5b0dd5127b5p-7, 0x1.371374acf777fp-61 }, + { 0x1.1843ded6af0f6p-9, -0x1.7779056d714p-64 } }, + { { 0x1.3ef64cb5ced7bp+1, 0x1.c3c21b0562715p-54 }, + { 0x1.654a3f497c726p+0, 0x1.3331f28ee09bbp-54 }, + { 0x1.f9f5117f295a1p-4, -0x1.ee3d2bb334106p-58 }, + { -0x1.4bb07b47ebf8dp-7, -0x1.64c2c019b90b5p-61 }, + { 0x1.449d9854bac59p-10, -0x1.d0a2827bf227p-64 } }, + { { 0x1.de185c1178ad9p+1, 0x1.d477f1a273bfcp-55 }, + { 0x1.9539397e34b21p+0, 0x1.9743cc0cd10f2p-54 }, + { 0x1.a3e2c09f7886dp-4, -0x1.17f6c25e05338p-59 }, + { -0x1.c98eb5fc97ce2p-8, 0x1.0a5104a9f402dp-63 }, + { 0x1.74b50213890abp-11, 0x1.ff0ae56647adp-65 } }, + { { 0x1.5c2be39a4c6fdp+2, 0x1.ff2814687494cp-52 }, + { 0x1.c59e5d40889c7p+0, 0x1.299ee0827992ap-55 }, + { 0x1.5bbf97b18270ep-4, -0x1.d04ddc6346897p-60 }, + { -0x1.3a2d0322cf70ep-8, 0x1.53fe131154027p-65 }, + { 0x1.a8c6d657c0cfdp-12, -0x1.b402fb82b45efp-66 } }, + { { 0x1.f07834a362b11p+2, -0x1.738a86a953af8p-52 }, + { 0x1.f68034cafc0d3p+0, 0x1.b8d6c9e2cd7d4p-56 }, + { 0x1.1f68e6efd00fap-4, -0x1.6083738e28e87p-61 }, + { -0x1.ad889b8da1552p-9, 0x1.1325e8a48689dp-64 }, + { 0x1.e0ae44f526429p-13, -0x1.997df9412e4aap-67 } }, + }; + + uint64_t t = asuint64 (x); + uint64_t nx = t << 1; + if (__glibc_unlikely (nx >= UINT64_C(0xfeaea9b24f16a34c))) + { + // |x| >= 0x1.006df1bfac84ep+1015 + *signgamp = 1; + if (t == UINT64_C(0x7f5754d9278b51a6)) + return 0x1.ffffffffffffep+1023 - 0x1p+969; + if (t == UINT64_C(0x7f5754d9278b51a7)) + return 0x1.fffffffffffffp+1023 - 0x1p+969; + if (__glibc_unlikely (nx >= (UINT64_C(0x7ff) << 53))) + { /* x=NaN or +/-Inf */ + if (nx == (UINT64_C(0x7ff) << 53)) /* x=+/-Inf */ + return fabs (x); /* +Inf */ + return x + + x; /* x=NaN, where x+x ensures the "Invalid operation" + exception is set if x is sNaN, and it yields a qNaN */ + } + /* The C standard says that if the function overflows, + errno is set to ERANGE. */ + if (t >> 63) + return __math_divzero (0); + else + return 0x1.fp1023 * 0x1.fp1023; // huge positive integer + } + double fx = floor (x); + if (__glibc_unlikely (fx == x)) + { /* x is integer */ + if (x <= 0.0) + { + *signgamp = 1 - 2 * (t >> 63); + return __math_divzero (0); + } + if (x == 1.0 || x == 2.0) + { + *signgamp = 1; + return 0.0; + } + } + unsigned int au = nx >> 38; + double fh, fl, eps; + if (au < ubrd[0]) + { // |x|<0.5 + *signgamp = 1 - 2 * (t >> 63); + double ll, lh = as_logd (fabs (x), &ll); + if (au < 0x1da0000) + { // |x|<0x1p-75 + fh = -lh; + fl = -ll; + eps = 1.5e-22; + } + else if (au < 0x1fd0000) + { // |x|<0.03125 + static const double c0[][2] + = { { -0x1.2788cfc6fb619p-1, 0x1.6cb9a4ff7c53bp-58 }, + { 0x1.a51a6625307d3p-1, 0x1.18722054895e9p-56 }, + { -0x1.9a4d55beab2d7p-2, -0x1.74ded0474fe66p-63 }, + { 0x1.151322ac7d848p-2, 0x1.825b3df1d5722p-56 } }; + static const double q[] + = { -0x1.a8b9c17aa5d3dp-3, 0x1.5b40cb100b9bfp-3, + -0x1.2703a1e13bcbcp-3, 0x1.010b36b6afdc1p-3, + -0x1.c8062dd09ec62p-4, 0x1.9a018c7345316p-4, + -0x1.7578ea8068cc4p-4, 0x1.566b51c990008p-4 }; + double z = x, z2 = z * z, z4 = z2 * z2; + double q0 = q[0] + z * q[1], q2 = q[2] + z * q[3], + q4 = q[4] + z * q[5], q6 = q[6] + z * q[7]; + fl = z * ((q0 + z2 * q2) + z4 * (q4 + z2 * q6)); + fh = polydddfst (z, 4, c0, &fl); + fh = mulddd (x, fh, fl, &fl); + fh = sumdd (-lh, -ll, fh, fl, &fl); + eps = 1.5e-22; + } + else + { + // -0.5> 37; + unsigned int ou = au - ubrd[0]; + int j = ((0x157ced865ul - ou * 0x150d) * ou + 0x128000000000) >> 45; + j -= au < ubrd[j]; + double z = (tf - offs[j]) + xl, z2 = z * z, z4 = z2 * z2; + const double *q = cl[j]; + double q0 = q[0] + z * q[1], q2 = q[2] + z * q[3], + q4 = q[4] + z * q[5], q6 = q[6] + z * q[7]; + fl = z * ((q0 + z2 * q2) + z4 * (q4 + z2 * q6)); + fh = polydddfst (z, 5, ch[j], &fl); + if (__glibc_unlikely (j == 4)) + { // treat the region around the root at 1 + z = -x; + fh = mulddd (z, fh, fl, &fl); + } + eps = fabs (fh) * 8.3e-20; + fh = sumdd (-lh, -ll, fh, fl, &fl); + eps += fabs (lh) * 5e-22; + } + } + else + { + double ax = fabs (x); + if (au >= ubrd[19]) + { // |x|>=8.29541 we use asymptotic expansion or Stirling's formula + double ll, lh = as_logd (ax, &ll); + lh -= 1; + // (x-0.5)*ln(x) = (x-0.5)*(lh + ll) + if (__glibc_unlikely (au >= 0x2198000)) + { // x >= 0x1p52 + // for large |x| use expansion + if (__glibc_unlikely (au >= 0x3fabaa6)) + lh = fasttwosum (lh, ll, &ll); // x>=0x1.754cp+1014 + double hlh = lh * 0.5; + lh = mulddd (ax, lh, ll, &ll); + ll -= hlh; + } + else + { + // for other |x| use a simple product + lh = mulddd (ax - 0.5, lh, ll, &ll); } - switch(i) { - case 0: - z = y*y; - p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10)))); - p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11))))); - p = y*p1+p2; - r += (p-0.5*y); break; - case 1: - z = y*y; - w = z*y; - p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */ - p2 = t1+w*(t4+w*(t7+w*(t10+w*t13))); - p3 = t2+w*(t5+w*(t8+w*(t11+w*t14))); - p = z*p1-(tt-w*(p2+y*p3)); - r += (tf + p); break; - case 2: - p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5))))); - p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5)))); - r += (-0.5*y + p1/p2); + static const double c[][2] + = { { 0x1.acfe390c97d6ap-2, -0x1.1d9792ced423ap-58 }, + { 0x1.55555555554c1p-4, -0x1.0143af34001bdp-59 } }; + static const double q[] + = { -0x1.6c16c1697de08p-9, 0x1.a019f47b230fdp-11, + -0x1.380aab821e42ep-11, 0x1.b617d2c5b5b66p-11, + -0x1.a7fd66a05ccfcp-10 }; + lh = fastsum (lh, ll, c[0][0], c[0][1], &ll); + if (ax < 0x1p100) + { + double zh = 1.0 / ax, zl = fma (zh, -ax, 1.0) * zh; + double z2h = zh * zh, z4h = z2h * z2h; + double q0 = q[0] + z2h * q[1], q2 = q[2] + z2h * q[3], q4 = q[4]; + fl = z2h * (q0 + z4h * (q2 + z4h * q4)); + fh = fasttwosum (c[1][0], fl, &fl); + fl += c[1][1]; + fh = muldd (fh, fl, zh, zl, &fl); } + else + { + fh = 0; + fl = 0; + } + fh = fastsum (lh, ll, fh, fl, &fl); + eps = fabs (fh) * 4.5e-20; } - else if(ix<0x40200000) { /* x < 8.0 */ - i = (int)x; - t = zero; - y = x-(double)i; - p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); - q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); - r = half*y+p/q; - z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ - switch(i) { - case 7: z *= (y+6.0); /* FALLTHRU */ - case 6: z *= (y+5.0); /* FALLTHRU */ - case 5: z *= (y+4.0); /* FALLTHRU */ - case 4: z *= (y+3.0); /* FALLTHRU */ - case 3: z *= (y+2.0); /* FALLTHRU */ - r += __ieee754_log(z); break; + else + { // x in [0.5, 8.29541] range + unsigned ou = au - ubrd[0]; + int j = ((0x157ced865ul - ou * 0x150d) * ou + 0x128000000000) >> 45; + j -= au < ubrd[j]; + double z = ax - offs[j], z2 = z * z, z4 = z2 * z2; + const double *q = cl[j]; + double q0 = q[0] + z * q[1], q2 = q[2] + z * q[3], + q4 = q[4] + z * q[5], q6 = q[6] + z * q[7]; + fl = z * ((q0 + z2 * q2) + z4 * (q4 + z2 * q6)); + fh = polydddfst (z, 5, ch[j], &fl); + if (__glibc_unlikely (j == 4)) + { // treat the region around the root at 1 + z = 1 - ax; + fh = mulddd (z, fh, fl, &fl); + } + if (__glibc_unlikely (j == 10)) + { // treat the region around the root at 2 + z = ax - 2; + fh = mulddd (z, fh, fl, &fl); } - /* 8.0 <= x < 2**58 */ - } else if (ix < 0x43900000) { - t = __ieee754_log(x); - z = one/x; - y = z*z; - w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); - r = (x-half)*(t-one)+w; - } else - /* 2**58 <= x <= inf */ - r = math_narrow_eval (x*(__ieee754_log(x)-one)); - /* NADJ is set for negative arguments but not otherwise, - resulting in warnings that it may be used uninitialized - although in the cases where it is used it has always been - set. */ - DIAG_PUSH_NEEDS_COMMENT; - DIAG_IGNORE_NEEDS_COMMENT_GCC (4.9, "-Wmaybe-uninitialized"); - if(hx<0) r = nadj - r; - DIAG_POP_NEEDS_COMMENT; - return r; + eps = fabs (fh) * 8.3e-20 + 1e-24; + } + if (t >> 63) + { // x<0 so use reflection formula + double sl, sh = as_sinpipid (x - floor (x), &sl); + sh = mulddd (-x, sh, sl, &sl); + double ll, lh = as_logd (sh, &ll); + ll += sl / sh; + fh = -sumdd (fh, fl, lh, ll, &fl); + fl = -fl; + eps += fabs (lh) * 4e-22; + int64_t k = fx; + *signgamp = 1 - 2 * (k & 1); + } + else + { + *signgamp = 1; + } + } + double ub = fh + (fl + eps), lb = fh + (fl - eps); + if (ub != lb) + { // rounding test + return as_lgamma_accurate (x); + } + return ub; } libm_alias_finite (__ieee754_lgamma_r, __lgamma_r) + +double as_logd (double x, double *l) +{ + static const struct + { + uint16_t c0; + short c1; + } B[] = { + { 301, 27565 }, { 7189, 24786 }, { 13383, 22167 }, { 18923, 19696 }, + { 23845, 17361 }, { 28184, 15150 }, { 31969, 13054 }, { 35231, 11064 }, + { 37996, 9173 }, { 40288, 7372 }, { 42129, 5657 }, { 43542, 4020 }, + { 44546, 2457 }, { 45160, 962 }, { 45399, -468 }, { 45281, -1838 }, + { 44821, -3151 }, { 44032, -4412 }, { 42929, -5622 }, { 41522, -6786 }, + { 39825, -7905 }, { 37848, -8982 }, { 35602, -10020 }, { 33097, -11020 }, + { 30341, -11985 }, { 27345, -12916 }, { 24115, -13816 }, { 20661, -14685 }, + { 16989, -15526 }, { 13107, -16339 }, { 9022, -17126 }, { 4740, -17889 } + }; + static const double r1[] = { + 0x1p+0, 0x1.f508p-1, 0x1.ea4ap-1, 0x1.dfcap-1, 0x1.d582p-1, + 0x1.cb72p-1, 0x1.c19ap-1, 0x1.b7f8p-1, 0x1.ae8ap-1, 0x1.a55p-1, + 0x1.9c4ap-1, 0x1.9374p-1, 0x1.8acep-1, 0x1.8258p-1, 0x1.7a12p-1, + 0x1.71f8p-1, 0x1.6a0ap-1, 0x1.6248p-1, 0x1.5abp-1, 0x1.5342p-1, + 0x1.4bfep-1, 0x1.44ep-1, 0x1.3deap-1, 0x1.371ap-1, 0x1.307p-1, + 0x1.29eap-1, 0x1.2388p-1, 0x1.1d48p-1, 0x1.172cp-1, 0x1.113p-1, + 0x1.0b56p-1, 0x1.059cp-1, 0x1p-1, + }; + static const double r2[] + = { 0x1p+0, 0x1.ffa7p-1, 0x1.ff4fp-1, 0x1.fef6p-1, 0x1.fe9ep-1, + 0x1.fe45p-1, 0x1.fdedp-1, 0x1.fd94p-1, 0x1.fd3cp-1, 0x1.fce4p-1, + 0x1.fc8cp-1, 0x1.fc34p-1, 0x1.fbdcp-1, 0x1.fb84p-1, 0x1.fb2cp-1, + 0x1.fad4p-1, 0x1.fa7cp-1, 0x1.fa24p-1, 0x1.f9cdp-1, 0x1.f975p-1, + 0x1.f91ep-1, 0x1.f8c6p-1, 0x1.f86fp-1, 0x1.f817p-1, 0x1.f7cp-1, + 0x1.f769p-1, 0x1.f711p-1, 0x1.f6bap-1, 0x1.f663p-1, 0x1.f60cp-1, + 0x1.f5b5p-1, 0x1.f55ep-1, 0x1.f507p-1 }; + static const double l1[][2] + = { { 0x0p+0, 0x0p+0 }, + { 0x1.9f5e440f128dbp-37, 0x1.62d07abp-6 }, + { -0x1.527d64b444fa3p-37, 0x1.62f483dp-5 }, + { 0x1.3aff57187d0cfp-39, 0x1.0a267214p-4 }, + { -0x1.4634c201e2b9cp-41, 0x1.62e04bcp-4 }, + { -0x1.d46364a8017c7p-36, 0x1.bb9db708p-4 }, + { -0x1.882b6acb3f696p-36, 0x1.0a29f69cp-3 }, + { 0x1.5a5833aeff542p-37, 0x1.368507dap-3 }, + { -0x1.3876d32b0cbf5p-36, 0x1.62e4116cp-3 }, + { 0x1.f5712171380e6p-37, 0x1.8f41d568p-3 }, + { 0x1.fc0b2e87a92c1p-36, 0x1.bb98bc4cp-3 }, + { 0x1.44c7ceb2f93f2p-36, 0x1.e7f71f08p-3 }, + { 0x1.a147c39e44ebap-37, 0x1.0a2bfe2cp-2 }, + { 0x1.36d8fc46707d1p-37, 0x1.205afe03p-2 }, + { -0x1.0fd8155ea585p-37, 0x1.3685b589p-2 }, + { 0x1.8954f1c1b010fp-37, 0x1.4cb42e19p-2 }, + { -0x1.5d0bcd7fa4afap-36, 0x1.62e3e78cp-2 }, + { -0x1.b0a96458bf187p-36, 0x1.79123647p-2 }, + { 0x1.c543eab5348b9p-36, 0x1.8f422996p-2 }, + { -0x1.15143e5c177e1p-37, 0x1.a5711c7ep-2 }, + { 0x1.3be09bf52475cp-38, 0x1.bb9c3cebp-2 }, + { -0x1.9b3b32e71e21dp-40, 0x1.d1cd255bp-2 }, + { -0x1.8f02175f93786p-38, 0x1.e7fb0671p-2 }, + { -0x1.c5fb374b7ddcfp-36, 0x1.fe2980ecp-2 }, + { -0x1.8e174c5571bbdp-36, 0x1.0a2aef35p-1 }, + { 0x1.fa33ff819b3ecp-36, 0x1.15420d49p-1 }, + { 0x1.23d2634096ca6p-38, 0x1.2058ca79p-1 }, + { 0x1.c8afc264146b2p-38, 0x1.2b7156ffp-1 }, + { 0x1.e21780abaa301p-37, 0x1.3686c62p-1 }, + { 0x1.3d67aee28cdc4p-36, 0x1.419f01cdp-1 }, + { 0x1.ccd8a77731be8p-36, 0x1.4cb504d68p-1 }, + { 0x1.0cc7dc4dbbcfdp-37, 0x1.57cb333b8p-1 }, + { 0x1.1cf79abc9e3b4p-36, 0x1.62e42fef8p-1 } }; + static const double l2[][2] = { { 0x0p+0, 0x0p+0 }, + { 0x1.2ccace5b018a7p-36, 0x1.641ef4p-11 }, + { -0x1.88a5cd275513ap-36, 0x1.623d3fp-10 }, + { -0x1.0006a77b80a2dp-38, 0x1.0a4531p-9 }, + { 0x1.81a0ebe451ddp-39, 0x1.627a998p-9 }, + { 0x1.4297627f3b4acp-37, 0x1.bbc015p-9 }, + { 0x1.afb8521676db1p-36, 0x1.0a0a0c8p-8 }, + { 0x1.080b4c8bf43cap-36, 0x1.36bc4ap-8 }, + { 0x1.bbcdc4ef244f5p-37, 0x1.62f5a48p-8 }, + { -0x1.eb4354215c794p-36, 0x1.8f36a44p-8 }, + { -0x1.020c23741371bp-36, 0x1.bb7f4b8p-8 }, + { -0x1.ff1289819f095p-36, 0x1.e7cf9d4p-8 }, + { 0x1.9b1383de1a3f4p-36, 0x1.0a13cdep-7 }, + { 0x1.7fa1788e44213p-38, 0x1.2043a52p-7 }, + { 0x1.f22c04fdaaa38p-37, 0x1.3677558p-7 }, + { -0x1.d9828c736de23p-36, 0x1.4caee08p-7 }, + { -0x1.41b9d7994644p-36, 0x1.62ea474p-7 }, + { 0x1.aced891ec8e07p-36, 0x1.79298b2p-7 }, + { -0x1.e4c5365e893ffp-36, 0x1.8f2be4ep-7 }, + { 0x1.e8d32fe7dc6f5p-36, 0x1.a572dbep-7 }, + { -0x1.3bf707f8ee0b7p-36, 0x1.bb7cd5p-7 }, + { -0x1.c60b95a619b91p-36, 0x1.d1cb84ap-7 }, + { 0x1.64cbc2b83b45cp-37, 0x1.e7dd224p-7 }, + { -0x1.68543f75f32c6p-37, 0x1.fe338fcp-7 }, + { 0x1.421a5be17c2ecp-36, 0x1.0a266bbp-6 }, + { -0x1.f6b329a1da537p-37, 0x1.1534f84p-6 }, + { -0x1.9a217f361c264p-37, 0x1.2065ff9p-6 }, + { 0x1.3856e49eac8ddp-40, 0x1.2b78651p-6 }, + { 0x1.c6c3e72a945a1p-39, 0x1.368cb54p-6 }, + { -0x1.dcfdc96d3a1c6p-36, 0x1.41a2f0dp-6 }, + { 0x1.a0da4813daf37p-38, 0x1.4cbb185p-6 }, + { 0x1.b4d3084ac1ad1p-36, 0x1.57d52c8p-6 } }; + uint64_t t = asuint64 (x); + int ex = t >> 52; + if (__glibc_unlikely (ex == 0)) + { + int k = stdc_leading_zeros (t); + t <<= k - 11; + ex -= k - 12; + } + int e = ex - 0x3ff; + t &= ~UINT64_C(0) >> 12; + double ed = e; + uint64_t i = t >> (52 - 5); + int64_t d = t & (~UINT64_C(0) >> 17); + uint64_t j + = (t + ((uint64_t) B[i].c0 << 33) + ((int64_t) B[i].c1 * (d >> 16))) + >> (52 - 10); + t |= (int64_t) 0x3ff << 52; + int i1 = j >> 5, i2 = j & 0x1f; + double r = r1[i1] * r2[i2]; + double tf = asdouble (t); + double o = r * tf, dxl = fma (r, tf, -o), dxh = o - 1; + static const double c[] = { -0x1.fffffffffffd3p-2, 0x1.55555555543d5p-2, + -0x1.000002bb2d74ep-2, 0x1.999a692c56e4ep-3 }; + double dx = fma (r, tf, -1), dx2 = dx * dx; + double f = dx2 * ((c[0] + dx * c[1]) + dx2 * (c[2] + dx * c[3])); + double lt = (l1[i1][1] + l2[i2][1]) + ed * 0x1.62e42fef8p-1; + double lh = lt + dxh, ll = (lt - lh) + dxh; + ll += ((l1[i1][0] + l2[i2][0]) + 0x1.1cf79abc9e3b4p-36 * ed) + dxl; + ll += f; + *l = ll; + return lh; +} + +double +as_logd_accurate (double x, double *l, double *l_) +{ + static const struct + { + uint16_t c0; + short c1; + } B[] = { + { 301, 27565 }, { 7189, 24786 }, { 13383, 22167 }, { 18923, 19696 }, + { 23845, 17361 }, { 28184, 15150 }, { 31969, 13054 }, { 35231, 11064 }, + { 37996, 9173 }, { 40288, 7372 }, { 42129, 5657 }, { 43542, 4020 }, + { 44546, 2457 }, { 45160, 962 }, { 45399, -468 }, { 45281, -1838 }, + { 44821, -3151 }, { 44032, -4412 }, { 42929, -5622 }, { 41522, -6786 }, + { 39825, -7905 }, { 37848, -8982 }, { 35602, -10020 }, { 33097, -11020 }, + { 30341, -11985 }, { 27345, -12916 }, { 24115, -13816 }, { 20661, -14685 }, + { 16989, -15526 }, { 13107, -16339 }, { 9022, -17126 }, { 4740, -17889 } + }; + static const double r1[] = { + 0x1p+0, 0x1.f508p-1, 0x1.ea4ap-1, 0x1.dfcap-1, 0x1.d582p-1, + 0x1.cb72p-1, 0x1.c19ap-1, 0x1.b7f8p-1, 0x1.ae8ap-1, 0x1.a55p-1, + 0x1.9c4ap-1, 0x1.9374p-1, 0x1.8acep-1, 0x1.8258p-1, 0x1.7a12p-1, + 0x1.71f8p-1, 0x1.6a0ap-1, 0x1.6248p-1, 0x1.5abp-1, 0x1.5342p-1, + 0x1.4bfep-1, 0x1.44ep-1, 0x1.3deap-1, 0x1.371ap-1, 0x1.307p-1, + 0x1.29eap-1, 0x1.2388p-1, 0x1.1d48p-1, 0x1.172cp-1, 0x1.113p-1, + 0x1.0b56p-1, 0x1.059cp-1, 0x1p-1, + }; + static const double r2[] + = { 0x1p+0, 0x1.ffa7p-1, 0x1.ff4fp-1, 0x1.fef6p-1, 0x1.fe9ep-1, + 0x1.fe45p-1, 0x1.fdedp-1, 0x1.fd94p-1, 0x1.fd3cp-1, 0x1.fce4p-1, + 0x1.fc8cp-1, 0x1.fc34p-1, 0x1.fbdcp-1, 0x1.fb84p-1, 0x1.fb2cp-1, + 0x1.fad4p-1, 0x1.fa7cp-1, 0x1.fa24p-1, 0x1.f9cdp-1, 0x1.f975p-1, + 0x1.f91ep-1, 0x1.f8c6p-1, 0x1.f86fp-1, 0x1.f817p-1, 0x1.f7cp-1, + 0x1.f769p-1, 0x1.f711p-1, 0x1.f6bap-1, 0x1.f663p-1, 0x1.f60cp-1, + 0x1.f5b5p-1, 0x1.f55ep-1, 0x1.f507p-1 }; + static const double h1[][3] = { + { 0x0p+0, 0x0p+0, 0x0p+0 }, + { 0x1.7052459b95dfcp-95, 0x1.d79103c4a36a4p-43, 0x1.62d07ab33p-6 }, + { 0x1.6eaee7553baap-99, 0x1.60a6d2eec1732p-43, 0x1.62f483cea8p-5 }, + { 0x1.a79804e7ea455p-95, 0x1.aff57187d0cecp-43, 0x1.0a26721424p-4 }, + { 0x1.3d01ad8f4bdb1p-95, 0x1.ce59eff0ea32p-44, 0x1.62e04bbff4p-4 }, + { 0x1.28aff4e83414bp-97, 0x1.ce4dabff41cbp-43, 0x1.bb9db70628p-4 }, + { 0x1.947079fb56ea3p-96, 0x1.ea4a9a604b4f2p-43, 0x1.0a29f69b3ap-3 }, + { 0x1.071b8f347abdcp-95, 0x1.2c19d77faa108p-44, 0x1.368507da56p-3 }, + { 0x1.60dc9eae5e8c7p-95, 0x1.c4966a79a05a8p-43, 0x1.62e4116b62p-3 }, + { 0x1.19de1930ac389p-95, 0x1.5c485c4e03992p-43, 0x1.8f41d5687cp-3 }, + { 0x1.38f75cf02d618p-95, 0x1.65d0f5258168p-49, 0x1.bb98bc4cfep-3 }, + { 0x1.66d0ece8cc055p-99, 0x1.8f9d65f27e3cp-45, 0x1.e7f71f08a2p-3 }, + { 0x1.110c0a2d0d02ap-97, 0x1.47c39e44eb9ap-45, 0x1.0a2bfe2c34p-2 }, + { 0x1.c2e1c7624937cp-95, 0x1.b63f119c1f42ap-43, 0x1.205afe0326p-2 }, + { 0x1.80a68afc2eecp-98, 0x1.3f550ad3d7d8p-48, 0x1.3685b588dep-2 }, + { 0x1.2cb8203dc0bb6p-95, 0x1.54f1c1b010eb8p-45, 0x1.4cb42e1931p-2 }, + { 0x1.e234ee6f272d2p-95, 0x1.7a19402da833cp-43, 0x1.62e3e78ba8p-2 }, + { 0x1.3d19c1b2c67cbp-96, 0x1.ab4dd3a073c68p-43, 0x1.7912364693p-2 }, + { 0x1.62e1c0308ab55p-95, 0x1.43eab5348b948p-44, 0x1.8f42299671p-2 }, + { 0x1.9b5e4e30c5f0fp-95, 0x1.75e0d1f440f44p-44, 0x1.a5711c7dddp-2 }, + { 0x1.32ba52f147867p-99, 0x1.7c137ea48eb8p-43, 0x1.bb9c3ceb13p-2 }, + { 0x1.32f3b1a0f392dp-95, 0x1.262668c70ef16p-43, 0x1.d1cd255af9p-2 }, + { 0x1.078661f35a007p-96, 0x1.fbd140d90f48p-47, 0x1.e7fb0670e7p-2 }, + { 0x1.0a4b6082eb02fp-96, 0x1.02645a4111858p-43, 0x1.fe2980eb8ep-2 }, + { 0x1.8f3673c008d9dp-95, 0x1.e8b3aa8e4433cp-44, 0x1.0a2aef34cep-1 }, + { 0x1.e4109a84a717ep-96, 0x1.19ffc0cd9f62p-43, 0x1.15420d493fp-1 }, + { 0x1.b400798743957p-96, 0x1.e931a04b652cp-45, 0x1.2058ca7909p-1 }, + { 0x1.3046edc6794ap-96, 0x1.15f84c828d642p-43, 0x1.2b7156ff0ep-1 }, + { 0x1.5cb5072b4b776p-95, 0x1.0bc055d518094p-44, 0x1.3686c6201ep-1 }, + { 0x1.051804915de87p-97, 0x1.67aee28cdc41cp-44, 0x1.419f01cd278p-1 }, + { 0x1.748ebb6888301p-98, 0x1.b14eee637d048p-45, 0x1.4cb504d6b98p-1 }, + { 0x1.55098cdd851edp-97, 0x1.31f7136ef3f3p-43, 0x1.57cb333b908p-1 }, + { 0x1.f97b57a079a19p-103, 0x1.ef35793c7673p-45, 0x1.62e42fefa38p-1 }, + }; + static const double h2[][3] = { + { 0x0p+0, 0x0p+0, 0x0p+0 }, + { 0x1.4826401258afcp-95, 0x1.959cb60314ep-45, 0x1.641ef496p-11 }, + { 0x1.9c71d4cdcb54bp-95, 0x1.ad196c55762e2p-43, 0x1.623d3e9dp-10 }, + { 0x1.012c2a9a87c83p-96, 0x1.ff2b108feba5p-43, 0x1.0a4530f78p-9 }, + { 0x1.437bb5bf9fc34p-95, 0x1.a0ebe451ddp-47, 0x1.627a9986p-9 }, + { 0x1.1c54e00c02c66p-96, 0x1.4bb13f9da5614p-44, 0x1.bbc01514p-9 }, + { 0x1.b959e749a5ea5p-95, 0x1.dc290b3b6d866p-43, 0x1.0a0a0c9acp-8 }, + { 0x1.9d74597e3c084p-95, 0x1.69917e8794bp-49, 0x1.36bc4a108p-8 }, + { 0x1.ee2ae1f36e35ap-98, 0x1.e6e2779227a9cp-44, 0x1.62f5a48dcp-8 }, + { 0x1.6236c8dcf0cfdp-95, 0x1.7957bd470d7fp-45, 0x1.8f36a4214p-8 }, + { 0x1.a291a656020e5p-95, 0x1.f3dc8bec8e52p-44, 0x1.bb7f4b6fcp-8 }, + { 0x1.a34920afb6399p-95, 0x1.daecfcc1ed628p-45, 0x1.e7cf9d2p-8 }, + { 0x1.7b470c6ab05c6p-95, 0x1.89c1ef0d1fa32p-43, 0x1.0a13cdeccp-7 }, + { 0x1.2c071f819920fp-95, 0x1.f42f11c88426ep-43, 0x1.2043a522ep-7 }, + { 0x1.3d45aed9f22ap-95, 0x1.16027ed551c38p-44, 0x1.36775587cp-7 }, + { 0x1.c4ec608241822p-95, 0x1.3eb9c6490ee9ep-43, 0x1.4caee0712p-7 }, + { 0x1.5c705f1ca65efp-95, 0x1.2314335cde02p-43, 0x1.62ea4735ep-7 }, + { 0x1.c6e86fe6c904cp-96, 0x1.db123d91c0e6p-45, 0x1.79298b2d6p-7 }, + { 0x1.cf83c85a09c0dp-95, 0x1.9d64d0bb600a6p-43, 0x1.8f2be4d0cp-7 }, + { 0x1.8a899399986f8p-96, 0x1.a65fcfb8de998p-45, 0x1.a572dbef4p-7 }, + { 0x1.75dd2d25f587cp-101, 0x1.1f00e23e91p-49, 0x1.bb7cd4f62p-7 }, + { 0x1.66dc5826a7f9dp-95, 0x1.f46a59e646ed8p-44, 0x1.d1cb8491cp-7 }, + { 0x1.0bdd376e8742ap-95, 0x1.32f0ae0ed1704p-43, 0x1.e7dd22458p-7 }, + { 0x1.3bcc642260db7p-95, 0x1.eaf0228334e7ep-43, 0x1.fe338fba4p-7 }, + { 0x1.fe76520a42aecp-95, 0x1.0d2df0be17626p-43, 0x1.0a266bb5p-6 }, + { 0x1.3f0eff67b7e5p-97, 0x1.4cd65e25ac96p-45, 0x1.1534f83c1p-6 }, + { 0x1.824751a1e6c5bp-99, 0x1.77a03278f66f6p-43, 0x1.2065ff8ccp-6 }, + { 0x1.506cfc4639accp-95, 0x1.c2b724f5646e6p-43, 0x1.2b7865104p-6 }, + { 0x1.a7ac1fdfa37ep-99, 0x1.b0f9caa516828p-45, 0x1.368cb540ep-6 }, + { 0x1.b61c029607481p-96, 0x1.811b4962f1cf2p-43, 0x1.41a2f0c88p-6 }, + { 0x1.2879f319383adp-97, 0x1.b49027b5e6dep-47, 0x1.4cbb1851ap-6 }, + { 0x1.f27f4e668c317p-97, 0x1.a61095835a298p-45, 0x1.57d52c86dp-6 }, + { 0x1.c301f232c0e74p-96, 0x1.b6fc11defa4a8p-43, 0x1.62f12e132p-6 }, + }; + + uint64_t t = asuint64 (x); + int ex = t >> 52; + if (__glibc_unlikely (ex == 0)) + { + int k = stdc_leading_zeros (t); + t <<= k - 11; + ex -= k - 12; + } + int e = ex - 0x3ff; + t &= ~UINT64_C(0) >> 12; + double ed = e; + uint64_t i = t >> (52 - 5); + int64_t d = t & (~UINT64_C(0) >> 17); + uint64_t j + = (t + ((uint64_t) B[i].c0 << 33) + ((int64_t) B[i].c1 * (d >> 16))) + >> (52 - 10); + t |= (int64_t) 0x3ff << 52; + int i1 = j >> 5, i2 = j & 0x1f; + double r = r1[i1] * r2[i2]; + double tf = asdouble (t); + double o = r * tf, dxl = fma (r, tf, -o), dxh = o - 1; + static const double c[][2] + = { { 0x1p+0, 0x1.a193d7f59d80ap-118 }, + { -0x1p-1, 0x1.8c7d7a8733406p-99 }, + { 0x1.5555555555555p-2, 0x1.55555554f571dp-56 }, + { -0x1p-2, -0x1.e516b5d7b8c15p-73 }, + { 0x1.999999999999ap-3, -0x1.97f2898534175p-57 }, + { -0x1.55555555554b5p-3, -0x1.3834d62d64ec4p-59 }, + { 0x1.249249248dbdcp-3, -0x1.8aa032979ebedp-58 }, + { -0x1.000004e71581bp-3, 0x1.2e7f17d1c0e63p-57 }, + { 0x1.c71e5ec7051f6p-4, 0x1.217ec3dcb2f03p-58 } }; + dxh = fasttwosum (dxh, dxl, &dxl); + double fl = dxh * (c[6][0] + dxh * (c[7][0] + dxh * (c[8][0]))), + fh = polydd (dxh, dxl, 6, c, &fl); + fh = muldd (dxh, dxl, fh, fl, &fl); + double s2 = h1[i1][2] + h2[i2][2], s1 = h1[i1][1] + h2[i2][1], + s0 = h1[i1][0] + h2[i2][0]; + double L0 = 0x1.62e42fefa38p-1 * ed, L1 = 0x1.ef35793c76p-45 * ed, + L2 = 0x1.cc01f97b57a08p-87 * ed; + L0 += s2; + L1 = sumdd (L1, L2, s1, s0, &L2); + L1 = sumdd (L1, L2, fh, fl, &L2); + + L0 = fasttwosum (L0, L1, &L1); + L1 = fasttwosum (L1, L2, &L2); + + *l = L1; + *l_ = L2; + + return L0; +} + +static const double stpi[][2] + = { { 0x0p+0, 0x0p+0 }, + { 0x1.c14eff99a3ff1p-64, 0x1.fff2d746c8895p-8 }, + { -0x1.8c4d4c1bbe38bp-62, 0x1.ffcb5e52d1f36p-7 }, + { -0x1.08ef2408930ebp-61, 0x1.7fa7329846febp-6 }, + { -0x1.14daa07929354p-60, 0x1.ff2d8cc5320c7p-6 }, + { 0x1.d845cf264d016p-60, 0x1.3f3289bb44643p-5 }, + { -0x1.43aa63f69aceap-60, 0x1.7e9d144d37f33p-5 }, + { -0x1.bc90382ed68a4p-59, 0x1.bdcc9ea69fc93p-5 }, + { 0x1.0fbc215a3c756p-60, 0x1.fcb76a6ecccabp-5 }, + { 0x1.72b75e84ab5e2p-58, 0x1.1da9e1f36c497p-4 }, + { -0x1.20d100fccf991p-59, 0x1.3ccc01b453709p-4 }, + { -0x1.f7aac846eccfdp-63, 0x1.5bbd477204bep-4 }, + { -0x1.17799578a6651p-59, 0x1.7a78edace5e27p-4 }, + { 0x1.0c85deb5bb812p-58, 0x1.98fa372a35c37p-4 }, + { -0x1.67d2eb81bbf36p-60, 0x1.b73c6faf2275cp-4 }, + { -0x1.14b2141507a9dp-63, 0x1.d53aecba7bfp-4 }, + { -0x1.8939cffeb036cp-58, 0x1.f2f10e3ce6d42p-4 }, + { -0x1.2f3fbb178d1c5p-57, 0x1.082d1fa7b9738p-3 }, + { 0x1.08479c62d3d77p-57, 0x1.16b8fb743c879p-3 }, + { -0x1.894149dc3b5f7p-57, 0x1.2519dc47527b3p-3 }, + { -0x1.44dad213ab344p-60, 0x1.334d8a850758dp-3 }, + { -0x1.2d415416bae28p-58, 0x1.4151d589a490fp-3 }, + { -0x1.2e0d0b51ed237p-57, 0x1.4f24940025067p-3 }, + { 0x1.e8045a3cf3213p-57, 0x1.5cc3a43788a3p-3 }, + { 0x1.be4e50e1bf91fp-57, 0x1.6a2cec76fa4bp-3 }, + { 0x1.b1e18c1f7f635p-62, 0x1.775e5b50bb365p-3 }, + { -0x1.10946c1f6f484p-63, 0x1.8455e7f3c6e5ap-3 }, + { 0x1.291a88889a4e6p-59, 0x1.9111927c231cfp-3 }, + { -0x1.bedd6f9a25da4p-57, 0x1.9d8f6441cf80bp-3 }, + { -0x1.ce108006670c7p-57, 0x1.a9cd702648a97p-3 }, + { 0x1.4c65624119572p-61, 0x1.b5c9d2e092baap-3 }, + { -0x1.a26e2a2682111p-57, 0x1.c182b347bfc21p-3 }, + { 0x1.fce159c2bb59bp-59, 0x1.ccf6429be6621p-3 }, + { -0x1.59b1cffa69603p-58, 0x1.d822bccd7d86ep-3 }, + { 0x1.677083288397ap-57, 0x1.e30668c31224ep-3 }, + { 0x1.9a49696faa0ecp-57, 0x1.ed9f989d4c415p-3 }, + { 0x1.ca323e77a3345p-58, 0x1.f7eca9f938c6fp-3 }, + { -0x1.c702625d3863bp-57, 0x1.00f6031866f76p-2 }, + { 0x1.180cf0e52237dp-56, 0x1.05ce114cd024ap-2 }, + { -0x1.4be56fec860b9p-56, 0x1.0a7dc060df5eep-2 }, + { 0x1.b5d970e5d9d07p-58, 0x1.0f045755560d9p-2 }, + { 0x1.4c32e06c67499p-58, 0x1.1361238136929p-2 }, + { -0x1.b512d49aedaa1p-56, 0x1.179378ad51274p-2 }, + { -0x1.161478130996dp-58, 0x1.1b9ab12ed2518p-2 }, + { -0x1.25feb091e921fp-59, 0x1.1f762e00ced83p-2 }, + { 0x1.3750bc95dae67p-56, 0x1.232556dcc945fp-2 }, + { -0x1.257966a1044c5p-56, 0x1.26a79a522d332p-2 }, + { -0x1.ac6af78c05e44p-57, 0x1.29fc6ddcbcb72p-2 }, + { 0x1.71dbd64ba4f95p-56, 0x1.2d234df9ec8c9p-2 }, + { 0x1.020107d2c17bp-57, 0x1.301bbe3d2b9c7p-2 }, + { -0x1.9d4016f0b15c4p-56, 0x1.32e5496312cfcp-2 }, + { 0x1.f557b51b587ccp-56, 0x1.357f81637a329p-2 }, + { 0x1.c88cee9bad9f9p-57, 0x1.37e9ff82709ecp-2 }, + { 0x1.ebde6bb284e87p-56, 0x1.3a246460134f7p-2 }, + { -0x1.8b1d8c40ffea3p-56, 0x1.3c2e580742edap-2 }, + { 0x1.de48797b477f2p-56, 0x1.3e0789fb33cf7p-2 }, + { -0x1.8bc6105a80fa5p-56, 0x1.3fafb143d754bp-2 }, + { -0x1.9f4cc680744f3p-56, 0x1.41268c791c743p-2 }, + { 0x1.2ed295e9d0ef2p-60, 0x1.426be1cd05c06p-2 }, + { -0x1.4a98b72ed3789p-60, 0x1.437f7f1493531p-2 }, + { -0x1.5c080cdd72ddfp-56, 0x1.446139cf7f413p-2 }, + { -0x1.b00c622ae015ep-57, 0x1.4510ef2ecb654p-2 }, + { 0x1.8dd5ec4960646p-56, 0x1.458e841a1f7dap-2 }, + { 0x1.e1f89d1adcbc6p-56, 0x1.45d9e533f6cacp-2 }, + { -0x1.6b01ec5417056p-56, 0x1.45f306dc9c883p-2 } }; + +double +as_sinpipid (double x, double *l) +{ + x -= 0.5; + double ax = fabs (x); + double sx = ax * 128; + double ix = roundeven_finite (sx); + int ky = ix, kx = 64 - ky; + if (__glibc_unlikely (kx < 2)) + { + static const double c[2] + = { -0x1.a51a6625307d3p+0, -0x1.16cc8f2044a4ap-55 }; + static const double cl[] = { 0x1.9f9cb402bc42ap-1, -0x1.86a8e46ddf78dp-3, + 0x1.ac644e7aa33e6p-6 }; + double z = 0.5 - ax, z2 = z * z, z2l = fma (z, z, -z2); + double fl = z2 * (cl[0] + z2 * (cl[1] + z2 * (cl[2]))), + fh = fasttwosum (c[0], fl, &fl), e; + fl += c[1]; + fh = muldd (z2, z2l, fh, fl, &fl); + fh = mulddd (z, fh, fl, &fl); + fh = fasttwosum (z, fh, &e); + fl += e; + *l = fl; + return fh; + } + double d = ix - sx, d2 = d * d; + + double sh = stpi[kx][1], sl = stpi[kx][0]; + double ch = stpi[ky][1], cl = stpi[ky][0]; + // sin(a + d) = sin(a)*(1-d^2*P(d^2)) + cos(a)*d*Q(d^2) + // sin(a_i + d) = s[i] + d*(c[i]*Q(d^2) - s[i]*d*P(d^2)) + // sin(a_i + d)/pi = s[i]/pi + d*(c[i]/pi*Q(d^2) - s[i]/pi*d*P(d^2)) + + static const double c[] = { -0x1.3bd3cc9be45dep-12, 0x1.03c1f081b5ac4p-26, + -0x1.55d3c7e3bd8bfp-42, 0x1.e1f4826790653p-59 }; + double c0 = -0x1.692b66e3cf6e8p-66; + static const double s[] = { 0x1.921fb54442d18p-6, -0x1.4abbce625be53p-19, + 0x1.466bc67748efcp-34, -0x1.32d26e446373ap-50 }; + double s0 = 0x1.1a624b88c9448p-60; + + double P = d2 * (c[1] + d2 * (c[2] + d2 * c[3])); + double Q = d2 * (s[1] + d2 * (s[2] + d2 * s[3])); + + double ql, qh = fasttwosum (s[0], Q, &ql); + ql += s0; + ch = muldd (qh, ql, ch, cl, &cl); + double tl, th = fasttwosum (c[0], P, &tl); + tl += c0; + th = mulddd (d, th, tl, &tl); + double pl, ph = muldd (th, tl, sh, sl, &pl); + ch = fastsum (ch, cl, ph, pl, &cl); + ch = mulddd (d, ch, cl, &cl); + sh = fastsum (sh, sl, ch, cl, l); + return sh; +} + +double +as_sinpipid_accurate (double x, double *l) +{ + x -= 0.5; + x = fabs (x); + x *= 128; + double ix = roundeven_finite (x), d = ix - x; + int ky = ix, kx = 64 - ky; + + double sh = stpi[kx][1], sl = stpi[kx][0]; + double ch = stpi[ky][1], cl = stpi[ky][0]; + // sin(a + d) = sin(a)*(1-d^2*P(d^2)) + cos(a)*d*Q(d^2) + // sin(a_i + d) = s[i] + d*(c[i]*Q(d^2) - s[i]*d*P(d^2)) + + static const double c[][2] + = { { -0x1.3bd3cc9be45dep-12, -0x1.692b71366c792p-66 }, + { 0x1.03c1f081b5ac4p-26, -0x1.32b342c5da62ep-80 }, + { -0x1.55d3c7e3cbffap-42, 0x1.17fcb68af9aaep-98 }, + { 0x1.e1f506890e556p-59, -0x1.56c9658a7cfd5p-114 }, + { -0x1.a6d193ca649ap-76, 0x1.61669b56f0275p-134 } }; + static const double s[][2] + = { { 0x1.921fb54442d18p-6, 0x1.1a62633145c07p-60 }, + { -0x1.4abbce625be53p-19, 0x1.05511c684796p-73 }, + { 0x1.466bc6775aae2p-34, -0x1.6dc0d14c26b21p-89 }, + { -0x1.32d2cce62bd86p-50, 0x1.2ba6dbc4b37cp-104 }, + { 0x1.50783487e6b5cp-67, -0x1.b77efed0f9c1fp-122 }, + { -0x1.e306ec8cf7c02p-85, 0x1.5d2601f85289ap-139 } }; + + double d2h = d * d, d2l = fma (d, d, -d2h); + double Pl = 0, Ph = polydd (d2h, d2l, 5, c, &Pl); + double Ql = 0, Qh = polydd (d2h, d2l, 6, s, &Ql); + + Ph = mulddd (d, Ph, Pl, &Pl); + Ph = muldd (sh, sl, Ph, Pl, &Pl); + Qh = muldd (ch, cl, Qh, Ql, &Ql); + + ch = fastsum (Qh, Ql, Ph, Pl, &cl); + ch = mulddd (d, ch, cl, &cl); + sh = fastsum (sh, sl, ch, cl, l); + return sh; +} + +double +as_lgamma_asym_accurate (double xh, double *xl, double *e) +{ + // x*ln(x) - x - 0.5*ln(x) + 0.5*ln(2*pi) + // (x-0.5)*ln(x) - (x - 0.5) - 0.5 + 0.5*ln(2*pi) + // (x-0.5)*(ln(x)-1) + 0.5*(ln(2*pi)-1) + + double l2, l1, l0 = as_logd_accurate (xh, &l1, &l2), l0x, l1x, l2x; + if (xh < 0x1p120) + { + double zh = 1.0 / xh, dz = *xl * zh, + zl = (fma (zh, -xh, 1.0) - dz) * zh; + if (*xl != 0) + { + double dl2, dl1 = mulddd (*xl, zh, fma (zh, -xh, 1.0) * zh, &dl2); + dl2 -= dl1 * dl1 / 2; + l1 = sumdd (l1, l2, dl1, dl2, &l2); + } + + double wl, wh; + if (asuint64 (xh) >> 52 > 0x3ff + 51) + { + wh = xh; + wl = *xl - 0.5; + } + else + { + wh = xh - 0.5; + wl = *xl; + } + l0 -= 1; + + l0x = l0 * wh; + double l0xl = fma (l0, wh, -l0x); + l1x = l1 * wh; + double l1xl = fma (l1, wh, -l1x); + l2x = l2 * wh; + + l1x = sumdd (l1x, l2x, l0xl, l1xl, &l2x); + l1x = sumdd (l1x, l2x, l0 * wl, l1 * wl, &l2x); + double z2l, z2h = muldd (zh, zl, zh, zl, &z2l); + double fh, fl; + if (xh >= 48) + { + static const double c[][2] + = { { 0x1.acfe390c97d69p-2, 0x1.3494bc9001766p-56 }, + { 0x1.5555555555555p-4, 0x1.55555554133c7p-58 }, + { -0x1.6c16c16c16c17p-9, 0x1.f4c03199d8517p-64 }, + { 0x1.a01a01a01a016p-11, 0x1.fb315e77b4883p-66 }, + { -0x1.381381380c1bp-11, -0x1.c4cc418316ed1p-65 }, + { 0x1.b951e22dd8dfcp-11, 0x1.e8da392824ecfp-65 }, + { -0x1.f6a875bb1ab7bp-10, 0x1.27c5fcbab6b5dp-64 }, + { 0x1.a0a6926f4992p-8, -0x1.1f355cbf82229p-63 } }; + l1x = sumdd (l1x, l2x, c[0][0], c[0][1], &l2x); + fl = 0; + fh = polydd (z2h, z2l, 7, c + 1, &fl); + } + else if (xh >= 14.5) + { + static const double c[][2] + = { { 0x1.acfe390c97d69p-2, 0x1.3494bc9007f28p-56 }, + { 0x1.5555555555555p-4, 0x1.5555554ad7655p-58 }, + { -0x1.6c16c16c16c17p-9, 0x1.f4b73ea546bd4p-64 }, + { 0x1.a01a01a01a01ap-11, -0x1.464e31b1a609ap-65 }, + { -0x1.3813813813692p-11, 0x1.63676fd3c851ep-66 }, + { 0x1.b951e2b143b5ep-11, 0x1.cfe48021143d1p-65 }, + { -0x1.f6ab0d459cadbp-10, 0x1.770df5ee0beeap-64 }, + { 0x1.a41a211f4e098p-8, -0x1.e00b3c1619519p-63 }, + { -0x1.e41f97ad4634dp-6, 0x1.05ffdbc72560fp-60 }, + { 0x1.6f15ef2b47719p-3, -0x1.5665cba69dbaap-57 }, + { -0x1.5762c9f49fe25p+0, 0x1.afeff5294ad13p-54 }, + { 0x1.2ea102098f818p+3, 0x1.609db97f1bc89p-51 } }; + l1x = sumdd (l1x, l2x, c[0][0], c[0][1], &l2x); + fl = 0; + fh = polydd (z2h, z2l, 11, c + 1, &fl); + } + else + { + static const double c[][2] + = { { 0x1.acfe390c97d69p-2, 0x1.3494bce9b5c5p-56 }, + { 0x1.5555555555555p-4, 0x1.55551a0d18a1dp-58 }, + { -0x1.6c16c16c16c17p-9, 0x1.07171d4a61bb9p-63 }, + { 0x1.a01a01a01a008p-11, -0x1.49af39145d6ep-65 }, + { -0x1.38138138124ccp-11, -0x1.b8f71068f5292p-66 }, + { 0x1.b951e2b09b07p-11, -0x1.3c8a4c099ae72p-66 }, + { -0x1.f6ab0d4de4a5cp-10, -0x1.fbb4ed542d17ep-65 }, + { 0x1.a41a384fafd09p-8, 0x1.1510cc5dec148p-64 }, + { -0x1.e4277d1a2b2e4p-6, 0x1.33cdd66b223e7p-61 }, + { 0x1.6fdf731005805p-3, -0x1.a5a306692214p-57 }, + { -0x1.641de6a9b2f1cp+0, 0x1.27cfac68728a3p-56 }, + { 0x1.aa463d5553a9fp+3, 0x1.aeba9a5b525f3p-52 }, + { -0x1.312d3aa56b6a2p+7, -0x1.2b2c4e643e5b9p-48 }, + { 0x1.f3e8a4cd3d268p+10, -0x1.ec8b086b0215ep-45 }, + { -0x1.ba9362f228307p+14, 0x1.b8ecae7dba56fp-40 }, + { 0x1.8ed4df00421e4p+18, 0x1.c620cd11f44a5p-36 }, + { -0x1.5af889993596dp+22, -0x1.23a78dbbfe515p-32 }, + { 0x1.1789aff6ddc93p+26, 0x1.fe7e161d8bdb1p-29 }, + { -0x1.94353def0e5f8p+29, -0x1.64a798491c211p-25 }, + { 0x1.ffb9253861b9cp+32, -0x1.469d6acfea485p-23 }, + { -0x1.159c1244c9ee4p+36, -0x1.f7491d4abcc7bp-19 }, + { 0x1.f991ba9776c6ap+38, 0x1.cf2c83cfe2e55p-19 }, + { -0x1.795f45ae3fa72p+41, -0x1.a777148d59de4p-13 }, + { 0x1.c05845a910b48p+43, 0x1.433e08525d4bdp-11 }, + { -0x1.96c399c145d71p+45, 0x1.d495c9370dc8ap-11 }, + { 0x1.083bb9c79529bp+47, -0x1.1640e2de8cf9cp-8 }, + { -0x1.b53232e707f4bp+47, 0x1.6d239acd71b5cp-8 }, + { 0x1.59bad61bd81d5p+47, 0x1.0e3f2ea42a0ep-7 } }; + l1x = sumdd (l1x, l2x, c[0][0], c[0][1], &l2x); + fl = 0; + fh = polydd (z2h, z2l, array_length (c) - 1, c + 1, &fl); + } + fh = muldd (zh, zl, fh, fl, &fl); + l1x = sumdd (l1x, l2x, fh, fl, &l2x); + l0x = fasttwosum (l0x, l1x, &l1x); + l1x = fasttwosum (l1x, l2x, &l2x); + } + else + { + double wl = *xl - 0.5; + l0 -= 1; + l0x = l0 * xh; + double l0xl = fma (l0, xh, -l0x); + l1x = l1 * xh; + double l1xl = fma (l1, xh, -l1x); + l2x = l2 * xh; + l1x = sumdd (l1x, l2x, l0xl, l1xl, &l2x); + l1x = sumdd (l1x, l2x, l0 * wl, l1 * wl, &l2x); + } + *xl = l1x; + *e = l2x; + return l0x; +} diff --git a/sysdeps/ieee754/dbl-64/lgamma_neg.c b/sysdeps/ieee754/dbl-64/lgamma_neg.c deleted file mode 100644 index ec76005e7a..0000000000 --- a/sysdeps/ieee754/dbl-64/lgamma_neg.c +++ /dev/null @@ -1,385 +0,0 @@ -/* lgamma expanding around zeros. - Copyright (C) 2015-2025 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 - . */ - -#include -#include -#include -#include -#include - -static const double lgamma_zeros[][2] = - { - { -0x2.74ff92c01f0d8p+0, -0x2.abec9f315f1ap-56 }, - { -0x2.bf6821437b202p+0, 0x6.866a5b4b9be14p-56 }, - { -0x3.24c1b793cb35ep+0, -0xf.b8be699ad3d98p-56 }, - { -0x3.f48e2a8f85fcap+0, -0x1.70d4561291237p-56 }, - { -0x4.0a139e1665604p+0, 0xf.3c60f4f21e7fp-56 }, - { -0x4.fdd5de9bbabf4p+0, 0xa.ef2f55bf89678p-56 }, - { -0x5.021a95fc2db64p+0, -0x3.2a4c56e595394p-56 }, - { -0x5.ffa4bd647d034p+0, -0x1.7dd4ed62cbd32p-52 }, - { -0x6.005ac9625f234p+0, 0x4.9f83d2692e9c8p-56 }, - { -0x6.fff2fddae1bcp+0, 0xc.29d949a3dc03p-60 }, - { -0x7.000cff7b7f87cp+0, 0x1.20bb7d2324678p-52 }, - { -0x7.fffe5fe05673cp+0, -0x3.ca9e82b522b0cp-56 }, - { -0x8.0001a01459fc8p+0, -0x1.f60cb3cec1cedp-52 }, - { -0x8.ffffd1c425e8p+0, -0xf.fc864e9574928p-56 }, - { -0x9.00002e3bb47d8p+0, -0x6.d6d843fedc35p-56 }, - { -0x9.fffffb606bep+0, 0x2.32f9d51885afap-52 }, - { -0xa.0000049f93bb8p+0, -0x1.927b45d95e154p-52 }, - { -0xa.ffffff9466eap+0, 0xe.4c92532d5243p-56 }, - { -0xb.0000006b9915p+0, -0x3.15d965a6ffea4p-52 }, - { -0xb.fffffff708938p+0, -0x7.387de41acc3d4p-56 }, - { -0xc.00000008f76c8p+0, 0x8.cea983f0fdafp-56 }, - { -0xc.ffffffff4f6ep+0, 0x3.09e80685a0038p-52 }, - { -0xd.00000000b092p+0, -0x3.09c06683dd1bap-52 }, - { -0xd.fffffffff3638p+0, 0x3.a5461e7b5c1f6p-52 }, - { -0xe.000000000c9c8p+0, -0x3.a545e94e75ec6p-52 }, - { -0xe.ffffffffff29p+0, 0x3.f9f399fb10cfcp-52 }, - { -0xf.0000000000d7p+0, -0x3.f9f399bd0e42p-52 }, - { -0xf.fffffffffff28p+0, -0xc.060c6621f513p-56 }, - { -0x1.000000000000dp+4, -0x7.3f9f399da1424p-52 }, - { -0x1.0ffffffffffffp+4, -0x3.569c47e7a93e2p-52 }, - { -0x1.1000000000001p+4, 0x3.569c47e7a9778p-52 }, - { -0x1.2p+4, 0xb.413c31dcbecdp-56 }, - { -0x1.2p+4, -0xb.413c31dcbeca8p-56 }, - { -0x1.3p+4, 0x9.7a4da340a0ab8p-60 }, - { -0x1.3p+4, -0x9.7a4da340a0ab8p-60 }, - { -0x1.4p+4, 0x7.950ae90080894p-64 }, - { -0x1.4p+4, -0x7.950ae90080894p-64 }, - { -0x1.5p+4, 0x5.c6e3bdb73d5c8p-68 }, - { -0x1.5p+4, -0x5.c6e3bdb73d5c8p-68 }, - { -0x1.6p+4, 0x4.338e5b6dfe14cp-72 }, - { -0x1.6p+4, -0x4.338e5b6dfe14cp-72 }, - { -0x1.7p+4, 0x2.ec368262c7034p-76 }, - { -0x1.7p+4, -0x2.ec368262c7034p-76 }, - { -0x1.8p+4, 0x1.f2cf01972f578p-80 }, - { -0x1.8p+4, -0x1.f2cf01972f578p-80 }, - { -0x1.9p+4, 0x1.3f3ccdd165fa9p-84 }, - { -0x1.9p+4, -0x1.3f3ccdd165fa9p-84 }, - { -0x1.ap+4, 0xc.4742fe35272dp-92 }, - { -0x1.ap+4, -0xc.4742fe35272dp-92 }, - { -0x1.bp+4, 0x7.46ac70b733a8cp-96 }, - { -0x1.bp+4, -0x7.46ac70b733a8cp-96 }, - { -0x1.cp+4, 0x4.2862898d42174p-100 }, - }; - -static const double e_hi = 0x2.b7e151628aed2p+0, e_lo = 0xa.6abf7158809dp-56; - -/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's - approximation to lgamma function. */ - -static const double lgamma_coeff[] = - { - 0x1.5555555555555p-4, - -0xb.60b60b60b60b8p-12, - 0x3.4034034034034p-12, - -0x2.7027027027028p-12, - 0x3.72a3c5631fe46p-12, - -0x7.daac36664f1f4p-12, - 0x1.a41a41a41a41ap-8, - -0x7.90a1b2c3d4e6p-8, - 0x2.dfd2c703c0dp-4, - -0x1.6476701181f3ap+0, - 0xd.672219167003p+0, - -0x9.cd9292e6660d8p+4, - }; - -#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 double poly_coeff[] = - { - /* Interval [-2.125, -2] (polynomial degree 10). */ - -0x1.0b71c5c54d42fp+0, - -0xc.73a1dc05f3758p-4, - -0x1.ec84140851911p-4, - -0xe.37c9da23847e8p-4, - -0x1.03cd87cdc0ac6p-4, - -0xe.ae9aedce12eep-4, - 0x9.b11a1780cfd48p-8, - -0xe.f25fc460bdebp-4, - 0x2.6e984c61ca912p-4, - -0xf.83fea1c6d35p-4, - 0x4.760c8c8909758p-4, - /* Interval [-2.25, -2.125] (polynomial degree 11). */ - -0xf.2930890d7d678p-4, - -0xc.a5cfde054eaa8p-4, - 0x3.9c9e0fdebd99cp-4, - -0x1.02a5ad35619d9p+0, - 0x9.6e9b1167c164p-4, - -0x1.4d8332eba090ap+0, - 0x1.1c0c94b1b2b6p+0, - -0x1.c9a70d138c74ep+0, - 0x1.d7d9cf1d4c196p+0, - -0x2.91fbf4cd6abacp+0, - 0x2.f6751f74b8ff8p+0, - -0x3.e1bb7b09e3e76p+0, - /* Interval [-2.375, -2.25] (polynomial degree 12). */ - -0xd.7d28d505d618p-4, - -0xe.69649a3040958p-4, - 0xb.0d74a2827cd6p-4, - -0x1.924b09228a86ep+0, - 0x1.d49b12bcf6175p+0, - -0x3.0898bb530d314p+0, - 0x4.207a6be8fda4cp+0, - -0x6.39eef56d4e9p+0, - 0x8.e2e42acbccec8p+0, - -0xd.0d91c1e596a68p+0, - 0x1.2e20d7099c585p+4, - -0x1.c4eb6691b4ca9p+4, - 0x2.96a1a11fd85fep+4, - /* Interval [-2.5, -2.375] (polynomial degree 13). */ - -0xb.74ea1bcfff948p-4, - -0x1.2a82bd590c376p+0, - 0x1.88020f828b81p+0, - -0x3.32279f040d7aep+0, - 0x5.57ac8252ce868p+0, - -0x9.c2aedd093125p+0, - 0x1.12c132716e94cp+4, - -0x1.ea94dfa5c0a6dp+4, - 0x3.66b61abfe858cp+4, - -0x6.0cfceb62a26e4p+4, - 0xa.beeba09403bd8p+4, - -0x1.3188d9b1b288cp+8, - 0x2.37f774dd14c44p+8, - -0x3.fdf0a64cd7136p+8, - /* Interval [-2.625, -2.5] (polynomial degree 13). */ - -0x3.d10108c27ebbp-4, - 0x1.cd557caff7d2fp+0, - 0x3.819b4856d36cep+0, - 0x6.8505cbacfc42p+0, - 0xb.c1b2e6567a4dp+0, - 0x1.50a53a3ce6c73p+4, - 0x2.57adffbb1ec0cp+4, - 0x4.2b15549cf400cp+4, - 0x7.698cfd82b3e18p+4, - 0xd.2decde217755p+4, - 0x1.7699a624d07b9p+8, - 0x2.98ecf617abbfcp+8, - 0x4.d5244d44d60b4p+8, - 0x8.e962bf7395988p+8, - /* Interval [-2.75, -2.625] (polynomial degree 12). */ - -0x6.b5d252a56e8a8p-4, - 0x1.28d60383da3a6p+0, - 0x1.db6513ada89bep+0, - 0x2.e217118fa8c02p+0, - 0x4.450112c651348p+0, - 0x6.4af990f589b8cp+0, - 0x9.2db5963d7a238p+0, - 0xd.62c03647da19p+0, - 0x1.379f81f6416afp+4, - 0x1.c5618b4fdb96p+4, - 0x2.9342d0af2ac4ep+4, - 0x3.d9cdf56d2b186p+4, - 0x5.ab9f91d5a27a4p+4, - /* Interval [-2.875, -2.75] (polynomial degree 11). */ - -0x8.a41b1e4f36ff8p-4, - 0xc.da87d3b69dbe8p-4, - 0x1.1474ad5c36709p+0, - 0x1.761ecb90c8c5cp+0, - 0x1.d279bff588826p+0, - 0x2.4e5d003fb36a8p+0, - 0x2.d575575566842p+0, - 0x3.85152b0d17756p+0, - 0x4.5213d921ca13p+0, - 0x5.55da7dfcf69c4p+0, - 0x6.acef729b9404p+0, - 0x8.483cc21dd0668p+0, - /* Interval [-3, -2.875] (polynomial degree 11). */ - -0xa.046d667e468f8p-4, - 0x9.70b88dcc006cp-4, - 0xa.a8a39421c94dp-4, - 0xd.2f4d1363f98ep-4, - 0xd.ca9aa19975b7p-4, - 0xf.cf09c2f54404p-4, - 0x1.04b1365a9adfcp+0, - 0x1.22b54ef213798p+0, - 0x1.2c52c25206bf5p+0, - 0x1.4aa3d798aace4p+0, - 0x1.5c3f278b504e3p+0, - 0x1.7e08292cc347bp+0, - }; - -static const size_t poly_deg[] = - { - 10, - 11, - 12, - 13, - 13, - 12, - 11, - 11, - }; - -static const size_t poly_end[] = - { - 10, - 22, - 35, - 49, - 63, - 76, - 88, - 100, - }; - -/* Compute sin (pi * X) for -0.25 <= X <= 0.5. */ - -static double -lg_sinpi (double x) -{ - if (x <= 0.25) - return __sin (M_PI * x); - else - return __cos (M_PI * (0.5 - x)); -} - -/* Compute cos (pi * X) for -0.25 <= X <= 0.5. */ - -static double -lg_cospi (double x) -{ - if (x <= 0.25) - return __cos (M_PI * x); - else - return __sin (M_PI * (0.5 - x)); -} - -/* Compute cot (pi * X) for -0.25 <= X <= 0.5. */ - -static double -lg_cotpi (double x) -{ - return lg_cospi (x) / lg_sinpi (x); -} - -/* Compute lgamma of a negative argument -28 < X < -2, setting - *SIGNGAMP accordingly. */ - -double -__lgamma_neg (double x, int *signgamp) -{ - /* Determine the half-integer region X lies in, handle exact - integers and determine the sign of the result. */ - int i = floor (-2 * x); - if ((i & 1) == 0 && i == -2 * x) - return 1.0 / 0.0; - double xn = ((i & 1) == 0 ? -i / 2 : (-i - 1) / 2); - i -= 4; - *signgamp = ((i & 2) == 0 ? -1 : 1); - - SET_RESTORE_ROUND (FE_TONEAREST); - - /* Expand around the zero X0 = X0_HI + X0_LO. */ - double x0_hi = lgamma_zeros[i][0], x0_lo = lgamma_zeros[i][1]; - double 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 = floor (-8 * x) - 16; - double xm = (-33 - 2 * j) * 0.0625; - double x_adj = x - xm; - size_t deg = poly_deg[j]; - size_t end = poly_end[j]; - double g = poly_coeff[end]; - for (size_t j = 1; j <= deg; j++) - g = g * x_adj + poly_coeff[end - j]; - return __log1p (g * xdiff / (x - xn)); - } - - /* The result we want is log (sinpi (X0) / sinpi (X)) - + log (gamma (1 - X0) / gamma (1 - X)). */ - double x_idiff = fabs (xn - x), x0_idiff = fabs (xn - x0_hi - x0_lo); - double log_sinpi_ratio; - if (x0_idiff < x_idiff * 0.5) - /* Use log not log1p to avoid inaccuracy from log1p of arguments - close to -1. */ - log_sinpi_ratio = __ieee754_log (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. */ - double x0diff2 = ((i & 1) == 0 ? xdiff : -xdiff) * 0.5; - double sx0d2 = lg_sinpi (x0diff2); - double cx0d2 = lg_cospi (x0diff2); - log_sinpi_ratio = __log1p (2 * sx0d2 - * (-sx0d2 + cx0d2 * lg_cotpi (x_idiff))); - } - - double log_gamma_ratio; - double y0 = math_narrow_eval (1 - x0_hi); - double y0_eps = -x0_hi + (1 - y0) - x0_lo; - double y = math_narrow_eval (1 - x); - double 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. First, we may need to - adjust into the range where Stirling's approximation is - sufficiently accurate. */ - double log_gamma_adj = 0; - if (i < 6) - { - int n_up = (7 - i) / 2; - double ny0, ny0_eps, ny, ny_eps; - ny0 = math_narrow_eval (y0 + n_up); - ny0_eps = y0 - (ny0 - n_up) + y0_eps; - y0 = ny0; - y0_eps = ny0_eps; - ny = math_narrow_eval (y + n_up); - ny_eps = y - (ny - n_up) + y_eps; - y = ny; - y_eps = ny_eps; - double prodm1 = __lgamma_product (xdiff, y - n_up, y_eps, n_up); - log_gamma_adj = -__log1p (prodm1); - } - double log_gamma_high - = (xdiff * __log1p ((y0 - e_hi - e_lo + y0_eps) / e_hi) - + (y - 0.5 + y_eps) * __log1p (xdiff / y) + log_gamma_adj); - /* Compute the sum of (B_2k / 2k(2k-1))(Y0^-(2k-1) - Y^-(2k-1)). */ - double y0r = 1 / y0, yr = 1 / y; - double y0r2 = y0r * y0r, yr2 = yr * yr; - double rdiff = -xdiff / (y * y0); - double bterm[NCOEFF]; - double dlast = rdiff, elast = rdiff * yr * (yr + y0r); - bterm[0] = dlast * lgamma_coeff[0]; - for (size_t j = 1; j < NCOEFF; j++) - { - double dnext = dlast * y0r2 + elast; - double enext = elast * yr2; - bterm[j] = dnext * lgamma_coeff[j]; - dlast = dnext; - elast = enext; - } - double 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; -} diff --git a/sysdeps/ieee754/dbl-64/lgamma_product.c b/sysdeps/ieee754/dbl-64/lgamma_product.c deleted file mode 100644 index 549f9688c7..0000000000 --- a/sysdeps/ieee754/dbl-64/lgamma_product.c +++ /dev/null @@ -1,52 +0,0 @@ -/* Compute a product of 1 + (T/X), 1 + (T/(X+1)), .... - Copyright (C) 2015-2025 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 - . */ - -#include -#include -#include - -/* Compute the product of 1 + (T / (X + X_EPS)), 1 + (T / (X + X_EPS + - 1)), ..., 1 + (T / (X + X_EPS + N - 1)), minus 1. 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. */ - -double -__lgamma_product (double t, double x, double x_eps, int n) -{ - double ret = 0, ret_eps = 0; - for (int i = 0; i < n; i++) - { - double xi = x + i; - double quot = t / xi; - double mhi, mlo; - mul_split (&mhi, &mlo, quot, xi); - double quot_lo = (t - mhi - mlo) / xi - t * x_eps / (xi * xi); - /* We want (1 + RET + RET_EPS) * (1 + QUOT + QUOT_LO) - 1. */ - double rhi, rlo; - mul_split (&rhi, &rlo, ret, quot); - double rpq = ret + quot; - double rpq_eps = (ret - rpq) + quot; - double nret = rpq + rhi; - double nret_eps = (rpq - nret) + rhi; - ret_eps += (rpq_eps + nret_eps + rlo + ret_eps * quot - + quot_lo + quot_lo * (ret + ret_eps)); - ret = nret; - } - return ret + ret_eps; -} diff --git a/sysdeps/ieee754/dbl-64/libm-test-ulps b/sysdeps/ieee754/dbl-64/libm-test-ulps index a1cb2490fd..38bb2423d9 100644 --- a/sysdeps/ieee754/dbl-64/libm-test-ulps +++ b/sysdeps/ieee754/dbl-64/libm-test-ulps @@ -34,3 +34,15 @@ double: 0 Function: "atanh_upward": double: 0 + +Function: "lgamma": +double: 0 + +Function: "lgamma_downward": +double: 0 + +Function: "lgamma_towardzero": +double: 0 + +Function: "lgamma_upward": +double: 0 diff --git a/sysdeps/ieee754/dbl-64/math_config.h b/sysdeps/ieee754/dbl-64/math_config.h index da687541be..3e30835edd 100644 --- a/sysdeps/ieee754/dbl-64/math_config.h +++ b/sysdeps/ieee754/dbl-64/math_config.h @@ -38,6 +38,33 @@ # define WANT_ERRNO_UFLOW (WANT_ROUNDING && WANT_ERRNO) #endif +#ifndef ROUNDEVEN_INTRINSICS +/* When set, roundeven_finite will route to the internal roundeven function. */ +# define ROUNDEVEN_INTRINSICS 1 +#endif + +/* Round x to nearest integer value in floating-point format, rounding halfway + cases to even. If the input is non finite the result is unspecified. */ +static inline double +roundeven_finite (double x) +{ + if (!isfinite (x)) + __builtin_unreachable (); +#if ROUNDEVEN_INTRINSICS + return roundeven (x); +#else + double y = round (x); + if (fabs (x - y) == 0.5) + { + union { double f; uint64_t i; } u = {y}; + union { double f; uint64_t i; } v = {y - copysign (1.0, x)}; + if (__builtin_ctzll (v.i) > __builtin_ctzll (u.i)) + y = v.f; + } + return y; +#endif +} + #ifndef TOINT_INTRINSICS /* When set, the roundtoint and converttoint functions are provided with the semantics documented below. */ diff --git a/sysdeps/ieee754/flt-32/lgamma_negf.c b/sysdeps/ieee754/flt-32/lgamma_negf.c deleted file mode 100644 index 1cc8931700..0000000000 --- a/sysdeps/ieee754/flt-32/lgamma_negf.c +++ /dev/null @@ -1 +0,0 @@ -/* Not needed. */ diff --git a/sysdeps/ieee754/flt-32/lgamma_productf.c b/sysdeps/ieee754/flt-32/lgamma_productf.c deleted file mode 100644 index 1cc8931700..0000000000 --- a/sysdeps/ieee754/flt-32/lgamma_productf.c +++ /dev/null @@ -1 +0,0 @@ -/* Not needed. */ diff --git a/sysdeps/ieee754/ldbl-96/lgamma_product.c b/sysdeps/ieee754/ldbl-96/lgamma_product.c deleted file mode 100644 index 06ea3430a7..0000000000 --- a/sysdeps/ieee754/ldbl-96/lgamma_product.c +++ /dev/null @@ -1,37 +0,0 @@ -/* Compute a product of 1 + (T/X), 1 + (T/(X+1)), .... - Copyright (C) 2015-2025 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 - . */ - -#include -#include -#include - -/* Compute the product of 1 + (T / (X + X_EPS)), 1 + (T / (X + X_EPS + - 1)), ..., 1 + (T / (X + X_EPS + N - 1)), minus 1. 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. */ - -double -__lgamma_product (double t, double x, double x_eps, int n) -{ - long double x_full = (long double) x + (long double) x_eps; - long double ret = 0; - for (int i = 0; i < n; i++) - ret += (t / (x_full + i)) * (1 + ret); - return ret; -}