diff options
Diffstat (limited to 'sysdeps/ieee754/ldbl-128ibm')
62 files changed, 839 insertions, 203 deletions
diff --git a/sysdeps/ieee754/ldbl-128ibm/Makefile b/sysdeps/ieee754/ldbl-128ibm/Makefile index 8e1e2e6ae5..5591814824 100644 --- a/sysdeps/ieee754/ldbl-128ibm/Makefile +++ b/sysdeps/ieee754/ldbl-128ibm/Makefile @@ -3,3 +3,8 @@ # when -mlong-double-64 is not used). long-double-fcts = yes sysdep-CFLAGS += -mlong-double-128 + +ifeq ($(subdir),stdlib) +tests += tst-strtold-ldbl-128ibm +$(objpfx)tst-strtold-ldbl-128ibm: $(libm) +endif diff --git a/sysdeps/ieee754/ldbl-128ibm/e_asinl.c b/sysdeps/ieee754/ldbl-128ibm/e_asinl.c index 57c3ac09c8..6ed5e8d68d 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_asinl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_asinl.c @@ -147,11 +147,7 @@ __ieee754_asinl (long double x) { if (a < 6.938893903907228e-18L) /* |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-128ibm/e_atanhl.c b/sysdeps/ieee754/ldbl-128ibm/e_atanhl.c index bcd1fce4ed..b576f42030 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_atanhl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_atanhl.c @@ -57,11 +57,7 @@ __ieee754_atanhl(long double x) } if(ix<0x3c70000000000000LL&&(huge+x)>zero) /* x<2**-56 */ { - if (fabsl (x) < LDBL_MIN) - { - long double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); return x; } x = fabsl (x); diff --git a/sysdeps/ieee754/ldbl-128ibm/e_exp10l.c b/sysdeps/ieee754/ldbl-128ibm/e_exp10l.c index deefe7f54f..5699b8e53d 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_exp10l.c +++ b/sysdeps/ieee754/ldbl-128ibm/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 @@ -19,8 +19,8 @@ #include <math_private.h> #include <float.h> -static const long double log10_high = 0x2.4d763776aaa2cp0L; -static const long double log10_low = -0xf.a456a4a751f4b3d75c75c04c18p-56L; +static const long double log10_high = 0x2.4d763776aaap+0L; +static const long double log10_low = 0x2.b05ba95b58ae0b4c28a38a3fb4p-48L; long double __ieee754_exp10l (long double arg) diff --git a/sysdeps/ieee754/ldbl-128ibm/e_expl.c b/sysdeps/ieee754/ldbl-128ibm/e_expl.c index 15ccc454e1..ca3cbb53af 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_expl.c +++ b/sysdeps/ieee754/ldbl-128ibm/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 diff --git a/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c b/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c index 48098c18f6..8dbb131f93 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c +++ b/sysdeps/ieee754/ldbl-128ibm/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-128ibm/e_hypotl.c b/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c index 3b07a47b40..c68dac03b0 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c @@ -125,7 +125,11 @@ __ieee754_hypotl(long double x, long double y) w = __ieee754_sqrtl(a1*b1-(w*(-w)-(a1*b2+a2*b))); } if(k!=0) - return w*kld; + { + w *= kld; + math_check_force_underflow_nonneg (w); + return w; + } else return w; } diff --git a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c b/sysdeps/ieee754/ldbl-128ibm/e_jnl.c index 7594196a1f..4a8ccb044e 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c +++ b/sysdeps/ieee754/ldbl-128ibm/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-128ibm/e_log2l.c b/sysdeps/ieee754/ldbl-128ibm/e_log2l.c index 442ad97254..e39eaba72a 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_log2l.c +++ b/sysdeps/ieee754/ldbl-128ibm/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-128ibm/e_logl.c b/sysdeps/ieee754/ldbl-128ibm/e_logl.c index 58d6bc6972..14acfc2db7 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_logl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_logl.c @@ -219,6 +219,8 @@ __ieee754_logl(long double x) /* On this interval the table is not used due to cancellation error. */ if ((x <= 1.0078125L) && (x >= 0.9921875L)) { + if (x == 1.0L) + return 0.0L; z = x - 1.0L; k = 64; t = 1.0L; @@ -268,7 +270,7 @@ __ieee754_logl(long double x) /* Series expansion of log(1+z). */ w = z * z; /* Avoid spurious underflows. */ - if (__glibc_unlikely(w <= ldbl_epsilon)) + if (__glibc_unlikely (fabsl (z) <= ldbl_epsilon)) y = 0.0L; else { diff --git a/sysdeps/ieee754/ldbl-128ibm/e_powl.c b/sysdeps/ieee754/ldbl-128ibm/e_powl.c index c942f2f249..90340e890e 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_powl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_powl.c @@ -409,6 +409,7 @@ __ieee754_powl (long double x, long double y) r = (z * t1) / (t1 - two) - (w + z * w); z = one - (r - z); z = __scalbnl (z, n); + math_check_force_underflow_nonneg (z); return s * z; } strong_alias (__ieee754_powl, __powl_finite) diff --git a/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c b/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c index 9c6ebbdbff..cc2b235534 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c +++ b/sysdeps/ieee754/ldbl-128ibm/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-128ibm/e_sinhl.c b/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c index 1790bef87e..67d9d24ce7 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c @@ -28,6 +28,7 @@ * only sinh(0)=0 is exact for finite x. */ +#include <float.h> #include <math.h> #include <math_private.h> @@ -52,8 +53,10 @@ __ieee754_sinhl(long double x) if (jx<0) h = -h; /* |x| in [0,40], return sign(x)*0.5*(E+E/(E+1))) */ if (ix < 0x4044000000000000LL) { /* |x|<40 */ - if (ix<0x3e20000000000000LL) /* |x|<2**-29 */ + if (ix<0x3c90000000000000LL) { /* |x|<2**-54 */ + math_check_force_underflow (x); if(shuge+x>one) return x;/* sinhl(tiny) = tiny with inexact */ + } t = __expm1l(fabsl(x)); if(ix<0x3ff0000000000000LL) return h*(2.0*t-t*t/(t+one)); w = t/(t+one); @@ -64,7 +67,7 @@ __ieee754_sinhl(long double x) if (ix < 0x40862e42fefa39efLL) return h*__ieee754_expl(fabsl(x)); /* |x| in [log(maxdouble), overflowthresold] */ - if (ix <= 0x408633ce8fb9f87dLL) { + if (ix <= 0x408633ce8fb9f87eLL) { w = __ieee754_expl(0.5*fabsl(x)); t = h*w; return t*w; diff --git a/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c b/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c index ac1b22175e..96845fe5f8 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2015 Free Software Foundation, Inc. + * Copyright (C) 2001-2016 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/ldbl-128ibm/gamma_productl.c b/sysdeps/ieee754/ldbl-128ibm/gamma_productl.c index 905a78025a..b631f90a44 100644 --- a/sysdeps/ieee754/ldbl-128ibm/gamma_productl.c +++ b/sysdeps/ieee754/ldbl-128ibm/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-128ibm/ieee754.h b/sysdeps/ieee754/ldbl-128ibm/ieee754.h index d11e040469..c07a6def4b 100644 --- a/sysdeps/ieee754/ldbl-128ibm/ieee754.h +++ b/sysdeps/ieee754/ldbl-128ibm/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-128ibm/k_cosl.c b/sysdeps/ieee754/ldbl-128ibm/k_cosl.c index e812786bcc..2a3189ad36 100644 --- a/sysdeps/ieee754/ldbl-128ibm/k_cosl.c +++ b/sysdeps/ieee754/ldbl-128ibm/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-128ibm/k_sincosl.c b/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c index 0a76e1c7e7..4e43c8622e 100644 --- a/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c +++ b/sysdeps/ieee754/ldbl-128ibm/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> @@ -114,11 +114,7 @@ __kernel_sincosl(long double x, long double y, long double *sinx, long double *c polynomial of degree 16(17). */ if (tix < 0x3c600000) /* |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-128ibm/k_sinl.c b/sysdeps/ieee754/ldbl-128ibm/k_sinl.c index 2050cd2b49..44da02b0f3 100644 --- a/sysdeps/ieee754/ldbl-128ibm/k_sinl.c +++ b/sysdeps/ieee754/ldbl-128ibm/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> @@ -95,11 +95,7 @@ __kernel_sinl(long double x, long double y, int iy) polynomial of degree 17. */ if (tix < 0x3c600000) /* |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-128ibm/k_tanl.c b/sysdeps/ieee754/ldbl-128ibm/k_tanl.c index 7f1caeebdf..3c1bf32af9 100644 --- a/sysdeps/ieee754/ldbl-128ibm/k_tanl.c +++ b/sysdeps/ieee754/ldbl-128ibm/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 | lx | (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 >= 0x3fe59420) /* |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-128ibm/ldbl2mpn.c b/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c index e9b5803de3..4f550ef47c 100644 --- a/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c +++ b/sysdeps/ieee754/ldbl-128ibm/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-128ibm/lgamma_negl.c b/sysdeps/ieee754/ldbl-128ibm/lgamma_negl.c new file mode 100644 index 0000000000..c2a5cd29b6 --- /dev/null +++ b/sysdeps/ieee754/ldbl-128ibm/lgamma_negl.c @@ -0,0 +1,532 @@ +/* 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.74ff92c01f0d82abec9f315f1ap+0L, -0x7.12c334804d9a79cb5d46094d46p-112L }, + { -0x2.bf6821437b20197995a4b4641fp+0L, 0x5.140b4ff4b7d6069e1bd7acc196p-108L }, + { -0x3.24c1b793cb35efb8be699ad3dap+0L, 0x4.59abab3480539f1c0e926287cp-108L }, + { -0x3.f48e2a8f85fca170d456129123p+0L, -0x6.cc320a4887d1cb4c711828a75ep-108L }, + { -0x4.0a139e16656030c39f0b0de182p+0L, 0xe.d53e84029416e1242006b2b3dp-108L }, + { -0x4.fdd5de9bbabf3510d0aa407698p+0L, -0x8.501d7d78125286f78d1e501f14p-108L }, + { -0x5.021a95fc2db6432a4c56e5953ap+0L, 0xb.2133950fbcf2b01a8b9058dcccp-108L }, + { -0x5.ffa4bd647d0357dd4ed62cbd32p+0L, 0x1.2071c071a2145d2982428f2269p-108L }, + { -0x6.005ac9625f233b607c2d96d164p+0L, 0x7.a347953a96cbf30e1a0db20856p-108L }, + { -0x6.fff2fddae1bbff3d626b65c24p+0L, 0x2.de0bfcff5c457ebcf4d3ad9674p-108L }, + { -0x7.000cff7b7f87adf4482dcdb988p+0L, 0x7.d54d99e35a74d6407b80292df2p-108L }, + { -0x7.fffe5fe05673c3ca9e82b522bp+0L, -0xc.a9d2e8837cd1f14bd3d05002e4p-108L }, + { -0x8.0001a01459fc9f60cb3cec1cecp+0L, -0x8.576677ca538d88084310983b8p-108L }, + { -0x8.ffffd1c425e80ffc864e957494p+0L, 0x1.a6181dfdef1807e3087e4bb163p-104L }, + { -0x9.00002e3bb47d86d6d843fedc34p+0L, -0x1.1deb7ad09ec5e9d6e8ae2d548bp-104L }, + { -0x9.fffffb606bdfdcd062ae77a504p+0L, -0x1.47c69d2eb6f33d170fce38ff818p-104L }, + { -0xa.0000049f93bb9927b45d95e154p+0L, -0x4.1e03086db9146a9287bd4f2172p-108L }, + { -0xa.ffffff9466e9f1b36dacd2adbcp+0L, -0x1.18d05a4e458062f3f95345a4dap-104L }, + { -0xb.0000006b9915315d965a6ffea4p+0L, -0xe.4bea39000dcc1848023c5f6bdcp-112L }, + { -0xb.fffffff7089387387de41acc3cp+0L, -0x1.3c978bd839c8c428b5efcf91ef8p-104L }, + { -0xc.00000008f76c7731567c0f025p+0L, -0xf.387920df5675833859190eb128p-108L }, + { -0xc.ffffffff4f6dcf617f97a5ffc8p+0L, 0xa.82ab72d76f32eaee2d1a42ed5p-108L }, + { -0xd.00000000b092309c06683dd1b8p+0L, -0x1.03e3700857a15c19ac5a611de98p-104L }, + { -0xd.fffffffff36345ab9e184a3e08p+0L, -0x1.d1176dc48e47f62d917973dd45p-104L }, + { -0xe.000000000c9cba545e94e75ec4p+0L, -0x1.718f753e2501e757a17cf2ecbfp-104L }, + { -0xe.ffffffffff28c060c6604ef304p+0L, 0x8.e0762c8ca8361c23e8393919c4p-108L }, + { -0xf.0000000000d73f9f399bd0e42p+0L, -0xf.85e9ee31b0b890744fc0e3fbcp-108L }, + { -0xf.fffffffffff28c060c6621f514p+0L, 0x1.18d1b2eec9d960bd9adc5be5f6p-104L }, + { -0x1.000000000000d73f9f399da1428p+4L, 0x3.406c46e0e88305d2800f0e414cp-104L }, + { -0x1.0ffffffffffff3569c47e7a93ep+4L, -0x1.c46a08a2e008a998ebabb8087fp-104L }, + { -0x1.1000000000000ca963b81856888p+4L, -0x7.6ca5a3a64ec15db0a95caf2cap-108L }, + { -0x1.1fffffffffffff4bec3ce23413p+4L, -0x2.d08b2b726187c841cb92cd5222p-104L }, + { -0x1.20000000000000b413c31dcbec8p+4L, -0x2.4c3b2ffacbb4932f18dceedfd7p-104L }, + { -0x1.2ffffffffffffff685b25cbf5f8p+4L, 0x2.ba3126cd1c7b7a0822d694705cp-104L }, + { -0x1.30000000000000097a4da340a08p+4L, -0x2.b81b7b1f1f001c72bf914141efp-104L }, + { -0x1.3fffffffffffffff86af516ff8p+4L, 0x8.9429818df2a87abafd48248a2p-108L }, + { -0x1.40000000000000007950ae9008p+4L, -0x8.9413ccc8a353fda263f8ce973cp-108L }, + { -0x1.4ffffffffffffffffa391c4249p+4L, 0x3.d5c63022b62b5484ba346524dbp-104L }, + { -0x1.500000000000000005c6e3bdb7p+4L, -0x3.d5c62f55ed5322b2685c5e9a52p-104L }, + { -0x1.5fffffffffffffffffbcc71a49p+4L, -0x2.01eb5aeb96c74d7ad25e060529p-104L }, + { -0x1.6000000000000000004338e5b7p+4L, 0x2.01eb5aec04b2f2eb663e4e3d8ap-104L }, + { -0x1.6ffffffffffffffffffd13c97d8p+4L, -0x1.d38fcc4d08d6fe5aa56ab04308p-104L }, + { -0x1.70000000000000000002ec36828p+4L, 0x1.d38fcc4d090cee2f5d0b69a99cp-104L }, + { -0x1.7fffffffffffffffffffe0d31p+4L, 0x1.972f577cca4b4c8cb1dc14001bp-104L }, + { -0x1.800000000000000000001f2cfp+4L, -0x1.972f577cca4b3442e35f0040b38p-104L }, + { -0x1.8ffffffffffffffffffffec0c3p+4L, -0x3.22e9a0572b1bb5b95f346a92d6p-104L }, + { -0x1.90000000000000000000013f3dp+4L, 0x3.22e9a0572b1bb5c371ddb35617p-104L }, + { -0x1.9ffffffffffffffffffffff3b88p+4L, -0x3.d01cad8d32e386fd783e97296dp-104L }, + { -0x1.a0000000000000000000000c478p+4L, 0x3.d01cad8d32e386fd7c1ab8c1fep-104L }, + { -0x1.afffffffffffffffffffffff8b8p+4L, -0x1.538f48cc5737d5979c39db806c8p-104L }, + { -0x1.b00000000000000000000000748p+4L, 0x1.538f48cc5737d5979c3b3a6bdap-104L }, + { -0x1.bffffffffffffffffffffffffcp+4L, 0x2.862898d42174dcf171470d8c8cp-104L }, + { -0x1.c0000000000000000000000004p+4L, -0x2.862898d42174dcf171470d18bap-104L }, + { -0x1.dp+4L, 0x2.4b3f31686b15af57c61ceecdf4p-104L }, + { -0x1.dp+4L, -0x2.4b3f31686b15af57c61ceecdd1p-104L }, + { -0x1.ep+4L, 0x1.3932c5047d60e60caded4c298ap-108L }, + { -0x1.ep+4L, -0x1.3932c5047d60e60caded4c29898p-108L }, + { -0x1.fp+4L, 0xa.1a6973c1fade2170f7237d36p-116L }, + { -0x1.fp+4L, -0xa.1a6973c1fade2170f7237d36p-116L }, + { -0x2p+4L, 0x5.0d34b9e0fd6f10b87b91be9bp-120L }, + { -0x2p+4L, -0x5.0d34b9e0fd6f10b87b91be9bp-120L }, + { -0x2.1p+4L, 0x2.73024a9ba1aa36a7059bff52e8p-124L }, + { -0x2.1p+4L, -0x2.73024a9ba1aa36a7059bff52e8p-124L }, + { -0x2.2p+4L, 0x1.2710231c0fd7a13f8a2b4af9d68p-128L }, + { -0x2.2p+4L, -0x1.2710231c0fd7a13f8a2b4af9d68p-128L }, + { -0x2.3p+4L, 0x8.6e2ce38b6c8f9419e3fad3f03p-136L }, + { -0x2.3p+4L, -0x8.6e2ce38b6c8f9419e3fad3f03p-136L }, + { -0x2.4p+4L, 0x3.bf30652185952560d71a254e4fp-140L }, + { -0x2.4p+4L, -0x3.bf30652185952560d71a254e4fp-140L }, + { -0x2.5p+4L, 0x1.9ec8d1c94e85af4c78b15c3d8ap-144L }, + { -0x2.5p+4L, -0x1.9ec8d1c94e85af4c78b15c3d8ap-144L }, + { -0x2.6p+4L, 0xa.ea565ce061d57489e9b8527628p-152L }, + { -0x2.6p+4L, -0xa.ea565ce061d57489e9b8527628p-152L }, + { -0x2.7p+4L, 0x4.7a6512692eb37804111dabad3p-156L }, + { -0x2.7p+4L, -0x4.7a6512692eb37804111dabad3p-156L }, + { -0x2.8p+4L, 0x1.ca8ed42a12ae3001a07244abadp-160L }, + { -0x2.8p+4L, -0x1.ca8ed42a12ae3001a07244abadp-160L }, + { -0x2.9p+4L, 0xb.2f30e1ce812063f12e7e8d8d98p-168L }, + { -0x2.9p+4L, -0xb.2f30e1ce812063f12e7e8d8d98p-168L }, + { -0x2.ap+4L, 0x4.42bd49d4c37a0db136489772e4p-172L }, + { -0x2.ap+4L, -0x4.42bd49d4c37a0db136489772e4p-172L }, + { -0x2.bp+4L, 0x1.95db45257e5122dcbae56def37p-176L }, + { -0x2.bp+4L, -0x1.95db45257e5122dcbae56def37p-176L }, + { -0x2.cp+4L, 0x9.3958d81ff63527ecf993f3fb7p-184L }, + { -0x2.cp+4L, -0x9.3958d81ff63527ecf993f3fb7p-184L }, + { -0x2.dp+4L, 0x3.47970e4440c8f1c058bd238c99p-188L }, + { -0x2.dp+4L, -0x3.47970e4440c8f1c058bd238c99p-188L }, + { -0x2.ep+4L, 0x1.240804f65951062ca46e4f25c6p-192L }, + { -0x2.ep+4L, -0x1.240804f65951062ca46e4f25c6p-192L }, + { -0x2.fp+4L, 0x6.36a382849fae6de2d15362d8a4p-200L }, + { -0x2.fp+4L, -0x6.36a382849fae6de2d15362d8a4p-200L }, + { -0x3p+4L, 0x2.123680d6dfe4cf4b9b1bcb9d8cp-204L }, + }; + +static const long double e_hi = 0x2.b7e151628aed2a6abf7158809dp+0L; +static const long double e_lo = -0xb.0c389d18e9f0c74b25a9587b28p-112L; + +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's + approximation to lgamma function. */ + +static const long double lgamma_coeff[] = + { + 0x1.555555555555555555555555558p-4L, + -0xb.60b60b60b60b60b60b60b60b6p-12L, + 0x3.4034034034034034034034034p-12L, + -0x2.7027027027027027027027027p-12L, + 0x3.72a3c5631fe46ae1d4e700dca9p-12L, + -0x7.daac36664f1f207daac36664f2p-12L, + 0x1.a41a41a41a41a41a41a41a41a4p-8L, + -0x7.90a1b2c3d4e5f708192a3b4c5ep-8L, + 0x2.dfd2c703c0cfff430edfd2c704p-4L, + -0x1.6476701181f39edbdb9ce625988p+0L, + 0xd.672219167002d3a7a9c886459cp+0L, + -0x9.cd9292e6660d55b3f712eb9e08p+4L, + 0x8.911a740da740da740da740da74p+8L, + -0x8.d0cc570e255bf59ff6eec24b48p+12L, + 0xa.8d1044d3708d1c219ee4fdc448p+16L, + -0xe.8844d8a169abbc406169abbc4p+20L, + 0x1.6d29a0f6433b79890cede624338p+28L, + -0x2.88a233b3c8cddaba9809357126p+32L, + 0x5.0dde6f27500939a85c40939a86p+36L, + -0xb.4005bde03d4642a243581714bp+40L, + 0x1.bc8cd6f8f1f755c78753cdb5d6p+48L, + -0x4.bbebb143bb94de5a0284fa7ec4p+52L, + 0xe.2e1337f5af0bed90b6b0a352d4p+56L, + -0x2.e78250162b62405ad3e4bfe61bp+64L, + 0xa.5f7eef9e71ac7c80326ab4cc8cp+68L, + -0x2.83be0395e550213369924971b2p+76L, + }; + +#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 21). */ + -0x1.0b71c5c54d42eb6c17f30b7aa9p+0L, + -0xc.73a1dc05f34951602554c6d76cp-4L, + -0x1.ec841408528b51473e6c42f1c58p-4L, + -0xe.37c9da26fc3c9a3c1844c04b84p-4L, + -0x1.03cd87c519305703b00b046ce4p-4L, + -0xe.ae9ada65e09aa7f1c817c91048p-4L, + 0x9.b11855a4864b571b6a4f571c88p-8L, + -0xe.f28c133e697a95ba2dabb97584p-4L, + 0x2.6ec14a1c586a7ddb6c4be90fe1p-4L, + -0xf.57cab973e14496f0900851c0d4p-4L, + 0x4.5b0fc25f16b0df37175495c70cp-4L, + -0xf.f50e59f1a8fb8c402091e3cd3cp-4L, + 0x6.5f5eae1681d1e50e575c3d4d36p-4L, + -0x1.0d2422dac7ea8a52db6bf0d14fp+0L, + 0x8.820008f221eae5a36e15913bacp-4L, + -0x1.1f492eec53b9481ea23a7e944ep+0L, + 0xa.cb55b4d662945e8cf1f81ee5b4p-4L, + -0x1.3616863983e131d7935700ccd48p+0L, + 0xd.43c783ebab66074d18709d5cap-4L, + -0x1.51d5dbc56bc85976871c6e51f78p+0L, + 0x1.06253af656eb6b2ed998387aabp+0L, + -0x1.7d910a0aadc63d7a1ef7690dbb8p+0L, + /* Interval [-2.25, -2.125] (polynomial degree 22). */ + -0xf.2930890d7d675a80c36afb0fd4p-4L, + -0xc.a5cfde054eab5c6770daeca684p-4L, + 0x3.9c9e0fdebb07cdf89c61d434adp-4L, + -0x1.02a5ad35605fcf4af65a67fe8a8p+0L, + 0x9.6e9b1185bb48be9de18d8bbeb8p-4L, + -0x1.4d8332f3cfbfa116fdf648372cp+0L, + 0x1.1c0c8cb4d9f4b1d495142b53ebp+0L, + -0x1.c9a6f5ae9130ccfb9b7e39136f8p+0L, + 0x1.d7e9307fd58a2e85209d0e83eap+0L, + -0x2.921cb3473d96462f22c171712fp+0L, + 0x2.e8d59113b6f3fc1ed3b556b62cp+0L, + -0x3.cbab931624e3b6cf299cea1213p+0L, + 0x4.7d9f0f05d2c4cf91e41ea1f048p+0L, + -0x5.ade9cba31affa276fe516135eep+0L, + 0x6.dc983a62cf6ddc935ae3c5b9ap+0L, + -0x8.8d9ed100b2a7813f82cbd83e3cp+0L, + 0xa.6fa0926892835a9a29c9b8db8p+0L, + -0xc.ebc90aff4ffe319d70bef0d61p+0L, + 0xf.d69cf50ab226bacece014c0b44p+0L, + -0x1.389964ac7cfef4578eec028e5c8p+4L, + 0x1.7ff0d2090164e25901f97cab3bp+4L, + -0x1.e9e6d282da6bd004619d073071p+4L, + 0x2.5d719ab6ad4be8b5c32b0fba2ap+4L, + /* Interval [-2.375, -2.25] (polynomial degree 24). */ + -0xd.7d28d505d6181218a25f31d5e4p-4L, + -0xe.69649a3040985140cdf946827cp-4L, + 0xb.0d74a2827d053a8d4459500f88p-4L, + -0x1.924b0922853617cac181b097e48p+0L, + 0x1.d49b12bccf0a568582e2dbf8ep+0L, + -0x3.0898bb7d8c4093e6360d26bbc5p+0L, + 0x4.207a6cac711cb538684f74619ep+0L, + -0x6.39ee63ea4fb1dcac86ab337e3cp+0L, + 0x8.e2e2556a797b64a1b9328a3978p+0L, + -0xd.0e83ac82552ee5596df1706ff4p+0L, + 0x1.2e4525e0ce666e48fac68ddcdep+4L, + -0x1.b8e350d6a8f6597ed2eb3c2eff8p+4L, + 0x2.805cd69b9197ee0089dd1b1c46p+4L, + -0x3.a42585423e4d00db075f2d687ep+4L, + 0x5.4b4f409f874e2a7dcd8aa4a62ap+4L, + -0x7.b3c5829962ca1b95535db9cc4ep+4L, + 0xb.33b7b928986ec6b219e2e15a98p+4L, + -0x1.04b76dec4115106bb16316d9cd8p+8L, + 0x1.7b366d8d46f179d5c5302d6534p+8L, + -0x2.2799846ddc54813d40da622b99p+8L, + 0x3.2253a862c1078a3ccabac65bebp+8L, + -0x4.8d92cebc90a4a29816f4952f4ep+8L, + 0x6.9ebb8f9d72c66c80c4f4492e7ap+8L, + -0xa.2850a483f9ba0e43f5848b5cd8p+8L, + 0xe.e1b6bdce83b27944edab8c428p+8L, + /* Interval [-2.5, -2.375] (polynomial degree 25). */ + -0xb.74ea1bcfff94b2c01afba9daa8p-4L, + -0x1.2a82bd590c37538cab143308e3p+0L, + 0x1.88020f828b966fec66b8648d16p+0L, + -0x3.32279f040eb694970e9db0308bp+0L, + 0x5.57ac82517767e68a72142041b4p+0L, + -0x9.c2aedcfe22833de438786dc658p+0L, + 0x1.12c132f1f5577f99dbfb7ecb408p+4L, + -0x1.ea94e26628a3de3557dc349db8p+4L, + 0x3.66b4ac4fa582f5cbe7e19d10c6p+4L, + -0x6.0cf746a9cf4cbcb0004cb01f66p+4L, + 0xa.c102ef2c20d5a313cbfd37f5b8p+4L, + -0x1.31ebff06e8f08f58d1c35eacfdp+8L, + 0x2.1fd6f0c0e788660ba1f1573722p+8L, + -0x3.c6d760404305e75356a86a11d6p+8L, + 0x6.b6d18e0c31a2ba4d5b5ac78676p+8L, + -0xb.efaf5426343e6b41a823ed6c44p+8L, + 0x1.53852db2fe01305b9f336d132d8p+12L, + -0x2.5b977cb2b568382e71ca93a36bp+12L, + 0x4.310d090a6119c7d85a2786a616p+12L, + -0x7.73a518387ef1d4d04917dfb25cp+12L, + 0xd.3f965798601aabd24bdaa6e68cp+12L, + -0x1.78db20b0b166480c93cf0031198p+16L, + 0x2.9be0068b65cf13bd1cf71f0eccp+16L, + -0x4.a221230466b9cd51d5b811d6b6p+16L, + 0x8.f6f8c13e2b52aa3e30a4ce6898p+16L, + -0x1.02145337ff16b44fa7c2adf7f28p+20L, + /* Interval [-2.625, -2.5] (polynomial degree 26). */ + -0x3.d10108c27ebafad533c20eac33p-4L, + 0x1.cd557caff7d2b2085f41dbec538p+0L, + 0x3.819b4856d399520dad9776ebb9p+0L, + 0x6.8505cbad03dc34c5e42e89c4b4p+0L, + 0xb.c1b2e653a9e38f82b3997134a8p+0L, + 0x1.50a53a38f1481381051544750ep+4L, + 0x2.57ae00cbe5232cbeef4e94eb2cp+4L, + 0x4.2b156301b8604db82856d5767p+4L, + 0x7.6989ed23ca3ca751fc9c32eb88p+4L, + 0xd.2dd29765579396f3a456772c44p+4L, + 0x1.76e1c3430eb8630991d1aa8a248p+8L, + 0x2.9a77bf548873743fe65d025f56p+8L, + 0x4.a0d62ed7266389753842d7be74p+8L, + 0x8.3a6184dd32d31ec73fc6f2d37cp+8L, + 0xe.a0ade153a3bf0247db49e11ae8p+8L, + 0x1.a01359fa74d4eaf8858bbc35f68p+12L, + 0x2.e3b0a32845cbc135bae4a5216cp+12L, + 0x5.23012653815fe88456170a7dc6p+12L, + 0x9.21c92dcde748ec199bc9c65738p+12L, + 0x1.03c0f3621b4c67d2d86e5e813d8p+16L, + 0x1.cdc884edcc9f5404f2708551cb8p+16L, + 0x3.35025f0b1624d1ffc86688bf03p+16L, + 0x5.b3bd9562ebf2409c5ce99929ep+16L, + 0xa.1a229b1986d9f89cb80abccfdp+16L, + 0x1.1e69136ebd520146d51837f3308p+20L, + 0x2.2d2738c72449db2524171b9271p+20L, + 0x4.036e80cc6621b836f94f426834p+20L, + /* Interval [-2.75, -2.625] (polynomial degree 24). */ + -0x6.b5d252a56e8a75458a27ed1c2ep-4L, + 0x1.28d60383da3ac721aed3c57949p+0L, + 0x1.db6513ada8a66ea77d87d9a796p+0L, + 0x2.e217118f9d348a27f7506c4b4fp+0L, + 0x4.450112c5cbf725a0fb982fc44cp+0L, + 0x6.4af99151eae7810a75a5fceac8p+0L, + 0x9.2db598b4a97a7f69ab7be31128p+0L, + 0xd.62bef9c22471f5f17955733c6p+0L, + 0x1.379f294e412bd6255506135f4a8p+4L, + 0x1.c5827349d8865d858d4f85f3c38p+4L, + 0x2.93a7e7a75b755bbea1785a1349p+4L, + 0x3.bf9bb882afed66a08b22ed7a45p+4L, + 0x5.73c737828d2044aca95fdef33ep+4L, + 0x7.ee46534920f1c81574db260f0ep+4L, + 0xb.891c6b837b513eaf1592fe78ccp+4L, + 0x1.0c775d815bf741526a3dd66ded8p+8L, + 0x1.867ee44cf11f26455a8924a56bp+8L, + 0x2.37fe968baa1018e55cae680f1dp+8L, + 0x3.3a2c557f686679eb5d8e960fd1p+8L, + 0x4.b1ba0539d4d80cc9174738b992p+8L, + 0x6.d3fd80155b6d2211956cb6bc5ap+8L, + 0x9.eb5a96b0ee3d9ca523f5fbc1fp+8L, + 0xe.6b37429c1acc7dc19ef312dda4p+8L, + 0x1.621132d6aa138b203a28e4792fp+12L, + 0x2.09610219270e2ce11a985d4d36p+12L, + /* Interval [-2.875, -2.75] (polynomial degree 23). */ + -0x8.a41b1e4f36ff88dc820815607cp-4L, + 0xc.da87d3b69dc0f2f9c6f368b8c8p-4L, + 0x1.1474ad5c36158a7bea04fd30b28p+0L, + 0x1.761ecb90c555df6555b7dbb9ce8p+0L, + 0x1.d279bff9ae291caf6c4b17497f8p+0L, + 0x2.4e5d00559a6e2b9b5d7e35b575p+0L, + 0x2.d57545a75cee8743b1ff6e22b8p+0L, + 0x3.8514eee3aac88b89d2d4ddef4ep+0L, + 0x4.5235e3b6e1891fd9c975383318p+0L, + 0x5.562acdb10eef3c14a780490e3cp+0L, + 0x6.8ec8965c76f0b261bc41b5e532p+0L, + 0x8.15251aca144a98a1e1c0981388p+0L, + 0x9.f08d56ab9e7eee9515a457214cp+0L, + 0xc.3dbbeda2620d5be4fe8621ce6p+0L, + 0xf.0f5bfd65b3feb6d745a2cdbf9cp+0L, + 0x1.28a6ccd8dd27fb90fcaa31d37dp+4L, + 0x1.6d0a3a3091c3d64cfd1a3c5769p+4L, + 0x1.c1570107e02d5ab0b8bea6d6c98p+4L, + 0x2.28fc9b295b583fa469de7acceap+4L, + 0x2.a8a4cac0217026bbdbce34f4adp+4L, + 0x3.4532c98bce75262ac0ede53edep+4L, + 0x4.062fd9ba18e00e55c25a4f0688p+4L, + 0x5.22e00e6d9846a3451fad5587f8p+4L, + 0x6.5d0f7ce92a0bf928d4a30e92c6p+4L, + /* Interval [-3, -2.875] (polynomial degree 22). */ + -0xa.046d667e468f3e44dcae1afcc8p-4L, + 0x9.70b88dcc006c214d8d996fdf7p-4L, + 0xa.a8a39421c86d3ff24931a093c4p-4L, + 0xd.2f4d1363f324da2b357c850124p-4L, + 0xd.ca9aa1a3a5c00de11bf5d7047p-4L, + 0xf.cf09c31eeb52a45dfb25e50ebcp-4L, + 0x1.04b133a39ed8a096914cc78812p+0L, + 0x1.22b547a06edda9447f516a2ee7p+0L, + 0x1.2c57fce7db86a91c8d0f12077b8p+0L, + 0x1.4aade4894708fb8b78365e9bf88p+0L, + 0x1.579c8b7b67ec5179ecc4e9c7dp+0L, + 0x1.776820e7fc7361c50e7ef40a88p+0L, + 0x1.883ab28c72ef238ada6c480ab18p+0L, + 0x1.aa2ef6e1d11b9fcea06a1dcab1p+0L, + 0x1.bf4ad50f2dd2aeb02395ea08648p+0L, + 0x1.e40206a5477615838e02279dfc8p+0L, + 0x1.fdcbcfd4b0777fb173b85d5b398p+0L, + 0x2.25e32b3b3c89e833029169a17bp+0L, + 0x2.44ce344ff0bda6570fe3d0a76dp+0L, + 0x2.70bfba6fa079faf2dbf31d2216p+0L, + 0x2.953e22a97725cc179ad21024fap+0L, + 0x2.d8ccc51524659a499eee0f267p+0L, + 0x3.080fbb09c14936c2171c8a51bcp+0L, + }; + +static const size_t poly_deg[] = + { + 21, + 22, + 24, + 25, + 26, + 24, + 23, + 22, + }; + +static const size_t poly_end[] = + { + 21, + 44, + 69, + 95, + 122, + 147, + 171, + 194, + }; + +/* 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 -48 < 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 < 18) + { + int n_up = (19 - 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-128ibm/lgamma_productl.c b/sysdeps/ieee754/ldbl-128ibm/lgamma_productl.c new file mode 100644 index 0000000000..c5fa81be8e --- /dev/null +++ b/sysdeps/ieee754/ldbl-128ibm/lgamma_productl.c @@ -0,0 +1,38 @@ +/* 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> + +/* 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 x_full = x + x_eps; + long double ret = 0; + for (int i = 0; i < n; i++) + /* FIXME: no extra precision used. */ + ret += (t / (x_full + i)) * (1 + ret); + return ret; +} diff --git a/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c b/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c index bacb2cc98b..42f5e6a02d 100644 --- a/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c +++ b/sysdeps/ieee754/ldbl-128ibm/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 @@ -15,12 +15,17 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ -#include "gmp.h" -#include "gmp-impl.h" #include <ieee754.h> +#include <errno.h> #include <float.h> #include <math.h> +/* Need to set this when including gmp headers after system headers. */ +#define HAVE_ALLOCA 1 + +#include "gmp.h" +#include "gmp-impl.h" + /* Convert a multi-precision integer of the needed number of bits (106 for long double) and an integral power of two to a `long double' in IBM extended format. */ @@ -111,6 +116,14 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign) { hi >>= 1; u.d[0].ieee.exponent++; + if (u.d[0].ieee.exponent == IEEE754_DOUBLE_BIAS + DBL_MAX_EXP) + { + /* Overflow. The appropriate overflowed result must + be produced (if an infinity, that means the low + part must be zero). */ + __set_errno (ERANGE); + return (sign ? -LDBL_MAX : LDBL_MAX) * LDBL_MAX; + } } u.d[1].ieee.negative = !sign; lo = (1LL << 53) - lo; diff --git a/sysdeps/ieee754/ldbl-128ibm/printf_fphex.c b/sysdeps/ieee754/ldbl-128ibm/printf_fphex.c index 1ceca89485..06be1c52d4 100644 --- a/sysdeps/ieee754/ldbl-128ibm/printf_fphex.c +++ b/sysdeps/ieee754/ldbl-128ibm/printf_fphex.c @@ -1,5 +1,5 @@ /* 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. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/ldbl-128ibm/s_asinhl.c b/sysdeps/ieee754/ldbl-128ibm/s_asinhl.c index dda7f780fc..aa9a9ba213 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_asinhl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_asinhl.c @@ -46,11 +46,7 @@ long double __asinhl(long double x) ix = hx&0x7fffffffffffffffLL; if(ix>=0x7ff0000000000000LL) return x+x; /* x is inf or NaN */ if(ix< 0x3c70000000000000LL) { /* |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 */ } if(ix>0x4370000000000000LL) { /* |x| > 2**56 */ diff --git a/sysdeps/ieee754/ldbl-128ibm/s_atanl.c b/sysdeps/ieee754/ldbl-128ibm/s_atanl.c index 6ddf4b1c5e..0560d820ae 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_atanl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_atanl.c @@ -199,11 +199,7 @@ __atanl (long double x) if (k <= 0x3c800000) /* |x| <= 2**-55. */ { - if (fabsl (x) < LDBL_MIN) - { - long double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); /* Raise inexact. */ if (1e300L + x > 0.0) return x; diff --git a/sysdeps/ieee754/ldbl-128ibm/s_ceill.c b/sysdeps/ieee754/ldbl-128ibm/s_ceill.c index dc4209e9a1..ac649b7215 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_ceill.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_ceill.c @@ -1,6 +1,6 @@ /* Ceil (round to +inf) long double floating-point values. IBM extended format long double version. - Copyright (C) 2006-2015 Free Software Foundation, Inc. + Copyright (C) 2006-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-128ibm/s_erfl.c b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c index f6fcf48cfa..7b761b0afa 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c @@ -106,6 +106,7 @@ #include <math.h> #include <math_private.h> #include <math_ldbl_opt.h> +#include <fix-int-fp-convert-zero.h> /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */ @@ -146,7 +147,6 @@ deval (long double x, const long double *p, int n) static const long double tiny = 1e-300L, - half = 0.5L, one = 1.0L, two = 2.0L, /* 2/sqrt(pi) - 1 */ @@ -805,11 +805,7 @@ __erfl (long double x) if (x == 0) return x; 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; @@ -843,7 +839,10 @@ __erfcl (long double x) if (ix >= 0x7ff00000) { /* erfc(nan)=nan */ /* erfc(+-inf)=0,2 */ - return (long double) ((hx >> 31) << 1) + one / x; + long double ret = (long double) ((hx >> 31) << 1) + one / x; + if (FIX_INT_FP_CONVERT_ZERO && ret == 0.0L) + return 0.0L; + return ret; } if (ix < 0x3fd00000) /* |x| <1/4 */ diff --git a/sysdeps/ieee754/ldbl-128ibm/s_expm1l.c b/sysdeps/ieee754/ldbl-128ibm/s_expm1l.c index 0464f79043..66f75e1c80 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_expm1l.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_expm1l.c @@ -82,8 +82,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 = 1e290L; @@ -105,14 +103,9 @@ __expm1l (long double x) return __expl (x); if (ix >= 0x7ff00000) { - /* Infinity. */ + /* Infinity (which must be negative infinity). */ if (((ix - 0x7ff00000) | lx) == 0) - { - if (sign) - return -1.0L; - else - return x; - } + return -1.0L; /* NaN. No invalid exception. */ return x; } @@ -121,13 +114,6 @@ __expm1l (long double x) if ((ix | lx) == 0) return x; - /* Overflow. */ - if (x > maxlog) - { - __set_errno (ERANGE); - return (big * big); - } - /* Minimum value. */ if (x < minarg) return (4.0/big - 1.0L); diff --git a/sysdeps/ieee754/ldbl-128ibm/s_finitel.c b/sysdeps/ieee754/ldbl-128ibm/s_finitel.c index b562ce60b6..3b9e3de292 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_finitel.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_finitel.c @@ -34,7 +34,7 @@ ___finitel (long double x) xhi = ldbl_high (x); EXTRACT_WORDS64 (hx, xhi); - hx &= 0x7fffffffffffffffLL; + hx &= 0x7ff0000000000000LL; hx -= 0x7ff0000000000000LL; return hx >> 63; } diff --git a/sysdeps/ieee754/ldbl-128ibm/s_floorl.c b/sysdeps/ieee754/ldbl-128ibm/s_floorl.c index c911f19247..912230870a 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_floorl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_floorl.c @@ -1,6 +1,6 @@ /* Round to int long double floating-point values. IBM extended format long double version. - Copyright (C) 2006-2015 Free Software Foundation, Inc. + Copyright (C) 2006-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-128ibm/s_fmal.c b/sysdeps/ieee754/ldbl-128ibm/s_fmal.c index 837444aa4c..eb3ee3cfb8 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_fmal.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_fmal.c @@ -1,5 +1,5 @@ /* Compute x * y + z as ternary operation. - Copyright (C) 2011-2015 Free Software Foundation, Inc. + Copyright (C) 2011-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by David Flaherty <flaherty@linux.vnet.ibm.com>. diff --git a/sysdeps/ieee754/ldbl-128ibm/s_fpclassifyl.c b/sysdeps/ieee754/ldbl-128ibm/s_fpclassifyl.c index a33c2cf73b..83c3a8dc51 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_fpclassifyl.c +++ b/sysdeps/ieee754/ldbl-128ibm/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-128ibm/s_isinf_nsl.c b/sysdeps/ieee754/ldbl-128ibm/s_isinf_nsl.c deleted file mode 100644 index 54e72c9166..0000000000 --- a/sysdeps/ieee754/ldbl-128ibm/s_isinf_nsl.c +++ /dev/null @@ -1,23 +0,0 @@ -/* - * __isinf_nsl(x) returns != 0 if x is ±inf, else 0; - * no branching! - * slightly dodgy in relying on signed shift right copying sign bit - */ - -#include <math.h> -#include <math_private.h> - -int -__isinf_nsl (long double x) -{ - double xhi; - int64_t hx, mask; - - xhi = ldbl_high (x); - EXTRACT_WORDS64 (hx, xhi); - - mask = (hx & 0x7fffffffffffffffLL) ^ 0x7ff0000000000000LL; - mask |= -mask; - mask >>= 63; - return ~mask; -} diff --git a/sysdeps/ieee754/ldbl-128ibm/s_issignalingl.c b/sysdeps/ieee754/ldbl-128ibm/s_issignalingl.c index 98f692d706..091513908b 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_issignalingl.c +++ b/sysdeps/ieee754/ldbl-128ibm/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-128ibm/s_llrintl.c b/sysdeps/ieee754/ldbl-128ibm/s_llrintl.c index 524340bbf4..860ede1e0b 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_llrintl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_llrintl.c @@ -1,6 +1,6 @@ /* Round to long long int long double floating-point values. IBM extended format long double version. - Copyright (C) 2006-2015 Free Software Foundation, Inc. + Copyright (C) 2006-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-128ibm/s_llroundl.c b/sysdeps/ieee754/ldbl-128ibm/s_llroundl.c index c8224b01a4..9fba087d5b 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_llroundl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_llroundl.c @@ -1,6 +1,6 @@ /* Round to long long int long double floating-point values. IBM extended format long double version. - Copyright (C) 2006-2015 Free Software Foundation, Inc. + Copyright (C) 2006-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-128ibm/s_log1pl.c b/sysdeps/ieee754/ldbl-128ibm/s_log1pl.c index a0e24d7f22..743693bfd6 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_log1pl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_log1pl.c @@ -117,8 +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 big = 2e300L; static const long double zero = 0.0L; @@ -149,7 +147,7 @@ __log1pl (long double xm1) if (x <= 0.0L) { if (x == 0.0L) - return (-1.0L / (x - x)); + return (-1.0L / 0.0L); else return (zero / (x - x)); } diff --git a/sysdeps/ieee754/ldbl-128ibm/s_logbl.c b/sysdeps/ieee754/ldbl-128ibm/s_logbl.c index 22e5fc24c0..3c07c5e8e2 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_logbl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_logbl.c @@ -22,6 +22,7 @@ #include <math.h> #include <math_private.h> #include <math_ldbl_opt.h> +#include <fix-int-fp-convert-zero.h> long double __logbl (long double x) @@ -53,6 +54,8 @@ __logbl (long double x) if ((hxs ^ lx) < 0 && (lx & 0x7fffffffffffffffLL) != 0) rhx--; } + if (FIX_INT_FP_CONVERT_ZERO && rhx == 1023) + return 0.0L; return (long double) (rhx - 1023); } #ifndef __logbl diff --git a/sysdeps/ieee754/ldbl-128ibm/s_lrintl.c b/sysdeps/ieee754/ldbl-128ibm/s_lrintl.c index c5756cae12..988de70c5a 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_lrintl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_lrintl.c @@ -1,6 +1,6 @@ /* Round to long int long double floating-point values. IBM extended format long double version. - Copyright (C) 2006-2015 Free Software Foundation, Inc. + Copyright (C) 2006-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-128ibm/s_lroundl.c b/sysdeps/ieee754/ldbl-128ibm/s_lroundl.c index 42790a5bab..aa48f680d4 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_lroundl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_lroundl.c @@ -1,6 +1,6 @@ /* Round to long int long double floating-point values. IBM extended format long double version. - Copyright (C) 2006-2015 Free Software Foundation, Inc. + Copyright (C) 2006-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-128ibm/s_nearbyintl.c b/sysdeps/ieee754/ldbl-128ibm/s_nearbyintl.c index ef1b3dca1e..08134edd10 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_nearbyintl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_nearbyintl.c @@ -1,6 +1,6 @@ /* Round to int long double floating-point values without raising inexact. IBM extended format long double version. - Copyright (C) 2006-2015 Free Software Foundation, Inc. + Copyright (C) 2006-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 @@ -36,7 +36,9 @@ __nearbyintl (long double x) union ibm_extended_long_double u; u.ld = x; - if (fabs (u.d[0].d) < TWO52) + if (!isfinite (u.d[0].d)) + return x; + else if (fabs (u.d[0].d) < TWO52) { double xh = u.d[0].d; double high = u.d[0].d; @@ -65,7 +67,7 @@ __nearbyintl (long double x) } else if (fabs (u.d[1].d) < TWO52 && u.d[1].d != 0.0) { - double high, low, tau; + double high = u.d[0].d, low = u.d[1].d, tau; /* In this case we have to round the low double and handle any adjustment to the high double that may be caused by rounding (up). This is complicated by the fact that the high double @@ -78,8 +80,6 @@ __nearbyintl (long double x) { /* If the high/low doubles are the same sign then simply round the low double. */ - high = u.d[0].d; - low = u.d[1].d; } else if (u.d[1].d < 0.0) { @@ -88,8 +88,8 @@ __nearbyintl (long double x) tau = __nextafter (u.d[0].d, 0.0); tau = (u.d[0].d - tau) * 2.0; - high = u.d[0].d - tau; - low = u.d[1].d + tau; + high -= tau; + low += tau; } low += TWO52; low -= TWO52; @@ -100,8 +100,6 @@ __nearbyintl (long double x) { /* If the high/low doubles are the same sign then simply round the low double. */ - high = u.d[0].d; - low = u.d[1].d; } else if (u.d[1].d > 0.0) { @@ -109,8 +107,8 @@ __nearbyintl (long double x) adjust for that. */ tau = __nextafter (u.d[0].d, 0.0); tau = (u.d[0].d - tau) * 2.0; - high = u.d[0].d - tau; - low = u.d[1].d + tau; + high -= tau; + low += tau; } low = TWO52 - low; low = -(low - TWO52); diff --git a/sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c b/sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c index bf57cb89d6..515aa1ef5b 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c @@ -24,6 +24,8 @@ static char rcsid[] = "$NetBSD: $"; * Special cases: */ +#include <errno.h> +#include <float.h> #include <math.h> #include <math_private.h> #include <math_ldbl_opt.h> @@ -68,6 +70,7 @@ long double __nextafterl(long double x, long double y) if((hx==0xffefffffffffffffLL)&&(lx==0xfc8ffffffffffffeLL)) { u = x+x; /* overflow, return -inf */ math_force_eval (u); + __set_errno (ERANGE); return y; } if (hx >= 0x7ff0000000000000LL) { @@ -76,12 +79,13 @@ long double __nextafterl(long double x, long double y) } if(ihx <= 0x0360000000000000LL) { /* x <= LDBL_MIN */ u = math_opt_barrier (x); - x -= __LDBL_DENORM_MIN__; + x -= LDBL_TRUE_MIN; if (ihx < 0x0360000000000000LL || (hx > 0 && lx <= 0) || (hx < 0 && lx > 1)) { u = u * u; math_force_eval (u); /* raise underflow flag */ + __set_errno (ERANGE); } return x; } @@ -107,6 +111,7 @@ long double __nextafterl(long double x, long double y) if((hx==0x7fefffffffffffffLL)&&(lx==0x7c8ffffffffffffeLL)) { u = x+x; /* overflow, return +inf */ math_force_eval (u); + __set_errno (ERANGE); return y; } if ((uint64_t) hx >= 0xfff0000000000000ULL) { @@ -115,14 +120,15 @@ long double __nextafterl(long double x, long double y) } if(ihx <= 0x0360000000000000LL) { /* x <= LDBL_MIN */ u = math_opt_barrier (x); - x += __LDBL_DENORM_MIN__; + x += LDBL_TRUE_MIN; if (ihx < 0x0360000000000000LL || (hx > 0 && lx < 0 && lx != 0x8000000000000001LL) || (hx < 0 && lx >= 0)) { u = u * u; math_force_eval (u); /* raise underflow flag */ + __set_errno (ERANGE); } - if (x == 0.0L) /* handle negative __LDBL_DENORM_MIN__ case */ + if (x == 0.0L) /* handle negative LDBL_TRUE_MIN case */ x = -0.0L; return x; } diff --git a/sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c b/sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c index b40cf167f3..d8f4fc6523 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c @@ -25,6 +25,7 @@ static char rcsid[] = "$NetBSD: $"; * Special cases: */ +#include <errno.h> #include <math.h> #include <math_private.h> #include <math_ldbl_opt.h> @@ -74,15 +75,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-128ibm/s_nexttowardf.c b/sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c index 19522f4762..7c5d1cc112 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c @@ -18,6 +18,7 @@ static char rcsid[] = "$NetBSD: $"; #endif +#include <errno.h> #include <math.h> #include <math_private.h> #include <math_ldbl_opt.h> @@ -63,15 +64,14 @@ float __nexttowardf(float x, long double y) } hy = hx&0x7f800000; if(hy>=0x7f800000) { - x = x+x; /* overflow */ - if (FLT_EVAL_METHOD != 0) - /* Force conversion to float. */ - asm ("" : "+m"(x)); - return x; + float u = x+x; /* overflow */ + math_force_eval (u); + __set_errno (ERANGE); } if(hy<0x00800000) { /* underflow */ float u = x*x; math_force_eval (u); /* raise underflow flag */ + __set_errno (ERANGE); } SET_FLOAT_WORD(x,hx); return x; diff --git a/sysdeps/ieee754/ldbl-128ibm/s_remquol.c b/sysdeps/ieee754/ldbl-128ibm/s_remquol.c index 57c50590c9..20e17cc823 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_remquol.c +++ b/sysdeps/ieee754/ldbl-128ibm/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-128ibm/s_rintl.c b/sysdeps/ieee754/ldbl-128ibm/s_rintl.c index d00abe860c..8c51ded1d6 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_rintl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_rintl.c @@ -1,6 +1,6 @@ /* Round to int long double floating-point values. IBM extended format long double version. - Copyright (C) 2006-2015 Free Software Foundation, Inc. + Copyright (C) 2006-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-128ibm/s_roundl.c b/sysdeps/ieee754/ldbl-128ibm/s_roundl.c index ed8e00ab71..20813ed366 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_roundl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_roundl.c @@ -1,6 +1,6 @@ /* Round to int long double floating-point values. IBM extended format long double version. - Copyright (C) 2006-2015 Free Software Foundation, Inc. + Copyright (C) 2006-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-128ibm/s_scalbnl.c b/sysdeps/ieee754/ldbl-128ibm/s_scalbnl.c index a0d7aa9d63..0c4508835e 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_scalbnl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_scalbnl.c @@ -102,8 +102,3 @@ long double __scalbnl (long double x, int n) x = ldbl_pack (xhi, xlo); return x*twolm54; } -#if IS_IN (libm) -long_double_symbol (libm, __scalbnl, scalbnl); -#else -long_double_symbol (libc, __scalbnl, scalbnl); -#endif diff --git a/sysdeps/ieee754/ldbl-128ibm/s_signbitl.c b/sysdeps/ieee754/ldbl-128ibm/s_signbitl.c index e95ad55e32..b4e8256329 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_signbitl.c +++ b/sysdeps/ieee754/ldbl-128ibm/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. @@ -18,18 +18,12 @@ <http://www.gnu.org/licenses/>. */ #include <math.h> -#include <math_private.h> #include <math_ldbl_opt.h> int ___signbitl (long double x) { - int64_t e; - double xhi; - - xhi = ldbl_high (x); - EXTRACT_WORDS64 (e, xhi); - return e < 0; + return __builtin_signbitl (x); } #if IS_IN (libm) long_double_symbol (libm, ___signbitl, __signbitl); diff --git a/sysdeps/ieee754/ldbl-128ibm/s_sincosl.c b/sysdeps/ieee754/ldbl-128ibm/s_sincosl.c index 8a49c54d6c..fae4020a7b 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_sincosl.c +++ b/sysdeps/ieee754/ldbl-128ibm/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>. @@ -42,7 +42,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-128ibm/s_tanhl.c b/sysdeps/ieee754/ldbl-128ibm/s_tanhl.c index c63e25345d..e6457a1c1c 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_tanhl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_tanhl.c @@ -29,15 +29,16 @@ static char rcsid[] = "$NetBSD: s_tanh.c,v 1.7 1995/05/10 20:48:22 jtc Exp $"; * 2**-57 < x <= 1 : tanh(x) := -----; t = expm1(-2x) * t + 2 * 2 - * 1 <= x <= 22.0 : tanh(x) := 1- ----- ; t=expm1(2x) + * 1 <= x <= 40.0 : tanh(x) := 1- ----- ; t=expm1(2x) * t + 2 - * 22.0 < x <= INF : tanh(x) := 1. + * 40.0 < x <= INF : tanh(x) := 1. * * Special cases: * tanh(NaN) is NaN; * only tanh(0)=0 is exact for finite argument. */ +#include <float.h> #include <math.h> #include <math_private.h> #include <math_ldbl_opt.h> @@ -61,12 +62,15 @@ long double __tanhl(long double x) else return one/x-one; /* tanh(NaN) = NaN */ } - /* |x| < 22 */ - if (ix < 0x4036000000000000LL) { /* |x|<22 */ + /* |x| < 40 */ + if (ix < 0x4044000000000000LL) { /* |x|<40 */ if (ix == 0) return x; /* x == +-0 */ if (ix<0x3c60000000000000LL) /* |x|<2**-57 */ - return x*(one+x); /* tanh(small) = small */ + { + math_check_force_underflow (x); + return x; /* tanh(small) = small */ + } if (ix>=0x3ff0000000000000LL) { /* |x|>=1 */ t = __expm1l(two*fabsl(x)); z = one - two/(t+two); @@ -74,7 +78,7 @@ long double __tanhl(long double x) t = __expm1l(-two*fabsl(x)); z= -t/(t+two); } - /* |x| > 22, return +-1 */ + /* |x| > 40, return +-1 */ } else { z = one - tiny; /* raised inexact flag */ } diff --git a/sysdeps/ieee754/ldbl-128ibm/s_truncl.c b/sysdeps/ieee754/ldbl-128ibm/s_truncl.c index 91dc9753f9..df58b64b53 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_truncl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_truncl.c @@ -1,6 +1,6 @@ /* Truncate (toward zero) long double floating-point values. IBM extended format long double version. - Copyright (C) 2006-2015 Free Software Foundation, Inc. + Copyright (C) 2006-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-128ibm/strtod_nan_ldouble.h b/sysdeps/ieee754/ldbl-128ibm/strtod_nan_ldouble.h new file mode 100644 index 0000000000..d827112446 --- /dev/null +++ b/sysdeps/ieee754/ldbl-128ibm/strtod_nan_ldouble.h @@ -0,0 +1,30 @@ +/* Convert string for NaN payload to corresponding NaN. For ldbl-128ibm. + 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 ibm_extended_long_double u; \ + u.ld = (flt); \ + u.d[0].ieee_nan.mantissa0 = (mant) >> 32; \ + u.d[0].ieee_nan.mantissa1 = (mant); \ + if ((u.d[0].ieee.mantissa0 | u.d[0].ieee.mantissa1) != 0) \ + (flt) = u.ld; \ + } \ + while (0) diff --git a/sysdeps/ieee754/ldbl-128ibm/strtold_l.c b/sysdeps/ieee754/ldbl-128ibm/strtold_l.c index 3e2f69e353..a8181740a8 100644 --- a/sysdeps/ieee754/ldbl-128ibm/strtold_l.c +++ b/sysdeps/ieee754/ldbl-128ibm/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 @@ -30,25 +30,19 @@ extern long double ____new_wcstold_l (const wchar_t *, wchar_t **, __locale_t); # define STRTOF __new_wcstold_l # define __STRTOF ____new_wcstold_l # define ____STRTOF_INTERNAL ____wcstold_l_internal +# define STRTOF_NAN __wcstold_nan #else extern long double ____new_strtold_l (const char *, char **, __locale_t); # define STRTOF __new_strtold_l # define __STRTOF ____new_strtold_l # define ____STRTOF_INTERNAL ____strtold_l_internal +# define STRTOF_NAN __strtold_nan #endif extern __typeof (__STRTOF) STRTOF; libc_hidden_proto (__STRTOF) libc_hidden_proto (STRTOF) #define MPN2FLOAT __mpn_construct_long_double #define FLOAT_HUGE_VAL HUGE_VALL -# define SET_MANTISSA(flt, mant) \ - do { union ibm_extended_long_double u; \ - u.ld = (flt); \ - u.d[0].ieee_nan.mantissa0 = (mant) >> 32; \ - u.d[0].ieee_nan.mantissa1 = (mant); \ - if ((u.d[0].ieee.mantissa0 | u.d[0].ieee.mantissa1) != 0) \ - (flt) = u.ld; \ - } while (0) #include <strtod_l.c> diff --git a/sysdeps/ieee754/ldbl-128ibm/t_sincosl.c b/sysdeps/ieee754/ldbl-128ibm/t_sincosl.c index a688577c2f..737c7c73fa 100644 --- a/sysdeps/ieee754/ldbl-128ibm/t_sincosl.c +++ b/sysdeps/ieee754/ldbl-128ibm/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-128ibm/tst-strtold-ldbl-128ibm.c b/sysdeps/ieee754/ldbl-128ibm/tst-strtold-ldbl-128ibm.c new file mode 100644 index 0000000000..14dc683619 --- /dev/null +++ b/sysdeps/ieee754/ldbl-128ibm/tst-strtold-ldbl-128ibm.c @@ -0,0 +1,85 @@ +/* Test for ldbl-128ibm strtold overflow to infinity (bug 14551). + 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 <errno.h> +#include <fenv.h> +#include <math.h> +#include <stdio.h> +#include <stdlib.h> + +static int +test_strtold_value (const char *s, double exp_hi, double exp_lo, int exp_exc, + int exp_errno) +{ + int result = 0; + union { long double ld; double d[2]; } x; + feclearexcept (FE_ALL_EXCEPT); + errno = 0; + x.ld = strtold (s, NULL); + int exc = fetestexcept (FE_ALL_EXCEPT); + int new_errno = errno; + printf ("strtold (\"%s\") returned (%a, %a), exceptions 0x%x, errno %d\n", + s, x.d[0], x.d[1], exc, new_errno); + if (x.d[0] == exp_hi) + printf ("PASS: strtold (\"%s\") high == %a\n", s, exp_hi); + else + { + printf ("FAIL: strtold (\"%s\") high == %a\n", s, exp_hi); + result = 1; + } + if (x.d[1] == exp_lo) + printf ("PASS: strtold (\"%s\") low == %a\n", s, exp_lo); + else + { + printf ("FAIL: strtold (\"%s\") low == %a\n", s, exp_lo); + result = 1; + } + if (exc == exp_exc) + printf ("PASS: strtold (\"%s\") exceptions 0x%x\n", s, exp_exc); + else + { + printf ("FAIL: strtold (\"%s\") exceptions 0x%x\n", s, exp_exc); + result = 1; + } + if (new_errno == exp_errno) + printf ("PASS: strtold (\"%s\") errno %d\n", s, exp_errno); + else + { + printf ("FAIL: strtold (\"%s\") errno %d\n", s, exp_errno); + result = 1; + } + return result; +} + +static int +do_test (void) +{ + int result = 0; + result |= test_strtold_value ("0x1.fffffffffffff8p+1023", INFINITY, 0, + FE_OVERFLOW | FE_INEXACT, ERANGE); + result |= test_strtold_value ("-0x1.fffffffffffff8p+1023", -INFINITY, 0, + FE_OVERFLOW | FE_INEXACT, ERANGE); + result |= test_strtold_value ("0x1.ffffffffffffffp+1023", INFINITY, 0, + FE_OVERFLOW | FE_INEXACT, ERANGE); + result |= test_strtold_value ("-0x1.ffffffffffffffp+1023", -INFINITY, 0, + FE_OVERFLOW | FE_INEXACT, ERANGE); + return result; +} + +#define TEST_FUNCTION do_test () +#include "../../../test-skeleton.c" diff --git a/sysdeps/ieee754/ldbl-128ibm/w_expl.c b/sysdeps/ieee754/ldbl-128ibm/w_expl.c index fb5c8d3629..c9d44b61dd 100644 --- a/sysdeps/ieee754/ldbl-128ibm/w_expl.c +++ b/sysdeps/ieee754/ldbl-128ibm/w_expl.c @@ -2,9 +2,6 @@ #include <math_private.h> #include <math_ldbl_opt.h> -static const long double o_thres = 709.78271289338399678773454114191496482L; -static const long double u_thres = -744.44007192138126231410729844608163411L; - long double __expl(long double x) /* wrapper exp */ { long double z; @@ -13,9 +10,9 @@ long double __expl(long double x) /* wrapper exp */ return z; if (isfinite(x)) { - if (x >= o_thres) + if (!isfinite (z)) return __kernel_standard_l(x,x,206); /* exp overflow */ - else if (x <= u_thres) + else if (z == 0.0L) return __kernel_standard_l(x,x,207); /* exp underflow */ } return z; diff --git a/sysdeps/ieee754/ldbl-128ibm/w_log1pl.c b/sysdeps/ieee754/ldbl-128ibm/w_log1pl.c index 279a62a0f4..969fadc205 100644 --- a/sysdeps/ieee754/ldbl-128ibm/w_log1pl.c +++ b/sysdeps/ieee754/ldbl-128ibm/w_log1pl.c @@ -1,5 +1,5 @@ /* Wrapper for __log1pl that handles setting errno. - Copyright (C) 2015 Free Software Foundation, Inc. + 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 diff --git a/sysdeps/ieee754/ldbl-128ibm/w_scalblnl.c b/sysdeps/ieee754/ldbl-128ibm/w_scalblnl.c index 32d4a14ca9..7e73c9abf8 100644 --- a/sysdeps/ieee754/ldbl-128ibm/w_scalblnl.c +++ b/sysdeps/ieee754/ldbl-128ibm/w_scalblnl.c @@ -1,5 +1,5 @@ /* Wrapper for __scalblnl handles setting errno. - Copyright (C) 2014-2015 Free Software Foundation, Inc. + Copyright (C) 2014-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-128ibm/x2y2m1l.c b/sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c index a001b58ca5..da2e929175 100644 --- a/sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c +++ b/sysdeps/ieee754/ldbl-128ibm/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,13 +80,13 @@ 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) { - double vals[12]; + double vals[13]; SET_RESTORE_ROUND (FE_TONEAREST); union ibm_extended_long_double xu, yu; xu.ld = x; @@ -105,25 +105,19 @@ __x2y2m1l (long double x, long double y) vals[8] *= 2.0; vals[9] *= 2.0; mul_split (&vals[11], &vals[10], yu.d[1].d, yu.d[1].d); - if (xu.d[0].d >= 0.75) - vals[1] -= 1.0; - else - { - vals[1] -= 0.5; - vals[7] -= 0.5; - } - qsort (vals, 12, sizeof (double), compare); + vals[12] = -1.0; + qsort (vals, 13, sizeof (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 <= 10; i++) + for (size_t i = 0; i <= 11; i++) { add_split (&vals[i + 1], &vals[i], vals[i + 1], vals[i]); - qsort (vals + i + 1, 11 - i, sizeof (double), compare); + qsort (vals + i + 1, 12 - i, sizeof (double), compare); } /* Now any error from this addition will be small. */ - long double retval = (long double) vals[11]; - for (size_t i = 10; i != (size_t) -1; i--) + long double retval = (long double) vals[12]; + for (size_t i = 11; i != (size_t) -1; i--) retval += (long double) vals[i]; return retval; } |