diff options
Diffstat (limited to 'sysdeps/ieee754/ldbl-128')
55 files changed, 976 insertions, 232 deletions
diff --git a/sysdeps/ieee754/ldbl-128/e_asinl.c b/sysdeps/ieee754/ldbl-128/e_asinl.c index 691ac267fb..5a0e473ef0 100644 --- a/sysdeps/ieee754/ldbl-128/e_asinl.c +++ b/sysdeps/ieee754/ldbl-128/e_asinl.c @@ -153,11 +153,7 @@ __ieee754_asinl (long double x) { if (ix < 0x3fc60000) /* |x| < 2**-57 */ { - if (fabsl (x) < LDBL_MIN) - { - long double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); long double force_inexact = huge + x; math_force_eval (force_inexact); return x; /* return x with inexact if x!=0 */ diff --git a/sysdeps/ieee754/ldbl-128/e_atanhl.c b/sysdeps/ieee754/ldbl-128/e_atanhl.c index a5a7ee0712..7fa53ef436 100644 --- a/sysdeps/ieee754/ldbl-128/e_atanhl.c +++ b/sysdeps/ieee754/ldbl-128/e_atanhl.c @@ -60,11 +60,7 @@ __ieee754_atanhl(long double x) } if(ix<0x3fc60000 && (huge+x)>zero) /* x < 2^-57 */ { - if (fabsl (x) < LDBL_MIN) - { - long double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); return x; } diff --git a/sysdeps/ieee754/ldbl-128/e_exp10l.c b/sysdeps/ieee754/ldbl-128/e_exp10l.c index c5b5cb7505..987002567a 100644 --- a/sysdeps/ieee754/ldbl-128/e_exp10l.c +++ b/sysdeps/ieee754/ldbl-128/e_exp10l.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2012-2015 Free Software Foundation, Inc. +/* Copyright (C) 2012-2016 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 diff --git a/sysdeps/ieee754/ldbl-128/e_expl.c b/sysdeps/ieee754/ldbl-128/e_expl.c index 1cd095cb77..7b71e644be 100644 --- a/sysdeps/ieee754/ldbl-128/e_expl.c +++ b/sysdeps/ieee754/ldbl-128/e_expl.c @@ -1,5 +1,5 @@ /* Quad-precision floating point e^x. - Copyright (C) 1999-2015 Free Software Foundation, Inc. + Copyright (C) 1999-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jj@ultra.linux.cz> Partly based on double-precision code @@ -232,11 +232,7 @@ __ieee754_expl (long double x) else { result *= scale_u.d; - if (result < LDBL_MIN) - { - long double force_underflow = result * result; - math_force_eval (force_underflow); - } + math_check_force_underflow_nonneg (result); return result; } } diff --git a/sysdeps/ieee754/ldbl-128/e_gammal_r.c b/sysdeps/ieee754/ldbl-128/e_gammal_r.c index c51b050e0e..d0286e31eb 100644 --- a/sysdeps/ieee754/ldbl-128/e_gammal_r.c +++ b/sysdeps/ieee754/ldbl-128/e_gammal_r.c @@ -1,5 +1,5 @@ /* Implementation of gamma function according to ISO C. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and Jakub Jelinek <jj@ultra.linux.cz, 1999. @@ -194,6 +194,7 @@ __ieee754_gammal_r (long double x, int *signgamp) ret = M_PIl / (-x * sinpix * gammal_positive (-x, &exp2_adj)); ret = __scalbnl (ret, -exp2_adj); + math_check_force_underflow_nonneg (ret); } } } diff --git a/sysdeps/ieee754/ldbl-128/e_hypotl.c b/sysdeps/ieee754/ldbl-128/e_hypotl.c index 01444cfb4e..80e5e38c72 100644 --- a/sysdeps/ieee754/ldbl-128/e_hypotl.c +++ b/sysdeps/ieee754/ldbl-128/e_hypotl.c @@ -130,7 +130,9 @@ __ieee754_hypotl(long double x, long double y) t1 = 1.0L; GET_LDOUBLE_MSW64(high,t1); SET_LDOUBLE_MSW64(t1,high+(k<<48)); - return t1*w; + w *= t1; + math_check_force_underflow_nonneg (w); + return w; } else return w; } strong_alias (__ieee754_hypotl, __hypotl_finite) diff --git a/sysdeps/ieee754/ldbl-128/e_j0l.c b/sysdeps/ieee754/ldbl-128/e_j0l.c index 1320de7044..c208916a79 100644 --- a/sysdeps/ieee754/ldbl-128/e_j0l.c +++ b/sysdeps/ieee754/ldbl-128/e_j0l.c @@ -693,6 +693,8 @@ __ieee754_j0l (long double x) xx = fabsl (x); if (xx <= 2.0L) { + if (xx < 0x1p-57L) + return 1.0L; /* 0 <= x <= 2 */ z = xx * xx; p = z * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D); diff --git a/sysdeps/ieee754/ldbl-128/e_j1l.c b/sysdeps/ieee754/ldbl-128/e_j1l.c index 591c38efd0..f5b04c073d 100644 --- a/sysdeps/ieee754/ldbl-128/e_j1l.c +++ b/sysdeps/ieee754/ldbl-128/e_j1l.c @@ -700,11 +700,9 @@ __ieee754_j1l (long double x) if (xx <= 0x1p-58L) { long double ret = x * 0.5L; - if (fabsl (ret) < LDBL_MIN) - { - long double force_underflow = ret * ret; - math_force_eval (force_underflow); - } + math_check_force_underflow (ret); + if (ret == 0) + __set_errno (ERANGE); return ret; } if (xx <= 2.0L) diff --git a/sysdeps/ieee754/ldbl-128/e_jnl.c b/sysdeps/ieee754/ldbl-128/e_jnl.c index a419994a86..98669e6e3e 100644 --- a/sysdeps/ieee754/ldbl-128/e_jnl.c +++ b/sysdeps/ieee754/ldbl-128/e_jnl.c @@ -296,12 +296,12 @@ __ieee754_jnl (int n, long double x) ret = b; } if (ret == 0) - ret = __copysignl (LDBL_MIN, ret) * LDBL_MIN; - else if (fabsl (ret) < LDBL_MIN) { - long double force_underflow = ret * ret; - math_force_eval (force_underflow); + ret = __copysignl (LDBL_MIN, ret) * LDBL_MIN; + __set_errno (ERANGE); } + else + math_check_force_underflow (ret); return ret; } strong_alias (__ieee754_jnl, __jnl_finite) diff --git a/sysdeps/ieee754/ldbl-128/e_lgammal_r.c b/sysdeps/ieee754/ldbl-128/e_lgammal_r.c index d8a5e5b9ec..5b513ea1df 100644 --- a/sysdeps/ieee754/ldbl-128/e_lgammal_r.c +++ b/sysdeps/ieee754/ldbl-128/e_lgammal_r.c @@ -70,19 +70,15 @@ #include <math.h> #include <math_private.h> -#include <libc-internal.h> #include <float.h> -/* BZ#16347: ldbl-128ibm uses this file as is, however the MAXLGM - definition overflows for IBM long double. This directive prevents the - overflow warnings until IBM long double version is fixed. */ static const long double PIL = 3.1415926535897932384626433832795028841972E0L; -DIAG_PUSH_NEEDS_COMMENT; -DIAG_IGNORE_NEEDS_COMMENT (4.6, "-Woverflow"); +#if LDBL_MANT_DIG == 106 +static const long double MAXLGM = 0x5.d53649e2d469dbc1f01e99fd66p+1012L; +#else static const long double MAXLGM = 1.0485738685148938358098967157129705071571E4928L; -DIAG_POP_NEEDS_COMMENT; +#endif static const long double one = 1.0L; -static const long double zero = 0.0L; static const long double huge = LDBL_MAX; /* log gamma(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x P(1/x^2) @@ -781,12 +777,14 @@ __ieee754_lgammal_r (long double x, int *signgamp) if (x < 0.0L) { + if (x < -2.0L && x > (LDBL_MANT_DIG == 106 ? -48.0L : -50.0L)) + return __lgamma_negl (x, signgamp); q = -x; p = __floorl (q); if (p == q) return (one / (p - p)); - i = p; - if ((i & 1) == 0) + long double halfp = p * 0.5L; + if (halfp == __floorl (halfp)) *signgamp = -1; else *signgamp = 1; @@ -1034,6 +1032,8 @@ __ieee754_lgammal_r (long double x, int *signgamp) if (x > MAXLGM) return (*signgamp * huge * huge); + if (x > 0x1p120L) + return x * (__logl (x) - 1.0L); q = ls2pi - x; q = (x - 0.5L) * __logl (x) + q; if (x > 1.0e18L) diff --git a/sysdeps/ieee754/ldbl-128/e_log2l.c b/sysdeps/ieee754/ldbl-128/e_log2l.c index 991a3b73e2..6b1faa0523 100644 --- a/sysdeps/ieee754/ldbl-128/e_log2l.c +++ b/sysdeps/ieee754/ldbl-128/e_log2l.c @@ -171,8 +171,7 @@ deval (long double x, const long double *p, int n) long double -__ieee754_log2l (x) - long double x; +__ieee754_log2l (long double x) { long double z; long double y; diff --git a/sysdeps/ieee754/ldbl-128/e_powl.c b/sysdeps/ieee754/ldbl-128/e_powl.c index f531385232..7f3037fb51 100644 --- a/sysdeps/ieee754/ldbl-128/e_powl.c +++ b/sysdeps/ieee754/ldbl-128/e_powl.c @@ -435,7 +435,11 @@ __ieee754_powl (long double x, long double y) j = o.parts32.w0; j += (n << 16); if ((j >> 16) <= 0) - z = __scalbnl (z, n); /* subnormal output */ + { + z = __scalbnl (z, n); /* subnormal output */ + long double force_underflow = z * z; + math_force_eval (force_underflow); + } else { o.parts32.w0 = j; diff --git a/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c b/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c index 2b60263119..101a4c9015 100644 --- a/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c +++ b/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c @@ -1,5 +1,5 @@ /* Quad-precision floating point argument reduction. - Copyright (C) 1999-2015 Free Software Foundation, Inc. + Copyright (C) 1999-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jj@ultra.linux.cz> diff --git a/sysdeps/ieee754/ldbl-128/e_sinhl.c b/sysdeps/ieee754/ldbl-128/e_sinhl.c index 1ca3c6e507..11974a39af 100644 --- a/sysdeps/ieee754/ldbl-128/e_sinhl.c +++ b/sysdeps/ieee754/ldbl-128/e_sinhl.c @@ -53,6 +53,7 @@ * only sinhl(0)=0 is exact for finite x. */ +#include <float.h> #include <math.h> #include <math_private.h> @@ -86,8 +87,11 @@ __ieee754_sinhl (long double x) if (ix <= 0x40044000) { if (ix < 0x3fc60000) /* |x| < 2^-57 */ - if (shuge + x > one) - return x; /* sinh(tiny) = tiny with inexact */ + { + math_check_force_underflow (x); + if (shuge + x > one) + return x; /* sinh(tiny) = tiny with inexact */ + } t = __expm1l (u.value); if (ix < 0x3fff0000) return h * (2.0 * t - t * t / (t + one)); diff --git a/sysdeps/ieee754/ldbl-128/gamma_productl.c b/sysdeps/ieee754/ldbl-128/gamma_productl.c index 32990b8f1c..849b57d95d 100644 --- a/sysdeps/ieee754/ldbl-128/gamma_productl.c +++ b/sysdeps/ieee754/ldbl-128/gamma_productl.c @@ -1,5 +1,5 @@ /* Compute a product of X, X+1, ..., with an error estimate. - Copyright (C) 2013-2015 Free Software Foundation, Inc. + Copyright (C) 2013-2016 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 diff --git a/sysdeps/ieee754/ldbl-128/ieee754.h b/sysdeps/ieee754/ldbl-128/ieee754.h index d8d7aa69a2..fb42140933 100644 --- a/sysdeps/ieee754/ldbl-128/ieee754.h +++ b/sysdeps/ieee754/ldbl-128/ieee754.h @@ -1,4 +1,4 @@ -/* Copyright (C) 1992-2015 Free Software Foundation, Inc. +/* Copyright (C) 1992-2016 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 diff --git a/sysdeps/ieee754/ldbl-128/k_cosl.c b/sysdeps/ieee754/ldbl-128/k_cosl.c index f2743ddf5f..3985b1225b 100644 --- a/sysdeps/ieee754/ldbl-128/k_cosl.c +++ b/sysdeps/ieee754/ldbl-128/k_cosl.c @@ -1,5 +1,5 @@ /* Quad-precision floating point cosine on <-pi/4,pi/4>. - Copyright (C) 1999-2015 Free Software Foundation, Inc. + Copyright (C) 1999-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jj@ultra.linux.cz> diff --git a/sysdeps/ieee754/ldbl-128/k_sincosl.c b/sysdeps/ieee754/ldbl-128/k_sincosl.c index 7b5c4b0e81..404df352d3 100644 --- a/sysdeps/ieee754/ldbl-128/k_sincosl.c +++ b/sysdeps/ieee754/ldbl-128/k_sincosl.c @@ -1,5 +1,5 @@ /* Quad-precision floating point sine and cosine on <-pi/4,pi/4>. - Copyright (C) 1999-2015 Free Software Foundation, Inc. + Copyright (C) 1999-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jj@ultra.linux.cz> @@ -111,11 +111,7 @@ __kernel_sincosl(long double x, long double y, long double *sinx, long double *c polynomial of degree 16(17). */ if (tix < 0x3fc60000) /* |x| < 2^-57 */ { - if (fabsl (x) < LDBL_MIN) - { - long double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); if (!((int)x)) /* generate inexact */ { *sinx = x; diff --git a/sysdeps/ieee754/ldbl-128/k_sinl.c b/sysdeps/ieee754/ldbl-128/k_sinl.c index 04d539ff1a..6290b173b4 100644 --- a/sysdeps/ieee754/ldbl-128/k_sinl.c +++ b/sysdeps/ieee754/ldbl-128/k_sinl.c @@ -1,5 +1,5 @@ /* Quad-precision floating point sine on <-pi/4,pi/4>. - Copyright (C) 1999-2015 Free Software Foundation, Inc. + Copyright (C) 1999-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jj@ultra.linux.cz> @@ -92,11 +92,7 @@ __kernel_sinl(long double x, long double y, int iy) polynomial of degree 17. */ if (tix < 0x3fc60000) /* |x| < 2^-57 */ { - if (fabsl (x) < LDBL_MIN) - { - long double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); if (!((int)x)) return x; /* generate inexact */ } z = x * x; diff --git a/sysdeps/ieee754/ldbl-128/k_tanl.c b/sysdeps/ieee754/ldbl-128/k_tanl.c index dfba2d9a76..6bb221e4a6 100644 --- a/sysdeps/ieee754/ldbl-128/k_tanl.c +++ b/sysdeps/ieee754/ldbl-128/k_tanl.c @@ -56,6 +56,7 @@ * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y))) */ +#include <float.h> #include <libc-internal.h> #include <math.h> #include <math_private.h> @@ -98,8 +99,13 @@ __kernel_tanl (long double x, long double y, int iy) if ((ix | u.parts32.w1 | u.parts32.w2 | u.parts32.w3 | (iy + 1)) == 0) return one / fabs (x); + else if (iy == 1) + { + math_check_force_underflow (x); + return x; + } else - return (iy == 1) ? x : -one / x; + return -one / x; } } if (ix >= 0x3ffe5942) /* |x| >= 0.6743316650390625 */ @@ -135,11 +141,7 @@ __kernel_tanl (long double x, long double y, int iy) uninitialized although in the cases where it is used it has always been set. */ DIAG_PUSH_NEEDS_COMMENT; -#if __GNUC_PREREQ (4, 7) DIAG_IGNORE_NEEDS_COMMENT (5, "-Wmaybe-uninitialized"); -#else - DIAG_IGNORE_NEEDS_COMMENT (5, "-Wuninitialized"); -#endif if (sign < 0) w = -w; DIAG_POP_NEEDS_COMMENT; diff --git a/sysdeps/ieee754/ldbl-128/ldbl2mpn.c b/sysdeps/ieee754/ldbl-128/ldbl2mpn.c index 08953af819..e3f97623ed 100644 --- a/sysdeps/ieee754/ldbl-128/ldbl2mpn.c +++ b/sysdeps/ieee754/ldbl-128/ldbl2mpn.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1995-2015 Free Software Foundation, Inc. +/* Copyright (C) 1995-2016 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 diff --git a/sysdeps/ieee754/ldbl-128/lgamma_negl.c b/sysdeps/ieee754/ldbl-128/lgamma_negl.c new file mode 100644 index 0000000000..df46199b82 --- /dev/null +++ b/sysdeps/ieee754/ldbl-128/lgamma_negl.c @@ -0,0 +1,551 @@ +/* lgammal expanding around zeros. + Copyright (C) 2015-2016 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <float.h> +#include <math.h> +#include <math_private.h> + +static const long double lgamma_zeros[][2] = + { + { -0x2.74ff92c01f0d82abec9f315f1a08p+0L, 0xe.d3ccb7fb2658634a2b9f6b2ba81p-116L }, + { -0x2.bf6821437b20197995a4b4641eaep+0L, -0xb.f4b00b4829f961e428533e6ad048p-116L }, + { -0x3.24c1b793cb35efb8be699ad3d9bap+0L, -0x6.5454cb7fac60e3f16d9d7840c2ep-116L }, + { -0x3.f48e2a8f85fca170d4561291236cp+0L, -0xc.320a4887d1cb4c711828a75d5758p-116L }, + { -0x4.0a139e16656030c39f0b0de18114p+0L, 0x1.53e84029416e1242006b2b3d1cfp-112L }, + { -0x4.fdd5de9bbabf3510d0aa40769884p+0L, -0x1.01d7d78125286f78d1e501f14966p-112L }, + { -0x5.021a95fc2db6432a4c56e595394cp+0L, -0x1.ecc6af0430d4fe5746fa7233356fp-112L }, + { -0x5.ffa4bd647d0357dd4ed62cbd31ecp+0L, -0x1.f8e3f8e5deba2d67dbd70dd96ce1p-112L }, + { -0x6.005ac9625f233b607c2d96d16384p+0L, -0x1.cb86ac569340cf1e5f24df7aab7bp-112L }, + { -0x6.fff2fddae1bbff3d626b65c23fd4p+0L, 0x1.e0bfcff5c457ebcf4d3ad9674167p-112L }, + { -0x7.000cff7b7f87adf4482dcdb98784p+0L, 0x1.54d99e35a74d6407b80292df199fp-112L }, + { -0x7.fffe5fe05673c3ca9e82b522b0ccp+0L, 0x1.62d177c832e0eb42c2faffd1b145p-112L }, + { -0x8.0001a01459fc9f60cb3cec1cec88p+0L, 0x2.8998835ac7277f7bcef67c47f188p-112L }, + { -0x8.ffffd1c425e80ffc864e95749258p+0L, -0x1.e7e20210e7f81cf781b44e9d2b02p-112L }, + { -0x9.00002e3bb47d86d6d843fedc352p+0L, 0x2.14852f613a16291751d2ab751f7ep-112L }, + { -0x9.fffffb606bdfdcd062ae77a50548p+0L, 0x3.962d1490cc2e8f031c7007eaa1ap-116L }, + { -0xa.0000049f93bb9927b45d95e1544p+0L, -0x1.e03086db9146a9287bd4f2172d5ap-112L }, + { -0xa.ffffff9466e9f1b36dacd2adbd18p+0L, -0xd.05a4e458062f3f95345a4d9c9b6p-116L }, + { -0xb.0000006b9915315d965a6ffea41p+0L, 0x1.b415c6fff233e7b7fdc3a094246fp-112L }, + { -0xb.fffffff7089387387de41acc3d4p+0L, 0x3.687427c6373bd74a10306e10a28ep-112L }, + { -0xc.00000008f76c7731567c0f0250fp+0L, -0x3.87920df5675833859190eb128ef6p-112L }, + { -0xc.ffffffff4f6dcf617f97a5ffc758p+0L, 0x2.ab72d76f32eaee2d1a42ed515d3ap-116L }, + { -0xd.00000000b092309c06683dd1b9p+0L, -0x3.e3700857a15c19ac5a611de9688ap-112L }, + { -0xd.fffffffff36345ab9e184a3e09dp+0L, -0x1.176dc48e47f62d917973dd44e553p-112L }, + { -0xe.000000000c9cba545e94e75ec57p+0L, -0x1.8f753e2501e757a17cf2ecbeeb89p-112L }, + { -0xe.ffffffffff28c060c6604ef3037p+0L, -0x1.f89d37357c9e3dc17c6c6e63becap-112L }, + { -0xf.0000000000d73f9f399bd0e420f8p+0L, -0x5.e9ee31b0b890744fc0e3fbc01048p-116L }, + { -0xf.fffffffffff28c060c6621f512e8p+0L, 0xd.1b2eec9d960bd9adc5be5f5fa5p-116L }, + { -0x1.000000000000d73f9f399da1424cp+4L, 0x6.c46e0e88305d2800f0e414c506a8p-116L }, + { -0x1.0ffffffffffff3569c47e7a93e1cp+4L, -0x4.6a08a2e008a998ebabb8087efa2cp-112L }, + { -0x1.1000000000000ca963b818568887p+4L, -0x6.ca5a3a64ec15db0a95caf2c9ffb4p-112L }, + { -0x1.1fffffffffffff4bec3ce234132dp+4L, -0x8.b2b726187c841cb92cd5221e444p-116L }, + { -0x1.20000000000000b413c31dcbeca5p+4L, 0x3.c4d005344b6cd0e7231120294abcp-112L }, + { -0x1.2ffffffffffffff685b25cbf5f54p+4L, -0x5.ced932e38485f7dd296b8fa41448p-112L }, + { -0x1.30000000000000097a4da340a0acp+4L, 0x7.e484e0e0ffe38d406ebebe112f88p-112L }, + { -0x1.3fffffffffffffff86af516ff7f7p+4L, -0x6.bd67e720d57854502b7db75e1718p-112L }, + { -0x1.40000000000000007950ae900809p+4L, 0x6.bec33375cac025d9c073168c5d9p-112L }, + { -0x1.4ffffffffffffffffa391c4248c3p+4L, 0x5.c63022b62b5484ba346524db607p-112L }, + { -0x1.500000000000000005c6e3bdb73dp+4L, -0x5.c62f55ed5322b2685c5e9a51e6a8p-112L }, + { -0x1.5fffffffffffffffffbcc71a492p+4L, -0x1.eb5aeb96c74d7ad25e060528fb5p-112L }, + { -0x1.6000000000000000004338e5b6ep+4L, 0x1.eb5aec04b2f2eb663e4e3d8a018cp-112L }, + { -0x1.6ffffffffffffffffffd13c97d9dp+4L, -0x3.8fcc4d08d6fe5aa56ab04307ce7ep-112L }, + { -0x1.70000000000000000002ec368263p+4L, 0x3.8fcc4d090cee2f5d0b69a99c353cp-112L }, + { -0x1.7fffffffffffffffffffe0d30fe7p+4L, 0x7.2f577cca4b4c8cb1dc14001ac5ecp-112L }, + { -0x1.800000000000000000001f2cf019p+4L, -0x7.2f577cca4b3442e35f0040b3b9e8p-112L }, + { -0x1.8ffffffffffffffffffffec0c332p+4L, -0x2.e9a0572b1bb5b95f346a92d67a6p-112L }, + { -0x1.90000000000000000000013f3ccep+4L, 0x2.e9a0572b1bb5c371ddb3561705ap-112L }, + { -0x1.9ffffffffffffffffffffff3b8bdp+4L, -0x1.cad8d32e386fd783e97296d63dcbp-116L }, + { -0x1.a0000000000000000000000c4743p+4L, 0x1.cad8d32e386fd7c1ab8c1fe34c0ep-116L }, + { -0x1.afffffffffffffffffffffff8b95p+4L, -0x3.8f48cc5737d5979c39db806c5406p-112L }, + { -0x1.b00000000000000000000000746bp+4L, 0x3.8f48cc5737d5979c3b3a6bda06f6p-112L }, + { -0x1.bffffffffffffffffffffffffbd8p+4L, 0x6.2898d42174dcf171470d8c8c6028p-112L }, + { -0x1.c000000000000000000000000428p+4L, -0x6.2898d42174dcf171470d18ba412cp-112L }, + { -0x1.cfffffffffffffffffffffffffdbp+4L, -0x4.c0ce9794ea50a839e311320bde94p-112L }, + { -0x1.d000000000000000000000000025p+4L, 0x4.c0ce9794ea50a839e311322f7cf8p-112L }, + { -0x1.dfffffffffffffffffffffffffffp+4L, 0x3.932c5047d60e60caded4c298a174p-112L }, + { -0x1.e000000000000000000000000001p+4L, -0x3.932c5047d60e60caded4c298973ap-112L }, + { -0x1.fp+4L, 0xa.1a6973c1fade2170f7237d35fe3p-116L }, + { -0x1.fp+4L, -0xa.1a6973c1fade2170f7237d35fe08p-116L }, + { -0x2p+4L, 0x5.0d34b9e0fd6f10b87b91be9aff1p-120L }, + { -0x2p+4L, -0x5.0d34b9e0fd6f10b87b91be9aff0cp-120L }, + { -0x2.1p+4L, 0x2.73024a9ba1aa36a7059bff52e844p-124L }, + { -0x2.1p+4L, -0x2.73024a9ba1aa36a7059bff52e844p-124L }, + { -0x2.2p+4L, 0x1.2710231c0fd7a13f8a2b4af9d6b7p-128L }, + { -0x2.2p+4L, -0x1.2710231c0fd7a13f8a2b4af9d6b7p-128L }, + { -0x2.3p+4L, 0x8.6e2ce38b6c8f9419e3fad3f0312p-136L }, + { -0x2.3p+4L, -0x8.6e2ce38b6c8f9419e3fad3f0312p-136L }, + { -0x2.4p+4L, 0x3.bf30652185952560d71a254e4eb8p-140L }, + { -0x2.4p+4L, -0x3.bf30652185952560d71a254e4eb8p-140L }, + { -0x2.5p+4L, 0x1.9ec8d1c94e85af4c78b15c3d89d3p-144L }, + { -0x2.5p+4L, -0x1.9ec8d1c94e85af4c78b15c3d89d3p-144L }, + { -0x2.6p+4L, 0xa.ea565ce061d57489e9b85276274p-152L }, + { -0x2.6p+4L, -0xa.ea565ce061d57489e9b85276274p-152L }, + { -0x2.7p+4L, 0x4.7a6512692eb37804111dabad30ecp-156L }, + { -0x2.7p+4L, -0x4.7a6512692eb37804111dabad30ecp-156L }, + { -0x2.8p+4L, 0x1.ca8ed42a12ae3001a07244abad2bp-160L }, + { -0x2.8p+4L, -0x1.ca8ed42a12ae3001a07244abad2bp-160L }, + { -0x2.9p+4L, 0xb.2f30e1ce812063f12e7e8d8d96e8p-168L }, + { -0x2.9p+4L, -0xb.2f30e1ce812063f12e7e8d8d96e8p-168L }, + { -0x2.ap+4L, 0x4.42bd49d4c37a0db136489772e428p-172L }, + { -0x2.ap+4L, -0x4.42bd49d4c37a0db136489772e428p-172L }, + { -0x2.bp+4L, 0x1.95db45257e5122dcbae56def372p-176L }, + { -0x2.bp+4L, -0x1.95db45257e5122dcbae56def372p-176L }, + { -0x2.cp+4L, 0x9.3958d81ff63527ecf993f3fb6f48p-184L }, + { -0x2.cp+4L, -0x9.3958d81ff63527ecf993f3fb6f48p-184L }, + { -0x2.dp+4L, 0x3.47970e4440c8f1c058bd238c9958p-188L }, + { -0x2.dp+4L, -0x3.47970e4440c8f1c058bd238c9958p-188L }, + { -0x2.ep+4L, 0x1.240804f65951062ca46e4f25c608p-192L }, + { -0x2.ep+4L, -0x1.240804f65951062ca46e4f25c608p-192L }, + { -0x2.fp+4L, 0x6.36a382849fae6de2d15362d8a394p-200L }, + { -0x2.fp+4L, -0x6.36a382849fae6de2d15362d8a394p-200L }, + { -0x3p+4L, 0x2.123680d6dfe4cf4b9b1bcb9d8bdcp-204L }, + { -0x3p+4L, -0x2.123680d6dfe4cf4b9b1bcb9d8bdcp-204L }, + { -0x3.1p+4L, 0xa.d21786ff5842eca51fea0870919p-212L }, + { -0x3.1p+4L, -0xa.d21786ff5842eca51fea0870919p-212L }, + { -0x3.2p+4L, 0x3.766dedc259af040be140a68a6c04p-216L }, + }; + +static const long double e_hi = 0x2.b7e151628aed2a6abf7158809cf4p+0L; +static const long double e_lo = 0xf.3c762e7160f38b4da56a784d9048p-116L; + + +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's + approximation to lgamma function. */ + +static const long double lgamma_coeff[] = + { + 0x1.5555555555555555555555555555p-4L, + -0xb.60b60b60b60b60b60b60b60b60b8p-12L, + 0x3.4034034034034034034034034034p-12L, + -0x2.7027027027027027027027027028p-12L, + 0x3.72a3c5631fe46ae1d4e700dca8f2p-12L, + -0x7.daac36664f1f207daac36664f1f4p-12L, + 0x1.a41a41a41a41a41a41a41a41a41ap-8L, + -0x7.90a1b2c3d4e5f708192a3b4c5d7p-8L, + 0x2.dfd2c703c0cfff430edfd2c703cp-4L, + -0x1.6476701181f39edbdb9ce625987dp+0L, + 0xd.672219167002d3a7a9c886459cp+0L, + -0x9.cd9292e6660d55b3f712eb9e07c8p+4L, + 0x8.911a740da740da740da740da741p+8L, + -0x8.d0cc570e255bf59ff6eec24b49p+12L, + 0xa.8d1044d3708d1c219ee4fdc446ap+16L, + -0xe.8844d8a169abbc406169abbc406p+20L, + 0x1.6d29a0f6433b79890cede62433b8p+28L, + -0x2.88a233b3c8cddaba9809357125d8p+32L, + 0x5.0dde6f27500939a85c40939a85c4p+36L, + -0xb.4005bde03d4642a243581714af68p+40L, + 0x1.bc8cd6f8f1f755c78753cdb5d5c9p+48L, + -0x4.bbebb143bb94de5a0284fa7ec424p+52L, + 0xe.2e1337f5af0bed90b6b0a352d4fp+56L, + -0x2.e78250162b62405ad3e4bfe61b38p+64L, + 0xa.5f7eef9e71ac7c80326ab4cc8bfp+68L, + -0x2.83be0395e550213369924971b21ap+76L, + 0xa.8ebfe48da17dd999790760b0cep+80L, + }; + +#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 long double poly_coeff[] = + { + /* Interval [-2.125, -2] (polynomial degree 23). */ + -0x1.0b71c5c54d42eb6c17f30b7aa8f5p+0L, + -0xc.73a1dc05f34951602554c6d7506p-4L, + -0x1.ec841408528b51473e6c425ee5ffp-4L, + -0xe.37c9da26fc3c9a3c1844c8c7f1cp-4L, + -0x1.03cd87c519305703b021fa33f827p-4L, + -0xe.ae9ada65e09aa7f1c75216128f58p-4L, + 0x9.b11855a4864b5731cf85736015a8p-8L, + -0xe.f28c133e697a95c28607c9701dep-4L, + 0x2.6ec14a1c586a72a7cc33ee569d6ap-4L, + -0xf.57cab973e14464a262fc24723c38p-4L, + 0x4.5b0fc25f16e52997b2886bbae808p-4L, + -0xf.f50e59f1a9b56e76e988dac9ccf8p-4L, + 0x6.5f5eae15e9a93369e1d85146c6fcp-4L, + -0x1.0d2422daac459e33e0994325ed23p+0L, + 0x8.82000a0e7401fb1117a0e6606928p-4L, + -0x1.1f492f178a3f1b19f58a2ca68e55p+0L, + 0xa.cb545f949899a04c160b19389abp-4L, + -0x1.36165a1b155ba3db3d1b77caf498p+0L, + 0xd.44c5d5576f74302e5cf79e183eep-4L, + -0x1.51f22e0cdd33d3d481e326c02f3ep+0L, + 0xf.f73a349c08244ac389c007779bfp-4L, + -0x1.73317bf626156ba716747c4ca866p+0L, + 0x1.379c3c97b9bc71e1c1c4802dd657p+0L, + -0x1.a72a351c54f902d483052000f5dfp+0L, + /* Interval [-2.25, -2.125] (polynomial degree 24). */ + -0xf.2930890d7d675a80c36afb0fd5e8p-4L, + -0xc.a5cfde054eab5c6770daeca577f8p-4L, + 0x3.9c9e0fdebb07cdf89c61d41c9238p-4L, + -0x1.02a5ad35605fcf4af65a6dbacb84p+0L, + 0x9.6e9b1185bb48be9de1918e00a2e8p-4L, + -0x1.4d8332f3cfbfa116fd611e9ce90dp+0L, + 0x1.1c0c8cb4d9f4b1d490e1a41fae4dp+0L, + -0x1.c9a6f5ae9130cd0299e293a42714p+0L, + 0x1.d7e9307fd58a2ea997f29573a112p+0L, + -0x2.921cb3473d96178ca2a11d2a8d46p+0L, + 0x2.e8d59113b6f3409ff8db226e9988p+0L, + -0x3.cbab931625a1ae2b26756817f264p+0L, + 0x4.7d9f0f05d5296d18663ca003912p+0L, + -0x5.ade9cba12a14ea485667b7135bbp+0L, + 0x6.dc983a5da74fb48e767b7fec0a3p+0L, + -0x8.8d9ed454ae31d9e138dd8ee0d1a8p+0L, + 0xa.6fa099d4e7c202e0c0fd6ed8492p+0L, + -0xc.ebc552a8090a0f0115e92d4ebbc8p+0L, + 0xf.d695e4772c0d829b53fba9ca5568p+0L, + -0x1.38c32ae38e5e9eb79b2a4c5570a9p+4L, + 0x1.8035145646cfab49306d0999a51bp+4L, + -0x1.d930adbb03dd342a4c2a8c4e1af6p+4L, + 0x2.45c2edb1b4943ddb3686cd9c6524p+4L, + -0x2.e818ebbfafe2f916fa21abf7756p+4L, + 0x3.9804ce51d0fb9a430a711fd7307p+4L, + /* Interval [-2.375, -2.25] (polynomial degree 25). */ + -0xd.7d28d505d6181218a25f31d5e45p-4L, + -0xe.69649a3040985140cdf946829fap-4L, + 0xb.0d74a2827d053a8d44595012484p-4L, + -0x1.924b0922853617cac181afbc08ddp+0L, + 0x1.d49b12bccf0a568582e2d3c410f3p+0L, + -0x3.0898bb7d8c4093e636279c791244p+0L, + 0x4.207a6cac711cb53868e8a5057eep+0L, + -0x6.39ee63ea4fb1dcab0c9144bf3ddcp+0L, + 0x8.e2e2556a797b649bf3f53bd26718p+0L, + -0xd.0e83ac82552ef12af508589e7a8p+0L, + 0x1.2e4525e0ce6670563c6484a82b05p+4L, + -0x1.b8e350d6a8f2b222fa390a57c23dp+4L, + 0x2.805cd69b919087d8a80295892c2cp+4L, + -0x3.a42585424a1b7e64c71743ab014p+4L, + 0x5.4b4f409f98de49f7bfb03c05f984p+4L, + -0x7.b3c5827fbe934bc820d6832fb9fcp+4L, + 0xb.33b7b90cc96c425526e0d0866e7p+4L, + -0x1.04b77047ac4f59ee3775ca10df0dp+8L, + 0x1.7b366f5e94a34f41386eac086313p+8L, + -0x2.2797338429385c9849ca6355bfc2p+8L, + 0x3.225273cf92a27c9aac1b35511256p+8L, + -0x4.8f078aa48afe6cb3a4e89690f898p+8L, + 0x6.9f311d7b6654fc1d0b5195141d04p+8L, + -0x9.a0c297b6b4621619ca9bacc48ed8p+8L, + 0xe.ce1f06b6f90d92138232a76e4cap+8L, + -0x1.5b0e6806fa064daf011613e43b17p+12L, + /* Interval [-2.5, -2.375] (polynomial degree 27). */ + -0xb.74ea1bcfff94b2c01afba9daa7d8p-4L, + -0x1.2a82bd590c37538cab143308de4dp+0L, + 0x1.88020f828b966fec66b8649fd6fcp+0L, + -0x3.32279f040eb694970e9db24863dcp+0L, + 0x5.57ac82517767e68a721005853864p+0L, + -0x9.c2aedcfe22833de43834a0a6cc4p+0L, + 0x1.12c132f1f5577f99e1a0ed3538e1p+4L, + -0x1.ea94e26628a3de3597f7bb55a948p+4L, + 0x3.66b4ac4fa582f58b59f96b2f7c7p+4L, + -0x6.0cf746a9cf4cba8c39afcc73fc84p+4L, + 0xa.c102ef2c20d75a342197df7fedf8p+4L, + -0x1.31ebff06e8f14626782df58db3b6p+8L, + 0x2.1fd6f0c0e710994e059b9dbdb1fep+8L, + -0x3.c6d76040407f447f8b5074f07706p+8L, + 0x6.b6d18e0d8feb4c2ef5af6a40ed18p+8L, + -0xb.efaf542c529f91e34217f24ae6a8p+8L, + 0x1.53852d873210e7070f5d9eb2296p+12L, + -0x2.5b977c0ddc6d540717173ac29fc8p+12L, + 0x4.310d452ae05100eff1e02343a724p+12L, + -0x7.73a5d8f20c4f986a7dd1912b2968p+12L, + 0xd.3f5ea2484f3fca15eab1f4d1a218p+12L, + -0x1.78d18aac156d1d93a2ffe7e08d3fp+16L, + 0x2.9df49ca75e5b567f5ea3e47106cp+16L, + -0x4.a7149af8961a08aa7c3233b5bb94p+16L, + 0x8.3db10ffa742c707c25197d989798p+16L, + -0xe.a26d6dd023cadd02041a049ec368p+16L, + 0x1.c825d90514e7c57c7fa5316f947cp+20L, + -0x3.34bb81e5a0952df8ca1abdc6684cp+20L, + /* Interval [-2.625, -2.5] (polynomial degree 28). */ + -0x3.d10108c27ebafad533c20eac32bp-4L, + 0x1.cd557caff7d2b2085f41dbec5106p+0L, + 0x3.819b4856d399520dad9776ea2cacp+0L, + 0x6.8505cbad03dc34c5e42e8b12eb78p+0L, + 0xb.c1b2e653a9e38f82b399c94e7f08p+0L, + 0x1.50a53a38f148138105124df65419p+4L, + 0x2.57ae00cbe5232cbeeed34d89727ap+4L, + 0x4.2b156301b8604db85a601544bfp+4L, + 0x7.6989ed23ca3ca7579b3462592b5cp+4L, + 0xd.2dd2976557939517f831f5552cc8p+4L, + 0x1.76e1c3430eb860969bce40cd494p+8L, + 0x2.9a77bf5488742466db3a2c7c1ec6p+8L, + 0x4.a0d62ed7266e8eb36f725a8ebcep+8L, + 0x8.3a6184dd3021067df2f8b91e99c8p+8L, + 0xe.a0ade1538245bf55d39d7e436b1p+8L, + 0x1.a01359fae8617b5826dd74428e9p+12L, + 0x2.e3b0a32caae77251169acaca1ad4p+12L, + 0x5.2301257c81589f62b38fb5993ee8p+12L, + 0x9.21c9275db253d4e719b73b18cb9p+12L, + 0x1.03c104bc96141cda3f3fa4b112bcp+16L, + 0x1.cdc8ed65119196a08b0c78f1445p+16L, + 0x3.34f31d2eaacf34382cdb0073572ap+16L, + 0x5.b37628cadf12bf0000907d0ef294p+16L, + 0xa.22d8b332c0b1e6a616f425dfe5ap+16L, + 0x1.205b01444804c3ff922cd78b4c42p+20L, + 0x1.fe8f0cea9d1e0ff25be2470b4318p+20L, + 0x3.8872aebeb368399aee02b39340aep+20L, + 0x6.ebd560d351e84e26a4381f5b293cp+20L, + 0xc.c3644d094b0dae2fbcbf682cd428p+20L, + /* Interval [-2.75, -2.625] (polynomial degree 26). */ + -0x6.b5d252a56e8a75458a27ed1c2dd4p-4L, + 0x1.28d60383da3ac721aed3c5794da9p+0L, + 0x1.db6513ada8a66ea77d87d9a8827bp+0L, + 0x2.e217118f9d348a27f7506a707e6ep+0L, + 0x4.450112c5cbf725a0fb9802396c9p+0L, + 0x6.4af99151eae7810a75df2a0303c4p+0L, + 0x9.2db598b4a97a7f69aeef32aec758p+0L, + 0xd.62bef9c22471f5ee47ea1b9c0b5p+0L, + 0x1.379f294e412bd62328326d4222f9p+4L, + 0x1.c5827349d8865f1e8825c37c31c6p+4L, + 0x2.93a7e7a75b7568cc8cbe8c016c12p+4L, + 0x3.bf9bb882afe57edb383d41879d3ap+4L, + 0x5.73c737828cee095c43a5566731c8p+4L, + 0x7.ee4653493a7f81e0442062b3823cp+4L, + 0xb.891c6b83fc8b55bd973b5d962d6p+4L, + 0x1.0c775d7de3bf9b246c0208e0207ep+8L, + 0x1.867ee43ec4bd4f4fd56abc05110ap+8L, + 0x2.37fe9ba6695821e9822d8c8af0a6p+8L, + 0x3.3a2c667e37c942f182cd3223a936p+8L, + 0x4.b1b500eb59f3f782c7ccec88754p+8L, + 0x6.d3efd3b65b3d0d8488d30b79fa4cp+8L, + 0x9.ee8224e65bed5ced8b75eaec609p+8L, + 0xe.72416e510cca77d53fc615c1f3dp+8L, + 0x1.4fb538b0a2dfe567a8904b7e0445p+12L, + 0x1.e7f56a9266cf525a5b8cf4cb76cep+12L, + 0x2.f0365c983f68c597ee49d099cce8p+12L, + 0x4.53aa229e1b9f5b5e59625265951p+12L, + /* Interval [-2.875, -2.75] (polynomial degree 24). */ + -0x8.a41b1e4f36ff88dc820815607d68p-4L, + 0xc.da87d3b69dc0f2f9c6f368b8ca1p-4L, + 0x1.1474ad5c36158a7bea04fd2f98c6p+0L, + 0x1.761ecb90c555df6555b7dba955b6p+0L, + 0x1.d279bff9ae291caf6c4b4bcb3202p+0L, + 0x2.4e5d00559a6e2b9b5d7fe1f6689cp+0L, + 0x2.d57545a75cee8743ae2b17bc8d24p+0L, + 0x3.8514eee3aac88b89bec2307021bap+0L, + 0x4.5235e3b6e1891ffeb87fed9f8a24p+0L, + 0x5.562acdb10eef3c9a773b3e27a864p+0L, + 0x6.8ec8965c76efe03c26bff60b1194p+0L, + 0x8.15251aca144877af32658399f9b8p+0L, + 0x9.f08d56aba174d844138af782c0f8p+0L, + 0xc.3dbbeda2679e8a1346ccc3f6da88p+0L, + 0xf.0f5bfd5eacc26db308ffa0556fa8p+0L, + 0x1.28a6ccd84476fbc713d6bab49ac9p+4L, + 0x1.6d0a3ae2a3b1c8ff400641a3a21fp+4L, + 0x1.c15701b28637f87acfb6a91d33b5p+4L, + 0x2.28fbe0eccf472089b017651ca55ep+4L, + 0x2.a8a453004f6e8ffaacd1603bc3dp+4L, + 0x3.45ae4d9e1e7cd1a5dba0e4ec7f6cp+4L, + 0x4.065fbfacb7fad3e473cb577a61e8p+4L, + 0x4.f3d1473020927acac1944734a39p+4L, + 0x6.54bb091245815a36fb74e314dd18p+4L, + 0x7.d7f445129f7fb6c055e582d3f6ep+4L, + /* Interval [-3, -2.875] (polynomial degree 23). */ + -0xa.046d667e468f3e44dcae1afcc648p-4L, + 0x9.70b88dcc006c214d8d996fdf5ccp-4L, + 0xa.a8a39421c86d3ff24931a0929fp-4L, + 0xd.2f4d1363f324da2b357c8b6ec94p-4L, + 0xd.ca9aa1a3a5c00de11bf60499a97p-4L, + 0xf.cf09c31eeb52a45dfa7ebe3778dp-4L, + 0x1.04b133a39ed8a09691205660468bp+0L, + 0x1.22b547a06edda944fcb12fd9b5ecp+0L, + 0x1.2c57fce7db86a91df09602d344b3p+0L, + 0x1.4aade4894708f84795212fe257eep+0L, + 0x1.579c8b7b67ec4afed5b28c8bf787p+0L, + 0x1.776820e7fc80ae5284239733078ap+0L, + 0x1.883ab28c7301fde4ca6b8ec26ec8p+0L, + 0x1.aa2ef6e1ae52eb42c9ee83b206e3p+0L, + 0x1.bf4ad50f0a9a9311300cf0c51ee7p+0L, + 0x1.e40206e0e96b1da463814dde0d09p+0L, + 0x1.fdcbcffef3a21b29719c2bd9feb1p+0L, + 0x2.25e2e8948939c4d42cf108fae4bep+0L, + 0x2.44ce14d2b59c1c0e6bf2cfa81018p+0L, + 0x2.70ee80bbd0387162be4861c43622p+0L, + 0x2.954b64d2c2ebf3489b949c74476p+0L, + 0x2.c616e133a811c1c9446105208656p+0L, + 0x3.05a69dfe1a9ba1079f90fcf26bd4p+0L, + 0x3.410d2ad16a0506de29736e6aafdap+0L, + }; + +static const size_t poly_deg[] = + { + 23, + 24, + 25, + 27, + 28, + 26, + 24, + 23, + }; + +static const size_t poly_end[] = + { + 23, + 48, + 74, + 102, + 131, + 158, + 183, + 207, + }; + +/* Compute sin (pi * X) for -0.25 <= X <= 0.5. */ + +static long double +lg_sinpi (long double x) +{ + if (x <= 0.25L) + return __sinl (M_PIl * x); + else + return __cosl (M_PIl * (0.5L - x)); +} + +/* Compute cos (pi * X) for -0.25 <= X <= 0.5. */ + +static long double +lg_cospi (long double x) +{ + if (x <= 0.25L) + return __cosl (M_PIl * x); + else + return __sinl (M_PIl * (0.5L - x)); +} + +/* Compute cot (pi * X) for -0.25 <= X <= 0.5. */ + +static long double +lg_cotpi (long double x) +{ + return lg_cospi (x) / lg_sinpi (x); +} + +/* Compute lgamma of a negative argument -50 < X < -2, setting + *SIGNGAMP accordingly. */ + +long double +__lgamma_negl (long double x, int *signgamp) +{ + /* Determine the half-integer region X lies in, handle exact + integers and determine the sign of the result. */ + int i = __floorl (-2 * x); + if ((i & 1) == 0 && i == -2 * x) + return 1.0L / 0.0L; + long double xn = ((i & 1) == 0 ? -i / 2 : (-i - 1) / 2); + i -= 4; + *signgamp = ((i & 2) == 0 ? -1 : 1); + + SET_RESTORE_ROUNDL (FE_TONEAREST); + + /* Expand around the zero X0 = X0_HI + X0_LO. */ + long double x0_hi = lgamma_zeros[i][0], x0_lo = lgamma_zeros[i][1]; + long 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 = __floorl (-8 * x) - 16; + long double xm = (-33 - 2 * j) * 0.0625L; + long double x_adj = x - xm; + size_t deg = poly_deg[j]; + size_t end = poly_end[j]; + long double g = poly_coeff[end]; + for (size_t j = 1; j <= deg; j++) + g = g * x_adj + poly_coeff[end - j]; + return __log1pl (g * xdiff / (x - xn)); + } + + /* The result we want is log (sinpi (X0) / sinpi (X)) + + log (gamma (1 - X0) / gamma (1 - X)). */ + long double x_idiff = fabsl (xn - x), x0_idiff = fabsl (xn - x0_hi - x0_lo); + long double log_sinpi_ratio; + if (x0_idiff < x_idiff * 0.5L) + /* Use log not log1p to avoid inaccuracy from log1p of arguments + close to -1. */ + log_sinpi_ratio = __ieee754_logl (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. */ + long double x0diff2 = ((i & 1) == 0 ? xdiff : -xdiff) * 0.5L; + long double sx0d2 = lg_sinpi (x0diff2); + long double cx0d2 = lg_cospi (x0diff2); + log_sinpi_ratio = __log1pl (2 * sx0d2 + * (-sx0d2 + cx0d2 * lg_cotpi (x_idiff))); + } + + long double log_gamma_ratio; + long double y0 = 1 - x0_hi; + long double y0_eps = -x0_hi + (1 - y0) - x0_lo; + long double y = 1 - x; + long 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. */ + long double log_gamma_adj = 0; + if (i < 20) + { + int n_up = (21 - i) / 2; + long double ny0, ny0_eps, ny, ny_eps; + ny0 = y0 + n_up; + ny0_eps = y0 - (ny0 - n_up) + y0_eps; + y0 = ny0; + y0_eps = ny0_eps; + ny = y + n_up; + ny_eps = y - (ny - n_up) + y_eps; + y = ny; + y_eps = ny_eps; + long double prodm1 = __lgamma_productl (xdiff, y - n_up, y_eps, n_up); + log_gamma_adj = -__log1pl (prodm1); + } + long double log_gamma_high + = (xdiff * __log1pl ((y0 - e_hi - e_lo + y0_eps) / e_hi) + + (y - 0.5L + y_eps) * __log1pl (xdiff / y) + log_gamma_adj); + /* Compute the sum of (B_2k / 2k(2k-1))(Y0^-(2k-1) - Y^-(2k-1)). */ + long double y0r = 1 / y0, yr = 1 / y; + long double y0r2 = y0r * y0r, yr2 = yr * yr; + long double rdiff = -xdiff / (y * y0); + long double bterm[NCOEFF]; + long double dlast = rdiff, elast = rdiff * yr * (yr + y0r); + bterm[0] = dlast * lgamma_coeff[0]; + for (size_t j = 1; j < NCOEFF; j++) + { + long double dnext = dlast * y0r2 + elast; + long double enext = elast * yr2; + bterm[j] = dnext * lgamma_coeff[j]; + dlast = dnext; + elast = enext; + } + long 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/ldbl-128/lgamma_productl.c b/sysdeps/ieee754/ldbl-128/lgamma_productl.c new file mode 100644 index 0000000000..de67cbe665 --- /dev/null +++ b/sysdeps/ieee754/ldbl-128/lgamma_productl.c @@ -0,0 +1,82 @@ +/* Compute a product of 1 + (T/X), 1 + (T/(X+1)), .... + Copyright (C) 2015-2016 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <float.h> + +/* Calculate X * Y exactly and store the result in *HI + *LO. It is + given that the values are small enough that no overflow occurs and + large enough (or zero) that no underflow occurs. */ + +static void +mul_split (long double *hi, long double *lo, long double x, long double y) +{ +#ifdef __FP_FAST_FMAL + /* Fast built-in fused multiply-add. */ + *hi = x * y; + *lo = __builtin_fmal (x, y, -*hi); +#elif defined FP_FAST_FMAL + /* Fast library fused multiply-add, compiler before GCC 4.6. */ + *hi = x * y; + *lo = __fmal (x, y, -*hi); +#else + /* Apply Dekker's algorithm. */ + *hi = x * y; +# define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1) + long double x1 = x * C; + long double y1 = y * C; +# undef C + x1 = (x - x1) + x1; + y1 = (y - y1) + y1; + long double x2 = x - x1; + long double y2 = y - y1; + *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2; +#endif +} + +/* Compute the product of 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. */ + +long double +__lgamma_productl (long double t, long double x, long double x_eps, int n) +{ + long double ret = 0, ret_eps = 0; + for (int i = 0; i < n; i++) + { + long double xi = x + i; + long double quot = t / xi; + long double mhi, mlo; + mul_split (&mhi, &mlo, quot, xi); + long double quot_lo = (t - mhi - mlo) / xi - t * x_eps / (xi * xi); + /* We want (1 + RET + RET_EPS) * (1 + QUOT + QUOT_LO) - 1. */ + long double rhi, rlo; + mul_split (&rhi, &rlo, ret, quot); + long double rpq = ret + quot; + long double rpq_eps = (ret - rpq) + quot; + long double nret = rpq + rhi; + long 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/ldbl-128/mpn2ldbl.c b/sysdeps/ieee754/ldbl-128/mpn2ldbl.c index b4cab91211..f9e4a39f59 100644 --- a/sysdeps/ieee754/ldbl-128/mpn2ldbl.c +++ b/sysdeps/ieee754/ldbl-128/mpn2ldbl.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1995-2015 Free Software Foundation, Inc. +/* Copyright (C) 1995-2016 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 diff --git a/sysdeps/ieee754/ldbl-128/printf_fphex.c b/sysdeps/ieee754/ldbl-128/printf_fphex.c index 59c276b295..83f8f39ae8 100644 --- a/sysdeps/ieee754/ldbl-128/printf_fphex.c +++ b/sysdeps/ieee754/ldbl-128/printf_fphex.c @@ -1,6 +1,6 @@ /* Print floating point number in hexadecimal notation according to ISO C99. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 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 diff --git a/sysdeps/ieee754/ldbl-128/s_asinhl.c b/sysdeps/ieee754/ldbl-128/s_asinhl.c index 8d5721a73c..5f3b9f2c76 100644 --- a/sysdeps/ieee754/ldbl-128/s_asinhl.c +++ b/sysdeps/ieee754/ldbl-128/s_asinhl.c @@ -52,11 +52,7 @@ __asinhl (long double x) return x + x; /* x is inf or NaN */ if (ix < 0x3fc70000) { /* |x| < 2^ -56 */ - if (fabsl (x) < LDBL_MIN) - { - long double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); if (huge + x > one) return x; /* return x inexact except 0 */ } diff --git a/sysdeps/ieee754/ldbl-128/s_atanl.c b/sysdeps/ieee754/ldbl-128/s_atanl.c index 1367b6b15d..aaae701551 100644 --- a/sysdeps/ieee754/ldbl-128/s_atanl.c +++ b/sysdeps/ieee754/ldbl-128/s_atanl.c @@ -202,11 +202,7 @@ __atanl (long double x) if (k <= 0x3fc50000) /* |x| < 2**-58 */ { - if (fabsl (x) < LDBL_MIN) - { - long double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); /* Raise inexact. */ if (huge + x > 0.0) return x; diff --git a/sysdeps/ieee754/ldbl-128/s_erfl.c b/sysdeps/ieee754/ldbl-128/s_erfl.c index fa4609f136..dd275a7dd9 100644 --- a/sysdeps/ieee754/ldbl-128/s_erfl.c +++ b/sysdeps/ieee754/ldbl-128/s_erfl.c @@ -140,7 +140,6 @@ deval (long double x, const long double *p, int n) static const long double tiny = 1e-4931L, - half = 0.5L, one = 1.0L, two = 2.0L, /* 2/sqrt(pi) - 1 */ @@ -786,11 +785,7 @@ __erfl (long double x) { /* Avoid spurious underflow. */ long double ret = 0.0625 * (16.0 * x + (16.0 * efx) * x); - if (fabsl (ret) < LDBL_MIN) - { - long double force_underflow = ret * ret; - math_force_eval (force_underflow); - } + math_check_force_underflow (ret); return ret; } return x + efx * x; diff --git a/sysdeps/ieee754/ldbl-128/s_expm1l.c b/sysdeps/ieee754/ldbl-128/s_expm1l.c index 573d00be57..b1100a9d0e 100644 --- a/sysdeps/ieee754/ldbl-128/s_expm1l.c +++ b/sysdeps/ieee754/ldbl-128/s_expm1l.c @@ -84,8 +84,6 @@ static const long double C1 = 6.93145751953125E-1L, C2 = 1.428606820309417232121458176568075500134E-6L, -/* ln (2^16384 * (1 - 2^-113)) */ - maxlog = 1.1356523406294143949491931077970764891253E4L, /* ln 2^-114 */ minarg = -7.9018778583833765273564461846232128760607E1L, big = 1e4932L; @@ -110,14 +108,9 @@ __expm1l (long double x) } if (ix >= 0x7fff0000) { - /* Infinity. */ + /* Infinity (which must be negative infinity). */ if (((ix & 0xffff) | u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0) - { - if (sign) - return -1.0L; - else - return x; - } + return -1.0L; /* NaN. No invalid exception. */ return x; } @@ -126,13 +119,6 @@ __expm1l (long double x) if ((ix == 0) && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0) return x; - /* Overflow. */ - if (x > maxlog) - { - __set_errno (ERANGE); - return (big * big); - } - /* Minimum value. */ if (x < minarg) return (4.0/big - 1.0L); @@ -142,11 +128,7 @@ __expm1l (long double x) when the result does underflow. */ if (fabsl (x) < 0x1p-113L) { - if (fabsl (x) < LDBL_MIN) - { - long double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); return x; } diff --git a/sysdeps/ieee754/ldbl-128/s_finitel.c b/sysdeps/ieee754/ldbl-128/s_finitel.c index f862a448ed..ea8a9ba379 100644 --- a/sysdeps/ieee754/ldbl-128/s_finitel.c +++ b/sysdeps/ieee754/ldbl-128/s_finitel.c @@ -29,7 +29,7 @@ int __finitel(long double x) { int64_t hx; GET_LDOUBLE_MSW64(hx,x); - return (int)((u_int64_t)((hx&0x7fffffffffffffffLL) + return (int)((u_int64_t)((hx&0x7fff000000000000LL) -0x7fff000000000000LL)>>63); } hidden_def (__finitel) diff --git a/sysdeps/ieee754/ldbl-128/s_fma.c b/sysdeps/ieee754/ldbl-128/s_fma.c index 813f54a9da..cbe708f4b8 100644 --- a/sysdeps/ieee754/ldbl-128/s_fma.c +++ b/sysdeps/ieee754/ldbl-128/s_fma.c @@ -1,5 +1,5 @@ /* Compute x * y + z as ternary operation. - Copyright (C) 2010-2015 Free Software Foundation, Inc. + Copyright (C) 2010-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jakub@redhat.com>, 2010. diff --git a/sysdeps/ieee754/ldbl-128/s_fmal.c b/sysdeps/ieee754/ldbl-128/s_fmal.c index b13178ffe7..728949c916 100644 --- a/sysdeps/ieee754/ldbl-128/s_fmal.c +++ b/sysdeps/ieee754/ldbl-128/s_fmal.c @@ -1,5 +1,5 @@ /* Compute x * y + z as ternary operation. - Copyright (C) 2010-2015 Free Software Foundation, Inc. + Copyright (C) 2010-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jakub@redhat.com>, 2010. @@ -68,7 +68,7 @@ __fmal (long double x, long double y, long double z) if (u.ieee.exponent + v.ieee.exponent > 0x7fff + IEEE854_LONG_DOUBLE_BIAS) return x * y; - /* If x * y is less than 1/4 of LDBL_DENORM_MIN, neither the + /* If x * y is less than 1/4 of LDBL_TRUE_MIN, neither the result nor whether there is underflow depends on its exact value, only on its sign. */ if (u.ieee.exponent + v.ieee.exponent @@ -94,8 +94,8 @@ __fmal (long double x, long double y, long double z) && w.ieee.mantissa1 == 0 && w.ieee.mantissa0 == 0))) { - volatile long double force_underflow = x * y; - (void) force_underflow; + long double force_underflow = x * y; + math_force_eval (force_underflow); } return v.d * 0x1p-114L; } @@ -121,7 +121,7 @@ __fmal (long double x, long double y, long double z) very small, adjust them up to avoid spurious underflows, rather than down. */ if (u.ieee.exponent + v.ieee.exponent - <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG) + <= IEEE854_LONG_DOUBLE_BIAS + 2 * LDBL_MANT_DIG) { if (u.ieee.exponent > v.ieee.exponent) u.ieee.exponent += 2 * LDBL_MANT_DIG + 2; @@ -179,7 +179,10 @@ __fmal (long double x, long double y, long double z) /* Ensure correct sign of exact 0 + 0. */ if (__glibc_unlikely ((x == 0 || y == 0) && z == 0)) - return x * y + z; + { + x = math_opt_barrier (x); + return x * y + z; + } fenv_t env; feholdexcept (&env); diff --git a/sysdeps/ieee754/ldbl-128/s_fpclassifyl.c b/sysdeps/ieee754/ldbl-128/s_fpclassifyl.c index 2d41185ccd..c89686f8f9 100644 --- a/sysdeps/ieee754/ldbl-128/s_fpclassifyl.c +++ b/sysdeps/ieee754/ldbl-128/s_fpclassifyl.c @@ -1,5 +1,5 @@ /* Return classification value corresponding to argument. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and Jakub Jelinek <jj@ultra.linux.cz>, 1999. diff --git a/sysdeps/ieee754/ldbl-128/s_isinf_nsl.c b/sysdeps/ieee754/ldbl-128/s_isinf_nsl.c deleted file mode 100644 index 7d6cfb9a62..0000000000 --- a/sysdeps/ieee754/ldbl-128/s_isinf_nsl.c +++ /dev/null @@ -1,19 +0,0 @@ -/* - * Written by Ulrich Drepper <drepper@gmail.com> - */ - -/* - * __isinf_nsl(x) returns != 0 if x is ±inf, else 0; - * no branching! - */ - -#include <math.h> -#include <math_private.h> - -int -__isinf_nsl (long double x) -{ - int64_t hx,lx; - GET_LDOUBLE_WORDS64(hx,lx,x); - return !(lx | ((hx & 0x7fffffffffffffffLL) ^ 0x7fff000000000000LL)); -} diff --git a/sysdeps/ieee754/ldbl-128/s_issignalingl.c b/sysdeps/ieee754/ldbl-128/s_issignalingl.c index 0db27c4d76..c06f14dd12 100644 --- a/sysdeps/ieee754/ldbl-128/s_issignalingl.c +++ b/sysdeps/ieee754/ldbl-128/s_issignalingl.c @@ -1,5 +1,5 @@ /* Test for signaling NaN. - Copyright (C) 2013-2015 Free Software Foundation, Inc. + Copyright (C) 2013-2016 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 diff --git a/sysdeps/ieee754/ldbl-128/s_llrintl.c b/sysdeps/ieee754/ldbl-128/s_llrintl.c index 27196c95a4..84fc576ab6 100644 --- a/sysdeps/ieee754/ldbl-128/s_llrintl.c +++ b/sysdeps/ieee754/ldbl-128/s_llrintl.c @@ -1,6 +1,6 @@ /* Round argument to nearest integral value according to current rounding direction. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and Jakub Jelinek <jj@ultra.linux.cz>, 1999. @@ -19,9 +19,12 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include <fenv.h> +#include <limits.h> #include <math.h> #include <math_private.h> +#include <fix-fp-int-convert-overflow.h> static const long double two112[2] = { @@ -34,7 +37,7 @@ __llrintl (long double x) { int32_t j0; u_int64_t i0,i1; - volatile long double w; + long double w; long double t; long long int result; int sx; @@ -47,8 +50,21 @@ __llrintl (long double x) if (j0 < (int32_t) (8 * sizeof (long long int)) - 1) { - w = two112[sx] + x; - t = w - two112[sx]; +#if defined FE_INVALID || defined FE_INEXACT + /* X < LLONG_MAX + 1 implied by J0 < 63. */ + if (x > (long double) LLONG_MAX) + { + /* In the event of overflow we must raise the "invalid" + exception, but not "inexact". */ + t = __nearbyintl (x); + feraiseexcept (t == LLONG_MAX ? FE_INEXACT : FE_INVALID); + } + else +#endif + { + w = two112[sx] + x; + t = w - two112[sx]; + } GET_LDOUBLE_WORDS64 (i0, i1, t); j0 = ((i0 >> 48) & 0x7fff) - 0x3fff; i0 &= 0x0000ffffffffffffLL; @@ -63,8 +79,26 @@ __llrintl (long double x) } else { - /* The number is too large. It is left implementation defined - what happens. */ + /* The number is too large. Unless it rounds to LLONG_MIN, + FE_INVALID must be raised and the return value is + unspecified. */ +#if defined FE_INVALID || defined FE_INEXACT + if (x < (long double) LLONG_MIN + && x > (long double) LLONG_MIN - 1.0L) + { + /* If truncation produces LLONG_MIN, the cast will not raise + the exception, but may raise "inexact". */ + t = __nearbyintl (x); + feraiseexcept (t == LLONG_MIN ? FE_INEXACT : FE_INVALID); + return LLONG_MIN; + } + else if (FIX_LDBL_LLONG_CONVERT_OVERFLOW && x != (long double) LLONG_MIN) + { + feraiseexcept (FE_INVALID); + return sx == 0 ? LLONG_MAX : LLONG_MIN; + } + +#endif return (long long int) x; } diff --git a/sysdeps/ieee754/ldbl-128/s_llroundl.c b/sysdeps/ieee754/ldbl-128/s_llroundl.c index 4adc50eaa4..bfc81cc534 100644 --- a/sysdeps/ieee754/ldbl-128/s_llroundl.c +++ b/sysdeps/ieee754/ldbl-128/s_llroundl.c @@ -1,5 +1,5 @@ /* Round long double value to long long int. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and Jakub Jelinek <jj@ultra.linux.cz>, 1999. @@ -18,10 +18,12 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include <fenv.h> +#include <limits.h> #include <math.h> #include <math_private.h> - +#include <fix-fp-int-convert-overflow.h> long long int __llroundl (long double x) @@ -60,13 +62,37 @@ __llroundl (long double x) if (j0 == 48) result = (long long int) i0; else - result = ((long long int) i0 << (j0 - 48)) | (j >> (112 - j0)); + { + result = ((long long int) i0 << (j0 - 48)) | (j >> (112 - j0)); +#ifdef FE_INVALID + if (sign == 1 && result == LLONG_MIN) + /* Rounding brought the value out of range. */ + feraiseexcept (FE_INVALID); +#endif + } } } else { - /* The number is too large. It is left implementation defined - what happens. */ + /* The number is too large. Unless it rounds to LLONG_MIN, + FE_INVALID must be raised and the return value is + unspecified. */ +#ifdef FE_INVALID + if (FIX_LDBL_LLONG_CONVERT_OVERFLOW + && !(sign == -1 && x > (long double) LLONG_MIN - 0.5L)) + { + feraiseexcept (FE_INVALID); + return sign == 1 ? LLONG_MAX : LLONG_MIN; + } + else if (!FIX_LDBL_LLONG_CONVERT_OVERFLOW + && x <= (long double) LLONG_MIN - 0.5L) + { + /* If truncation produces LLONG_MIN, the cast will not raise + the exception, but may raise "inexact". */ + feraiseexcept (FE_INVALID); + return LLONG_MIN; + } +#endif return (long long int) x; } diff --git a/sysdeps/ieee754/ldbl-128/s_log1pl.c b/sysdeps/ieee754/ldbl-128/s_log1pl.c index ff759bc000..b348f41e55 100644 --- a/sysdeps/ieee754/ldbl-128/s_log1pl.c +++ b/sysdeps/ieee754/ldbl-128/s_log1pl.c @@ -117,7 +117,6 @@ static const long double C2 = 1.428606820309417232121458176568075500134E-6L; static const long double sqrth = 0.7071067811865475244008443621048490392848L; /* ln (2^16384 * (1 - 2^-113)) */ -static const long double maxlog = 1.1356523406294143949491931077970764891253E4L; static const long double zero = 0.0L; long double @@ -131,8 +130,8 @@ __log1pl (long double xm1) /* Test for NaN or infinity input. */ u.value = xm1; hx = u.parts32.w0; - if (hx >= 0x7fff0000) - return xm1; + if ((hx & 0x7fffffff) >= 0x7fff0000) + return xm1 + fabsl (xm1); /* log1p(+- 0) = +- 0. */ if (((hx & 0x7fffffff) == 0) @@ -141,11 +140,7 @@ __log1pl (long double xm1) if ((hx & 0x7fffffff) < 0x3f8e0000) { - if (fabsl (xm1) < LDBL_MIN) - { - long double force_underflow = xm1 * xm1; - math_force_eval (force_underflow); - } + math_check_force_underflow (xm1); if ((int) xm1 == 0) return xm1; } diff --git a/sysdeps/ieee754/ldbl-128/s_lrintl.c b/sysdeps/ieee754/ldbl-128/s_lrintl.c index dc678d585d..23f828f862 100644 --- a/sysdeps/ieee754/ldbl-128/s_lrintl.c +++ b/sysdeps/ieee754/ldbl-128/s_lrintl.c @@ -1,6 +1,6 @@ /* Round argument to nearest integral value according to current rounding direction. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and Jakub Jelinek <jj@ultra.linux.cz>, 1999. @@ -19,9 +19,12 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include <fenv.h> +#include <limits.h> #include <math.h> #include <math_private.h> +#include <fix-fp-int-convert-overflow.h> static const long double two112[2] = { @@ -34,7 +37,7 @@ __lrintl (long double x) { int32_t j0; u_int64_t i0,i1; - volatile long double w; + long double w; long double t; long int result; int sx; @@ -45,25 +48,53 @@ __lrintl (long double x) i0 &= 0x0000ffffffffffffLL; i0 |= 0x0001000000000000LL; - if (j0 < 48) + if (j0 < (int32_t) (8 * sizeof (long int)) - 1) { - w = two112[sx] + x; - t = w - two112[sx]; - GET_LDOUBLE_WORDS64 (i0, i1, t); - j0 = ((i0 >> 48) & 0x7fff) - 0x3fff; - i0 &= 0x0000ffffffffffffLL; - i0 |= 0x0001000000000000LL; + if (j0 < 48) + { +#if defined FE_INVALID || defined FE_INEXACT + /* X < LONG_MAX + 1 implied by J0 < 31. */ + if (sizeof (long int) == 4 + && x > (long double) LONG_MAX) + { + /* In the event of overflow we must raise the "invalid" + exception, but not "inexact". */ + t = __nearbyintl (x); + feraiseexcept (t == LONG_MAX ? FE_INEXACT : FE_INVALID); + } + else +#endif + { + w = two112[sx] + x; + t = w - two112[sx]; + } + GET_LDOUBLE_WORDS64 (i0, i1, t); + j0 = ((i0 >> 48) & 0x7fff) - 0x3fff; + i0 &= 0x0000ffffffffffffLL; + i0 |= 0x0001000000000000LL; - result = (j0 < 0 ? 0 : i0 >> (48 - j0)); - } - else if (j0 < (int32_t) (8 * sizeof (long int)) - 1) - { - if (j0 >= 112) + result = (j0 < 0 ? 0 : i0 >> (48 - j0)); + } + else if (j0 >= 112) result = ((long int) i0 << (j0 - 48)) | (i1 << (j0 - 112)); else { - w = two112[sx] + x; - t = w - two112[sx]; +#if defined FE_INVALID || defined FE_INEXACT + /* X < LONG_MAX + 1 implied by J0 < 63. */ + if (sizeof (long int) == 8 + && x > (long double) LONG_MAX) + { + /* In the event of overflow we must raise the "invalid" + exception, but not "inexact". */ + t = __nearbyintl (x); + feraiseexcept (t == LONG_MAX ? FE_INEXACT : FE_INVALID); + } + else +#endif + { + w = two112[sx] + x; + t = w - two112[sx]; + } GET_LDOUBLE_WORDS64 (i0, i1, t); j0 = ((i0 >> 48) & 0x7fff) - 0x3fff; i0 &= 0x0000ffffffffffffLL; @@ -77,8 +108,26 @@ __lrintl (long double x) } else { - /* The number is too large. It is left implementation defined - what happens. */ + /* The number is too large. Unless it rounds to LONG_MIN, + FE_INVALID must be raised and the return value is + unspecified. */ +#if defined FE_INVALID || defined FE_INEXACT + if (x < (long double) LONG_MIN + && x > (long double) LONG_MIN - 1.0L) + { + /* If truncation produces LONG_MIN, the cast will not raise + the exception, but may raise "inexact". */ + t = __nearbyintl (x); + feraiseexcept (t == LONG_MIN ? FE_INEXACT : FE_INVALID); + return LONG_MIN; + } + else if (FIX_LDBL_LONG_CONVERT_OVERFLOW && x != (long double) LONG_MIN) + { + feraiseexcept (FE_INVALID); + return sx == 0 ? LONG_MAX : LONG_MIN; + } + +#endif return (long int) x; } diff --git a/sysdeps/ieee754/ldbl-128/s_lroundl.c b/sysdeps/ieee754/ldbl-128/s_lroundl.c index 8421609676..f03262543f 100644 --- a/sysdeps/ieee754/ldbl-128/s_lroundl.c +++ b/sysdeps/ieee754/ldbl-128/s_lroundl.c @@ -1,5 +1,5 @@ /* Round long double value to long int. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and Jakub Jelinek <jj@ultra.linux.cz>, 1999. @@ -18,10 +18,12 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include <fenv.h> +#include <limits.h> #include <math.h> #include <math_private.h> - +#include <fix-fp-int-convert-overflow.h> long int __lroundl (long double x) @@ -37,19 +39,26 @@ __lroundl (long double x) i0 &= 0x0000ffffffffffffLL; i0 |= 0x0001000000000000LL; - if (j0 < 48) + if (j0 < (int32_t) (8 * sizeof (long int)) - 1) { - if (j0 < 0) - return j0 < -1 ? 0 : sign; - else + if (j0 < 48) { - i0 += 0x0000800000000000LL >> j0; - result = i0 >> (48 - j0); + if (j0 < 0) + return j0 < -1 ? 0 : sign; + else + { + i0 += 0x0000800000000000LL >> j0; + result = i0 >> (48 - j0); +#ifdef FE_INVALID + if (sizeof (long int) == 4 + && sign == 1 + && result == LONG_MIN) + /* Rounding brought the value out of range. */ + feraiseexcept (FE_INVALID); +#endif + } } - } - else if (j0 < (int32_t) (8 * sizeof (long int)) - 1) - { - if (j0 >= 112) + else if (j0 >= 112) result = ((long int) i0 << (j0 - 48)) | (i1 << (j0 - 112)); else { @@ -60,11 +69,39 @@ __lroundl (long double x) if (j0 == 48) result = (long int) i0; else - result = ((long int) i0 << (j0 - 48)) | (j >> (112 - j0)); + { + result = ((long int) i0 << (j0 - 48)) | (j >> (112 - j0)); +#ifdef FE_INVALID + if (sizeof (long int) == 8 + && sign == 1 + && result == LONG_MIN) + /* Rounding brought the value out of range. */ + feraiseexcept (FE_INVALID); +#endif + } } } else { + /* The number is too large. Unless it rounds to LONG_MIN, + FE_INVALID must be raised and the return value is + unspecified. */ +#ifdef FE_INVALID + if (FIX_LDBL_LONG_CONVERT_OVERFLOW + && !(sign == -1 && x > (long double) LONG_MIN - 0.5L)) + { + feraiseexcept (FE_INVALID); + return sign == 1 ? LONG_MAX : LONG_MIN; + } + else if (!FIX_LDBL_LONG_CONVERT_OVERFLOW + && x <= (long double) LONG_MIN - 0.5L) + { + /* If truncation produces LONG_MIN, the cast will not raise + the exception, but may raise "inexact". */ + feraiseexcept (FE_INVALID); + return LONG_MIN; + } +#endif /* The number is too large. It is left implementation defined what happens. */ return (long int) x; diff --git a/sysdeps/ieee754/ldbl-128/s_nextafterl.c b/sysdeps/ieee754/ldbl-128/s_nextafterl.c index d5eaa1cc91..4e9a2ce520 100644 --- a/sysdeps/ieee754/ldbl-128/s_nextafterl.c +++ b/sysdeps/ieee754/ldbl-128/s_nextafterl.c @@ -24,6 +24,7 @@ static char rcsid[] = "$NetBSD: $"; * Special cases: */ +#include <errno.h> #include <math.h> #include <math_private.h> @@ -70,10 +71,12 @@ long double __nextafterl(long double x, long double y) if(hy==0x7fff000000000000LL) { long double u = x + x; /* overflow */ math_force_eval (u); + __set_errno (ERANGE); } if(hy==0) { long double u = x*x; /* underflow */ math_force_eval (u); /* raise underflow flag */ + __set_errno (ERANGE); } SET_LDOUBLE_WORDS64(x,hx,lx); return x; diff --git a/sysdeps/ieee754/ldbl-128/s_nexttoward.c b/sysdeps/ieee754/ldbl-128/s_nexttoward.c index 2cd2e58391..4343fe83f8 100644 --- a/sysdeps/ieee754/ldbl-128/s_nexttoward.c +++ b/sysdeps/ieee754/ldbl-128/s_nexttoward.c @@ -25,6 +25,7 @@ static char rcsid[] = "$NetBSD: $"; * Special cases: */ +#include <errno.h> #include <math.h> #include <math_private.h> #include <float.h> @@ -73,15 +74,14 @@ double __nexttoward(double x, long double y) } hy = hx&0x7ff00000; if(hy>=0x7ff00000) { - x = x+x; /* overflow */ - if (FLT_EVAL_METHOD != 0 && FLT_EVAL_METHOD != 1) - /* Force conversion to double. */ - asm ("" : "+m"(x)); - return x; + double u = x+x; /* overflow */ + math_force_eval (u); + __set_errno (ERANGE); } if(hy<0x00100000) { double u = x*x; /* underflow */ math_force_eval (u); /* raise underflow flag */ + __set_errno (ERANGE); } INSERT_WORDS(x,hx,lx); return x; diff --git a/sysdeps/ieee754/ldbl-128/s_nexttowardf.c b/sysdeps/ieee754/ldbl-128/s_nexttowardf.c index ccd49d235d..8703359d4f 100644 --- a/sysdeps/ieee754/ldbl-128/s_nexttowardf.c +++ b/sysdeps/ieee754/ldbl-128/s_nexttowardf.c @@ -18,6 +18,7 @@ static char rcsid[] = "$NetBSD: $"; #endif +#include <errno.h> #include <math.h> #include <math_private.h> @@ -59,10 +60,15 @@ float __nexttowardf(float x, long double y) } } hy = hx&0x7f800000; - if(hy>=0x7f800000) return x+x; /* overflow */ + if(hy>=0x7f800000) { + float u = x+x; /* overflow */ + math_force_eval (u); + __set_errno (ERANGE); + } if(hy<0x00800000) { float u = x*x; /* underflow */ math_force_eval (u); /* raise underflow flag */ + __set_errno (ERANGE); } SET_FLOAT_WORD(x,hx); return x; diff --git a/sysdeps/ieee754/ldbl-128/s_remquol.c b/sysdeps/ieee754/ldbl-128/s_remquol.c index 82cdf1a589..7356f5f4ef 100644 --- a/sysdeps/ieee754/ldbl-128/s_remquol.c +++ b/sysdeps/ieee754/ldbl-128/s_remquol.c @@ -1,5 +1,5 @@ /* Compute remainder and a congruent to the quotient. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and Jakub Jelinek <jj@ultra.linux.cz>, 1999. diff --git a/sysdeps/ieee754/ldbl-128/s_roundl.c b/sysdeps/ieee754/ldbl-128/s_roundl.c index 3ff6ebe7cc..98b448e4f3 100644 --- a/sysdeps/ieee754/ldbl-128/s_roundl.c +++ b/sysdeps/ieee754/ldbl-128/s_roundl.c @@ -1,5 +1,5 @@ /* Round long double to integer away from zero. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and Jakub Jelinek <jj@ultra.linux.cz>, 1999. diff --git a/sysdeps/ieee754/ldbl-128/s_scalbnl.c b/sysdeps/ieee754/ldbl-128/s_scalbnl.c index 8916dcbfdd..a5cbd0d9d2 100644 --- a/sysdeps/ieee754/ldbl-128/s_scalbnl.c +++ b/sysdeps/ieee754/ldbl-128/s_scalbnl.c @@ -60,4 +60,3 @@ long double __scalbnl (long double x, int n) SET_LDOUBLE_MSW64(x,(hx&0x8000ffffffffffffULL)|(k<<48)); return x*twom114; } -weak_alias (__scalbnl, scalbnl) diff --git a/sysdeps/ieee754/ldbl-128/s_signbitl.c b/sysdeps/ieee754/ldbl-128/s_signbitl.c index acfe859a2e..ee5d77e27a 100644 --- a/sysdeps/ieee754/ldbl-128/s_signbitl.c +++ b/sysdeps/ieee754/ldbl-128/s_signbitl.c @@ -1,5 +1,5 @@ /* Return nonzero value if number is negative. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -19,13 +19,8 @@ #include <math.h> -#include <math_private.h> - int __signbitl (long double x) { - int64_t e; - - GET_LDOUBLE_MSW64 (e, x); - return e < 0; + return __builtin_signbitl (x); } diff --git a/sysdeps/ieee754/ldbl-128/s_sincosl.c b/sysdeps/ieee754/ldbl-128/s_sincosl.c index 6c4cbd7237..1abdb4419e 100644 --- a/sysdeps/ieee754/ldbl-128/s_sincosl.c +++ b/sysdeps/ieee754/ldbl-128/s_sincosl.c @@ -1,5 +1,5 @@ /* Compute sine and cosine of argument. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and Jakub Jelinek <jj@ultra.linux.cz>. @@ -39,7 +39,7 @@ __sincosl (long double x, long double *sinx, long double *cosx) { /* sin(Inf or NaN) is NaN */ *sinx = *cosx = x - x; - if (__isinf_nsl (x)) + if (isinf (x)) __set_errno (EDOM); } else diff --git a/sysdeps/ieee754/ldbl-128/s_tanhl.c b/sysdeps/ieee754/ldbl-128/s_tanhl.c index 129735b1b5..f7a1d20f79 100644 --- a/sysdeps/ieee754/ldbl-128/s_tanhl.c +++ b/sysdeps/ieee754/ldbl-128/s_tanhl.c @@ -41,6 +41,7 @@ * only tanhl(0)=0 is exact for finite argument. */ +#include <float.h> #include <math.h> #include <math_private.h> @@ -73,7 +74,10 @@ __tanhl (long double x) if (u.value == 0) return x; /* x == +- 0 */ if (ix < 0x3fc60000) /* |x| < 2^-57 */ - return x * (one + tiny); /* tanh(small) = small */ + { + math_check_force_underflow (x); + return x * (one + tiny); /* tanh(small) = small */ + } u.parts32.w0 = ix; /* Absolute value of x. */ if (ix >= 0x3fff0000) { /* |x| >= 1 */ diff --git a/sysdeps/ieee754/ldbl-128/s_truncl.c b/sysdeps/ieee754/ldbl-128/s_truncl.c index 14205c76a4..c71f0aba92 100644 --- a/sysdeps/ieee754/ldbl-128/s_truncl.c +++ b/sysdeps/ieee754/ldbl-128/s_truncl.c @@ -1,5 +1,5 @@ /* Truncate argument to nearest integral value not larger than the argument. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and Jakub Jelinek <jj@ultra.linux.cz>, 1999. diff --git a/sysdeps/ieee754/ldbl-128/strtod_nan_ldouble.h b/sysdeps/ieee754/ldbl-128/strtod_nan_ldouble.h new file mode 100644 index 0000000000..aed8953421 --- /dev/null +++ b/sysdeps/ieee754/ldbl-128/strtod_nan_ldouble.h @@ -0,0 +1,33 @@ +/* Convert string for NaN payload to corresponding NaN. For ldbl-128. + Copyright (C) 1997-2016 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#define FLOAT long double +#define SET_MANTISSA(flt, mant) \ + do \ + { \ + union ieee854_long_double u; \ + u.d = (flt); \ + u.ieee_nan.mantissa0 = 0; \ + u.ieee_nan.mantissa1 = 0; \ + u.ieee_nan.mantissa2 = (mant) >> 32; \ + u.ieee_nan.mantissa3 = (mant); \ + if ((u.ieee.mantissa0 | u.ieee.mantissa1 \ + | u.ieee.mantissa2 | u.ieee.mantissa3) != 0) \ + (flt) = u.d; \ + } \ + while (0) diff --git a/sysdeps/ieee754/ldbl-128/strtold_l.c b/sysdeps/ieee754/ldbl-128/strtold_l.c index d1ae57e304..12ccd93ee5 100644 --- a/sysdeps/ieee754/ldbl-128/strtold_l.c +++ b/sysdeps/ieee754/ldbl-128/strtold_l.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1999-2015 Free Software Foundation, Inc. +/* Copyright (C) 1999-2016 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 @@ -25,22 +25,13 @@ #ifdef USE_WIDE_CHAR # define STRTOF wcstold_l # define __STRTOF __wcstold_l +# define STRTOF_NAN __wcstold_nan #else # define STRTOF strtold_l # define __STRTOF __strtold_l +# define STRTOF_NAN __strtold_nan #endif #define MPN2FLOAT __mpn_construct_long_double #define FLOAT_HUGE_VAL HUGE_VALL -#define SET_MANTISSA(flt, mant) \ - do { union ieee854_long_double u; \ - u.d = (flt); \ - u.ieee_nan.mantissa0 = 0; \ - u.ieee_nan.mantissa1 = 0; \ - u.ieee_nan.mantissa2 = (mant) >> 32; \ - u.ieee_nan.mantissa3 = (mant); \ - if ((u.ieee.mantissa0 | u.ieee.mantissa1 \ - | u.ieee.mantissa2 | u.ieee.mantissa3) != 0) \ - (flt) = u.d; \ - } while (0) #include <strtod_l.c> diff --git a/sysdeps/ieee754/ldbl-128/t_expl.h b/sysdeps/ieee754/ldbl-128/t_expl.h index a9bf752e8f..a23d3d9344 100644 --- a/sysdeps/ieee754/ldbl-128/t_expl.h +++ b/sysdeps/ieee754/ldbl-128/t_expl.h @@ -1,5 +1,5 @@ /* Accurate table for expl(). - Copyright (C) 1999-2015 Free Software Foundation, Inc. + Copyright (C) 1999-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jj@ultra.linux.cz> diff --git a/sysdeps/ieee754/ldbl-128/t_sincosl.c b/sysdeps/ieee754/ldbl-128/t_sincosl.c index b6125f5612..de5fe0d9e9 100644 --- a/sysdeps/ieee754/ldbl-128/t_sincosl.c +++ b/sysdeps/ieee754/ldbl-128/t_sincosl.c @@ -1,5 +1,5 @@ /* Quad-precision floating point sine and cosine tables. - Copyright (C) 1999-2015 Free Software Foundation, Inc. + Copyright (C) 1999-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jj@ultra.linux.cz> diff --git a/sysdeps/ieee754/ldbl-128/x2y2m1l.c b/sysdeps/ieee754/ldbl-128/x2y2m1l.c index 11757c6af8..733742da04 100644 --- a/sysdeps/ieee754/ldbl-128/x2y2m1l.c +++ b/sysdeps/ieee754/ldbl-128/x2y2m1l.c @@ -1,5 +1,5 @@ /* Compute x^2 + y^2 - 1, without large cancellation error. - Copyright (C) 2012-2015 Free Software Foundation, Inc. + Copyright (C) 2012-2016 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 @@ -80,32 +80,26 @@ compare (const void *p, const void *q) } /* Return X^2 + Y^2 - 1, computed without large cancellation error. - It is given that 1 > X >= Y >= epsilon / 2, and that either X >= - 0.75 or Y >= 0.5. */ + It is given that 1 > X >= Y >= epsilon / 2, and that X^2 + Y^2 >= + 0.5. */ long double __x2y2m1l (long double x, long double y) { - long double vals[4]; + long double vals[5]; SET_RESTORE_ROUNDL (FE_TONEAREST); mul_split (&vals[1], &vals[0], x, x); mul_split (&vals[3], &vals[2], y, y); - if (x >= 0.75L) - vals[1] -= 1.0L; - else - { - vals[1] -= 0.5L; - vals[3] -= 0.5L; - } - qsort (vals, 4, sizeof (long double), compare); + vals[4] = -1.0L; + qsort (vals, 5, sizeof (long double), compare); /* Add up the values so that each element of VALS has absolute value at most equal to the last set bit of the next nonzero element. */ - for (size_t i = 0; i <= 2; i++) + for (size_t i = 0; i <= 3; i++) { add_split (&vals[i + 1], &vals[i], vals[i + 1], vals[i]); - qsort (vals + i + 1, 3 - i, sizeof (long double), compare); + qsort (vals + i + 1, 4 - i, sizeof (long double), compare); } /* Now any error from this addition will be small. */ - return vals[3] + vals[2] + vals[1] + vals[0]; + return vals[4] + vals[3] + vals[2] + vals[1] + vals[0]; } |