diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64')
117 files changed, 1293 insertions, 953 deletions
diff --git a/sysdeps/ieee754/dbl-64/MathLib.h b/sysdeps/ieee754/dbl-64/MathLib.h index 8d8de6fa1b..94182e45c1 100644 --- a/sysdeps/ieee754/dbl-64/MathLib.h +++ b/sysdeps/ieee754/dbl-64/MathLib.h @@ -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/dbl-64/asincos.tbl b/sysdeps/ieee754/dbl-64/asincos.tbl index 84efa43a2d..8da9ff3c07 100644 --- a/sysdeps/ieee754/dbl-64/asincos.tbl +++ b/sysdeps/ieee754/dbl-64/asincos.tbl @@ -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/dbl-64/atnat.h b/sysdeps/ieee754/dbl-64/atnat.h index a9c265fe41..a9649c57ea 100644 --- a/sysdeps/ieee754/dbl-64/atnat.h +++ b/sysdeps/ieee754/dbl-64/atnat.h @@ -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/dbl-64/atnat2.h b/sysdeps/ieee754/dbl-64/atnat2.h index e0d65afd06..d5f4ab946d 100644 --- a/sysdeps/ieee754/dbl-64/atnat2.h +++ b/sysdeps/ieee754/dbl-64/atnat2.h @@ -2,7 +2,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 @@ -65,10 +65,8 @@ /**/ hpi1 = {{0x3c91a626, 0x33145c07} }, /* pi/2-hpi */ /**/ mhpi = {{0xbff921fb, 0x54442d18} }, /* -pi/2 */ /**/ qpi = {{0x3fe921fb, 0x54442d18} }, /* pi/4 */ -/**/ qpi1 = {{0x3c81a626, 0x33145c07} }, /* pi/4-qpi */ /**/ mqpi = {{0xbfe921fb, 0x54442d18} }, /* -pi/4 */ /**/ tqpi = {{0x4002d97c, 0x7f3321d2} }, /* 3pi/4 */ -/**/ tqpi1 = {{0x3c9a7939, 0x4c9e8a0a} }, /* 3pi/4-tqpi */ /**/ mtqpi = {{0xc002d97c, 0x7f3321d2} }, /* -3pi/4 */ /**/ u1 = {{0x3c314c2a, 0x00000000} }, /* 9.377e-19 */ /**/ u2 = {{0x3bf955e4, 0x00000000} }, /* 8.584e-20 */ @@ -129,10 +127,8 @@ /**/ hpi1 = {{0x33145c07, 0x3c91a626} }, /* pi/2-hpi */ /**/ mhpi = {{0x54442d18, 0xbff921fb} }, /* -pi/2 */ /**/ qpi = {{0x54442d18, 0x3fe921fb} }, /* pi/4 */ -/**/ qpi1 = {{0x33145c07, 0x3c81a626} }, /* pi/4-qpi */ /**/ mqpi = {{0x54442d18, 0xbfe921fb} }, /* -pi/4 */ /**/ tqpi = {{0x7f3321d2, 0x4002d97c} }, /* 3pi/4 */ -/**/ tqpi1 = {{0x4c9e8a0a, 0x3c9a7939} }, /* 3pi/4-tqpi */ /**/ mtqpi = {{0x7f3321d2, 0xc002d97c} }, /* -3pi/4 */ /**/ u1 = {{0x00000000, 0x3c314c2a} }, /* 9.377e-19 */ /**/ u2 = {{0x00000000, 0x3bf955e4} }, /* 8.584e-20 */ diff --git a/sysdeps/ieee754/dbl-64/branred.c b/sysdeps/ieee754/dbl-64/branred.c index 7757a9c833..ca9dd90aff 100644 --- a/sysdeps/ieee754/dbl-64/branred.c +++ b/sysdeps/ieee754/dbl-64/branred.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/dbl-64/branred.h b/sysdeps/ieee754/dbl-64/branred.h index 431bce0c84..7b6b5473b6 100644 --- a/sysdeps/ieee754/dbl-64/branred.h +++ b/sysdeps/ieee754/dbl-64/branred.h @@ -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/dbl-64/dbl2mpn.c b/sysdeps/ieee754/dbl-64/dbl2mpn.c index bb28192883..cb5a070914 100644 --- a/sysdeps/ieee754/dbl-64/dbl2mpn.c +++ b/sysdeps/ieee754/dbl-64/dbl2mpn.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1993-2015 Free Software Foundation, Inc. +/* Copyright (C) 1993-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/dbl-64/dla.h b/sysdeps/ieee754/dbl-64/dla.h index 641644eeb1..d21c47a68f 100644 --- a/sysdeps/ieee754/dbl-64/dla.h +++ b/sysdeps/ieee754/dbl-64/dla.h @@ -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/dbl-64/doasin.c b/sysdeps/ieee754/dbl-64/doasin.c index 2372e1bc53..903b5484ba 100644 --- a/sysdeps/ieee754/dbl-64/doasin.c +++ b/sysdeps/ieee754/dbl-64/doasin.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/dbl-64/doasin.h b/sysdeps/ieee754/dbl-64/doasin.h index ac4c84367e..044d597f9c 100644 --- a/sysdeps/ieee754/dbl-64/doasin.h +++ b/sysdeps/ieee754/dbl-64/doasin.h @@ -2,7 +2,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/dbl-64/dosincos.c b/sysdeps/ieee754/dbl-64/dosincos.c index ed5e2c9027..931975e6c3 100644 --- a/sysdeps/ieee754/dbl-64/dosincos.c +++ b/sysdeps/ieee754/dbl-64/dosincos.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/dbl-64/dosincos.h b/sysdeps/ieee754/dbl-64/dosincos.h index 7d9324f736..02772592e8 100644 --- a/sysdeps/ieee754/dbl-64/dosincos.h +++ b/sysdeps/ieee754/dbl-64/dosincos.h @@ -2,7 +2,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/dbl-64/e_asin.c b/sysdeps/ieee754/dbl-64/e_asin.c index a7684d1078..d8c012d1d3 100644 --- a/sysdeps/ieee754/dbl-64/e_asin.c +++ b/sysdeps/ieee754/dbl-64/e_asin.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 @@ -71,11 +71,7 @@ __ieee754_asin(double x){ if (k < 0x3e500000) { - if (fabs (x) < DBL_MIN) - { - double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); return x; /* for x->0 => sin(x)=x */ } /*----------------------2^-26 <= |x| < 2^ -3 -----------------*/ diff --git a/sysdeps/ieee754/dbl-64/e_atan2.c b/sysdeps/ieee754/dbl-64/e_atan2.c index 0c2902855e..22e8fb8fef 100644 --- a/sysdeps/ieee754/dbl-64/e_atan2.c +++ b/sysdeps/ieee754/dbl-64/e_atan2.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/dbl-64/e_atanh.c b/sysdeps/ieee754/dbl-64/e_atanh.c index 6b00b800f2..395118e92a 100644 --- a/sysdeps/ieee754/dbl-64/e_atanh.c +++ b/sysdeps/ieee754/dbl-64/e_atanh.c @@ -1,4 +1,4 @@ -/* 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 Ulrich Drepper <drepper@gmail.com>, 2011. @@ -52,11 +52,7 @@ __ieee754_atanh (double x) if (__glibc_unlikely (xa < 0x1.0p-28)) { math_force_eval (huge + x); - if (fabs (x) < DBL_MIN) - { - double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); return x; } diff --git a/sysdeps/ieee754/dbl-64/e_cosh.c b/sysdeps/ieee754/dbl-64/e_cosh.c index af3910dd6e..52a5d5007d 100644 --- a/sysdeps/ieee754/dbl-64/e_cosh.c +++ b/sysdeps/ieee754/dbl-64/e_cosh.c @@ -83,6 +83,6 @@ __ieee754_cosh (double x) return x * x; /* |x| > overflowthresold, cosh(x) overflow */ - return huge * huge; + return math_narrow_eval (huge * huge); } strong_alias (__ieee754_cosh, __cosh_finite) diff --git a/sysdeps/ieee754/dbl-64/e_exp.c b/sysdeps/ieee754/dbl-64/e_exp.c index bb76907f74..ad1bc84625 100644 --- a/sysdeps/ieee754/dbl-64/e_exp.c +++ b/sysdeps/ieee754/dbl-64/e_exp.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 @@ -200,10 +200,7 @@ __ieee754_exp (double x) check_uflow_ret: if (retval < DBL_MIN) { -#if FLT_EVAL_METHOD != 0 - volatile -#endif - double force_underflow = tiny * tiny; + double force_underflow = tiny * tiny; math_force_eval (force_underflow); } if (retval == 0) diff --git a/sysdeps/ieee754/dbl-64/e_exp10.c b/sysdeps/ieee754/dbl-64/e_exp10.c index b14eaa9833..e3b39f7e31 100644 --- a/sysdeps/ieee754/dbl-64/e_exp10.c +++ b/sysdeps/ieee754/dbl-64/e_exp10.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/dbl-64/e_exp2.c b/sysdeps/ieee754/dbl-64/e_exp2.c index 30f0a8f529..7402bd7d89 100644 --- a/sysdeps/ieee754/dbl-64/e_exp2.c +++ b/sysdeps/ieee754/dbl-64/e_exp2.c @@ -1,5 +1,5 @@ /* Double-precision floating point 2^x. - 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 Geoffrey Keating <geoffk@ozemail.com.au> @@ -120,7 +120,11 @@ __ieee754_exp2 (double x) if (!unsafe) return result; else - return result * scale_u.d; + { + result *= scale_u.d; + math_check_force_underflow_nonneg (result); + return result; + } } else /* Return x, if x is a NaN or Inf; or overflow, otherwise. */ diff --git a/sysdeps/ieee754/dbl-64/e_gamma_r.c b/sysdeps/ieee754/dbl-64/e_gamma_r.c index adeb61a248..2d156850d5 100644 --- a/sysdeps/ieee754/dbl-64/e_gamma_r.c +++ b/sysdeps/ieee754/dbl-64/e_gamma_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. @@ -76,11 +76,7 @@ gamma_positive (double x, int *exp2_adj) /* Adjust into the range for applying Stirling's approximation. */ double n = __ceil (12.0 - x); -#if FLT_EVAL_METHOD != 0 - volatile -#endif - double x_tmp = x + n; - x_adj = x_tmp; + x_adj = math_narrow_eval (x + n); x_eps = (x - (x_adj - n)); prod = __gamma_product (x_adj - n, x_eps, n, &eps); } @@ -119,9 +115,6 @@ __ieee754_gamma_r (double x, int *signgamp) { int32_t hx; u_int32_t lx; -#if FLT_EVAL_METHOD != 0 - volatile -#endif double ret; EXTRACT_WORDS (hx, lx, x); @@ -157,7 +150,7 @@ __ieee754_gamma_r (double x, int *signgamp) { /* Overflow. */ *signgamp = 0; - ret = DBL_MAX * DBL_MAX; + ret = math_narrow_eval (DBL_MAX * DBL_MAX); return ret; } else @@ -194,29 +187,31 @@ __ieee754_gamma_r (double x, int *signgamp) double tret = M_PI / (-x * sinpix * gamma_positive (-x, &exp2_adj)); ret = __scalbn (tret, -exp2_adj); + math_check_force_underflow_nonneg (ret); } } + ret = math_narrow_eval (ret); } if (isinf (ret) && x != 0) { if (*signgamp < 0) { - ret = -__copysign (DBL_MAX, ret) * DBL_MAX; + ret = math_narrow_eval (-__copysign (DBL_MAX, ret) * DBL_MAX); ret = -ret; } else - ret = __copysign (DBL_MAX, ret) * DBL_MAX; + ret = math_narrow_eval (__copysign (DBL_MAX, ret) * DBL_MAX); return ret; } else if (ret == 0) { if (*signgamp < 0) { - ret = -__copysign (DBL_MIN, ret) * DBL_MIN; + ret = math_narrow_eval (-__copysign (DBL_MIN, ret) * DBL_MIN); ret = -ret; } else - ret = __copysign (DBL_MIN, ret) * DBL_MIN; + ret = math_narrow_eval (__copysign (DBL_MIN, ret) * DBL_MIN); return ret; } else diff --git a/sysdeps/ieee754/dbl-64/e_hypot.c b/sysdeps/ieee754/dbl-64/e_hypot.c index 5cbfcbeb48..f142c450a2 100644 --- a/sysdeps/ieee754/dbl-64/e_hypot.c +++ b/sysdeps/ieee754/dbl-64/e_hypot.c @@ -149,7 +149,9 @@ __ieee754_hypot (double x, double y) t1 = 1.0; GET_HIGH_WORD (high, t1); SET_HIGH_WORD (t1, high + (k << 20)); - return t1 * w; + w *= t1; + math_check_force_underflow_nonneg (w); + return w; } else return w; diff --git a/sysdeps/ieee754/dbl-64/e_j1.c b/sysdeps/ieee754/dbl-64/e_j1.c index 26ffdfe282..4827fbf3d3 100644 --- a/sysdeps/ieee754/dbl-64/e_j1.c +++ b/sysdeps/ieee754/dbl-64/e_j1.c @@ -127,12 +127,10 @@ __ieee754_j1 (double x) { if (huge + x > one) /* inexact if x!=0 necessary */ { - double ret = 0.5 * x; - if (fabs (ret) < DBL_MIN) - { - double force_underflow = ret * ret; - math_force_eval (force_underflow); - } + double ret = math_narrow_eval (0.5 * x); + math_check_force_underflow (ret); + if (ret == 0 && x != 0) + __set_errno (ERANGE); return ret; } } diff --git a/sysdeps/ieee754/dbl-64/e_jn.c b/sysdeps/ieee754/dbl-64/e_jn.c index ccef2dcd80..3fecf82f10 100644 --- a/sysdeps/ieee754/dbl-64/e_jn.c +++ b/sysdeps/ieee754/dbl-64/e_jn.c @@ -243,14 +243,15 @@ __ieee754_jn (int n, double x) ret = -b; else ret = b; + ret = math_narrow_eval (ret); } if (ret == 0) - ret = __copysign (DBL_MIN, ret) * DBL_MIN; - else if (fabs (ret) < DBL_MIN) { - double force_underflow = ret * ret; - math_force_eval (force_underflow); + ret = math_narrow_eval (__copysign (DBL_MIN, ret) * DBL_MIN); + __set_errno (ERANGE); } + else + math_check_force_underflow (ret); return ret; } strong_alias (__ieee754_jn, __jn_finite) diff --git a/sysdeps/ieee754/dbl-64/e_lgamma_r.c b/sysdeps/ieee754/dbl-64/e_lgamma_r.c index fc6f594d62..15154c0f43 100644 --- a/sysdeps/ieee754/dbl-64/e_lgamma_r.c +++ b/sysdeps/ieee754/dbl-64/e_lgamma_r.c @@ -226,6 +226,8 @@ __ieee754_lgamma_r(double x, int *signgamp) if(__builtin_expect(ix>=0x43300000, 0)) /* |x|>=2**52, must be -integer */ return x/zero; + if (x < -2.0 && x > -28.0) + return __lgamma_neg (x, signgamp); t = sin_pi(x); if(t==zero) return one/fabsf(t); /* -integer */ nadj = __ieee754_log(pi/fabs(t*x)); @@ -294,17 +296,13 @@ __ieee754_lgamma_r(double x, int *signgamp) r = (x-half)*(t-one)+w; } else /* 2**58 <= x <= inf */ - r = x*(__ieee754_log(x)-one); + r = math_narrow_eval (x*(__ieee754_log(x)-one)); /* NADJ is set for negative arguments but not otherwise, resulting in warnings that it may be used uninitialized although in the cases where it is used it has always been set. */ DIAG_PUSH_NEEDS_COMMENT; -#if __GNUC_PREREQ (4, 7) DIAG_IGNORE_NEEDS_COMMENT (4.9, "-Wmaybe-uninitialized"); -#else - DIAG_IGNORE_NEEDS_COMMENT (4.9, "-Wuninitialized"); -#endif if(hx<0) r = nadj - r; DIAG_POP_NEEDS_COMMENT; return r; diff --git a/sysdeps/ieee754/dbl-64/e_log.c b/sysdeps/ieee754/dbl-64/e_log.c index a9117fb4c5..9917dc236f 100644 --- a/sysdeps/ieee754/dbl-64/e_log.c +++ b/sysdeps/ieee754/dbl-64/e_log.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/dbl-64/e_log10.c b/sysdeps/ieee754/dbl-64/e_log10.c index 8548ee3942..df59d9dce4 100644 --- a/sysdeps/ieee754/dbl-64/e_log10.c +++ b/sysdeps/ieee754/dbl-64/e_log10.c @@ -45,6 +45,7 @@ #include <math.h> #include <math_private.h> +#include <fix-int-fp-convert-zero.h> static const double two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */ static const double ivln10 = 4.34294481903251816668e-01; /* 0x3FDBCB7B, 0x1526E50E */ @@ -77,6 +78,8 @@ __ieee754_log10 (double x) i = ((u_int32_t) k & 0x80000000) >> 31; hx = (hx & 0x000fffff) | ((0x3ff - i) << 20); y = (double) (k + i); + if (FIX_INT_FP_CONVERT_ZERO && y == 0.0) + y = 0.0; SET_HIGH_WORD (x, hx); z = y * log10_2lo + ivln10 * __ieee754_log (x); return z + y * log10_2hi; diff --git a/sysdeps/ieee754/dbl-64/e_log2.c b/sysdeps/ieee754/dbl-64/e_log2.c index 997d7cefc8..bc6a34192a 100644 --- a/sysdeps/ieee754/dbl-64/e_log2.c +++ b/sysdeps/ieee754/dbl-64/e_log2.c @@ -56,6 +56,7 @@ #include <math.h> #include <math_private.h> +#include <fix-int-fp-convert-zero.h> static const double ln2 = 0.69314718055994530942; static const double two54 = 1.80143985094819840000e+16; /* 43500000 00000000 */ @@ -101,7 +102,11 @@ __ieee754_log2 (double x) if ((0x000fffff & (2 + hx)) < 3) { /* |f| < 2**-20 */ if (f == zero) - return dk; + { + if (FIX_INT_FP_CONVERT_ZERO && dk == 0.0) + dk = 0.0; + return dk; + } R = f * f * (0.5 - 0.33333333333333333 * f); return dk - (R - f) / ln2; } diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c index 3c027fe7f5..663fa392c2 100644 --- a/sysdeps/ieee754/dbl-64/e_pow.c +++ b/sysdeps/ieee754/dbl-64/e_pow.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 @@ -119,6 +119,8 @@ __ieee754_pow (double x, double y) retval = huge * huge; else if (retval == 0) retval = tiny * tiny; + else + math_check_force_underflow_nonneg (retval); return retval; } @@ -243,7 +245,8 @@ static double SECTION log1 (double x, double *delta, double *error) { - int i, j, m; + unsigned int i, j; + int m; double uu, vv, eps, nx, e, e1, e2, t, t1, t2, res, add = 0; mynumber u, v; #ifdef BIG_ENDI @@ -342,7 +345,8 @@ static double SECTION my_log2 (double x, double *delta, double *error) { - int i, j, m; + unsigned int i, j; + int m; double uu, vv, eps, nx, e, e1, e2, t, t1, t2, res, add = 0; double ou1, ou2, lu1, lu2, ov, lv1, lv2, a, a1, a2; double y, yy, z, zz, j1, j2, j7, j8; diff --git a/sysdeps/ieee754/dbl-64/e_remainder.c b/sysdeps/ieee754/dbl-64/e_remainder.c index 7f3a02b3c2..a445e74b5c 100644 --- a/sysdeps/ieee754/dbl-64/e_remainder.c +++ b/sysdeps/ieee754/dbl-64/e_remainder.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 @@ -130,6 +130,8 @@ __ieee754_remainder (double x, double y) d = fabs (z); if (d <= fabs (d - y)) return z; + else if (d == y) + return 0.0 * x; else return (z > 0) ? z - y : z + y; } diff --git a/sysdeps/ieee754/dbl-64/e_sinh.c b/sysdeps/ieee754/dbl-64/e_sinh.c index 4ff28bf85d..8479bdd9b8 100644 --- a/sysdeps/ieee754/dbl-64/e_sinh.c +++ b/sysdeps/ieee754/dbl-64/e_sinh.c @@ -32,6 +32,7 @@ static char rcsid[] = "$NetBSD: e_sinh.c,v 1.7 1995/05/10 20:46:13 jtc Exp $"; * only sinh(0)=0 is exact for finite x. */ +#include <float.h> #include <math.h> #include <math_private.h> @@ -58,10 +59,12 @@ __ieee754_sinh (double x) /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */ if (ix < 0x40360000) /* |x|<22 */ { - if (__glibc_unlikely (ix < 0x3e300000)) /* |x|<2**-28 */ + if (__glibc_unlikely (ix < 0x3e300000)) { /* |x|<2**-28 */ + math_check_force_underflow (x); if (shuge + x > one) return x; - /* sinh(tiny) = tiny with inexact */ + /* sinh(tiny) = tiny with inexact */ + } t = __expm1 (fabs (x)); if (ix < 0x3ff00000) return h * (2.0 * t - t * t / (t + one)); @@ -82,6 +85,6 @@ __ieee754_sinh (double x) } /* |x| > overflowthresold, sinh(x) overflow */ - return x * shuge; + return math_narrow_eval (x * shuge); } strong_alias (__ieee754_sinh, __sinh_finite) diff --git a/sysdeps/ieee754/dbl-64/e_sqrt.c b/sysdeps/ieee754/dbl-64/e_sqrt.c index fff6d148fe..8304a2bb63 100644 --- a/sysdeps/ieee754/dbl-64/e_sqrt.c +++ b/sysdeps/ieee754/dbl-64/e_sqrt.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 @@ -21,7 +21,7 @@ /* */ /* FUNCTION: usqrt */ /* */ -/* FILES NEEDED: dla.h endian.h mydefs.h uroot.h */ +/* FILES NEEDED: dla.h endian.h mydefs.h */ /* uroot.tbl */ /* */ /* An ultimate sqrt routine. Given an IEEE double machine number x */ @@ -47,7 +47,6 @@ double __ieee754_sqrt (double x) { -#include "uroot.h" static const double rt0 = 9.99999999859990725855365213134618E-01, rt1 = 4.99999999495955425917856814202739E-01, @@ -134,7 +133,7 @@ __ieee754_sqrt (double x) return x; /* sqrt(+0)=+0, sqrt(-0)=-0 */ if (k < 0) return (x - x) / (x - x); /* sqrt(-ve)=sNaN */ - return tm256.x * __ieee754_sqrt (x * t512.x); + return 0x1p-256 * __ieee754_sqrt (x * 0x1p512); } } strong_alias (__ieee754_sqrt, __sqrt_finite) diff --git a/sysdeps/ieee754/dbl-64/gamma_product.c b/sysdeps/ieee754/dbl-64/gamma_product.c index 6c0d6210d9..7ae144aeb3 100644 --- a/sysdeps/ieee754/dbl-64/gamma_product.c +++ b/sysdeps/ieee754/dbl-64/gamma_product.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/dbl-64/gamma_productf.c b/sysdeps/ieee754/dbl-64/gamma_productf.c index df59e1beec..58b4e761fe 100644 --- a/sysdeps/ieee754/dbl-64/gamma_productf.c +++ b/sysdeps/ieee754/dbl-64/gamma_productf.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 @@ -36,10 +36,7 @@ __gamma_productf (float x, float x_eps, int n, float *eps) for (int i = 1; i < n; i++) ret *= x_full + i; -#if FLT_EVAL_METHOD != 0 - volatile -#endif - float fret = ret; + float fret = math_narrow_eval ((float) ret); *eps = (ret - fret) / fret; return fret; diff --git a/sysdeps/ieee754/dbl-64/halfulp.c b/sysdeps/ieee754/dbl-64/halfulp.c index 028cb748f6..5e3e731754 100644 --- a/sysdeps/ieee754/dbl-64/halfulp.c +++ b/sysdeps/ieee754/dbl-64/halfulp.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/dbl-64/k_rem_pio2.c b/sysdeps/ieee754/dbl-64/k_rem_pio2.c index 047c6c2886..e58c9e854c 100644 --- a/sysdeps/ieee754/dbl-64/k_rem_pio2.c +++ b/sysdeps/ieee754/dbl-64/k_rem_pio2.c @@ -315,34 +315,25 @@ recompute: break; case 1: case 2:; -#if __FLT_EVAL_METHOD__ != 0 - volatile -#endif double fv = 0.0; for (i = jz; i >= 0; i--) - fv += fq[i]; + fv = math_narrow_eval (fv + fq[i]); y[0] = (ih == 0) ? fv : -fv; - fv = fq[0] - fv; + fv = math_narrow_eval (fq[0] - fv); for (i = 1; i <= jz; i++) - fv += fq[i]; + fv = math_narrow_eval (fv + fq[i]); y[1] = (ih == 0) ? fv : -fv; break; case 3: /* painful */ for (i = jz; i > 0; i--) { -#if __FLT_EVAL_METHOD__ != 0 - volatile -#endif - double fv = (double) (fq[i - 1] + fq[i]); + double fv = math_narrow_eval (fq[i - 1] + fq[i]); fq[i] += fq[i - 1] - fv; fq[i - 1] = fv; } for (i = jz; i > 1; i--) { -#if __FLT_EVAL_METHOD__ != 0 - volatile -#endif - double fv = (double) (fq[i - 1] + fq[i]); + double fv = math_narrow_eval (fq[i - 1] + fq[i]); fq[i] += fq[i - 1] - fv; fq[i - 1] = fv; } diff --git a/sysdeps/ieee754/dbl-64/lgamma_neg.c b/sysdeps/ieee754/dbl-64/lgamma_neg.c new file mode 100644 index 0000000000..fbf992203c --- /dev/null +++ b/sysdeps/ieee754/dbl-64/lgamma_neg.c @@ -0,0 +1,383 @@ +/* lgamma 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 double lgamma_zeros[][2] = + { + { -0x2.74ff92c01f0d8p+0, -0x2.abec9f315f1ap-56 }, + { -0x2.bf6821437b202p+0, 0x6.866a5b4b9be14p-56 }, + { -0x3.24c1b793cb35ep+0, -0xf.b8be699ad3d98p-56 }, + { -0x3.f48e2a8f85fcap+0, -0x1.70d4561291237p-56 }, + { -0x4.0a139e1665604p+0, 0xf.3c60f4f21e7fp-56 }, + { -0x4.fdd5de9bbabf4p+0, 0xa.ef2f55bf89678p-56 }, + { -0x5.021a95fc2db64p+0, -0x3.2a4c56e595394p-56 }, + { -0x5.ffa4bd647d034p+0, -0x1.7dd4ed62cbd32p-52 }, + { -0x6.005ac9625f234p+0, 0x4.9f83d2692e9c8p-56 }, + { -0x6.fff2fddae1bcp+0, 0xc.29d949a3dc03p-60 }, + { -0x7.000cff7b7f87cp+0, 0x1.20bb7d2324678p-52 }, + { -0x7.fffe5fe05673cp+0, -0x3.ca9e82b522b0cp-56 }, + { -0x8.0001a01459fc8p+0, -0x1.f60cb3cec1cedp-52 }, + { -0x8.ffffd1c425e8p+0, -0xf.fc864e9574928p-56 }, + { -0x9.00002e3bb47d8p+0, -0x6.d6d843fedc35p-56 }, + { -0x9.fffffb606bep+0, 0x2.32f9d51885afap-52 }, + { -0xa.0000049f93bb8p+0, -0x1.927b45d95e154p-52 }, + { -0xa.ffffff9466eap+0, 0xe.4c92532d5243p-56 }, + { -0xb.0000006b9915p+0, -0x3.15d965a6ffea4p-52 }, + { -0xb.fffffff708938p+0, -0x7.387de41acc3d4p-56 }, + { -0xc.00000008f76c8p+0, 0x8.cea983f0fdafp-56 }, + { -0xc.ffffffff4f6ep+0, 0x3.09e80685a0038p-52 }, + { -0xd.00000000b092p+0, -0x3.09c06683dd1bap-52 }, + { -0xd.fffffffff3638p+0, 0x3.a5461e7b5c1f6p-52 }, + { -0xe.000000000c9c8p+0, -0x3.a545e94e75ec6p-52 }, + { -0xe.ffffffffff29p+0, 0x3.f9f399fb10cfcp-52 }, + { -0xf.0000000000d7p+0, -0x3.f9f399bd0e42p-52 }, + { -0xf.fffffffffff28p+0, -0xc.060c6621f513p-56 }, + { -0x1.000000000000dp+4, -0x7.3f9f399da1424p-52 }, + { -0x1.0ffffffffffffp+4, -0x3.569c47e7a93e2p-52 }, + { -0x1.1000000000001p+4, 0x3.569c47e7a9778p-52 }, + { -0x1.2p+4, 0xb.413c31dcbecdp-56 }, + { -0x1.2p+4, -0xb.413c31dcbeca8p-56 }, + { -0x1.3p+4, 0x9.7a4da340a0ab8p-60 }, + { -0x1.3p+4, -0x9.7a4da340a0ab8p-60 }, + { -0x1.4p+4, 0x7.950ae90080894p-64 }, + { -0x1.4p+4, -0x7.950ae90080894p-64 }, + { -0x1.5p+4, 0x5.c6e3bdb73d5c8p-68 }, + { -0x1.5p+4, -0x5.c6e3bdb73d5c8p-68 }, + { -0x1.6p+4, 0x4.338e5b6dfe14cp-72 }, + { -0x1.6p+4, -0x4.338e5b6dfe14cp-72 }, + { -0x1.7p+4, 0x2.ec368262c7034p-76 }, + { -0x1.7p+4, -0x2.ec368262c7034p-76 }, + { -0x1.8p+4, 0x1.f2cf01972f578p-80 }, + { -0x1.8p+4, -0x1.f2cf01972f578p-80 }, + { -0x1.9p+4, 0x1.3f3ccdd165fa9p-84 }, + { -0x1.9p+4, -0x1.3f3ccdd165fa9p-84 }, + { -0x1.ap+4, 0xc.4742fe35272dp-92 }, + { -0x1.ap+4, -0xc.4742fe35272dp-92 }, + { -0x1.bp+4, 0x7.46ac70b733a8cp-96 }, + { -0x1.bp+4, -0x7.46ac70b733a8cp-96 }, + { -0x1.cp+4, 0x4.2862898d42174p-100 }, + }; + +static const double e_hi = 0x2.b7e151628aed2p+0, e_lo = 0xa.6abf7158809dp-56; + +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's + approximation to lgamma function. */ + +static const double lgamma_coeff[] = + { + 0x1.5555555555555p-4, + -0xb.60b60b60b60b8p-12, + 0x3.4034034034034p-12, + -0x2.7027027027028p-12, + 0x3.72a3c5631fe46p-12, + -0x7.daac36664f1f4p-12, + 0x1.a41a41a41a41ap-8, + -0x7.90a1b2c3d4e6p-8, + 0x2.dfd2c703c0dp-4, + -0x1.6476701181f3ap+0, + 0xd.672219167003p+0, + -0x9.cd9292e6660d8p+4, + }; + +#define NCOEFF (sizeof (lgamma_coeff) / sizeof (lgamma_coeff[0])) + +/* Polynomial approximations to (|gamma(x)|-1)(x-n)/(x-x0), where n is + the integer end-point of the half-integer interval containing x and + x0 is the zero of lgamma in that half-integer interval. Each + polynomial is expressed in terms of x-xm, where xm is the midpoint + of the interval for which the polynomial applies. */ + +static const double poly_coeff[] = + { + /* Interval [-2.125, -2] (polynomial degree 10). */ + -0x1.0b71c5c54d42fp+0, + -0xc.73a1dc05f3758p-4, + -0x1.ec84140851911p-4, + -0xe.37c9da23847e8p-4, + -0x1.03cd87cdc0ac6p-4, + -0xe.ae9aedce12eep-4, + 0x9.b11a1780cfd48p-8, + -0xe.f25fc460bdebp-4, + 0x2.6e984c61ca912p-4, + -0xf.83fea1c6d35p-4, + 0x4.760c8c8909758p-4, + /* Interval [-2.25, -2.125] (polynomial degree 11). */ + -0xf.2930890d7d678p-4, + -0xc.a5cfde054eaa8p-4, + 0x3.9c9e0fdebd99cp-4, + -0x1.02a5ad35619d9p+0, + 0x9.6e9b1167c164p-4, + -0x1.4d8332eba090ap+0, + 0x1.1c0c94b1b2b6p+0, + -0x1.c9a70d138c74ep+0, + 0x1.d7d9cf1d4c196p+0, + -0x2.91fbf4cd6abacp+0, + 0x2.f6751f74b8ff8p+0, + -0x3.e1bb7b09e3e76p+0, + /* Interval [-2.375, -2.25] (polynomial degree 12). */ + -0xd.7d28d505d618p-4, + -0xe.69649a3040958p-4, + 0xb.0d74a2827cd6p-4, + -0x1.924b09228a86ep+0, + 0x1.d49b12bcf6175p+0, + -0x3.0898bb530d314p+0, + 0x4.207a6be8fda4cp+0, + -0x6.39eef56d4e9p+0, + 0x8.e2e42acbccec8p+0, + -0xd.0d91c1e596a68p+0, + 0x1.2e20d7099c585p+4, + -0x1.c4eb6691b4ca9p+4, + 0x2.96a1a11fd85fep+4, + /* Interval [-2.5, -2.375] (polynomial degree 13). */ + -0xb.74ea1bcfff948p-4, + -0x1.2a82bd590c376p+0, + 0x1.88020f828b81p+0, + -0x3.32279f040d7aep+0, + 0x5.57ac8252ce868p+0, + -0x9.c2aedd093125p+0, + 0x1.12c132716e94cp+4, + -0x1.ea94dfa5c0a6dp+4, + 0x3.66b61abfe858cp+4, + -0x6.0cfceb62a26e4p+4, + 0xa.beeba09403bd8p+4, + -0x1.3188d9b1b288cp+8, + 0x2.37f774dd14c44p+8, + -0x3.fdf0a64cd7136p+8, + /* Interval [-2.625, -2.5] (polynomial degree 13). */ + -0x3.d10108c27ebbp-4, + 0x1.cd557caff7d2fp+0, + 0x3.819b4856d36cep+0, + 0x6.8505cbacfc42p+0, + 0xb.c1b2e6567a4dp+0, + 0x1.50a53a3ce6c73p+4, + 0x2.57adffbb1ec0cp+4, + 0x4.2b15549cf400cp+4, + 0x7.698cfd82b3e18p+4, + 0xd.2decde217755p+4, + 0x1.7699a624d07b9p+8, + 0x2.98ecf617abbfcp+8, + 0x4.d5244d44d60b4p+8, + 0x8.e962bf7395988p+8, + /* Interval [-2.75, -2.625] (polynomial degree 12). */ + -0x6.b5d252a56e8a8p-4, + 0x1.28d60383da3a6p+0, + 0x1.db6513ada89bep+0, + 0x2.e217118fa8c02p+0, + 0x4.450112c651348p+0, + 0x6.4af990f589b8cp+0, + 0x9.2db5963d7a238p+0, + 0xd.62c03647da19p+0, + 0x1.379f81f6416afp+4, + 0x1.c5618b4fdb96p+4, + 0x2.9342d0af2ac4ep+4, + 0x3.d9cdf56d2b186p+4, + 0x5.ab9f91d5a27a4p+4, + /* Interval [-2.875, -2.75] (polynomial degree 11). */ + -0x8.a41b1e4f36ff8p-4, + 0xc.da87d3b69dbe8p-4, + 0x1.1474ad5c36709p+0, + 0x1.761ecb90c8c5cp+0, + 0x1.d279bff588826p+0, + 0x2.4e5d003fb36a8p+0, + 0x2.d575575566842p+0, + 0x3.85152b0d17756p+0, + 0x4.5213d921ca13p+0, + 0x5.55da7dfcf69c4p+0, + 0x6.acef729b9404p+0, + 0x8.483cc21dd0668p+0, + /* Interval [-3, -2.875] (polynomial degree 11). */ + -0xa.046d667e468f8p-4, + 0x9.70b88dcc006cp-4, + 0xa.a8a39421c94dp-4, + 0xd.2f4d1363f98ep-4, + 0xd.ca9aa19975b7p-4, + 0xf.cf09c2f54404p-4, + 0x1.04b1365a9adfcp+0, + 0x1.22b54ef213798p+0, + 0x1.2c52c25206bf5p+0, + 0x1.4aa3d798aace4p+0, + 0x1.5c3f278b504e3p+0, + 0x1.7e08292cc347bp+0, + }; + +static const size_t poly_deg[] = + { + 10, + 11, + 12, + 13, + 13, + 12, + 11, + 11, + }; + +static const size_t poly_end[] = + { + 10, + 22, + 35, + 49, + 63, + 76, + 88, + 100, + }; + +/* Compute sin (pi * X) for -0.25 <= X <= 0.5. */ + +static double +lg_sinpi (double x) +{ + if (x <= 0.25) + return __sin (M_PI * x); + else + return __cos (M_PI * (0.5 - x)); +} + +/* Compute cos (pi * X) for -0.25 <= X <= 0.5. */ + +static double +lg_cospi (double x) +{ + if (x <= 0.25) + return __cos (M_PI * x); + else + return __sin (M_PI * (0.5 - x)); +} + +/* Compute cot (pi * X) for -0.25 <= X <= 0.5. */ + +static double +lg_cotpi (double x) +{ + return lg_cospi (x) / lg_sinpi (x); +} + +/* Compute lgamma of a negative argument -28 < X < -2, setting + *SIGNGAMP accordingly. */ + +double +__lgamma_neg (double x, int *signgamp) +{ + /* Determine the half-integer region X lies in, handle exact + integers and determine the sign of the result. */ + int i = __floor (-2 * x); + if ((i & 1) == 0 && i == -2 * x) + return 1.0 / 0.0; + double xn = ((i & 1) == 0 ? -i / 2 : (-i - 1) / 2); + i -= 4; + *signgamp = ((i & 2) == 0 ? -1 : 1); + + SET_RESTORE_ROUND (FE_TONEAREST); + + /* Expand around the zero X0 = X0_HI + X0_LO. */ + double x0_hi = lgamma_zeros[i][0], x0_lo = lgamma_zeros[i][1]; + double xdiff = x - x0_hi - x0_lo; + + /* For arguments in the range -3 to -2, use polynomial + approximations to an adjusted version of the gamma function. */ + if (i < 2) + { + int j = __floor (-8 * x) - 16; + double xm = (-33 - 2 * j) * 0.0625; + double x_adj = x - xm; + size_t deg = poly_deg[j]; + size_t end = poly_end[j]; + double g = poly_coeff[end]; + for (size_t j = 1; j <= deg; j++) + g = g * x_adj + poly_coeff[end - j]; + return __log1p (g * xdiff / (x - xn)); + } + + /* The result we want is log (sinpi (X0) / sinpi (X)) + + log (gamma (1 - X0) / gamma (1 - X)). */ + double x_idiff = fabs (xn - x), x0_idiff = fabs (xn - x0_hi - x0_lo); + double log_sinpi_ratio; + if (x0_idiff < x_idiff * 0.5) + /* Use log not log1p to avoid inaccuracy from log1p of arguments + close to -1. */ + log_sinpi_ratio = __ieee754_log (lg_sinpi (x0_idiff) + / lg_sinpi (x_idiff)); + else + { + /* Use log1p not log to avoid inaccuracy from log of arguments + close to 1. X0DIFF2 has positive sign if X0 is further from + XN than X is from XN, negative sign otherwise. */ + double x0diff2 = ((i & 1) == 0 ? xdiff : -xdiff) * 0.5; + double sx0d2 = lg_sinpi (x0diff2); + double cx0d2 = lg_cospi (x0diff2); + log_sinpi_ratio = __log1p (2 * sx0d2 + * (-sx0d2 + cx0d2 * lg_cotpi (x_idiff))); + } + + double log_gamma_ratio; + double y0 = math_narrow_eval (1 - x0_hi); + double y0_eps = -x0_hi + (1 - y0) - x0_lo; + double y = math_narrow_eval (1 - x); + double y_eps = -x + (1 - y); + /* We now wish to compute LOG_GAMMA_RATIO + = log (gamma (Y0 + Y0_EPS) / gamma (Y + Y_EPS)). XDIFF + accurately approximates the difference Y0 + Y0_EPS - Y - + Y_EPS. Use Stirling's approximation. First, we may need to + adjust into the range where Stirling's approximation is + sufficiently accurate. */ + double log_gamma_adj = 0; + if (i < 6) + { + int n_up = (7 - i) / 2; + double ny0, ny0_eps, ny, ny_eps; + ny0 = math_narrow_eval (y0 + n_up); + ny0_eps = y0 - (ny0 - n_up) + y0_eps; + y0 = ny0; + y0_eps = ny0_eps; + ny = math_narrow_eval (y + n_up); + ny_eps = y - (ny - n_up) + y_eps; + y = ny; + y_eps = ny_eps; + double prodm1 = __lgamma_product (xdiff, y - n_up, y_eps, n_up); + log_gamma_adj = -__log1p (prodm1); + } + double log_gamma_high + = (xdiff * __log1p ((y0 - e_hi - e_lo + y0_eps) / e_hi) + + (y - 0.5 + y_eps) * __log1p (xdiff / y) + log_gamma_adj); + /* Compute the sum of (B_2k / 2k(2k-1))(Y0^-(2k-1) - Y^-(2k-1)). */ + double y0r = 1 / y0, yr = 1 / y; + double y0r2 = y0r * y0r, yr2 = yr * yr; + double rdiff = -xdiff / (y * y0); + double bterm[NCOEFF]; + double dlast = rdiff, elast = rdiff * yr * (yr + y0r); + bterm[0] = dlast * lgamma_coeff[0]; + for (size_t j = 1; j < NCOEFF; j++) + { + double dnext = dlast * y0r2 + elast; + double enext = elast * yr2; + bterm[j] = dnext * lgamma_coeff[j]; + dlast = dnext; + elast = enext; + } + double log_gamma_low = 0; + for (size_t j = 0; j < NCOEFF; j++) + log_gamma_low += bterm[NCOEFF - 1 - j]; + log_gamma_ratio = log_gamma_high + log_gamma_low; + + return log_sinpi_ratio + log_gamma_ratio; +} diff --git a/sysdeps/ieee754/dbl-64/lgamma_product.c b/sysdeps/ieee754/dbl-64/lgamma_product.c new file mode 100644 index 0000000000..d956575bc7 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/lgamma_product.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 (double *hi, double *lo, double x, double y) +{ +#ifdef __FP_FAST_FMA + /* Fast built-in fused multiply-add. */ + *hi = x * y; + *lo = __builtin_fma (x, y, -*hi); +#elif defined FP_FAST_FMA + /* Fast library fused multiply-add, compiler before GCC 4.6. */ + *hi = x * y; + *lo = __fma (x, y, -*hi); +#else + /* Apply Dekker's algorithm. */ + *hi = x * y; +# define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1) + double x1 = x * C; + double y1 = y * C; +# undef C + x1 = (x - x1) + x1; + y1 = (y - y1) + y1; + double x2 = x - x1; + 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. */ + +double +__lgamma_product (double t, double x, double x_eps, int n) +{ + double ret = 0, ret_eps = 0; + for (int i = 0; i < n; i++) + { + double xi = x + i; + double quot = t / xi; + double mhi, mlo; + mul_split (&mhi, &mlo, quot, xi); + double quot_lo = (t - mhi - mlo) / xi - t * x_eps / (xi * xi); + /* We want (1 + RET + RET_EPS) * (1 + QUOT + QUOT_LO) - 1. */ + double rhi, rlo; + mul_split (&rhi, &rlo, ret, quot); + double rpq = ret + quot; + double rpq_eps = (ret - rpq) + quot; + double nret = rpq + rhi; + double nret_eps = (rpq - nret) + rhi; + ret_eps += (rpq_eps + nret_eps + rlo + ret_eps * quot + + quot_lo + quot_lo * (ret + ret_eps)); + ret = nret; + } + return ret + ret_eps; +} diff --git a/sysdeps/ieee754/dbl-64/mpa-arch.h b/sysdeps/ieee754/dbl-64/mpa-arch.h index 0940f0bb57..ea33a14219 100644 --- a/sysdeps/ieee754/dbl-64/mpa-arch.h +++ b/sysdeps/ieee754/dbl-64/mpa-arch.h @@ -1,5 +1,5 @@ /* Overridable constants and operations. - Copyright (C) 2013-2015 Free Software Foundation, Inc. + Copyright (C) 2013-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/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c index 7b52da91d5..4b21d1d2e6 100644 --- a/sysdeps/ieee754/dbl-64/mpa.c +++ b/sysdeps/ieee754/dbl-64/mpa.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/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h index 47dd6c457c..b8ca297eb7 100644 --- a/sysdeps/ieee754/dbl-64/mpa.h +++ b/sysdeps/ieee754/dbl-64/mpa.h @@ -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/dbl-64/mpatan.c b/sysdeps/ieee754/dbl-64/mpatan.c index d9ac3d62b1..dbd49b2bee 100644 --- a/sysdeps/ieee754/dbl-64/mpatan.c +++ b/sysdeps/ieee754/dbl-64/mpatan.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/dbl-64/mpatan.h b/sysdeps/ieee754/dbl-64/mpatan.h index 2ceb966998..36c02602ac 100644 --- a/sysdeps/ieee754/dbl-64/mpatan.h +++ b/sysdeps/ieee754/dbl-64/mpatan.h @@ -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/dbl-64/mpatan2.c b/sysdeps/ieee754/dbl-64/mpatan2.c index 78f01d9839..e7de6bc158 100644 --- a/sysdeps/ieee754/dbl-64/mpatan2.c +++ b/sysdeps/ieee754/dbl-64/mpatan2.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/dbl-64/mpexp.c b/sysdeps/ieee754/dbl-64/mpexp.c index b5a1def93d..f17baf2139 100644 --- a/sysdeps/ieee754/dbl-64/mpexp.c +++ b/sysdeps/ieee754/dbl-64/mpexp.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/dbl-64/mplog.c b/sysdeps/ieee754/dbl-64/mplog.c index 53e1b57621..b297153b10 100644 --- a/sysdeps/ieee754/dbl-64/mplog.c +++ b/sysdeps/ieee754/dbl-64/mplog.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/dbl-64/mpn2dbl.c b/sysdeps/ieee754/dbl-64/mpn2dbl.c index c0cf8e238d..1050ec3b0e 100644 --- a/sysdeps/ieee754/dbl-64/mpn2dbl.c +++ b/sysdeps/ieee754/dbl-64/mpn2dbl.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/dbl-64/mpsqrt.c b/sysdeps/ieee754/dbl-64/mpsqrt.c index f3ef55eb7f..a77d3938a9 100644 --- a/sysdeps/ieee754/dbl-64/mpsqrt.c +++ b/sysdeps/ieee754/dbl-64/mpsqrt.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/dbl-64/mpsqrt.h b/sysdeps/ieee754/dbl-64/mpsqrt.h index e322398a7a..0f1cf58fb2 100644 --- a/sysdeps/ieee754/dbl-64/mpsqrt.h +++ b/sysdeps/ieee754/dbl-64/mpsqrt.h @@ -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/dbl-64/mptan.c b/sysdeps/ieee754/dbl-64/mptan.c index 0d58d9c786..280c4909d6 100644 --- a/sysdeps/ieee754/dbl-64/mptan.c +++ b/sysdeps/ieee754/dbl-64/mptan.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/dbl-64/mydefs.h b/sysdeps/ieee754/dbl-64/mydefs.h index 31585142b7..b1d67d4445 100644 --- a/sysdeps/ieee754/dbl-64/mydefs.h +++ b/sysdeps/ieee754/dbl-64/mydefs.h @@ -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/dbl-64/powtwo.tbl b/sysdeps/ieee754/dbl-64/powtwo.tbl index 2406ee5d65..2d828f6022 100644 --- a/sysdeps/ieee754/dbl-64/powtwo.tbl +++ b/sysdeps/ieee754/dbl-64/powtwo.tbl @@ -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/dbl-64/root.tbl b/sysdeps/ieee754/dbl-64/root.tbl index cb0c8b11ec..61c859248c 100644 --- a/sysdeps/ieee754/dbl-64/root.tbl +++ b/sysdeps/ieee754/dbl-64/root.tbl @@ -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/dbl-64/s_asinh.c b/sysdeps/ieee754/dbl-64/s_asinh.c index ebe471015b..9193301b5e 100644 --- a/sysdeps/ieee754/dbl-64/s_asinh.c +++ b/sysdeps/ieee754/dbl-64/s_asinh.c @@ -39,11 +39,7 @@ __asinh (double x) ix = hx & 0x7fffffff; if (__glibc_unlikely (ix < 0x3e300000)) /* |x|<2**-28 */ { - if (fabs (x) < DBL_MIN) - { - 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/dbl-64/s_atan.c b/sysdeps/ieee754/dbl-64/s_atan.c index 5035ae87bc..780c3ff17a 100644 --- a/sysdeps/ieee754/dbl-64/s_atan.c +++ b/sysdeps/ieee754/dbl-64/s_atan.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 @@ -90,11 +90,7 @@ atan (double x) { if (u < A) { - if (u < DBL_MIN) - { - double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow_nonneg (u); return x; } else diff --git a/sysdeps/ieee754/dbl-64/s_cbrt.c b/sysdeps/ieee754/dbl-64/s_cbrt.c index a728968a1d..647f30bc71 100644 --- a/sysdeps/ieee754/dbl-64/s_cbrt.c +++ b/sysdeps/ieee754/dbl-64/s_cbrt.c @@ -1,5 +1,5 @@ /* Compute cubic root of double value. - 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 Dirk Alboth <dirka@uni-paderborn.de> and Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/dbl-64/s_erf.c b/sysdeps/ieee754/dbl-64/s_erf.c index ea0a73424e..b4975a8af8 100644 --- a/sysdeps/ieee754/dbl-64/s_erf.c +++ b/sysdeps/ieee754/dbl-64/s_erf.c @@ -116,6 +116,7 @@ static char rcsid[] = "$NetBSD: s_erf.c,v 1.8 1995/05/10 20:47:05 jtc Exp $"; #include <float.h> #include <math.h> #include <math_private.h> +#include <fix-int-fp-convert-zero.h> static const double tiny = 1e-300, @@ -213,11 +214,7 @@ __erf (double x) { /* Avoid spurious underflow. */ double ret = 0.0625 * (16.0 * x + (16.0 * efx) * x); - if (fabs (ret) < DBL_MIN) - { - double force_underflow = ret * ret; - math_force_eval (force_underflow); - } + math_check_force_underflow (ret); return ret; } return x + efx * x; @@ -312,7 +309,10 @@ __erfc (double x) ix = hx & 0x7fffffff; if (ix >= 0x7ff00000) /* erfc(nan)=nan */ { /* erfc(+-inf)=0,2 */ - return (double) (((u_int32_t) hx >> 31) << 1) + one / x; + double ret = (double) (((u_int32_t) hx >> 31) << 1) + one / x; + if (FIX_INT_FP_CONVERT_ZERO && ret == 0.0) + return 0.0; + return ret; } if (ix < 0x3feb0000) /* |x|<0.84375 */ @@ -402,10 +402,7 @@ __erfc (double x) __ieee754_exp ((z - x) * (z + x) + R / S); if (hx > 0) { -#if FLT_EVAL_METHOD != 0 - volatile -#endif - double ret = r / x; + double ret = math_narrow_eval (r / x); if (ret == 0) __set_errno (ERANGE); return ret; diff --git a/sysdeps/ieee754/dbl-64/s_expm1.c b/sysdeps/ieee754/dbl-64/s_expm1.c index 41ef63a786..54d771007a 100644 --- a/sysdeps/ieee754/dbl-64/s_expm1.c +++ b/sysdeps/ieee754/dbl-64/s_expm1.c @@ -195,11 +195,7 @@ __expm1 (double x) } else if (hx < 0x3c900000) /* when |x|<2**-54, return x */ { - if (fabs (x) < DBL_MIN) - { - double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); t = huge + x; /* return x with inexact flags when x!=0 */ return x - (t - (huge + x)); } diff --git a/sysdeps/ieee754/dbl-64/s_finite.c b/sysdeps/ieee754/dbl-64/s_finite.c index 49986bbde9..69141db75d 100644 --- a/sysdeps/ieee754/dbl-64/s_finite.c +++ b/sysdeps/ieee754/dbl-64/s_finite.c @@ -21,6 +21,7 @@ static char rcsid[] = "$NetBSD: s_finite.c,v 1.8 1995/05/10 20:47:17 jtc Exp $"; #include <math.h> #include <math_private.h> +#include <shlib-compat.h> #undef __finite @@ -32,11 +33,18 @@ int FINITE(double x) { int32_t hx; GET_HIGH_WORD (hx, x); - return (int) ((u_int32_t) ((hx & 0x7fffffff) - 0x7ff00000) >> 31); + return (int) ((u_int32_t) ((hx & 0x7ff00000) - 0x7ff00000) >> 31); } hidden_def (__finite) weak_alias (__finite, finite) #ifdef NO_LONG_DOUBLE -strong_alias (__finite, __finitel) +# ifdef LDBL_CLASSIFY_COMPAT +# if SHLIB_COMPAT (libc, GLIBC_2_0, GLIBC_2_23) +compat_symbol (libc, __finite, __finitel, GLIBC_2_0); +# endif +# if SHLIB_COMPAT (libm, GLIBC_2_1, GLIBC_2_23) +compat_symbol (libm, __finite, __finitel, GLIBC_2_1); +# endif +# endif weak_alias (__finite, finitel) #endif diff --git a/sysdeps/ieee754/dbl-64/s_fma.c b/sysdeps/ieee754/dbl-64/s_fma.c index 716b41273d..a3492434e4 100644 --- a/sysdeps/ieee754/dbl-64/s_fma.c +++ b/sysdeps/ieee754/dbl-64/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. @@ -66,7 +66,7 @@ __fma (double x, double y, double z) /* If fma will certainly overflow, compute as x * y. */ if (u.ieee.exponent + v.ieee.exponent > 0x7ff + IEEE754_DOUBLE_BIAS) return x * y; - /* If x * y is less than 1/4 of DBL_DENORM_MIN, neither the + /* If x * y is less than 1/4 of DBL_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 @@ -90,8 +90,8 @@ __fma (double x, double y, double z) && w.ieee.mantissa1 == 0 && w.ieee.mantissa0 == 0))) { - volatile double force_underflow = x * y; - (void) force_underflow; + double force_underflow = x * y; + math_force_eval (force_underflow); } return v.d * 0x1p-54; } @@ -117,7 +117,7 @@ __fma (double x, double y, double z) very small, adjust them up to avoid spurious underflows, rather than down. */ if (u.ieee.exponent + v.ieee.exponent - <= IEEE754_DOUBLE_BIAS + DBL_MANT_DIG) + <= IEEE754_DOUBLE_BIAS + 2 * DBL_MANT_DIG) { if (u.ieee.exponent > v.ieee.exponent) u.ieee.exponent += 2 * DBL_MANT_DIG + 2; @@ -175,7 +175,10 @@ __fma (double x, double y, 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; libc_feholdexcept_setround (&env, FE_TONEAREST); diff --git a/sysdeps/ieee754/dbl-64/s_fmaf.c b/sysdeps/ieee754/dbl-64/s_fmaf.c index 86b8662dbf..582c205572 100644 --- a/sysdeps/ieee754/dbl-64/s_fmaf.c +++ b/sysdeps/ieee754/dbl-64/s_fmaf.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/dbl-64/s_fpclassify.c b/sysdeps/ieee754/dbl-64/s_fpclassify.c index 2e8751277d..176e29fccb 100644 --- a/sysdeps/ieee754/dbl-64/s_fpclassify.c +++ b/sysdeps/ieee754/dbl-64/s_fpclassify.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. diff --git a/sysdeps/ieee754/dbl-64/s_isinf.c b/sysdeps/ieee754/dbl-64/s_isinf.c index 46a7266e7f..c0ad54538a 100644 --- a/sysdeps/ieee754/dbl-64/s_isinf.c +++ b/sysdeps/ieee754/dbl-64/s_isinf.c @@ -15,6 +15,7 @@ static char rcsid[] = "$NetBSD: s_isinf.c,v 1.3 1995/05/11 23:20:14 jtc Exp $"; #include <math.h> #include <math_private.h> +#include <shlib-compat.h> int __isinf (double x) @@ -28,6 +29,8 @@ __isinf (double x) hidden_def (__isinf) weak_alias (__isinf, isinf) #ifdef NO_LONG_DOUBLE -strong_alias (__isinf, __isinfl) +# if defined LDBL_CLASSIFY_COMPAT && SHLIB_COMPAT (libc, GLIBC_2_0, GLIBC_2_23) +compat_symbol (libc, __isinf, __isinfl, GLIBC_2_0); +# endif weak_alias (__isinf, isinfl) #endif diff --git a/sysdeps/ieee754/dbl-64/s_isinf_ns.c b/sysdeps/ieee754/dbl-64/s_isinf_ns.c deleted file mode 100644 index d8142bcf6e..0000000000 --- a/sysdeps/ieee754/dbl-64/s_isinf_ns.c +++ /dev/null @@ -1,20 +0,0 @@ -/* - * Written by Ulrich Drepper <drepper@gmail.com>. - */ - -/* - * __isinf_ns(x) returns != 0 if x is ±inf, else 0; - * no branching! - */ - -#include <math.h> -#include <math_private.h> - -#undef __isinf_ns -int -__isinf_ns (double x) -{ - int32_t hx, lx; - EXTRACT_WORDS (hx, lx, x); - return !(lx | ((hx & 0x7fffffff) ^ 0x7ff00000)); -} diff --git a/sysdeps/ieee754/dbl-64/s_isnan.c b/sysdeps/ieee754/dbl-64/s_isnan.c index 5d9f31be20..2174d988d8 100644 --- a/sysdeps/ieee754/dbl-64/s_isnan.c +++ b/sysdeps/ieee754/dbl-64/s_isnan.c @@ -21,6 +21,7 @@ static char rcsid[] = "$NetBSD: s_isnan.c,v 1.8 1995/05/10 20:47:36 jtc Exp $"; #include <math.h> #include <math_private.h> +#include <shlib-compat.h> #undef __isnan int @@ -36,6 +37,8 @@ __isnan (double x) hidden_def (__isnan) weak_alias (__isnan, isnan) #ifdef NO_LONG_DOUBLE -strong_alias (__isnan, __isnanl) +# if defined LDBL_CLASSIFY_COMPAT && SHLIB_COMPAT (libc, GLIBC_2_0, GLIBC_2_23) +compat_symbol (libc, __isnan, __isnanl, GLIBC_2_0); +# endif weak_alias (__isnan, isnanl) #endif diff --git a/sysdeps/ieee754/dbl-64/s_issignaling.c b/sysdeps/ieee754/dbl-64/s_issignaling.c index 925d85fb21..4b93d6ae60 100644 --- a/sysdeps/ieee754/dbl-64/s_issignaling.c +++ b/sysdeps/ieee754/dbl-64/s_issignaling.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/dbl-64/s_llrint.c b/sysdeps/ieee754/dbl-64/s_llrint.c index 62acb67c38..ff40ee626d 100644 --- a/sysdeps/ieee754/dbl-64/s_llrint.c +++ b/sysdeps/ieee754/dbl-64/s_llrint.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. @@ -18,9 +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> static const double two52[2] = { @@ -35,7 +38,7 @@ __llrint (double x) int32_t j0; u_int32_t i1, i0; long long int result; - volatile double w; + double w; double t; int sx; @@ -47,7 +50,7 @@ __llrint (double x) if (j0 < 20) { - w = two52[sx] + x; + w = math_narrow_eval (two52[sx] + x); t = w - two52[sx]; EXTRACT_WORDS (i0, i1, t); j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; @@ -62,7 +65,7 @@ __llrint (double x) result = (((long long int) i0 << 32) | i1) << (j0 - 52); else { - w = two52[sx] + x; + w = math_narrow_eval (two52[sx] + x); t = w - two52[sx]; EXTRACT_WORDS (i0, i1, t); j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; @@ -77,8 +80,16 @@ __llrint (double x) } else { - /* The number is too large. It is left implementation defined - what happens. */ +#ifdef FE_INVALID + /* The number is too large. Unless it rounds to LLONG_MIN, + FE_INVALID must be raised and the return value is + unspecified. */ + if (FIX_DBL_LLONG_CONVERT_OVERFLOW && x != (double) LLONG_MIN) + { + feraiseexcept (FE_INVALID); + return sx == 0 ? LLONG_MAX : LLONG_MIN; + } +#endif return (long long int) x; } diff --git a/sysdeps/ieee754/dbl-64/s_llround.c b/sysdeps/ieee754/dbl-64/s_llround.c index c8b3c190b5..1d6a5205f1 100644 --- a/sysdeps/ieee754/dbl-64/s_llround.c +++ b/sysdeps/ieee754/dbl-64/s_llround.c @@ -1,5 +1,5 @@ /* Round 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. @@ -17,9 +17,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 @@ -65,8 +68,16 @@ __llround (double x) } else { - /* The number is too large. It is left implementation defined - what happens. */ +#ifdef FE_INVALID + /* The number is too large. Unless it rounds to LLONG_MIN, + FE_INVALID must be raised and the return value is + unspecified. */ + if (FIX_DBL_LLONG_CONVERT_OVERFLOW && x != (double) LLONG_MIN) + { + feraiseexcept (FE_INVALID); + return sign == 1 ? LLONG_MAX : LLONG_MIN; + } +#endif return (long long int) x; } diff --git a/sysdeps/ieee754/dbl-64/s_log1p.c b/sysdeps/ieee754/dbl-64/s_log1p.c index cff555b0aa..340f6377f7 100644 --- a/sysdeps/ieee754/dbl-64/s_log1p.c +++ b/sysdeps/ieee754/dbl-64/s_log1p.c @@ -120,11 +120,7 @@ __log1p (double x) math_force_eval (two54 + x); /* raise inexact */ if (ax < 0x3c900000) /* |x| < 2**-54 */ { - if (fabs (x) < DBL_MIN) - { - double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); return x; } else diff --git a/sysdeps/ieee754/dbl-64/s_logb.c b/sysdeps/ieee754/dbl-64/s_logb.c index 7a6c49abf5..3a26b18f78 100644 --- a/sysdeps/ieee754/dbl-64/s_logb.c +++ b/sysdeps/ieee754/dbl-64/s_logb.c @@ -18,6 +18,7 @@ #include <math.h> #include <math_private.h> +#include <fix-int-fp-convert-zero.h> double __logb (double x) @@ -41,6 +42,8 @@ __logb (double x) ma = __builtin_clz (ix); rix -= ma - 12; } + if (FIX_INT_FP_CONVERT_ZERO && rix == 1023) + return 0.0; return (double) (rix - 1023); } weak_alias (__logb, logb) diff --git a/sysdeps/ieee754/dbl-64/s_lrint.c b/sysdeps/ieee754/dbl-64/s_lrint.c index 9cce37f3cf..4e7cbef97b 100644 --- a/sysdeps/ieee754/dbl-64/s_lrint.c +++ b/sysdeps/ieee754/dbl-64/s_lrint.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. @@ -18,9 +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> static const double two52[2] = { @@ -34,7 +37,7 @@ __lrint (double x) { int32_t j0; u_int32_t i0, i1; - volatile double w; + double w; double t; long int result; int sx; @@ -47,7 +50,7 @@ __lrint (double x) if (j0 < 20) { - w = two52[sx] + x; + w = math_narrow_eval (two52[sx] + x); t = w - two52[sx]; EXTRACT_WORDS (i0, i1, t); j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; @@ -59,11 +62,25 @@ __lrint (double x) else if (j0 < (int32_t) (8 * sizeof (long int)) - 1) { if (j0 >= 52) - result = ((long int) i0 << (j0 - 20)) | (i1 << (j0 - 52)); + result = ((long int) i0 << (j0 - 20)) | ((long int) i1 << (j0 - 52)); else { - w = two52[sx] + x; - t = w - two52[sx]; +#if defined FE_INVALID || defined FE_INEXACT + /* X < LONG_MAX + 1 implied by J0 < 31. */ + if (sizeof (long int) == 4 + && x > (double) LONG_MAX) + { + /* In the event of overflow we must raise the "invalid" + exception, but not "inexact". */ + t = __nearbyint (x); + feraiseexcept (t == LONG_MAX ? FE_INEXACT : FE_INVALID); + } + else +#endif + { + w = math_narrow_eval (two52[sx] + x); + t = w - two52[sx]; + } EXTRACT_WORDS (i0, i1, t); j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; i0 &= 0xfffff; @@ -77,8 +94,26 @@ __lrint (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 (sizeof (long int) == 4 + && x < (double) LONG_MIN + && x > (double) LONG_MIN - 1.0) + { + /* If truncation produces LONG_MIN, the cast will not raise + the exception, but may raise "inexact". */ + t = __nearbyint (x); + feraiseexcept (t == LONG_MIN ? FE_INEXACT : FE_INVALID); + return LONG_MIN; + } + else if (FIX_DBL_LONG_CONVERT_OVERFLOW && x != (double) LONG_MIN) + { + feraiseexcept (FE_INVALID); + return sx == 0 ? LONG_MAX : LONG_MIN; + } +#endif return (long int) x; } diff --git a/sysdeps/ieee754/dbl-64/s_lround.c b/sysdeps/ieee754/dbl-64/s_lround.c index bdc838a676..aa1275ddc3 100644 --- a/sysdeps/ieee754/dbl-64/s_lround.c +++ b/sysdeps/ieee754/dbl-64/s_lround.c @@ -1,5 +1,5 @@ /* Round 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. @@ -17,9 +17,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 @@ -60,13 +63,43 @@ __lround (double x) if (j0 == 20) result = (long int) i0; else - result = ((long int) i0 << (j0 - 20)) | (j >> (52 - j0)); + { + result = ((long int) i0 << (j0 - 20)) | (j >> (52 - 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 { - /* 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. */ +#ifdef FE_INVALID + if (FIX_DBL_LONG_CONVERT_OVERFLOW + && !(sign == -1 + && (sizeof (long int) == 4 + ? x > (double) LONG_MIN - 0.5 + : x >= (double) LONG_MIN))) + { + feraiseexcept (FE_INVALID); + return sign == 1 ? LONG_MAX : LONG_MIN; + } + else if (!FIX_DBL_LONG_CONVERT_OVERFLOW + && sizeof (long int) == 4 + && x <= (double) LONG_MIN - 0.5) + { + /* If truncation produces LONG_MIN, the cast will not raise + the exception, but may raise "inexact". */ + feraiseexcept (FE_INVALID); + return LONG_MIN; + } +#endif return (long int) x; } diff --git a/sysdeps/ieee754/dbl-64/s_remquo.c b/sysdeps/ieee754/dbl-64/s_remquo.c index b081d26b0d..16215e1e29 100644 --- a/sysdeps/ieee754/dbl-64/s_remquo.c +++ b/sysdeps/ieee754/dbl-64/s_remquo.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. diff --git a/sysdeps/ieee754/dbl-64/s_round.c b/sysdeps/ieee754/dbl-64/s_round.c index 770f8e6386..aeed6faf0e 100644 --- a/sysdeps/ieee754/dbl-64/s_round.c +++ b/sysdeps/ieee754/dbl-64/s_round.c @@ -1,5 +1,5 @@ /* Round 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. diff --git a/sysdeps/ieee754/dbl-64/s_scalbn.c b/sysdeps/ieee754/dbl-64/s_scalbn.c index 0f58034df5..58c7e1b33a 100644 --- a/sysdeps/ieee754/dbl-64/s_scalbn.c +++ b/sysdeps/ieee754/dbl-64/s_scalbn.c @@ -58,8 +58,6 @@ __scalbn (double x, int n) SET_HIGH_WORD (x, (hx & 0x800fffff) | (k << 20)); return x * twom54; } -weak_alias (__scalbn, scalbn) #ifdef NO_LONG_DOUBLE strong_alias (__scalbn, __scalbnl) -weak_alias (__scalbn, scalbnl) #endif diff --git a/sysdeps/ieee754/dbl-64/s_signbit.c b/sysdeps/ieee754/dbl-64/s_signbit.c index 764f11a2d2..f4ef4d238f 100644 --- a/sysdeps/ieee754/dbl-64/s_signbit.c +++ b/sysdeps/ieee754/dbl-64/s_signbit.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 __signbit (double x) { - int32_t hx; - - GET_HIGH_WORD (hx, x); - return hx & 0x80000000; + return __builtin_signbit (x); } diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c index eff120e88d..ca2532fb63 100644 --- a/sysdeps/ieee754/dbl-64/s_sin.c +++ b/sysdeps/ieee754/dbl-64/s_sin.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 @@ -32,9 +32,6 @@ /* bsloww1 */ /* bsloww2 */ /* cslow2 */ -/* csloww */ -/* csloww1 */ -/* csloww2 */ /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h usncs.h */ /* branred.c sincos32.c dosincos.c mpa.c */ /* sincos.tbl */ @@ -135,17 +132,14 @@ double __mpcos (double x, double dx, bool reduce_range); static double slow (double x); static double slow1 (double x); static double slow2 (double x); -static double sloww (double x, double dx, double orig); -static double sloww1 (double x, double dx, double orig, int m); +static double sloww (double x, double dx, double orig, int n); +static double sloww1 (double x, double dx, double orig, int m, int n); static double sloww2 (double x, double dx, double orig, int n); static double bsloww (double x, double dx, double orig, int n); static double bsloww1 (double x, double dx, double orig, int n); static double bsloww2 (double x, double dx, double orig, int n); int __branred (double x, double *a, double *aa); static double cslow2 (double x); -static double csloww (double x, double dx, double orig); -static double csloww1 (double x, double dx, double orig, int m); -static double csloww2 (double x, double dx, double orig, int n); /* Given a number partitioned into U and X such that U is an index into the sin/cos table, this macro computes the cosine of the number by combining @@ -276,32 +270,216 @@ reduce_and_compute (double x, unsigned int k) return retval; } +static inline int4 +__always_inline +reduce_sincos_1 (double x, double *a, double *da) +{ + mynumber v; + + double t = (x * hpinv + toint); + double xn = t - toint; + v.x = t; + double y = (x - xn * mp1) - xn * mp2; + int4 n = v.i[LOW_HALF] & 3; + double db = xn * mp3; + double b = y - db; + db = (y - b) - db; + + *a = b; + *da = db; + + return n; +} + +/* Compute sin (A + DA). cos can be computed by shifting the quadrant N + clockwise. */ +static double +__always_inline +do_sincos_1 (double a, double da, double x, int4 n, int4 k) +{ + double xx, retval, res, cor, y; + mynumber u; + int m; + double eps = fabs (x) * 1.2e-30; + + int k1 = (n + k) & 3; + switch (k1) + { /* quarter of unit circle */ + case 2: + a = -a; + da = -da; + case 0: + xx = a * a; + if (xx < 0.01588) + { + /* Taylor series. */ + res = TAYLOR_SIN (xx, a, da, cor); + cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; + retval = (res == res + cor) ? res : sloww (a, da, x, k); + } + else + { + if (a > 0) + m = 1; + else + { + m = 0; + a = -a; + da = -da; + } + u.x = big + a; + y = a - (u.x - big); + res = do_sin (u, y, da, &cor); + cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; + retval = ((res == res + cor) ? ((m) ? res : -res) + : sloww1 (a, da, x, m, k)); + } + break; + + case 1: + case 3: + if (a < 0) + { + a = -a; + da = -da; + } + u.x = big + a; + y = a - (u.x - big) + da; + res = do_cos (u, y, &cor); + cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; + retval = ((res == res + cor) ? ((k1 & 2) ? -res : res) + : sloww2 (a, da, x, n)); + break; + } + + return retval; +} + +static inline int4 +__always_inline +reduce_sincos_2 (double x, double *a, double *da) +{ + mynumber v; + + double t = (x * hpinv + toint); + double xn = t - toint; + v.x = t; + double xn1 = (xn + 8.0e22) - 8.0e22; + double xn2 = xn - xn1; + double y = ((((x - xn1 * mp1) - xn1 * mp2) - xn2 * mp1) - xn2 * mp2); + int4 n = v.i[LOW_HALF] & 3; + double db = xn1 * pp3; + t = y - db; + db = (y - t) - db; + db = (db - xn2 * pp3) - xn * pp4; + double b = t + db; + db = (t - b) + db; + + *a = b; + *da = db; + + return n; +} + +/* Compute sin (A + DA). cos can be computed by shifting the quadrant N + clockwise. */ +static double +__always_inline +do_sincos_2 (double a, double da, double x, int4 n, int4 k) +{ + double res, retval, cor, xx; + mynumber u; + + double eps = 1.0e-24; + + k = (n + k) & 3; + + switch (k) + { + case 2: + a = -a; + da = -da; + /* Fall through. */ + case 0: + xx = a * a; + if (xx < 0.01588) + { + /* Taylor series. */ + res = TAYLOR_SIN (xx, a, da, cor); + cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; + retval = (res == res + cor) ? res : bsloww (a, da, x, n); + } + else + { + double t, db, y; + int m; + if (a > 0) + { + m = 1; + t = a; + db = da; + } + else + { + m = 0; + t = -a; + db = -da; + } + u.x = big + t; + y = t - (u.x - big); + res = do_sin (u, y, db, &cor); + cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; + retval = ((res == res + cor) ? ((m) ? res : -res) + : bsloww1 (a, da, x, n)); + } + break; + + case 1: + case 3: + if (a < 0) + { + a = -a; + da = -da; + } + u.x = big + a; + double y = a - (u.x - big) + da; + res = do_cos (u, y, &cor); + cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; + retval = ((res == res + cor) ? ((n & 2) ? -res : res) + : bsloww2 (a, da, x, n)); + break; + } + + return retval; +} + /*******************************************************************/ /* An ultimate sin routine. Given an IEEE double machine number x */ /* it computes the correctly rounded (to nearest) value of sin(x) */ /*******************************************************************/ +#ifdef IN_SINCOS +static double +#else double SECTION +#endif __sin (double x) { - double xx, res, t, cor, y, s, c, sn, ssn, cs, ccs, xn, a, da, db, eps, xn1, - xn2; - mynumber u, v; - int4 k, m, n; + double xx, res, t, cor, y, s, c, sn, ssn, cs, ccs; + mynumber u; + int4 k, m; double retval = 0; +#ifndef IN_SINCOS SET_RESTORE_ROUND_53BIT (FE_TONEAREST); +#endif u.x = x; m = u.i[HIGH_HALF]; k = 0x7fffffff & m; /* no sign */ if (k < 0x3e500000) /* if x->0 =>sin(x)=x */ { - if (fabs (x) < DBL_MIN) - { - double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); retval = x; } /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/ @@ -353,146 +531,22 @@ __sin (double x) retval = (res == res + 1.020 * cor) ? ((m > 0) ? res : -res) : slow2 (x); } /* else if (k < 0x400368fd) */ +#ifndef IN_SINCOS /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/ else if (k < 0x419921FB) { - t = (x * hpinv + toint); - xn = t - toint; - v.x = t; - y = (x - xn * mp1) - xn * mp2; - n = v.i[LOW_HALF] & 3; - da = xn * mp3; - a = y - da; - da = (y - a) - da; - eps = fabs (x) * 1.2e-30; - - switch (n) - { /* quarter of unit circle */ - case 0: - case 2: - xx = a * a; - if (n) - { - a = -a; - da = -da; - } - if (xx < 0.01588) - { - /* Taylor series. */ - res = TAYLOR_SIN (xx, a, da, cor); - cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; - retval = (res == res + cor) ? res : sloww (a, da, x); - } - else - { - if (a > 0) - m = 1; - else - { - m = 0; - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big); - res = do_sin (u, y, da, &cor); - cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; - retval = ((res == res + cor) ? ((m) ? res : -res) - : sloww1 (a, da, x, m)); - } - break; - - case 1: - case 3: - if (a < 0) - { - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big) + da; - res = do_cos (u, y, &cor); - cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; - retval = ((res == res + cor) ? ((n & 2) ? -res : res) - : sloww2 (a, da, x, n)); - break; - } + double a, da; + int4 n = reduce_sincos_1 (x, &a, &da); + retval = do_sincos_1 (a, da, x, n, 0); } /* else if (k < 0x419921FB ) */ /*---------------------105414350 <|x|< 281474976710656 --------------------*/ else if (k < 0x42F00000) { - t = (x * hpinv + toint); - xn = t - toint; - v.x = t; - xn1 = (xn + 8.0e22) - 8.0e22; - xn2 = xn - xn1; - y = ((((x - xn1 * mp1) - xn1 * mp2) - xn2 * mp1) - xn2 * mp2); - n = v.i[LOW_HALF] & 3; - da = xn1 * pp3; - t = y - da; - da = (y - t) - da; - da = (da - xn2 * pp3) - xn * pp4; - a = t + da; - da = (t - a) + da; - eps = 1.0e-24; - - switch (n) - { - case 0: - case 2: - xx = a * a; - if (n) - { - a = -a; - da = -da; - } - if (xx < 0.01588) - { - /* Taylor series. */ - res = TAYLOR_SIN (xx, a, da, cor); - cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; - retval = (res == res + cor) ? res : bsloww (a, da, x, n); - } - else - { - double t; - if (a > 0) - { - m = 1; - t = a; - db = da; - } - else - { - m = 0; - t = -a; - db = -da; - } - u.x = big + t; - y = t - (u.x - big); - res = do_sin (u, y, db, &cor); - cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; - retval = ((res == res + cor) ? ((m) ? res : -res) - : bsloww1 (a, da, x, n)); - } - break; + double a, da; - case 1: - case 3: - if (a < 0) - { - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big) + da; - res = do_cos (u, y, &cor); - cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; - retval = ((res == res + cor) ? ((n & 2) ? -res : res) - : bsloww2 (a, da, x, n)); - break; - } + int4 n = reduce_sincos_2 (x, &a, &da); + retval = do_sincos_2 (a, da, x, n, 0); } /* else if (k < 0x42F00000 ) */ /* -----------------281474976710656 <|x| <2^1024----------------------------*/ @@ -506,6 +560,7 @@ __sin (double x) __set_errno (EDOM); retval = x / x; } +#endif return retval; } @@ -516,18 +571,23 @@ __sin (double x) /* it computes the correctly rounded (to nearest) value of cos(x) */ /*******************************************************************/ +#ifdef IN_SINCOS +static double +#else double SECTION +#endif __cos (double x) { - double y, xx, res, t, cor, xn, a, da, db, eps, xn1, - xn2; - mynumber u, v; - int4 k, m, n; + double y, xx, res, cor, a, da; + mynumber u; + int4 k, m; double retval = 0; +#ifndef IN_SINCOS SET_RESTORE_ROUND_53BIT (FE_TONEAREST); +#endif u.x = x; m = u.i[HIGH_HALF]; @@ -556,7 +616,7 @@ __cos (double x) { res = TAYLOR_SIN (xx, a, da, cor); cor = (cor > 0) ? 1.02 * cor + 1.0e-31 : 1.02 * cor - 1.0e-31; - retval = (res == res + cor) ? res : csloww (a, da, x); + retval = (res == res + cor) ? res : sloww (a, da, x, 1); } else { @@ -575,150 +635,26 @@ __cos (double x) res = do_sin (u, y, da, &cor); cor = (cor > 0) ? 1.035 * cor + 1.0e-31 : 1.035 * cor - 1.0e-31; retval = ((res == res + cor) ? ((m) ? res : -res) - : csloww1 (a, da, x, m)); + : sloww1 (a, da, x, m, 1)); } } /* else if (k < 0x400368fd) */ +#ifndef IN_SINCOS else if (k < 0x419921FB) { /* 2.426265<|x|< 105414350 */ - t = (x * hpinv + toint); - xn = t - toint; - v.x = t; - y = (x - xn * mp1) - xn * mp2; - n = v.i[LOW_HALF] & 3; - da = xn * mp3; - a = y - da; - da = (y - a) - da; - eps = fabs (x) * 1.2e-30; - - switch (n) - { - case 1: - case 3: - xx = a * a; - if (n == 1) - { - a = -a; - da = -da; - } - if (xx < 0.01588) - { - res = TAYLOR_SIN (xx, a, da, cor); - cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; - retval = (res == res + cor) ? res : csloww (a, da, x); - } - else - { - if (a > 0) - { - m = 1; - } - else - { - m = 0; - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big); - res = do_sin (u, y, da, &cor); - cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; - retval = ((res == res + cor) ? ((m) ? res : -res) - : csloww1 (a, da, x, m)); - } - break; - - case 0: - case 2: - if (a < 0) - { - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big) + da; - res = do_cos (u, y, &cor); - cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; - retval = ((res == res + cor) ? ((n) ? -res : res) - : csloww2 (a, da, x, n)); - break; - } + double a, da; + int4 n = reduce_sincos_1 (x, &a, &da); + retval = do_sincos_1 (a, da, x, n, 1); } /* else if (k < 0x419921FB ) */ else if (k < 0x42F00000) { - t = (x * hpinv + toint); - xn = t - toint; - v.x = t; - xn1 = (xn + 8.0e22) - 8.0e22; - xn2 = xn - xn1; - y = ((((x - xn1 * mp1) - xn1 * mp2) - xn2 * mp1) - xn2 * mp2); - n = v.i[LOW_HALF] & 3; - da = xn1 * pp3; - t = y - da; - da = (y - t) - da; - da = (da - xn2 * pp3) - xn * pp4; - a = t + da; - da = (t - a) + da; - eps = 1.0e-24; - - switch (n) - { - case 1: - case 3: - xx = a * a; - if (n == 1) - { - a = -a; - da = -da; - } - if (xx < 0.01588) - { - res = TAYLOR_SIN (xx, a, da, cor); - cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; - retval = (res == res + cor) ? res : bsloww (a, da, x, n); - } - else - { - double t; - if (a > 0) - { - m = 1; - t = a; - db = da; - } - else - { - m = 0; - t = -a; - db = -da; - } - u.x = big + t; - y = t - (u.x - big); - res = do_sin (u, y, db, &cor); - cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; - retval = ((res == res + cor) ? ((m) ? res : -res) - : bsloww1 (a, da, x, n)); - } - break; + double a, da; - case 0: - case 2: - if (a < 0) - { - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big) + da; - res = do_cos (u, y, &cor); - cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; - retval = ((res == res + cor) ? ((n) ? -res : res) - : bsloww2 (a, da, x, n)); - break; - } + int4 n = reduce_sincos_2 (x, &a, &da); + retval = do_sincos_2 (a, da, x, n, 1); } /* else if (k < 0x42F00000 ) */ /* 281474976710656 <|x| <2^1024 */ @@ -731,6 +667,7 @@ __cos (double x) __set_errno (EDOM); retval = x / x; /* |x| > 2^1024 */ } +#endif return retval; } @@ -748,14 +685,12 @@ slow (double x) res = TAYLOR_SLOW (x, 0, cor); if (res == res + 1.0007 * cor) return res; - else - { - __dubsin (fabs (x), 0, w); - if (w[0] == w[0] + 1.000000001 * w[1]) - return (x > 0) ? w[0] : -w[0]; - else - return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); - } + + __dubsin (fabs (x), 0, w); + if (w[0] == w[0] + 1.000000001 * w[1]) + return (x > 0) ? w[0] : -w[0]; + + return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); } /*******************************************************************************/ @@ -775,14 +710,12 @@ slow1 (double x) res = do_sin_slow (u, y, 0, 0, &cor); if (res == res + cor) return (x > 0) ? res : -res; - else - { - __dubsin (fabs (x), 0, w); - if (w[0] == w[0] + 1.000000005 * w[1]) - return (x > 0) ? w[0] : -w[0]; - else - return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); - } + + __dubsin (fabs (x), 0, w); + if (w[0] == w[0] + 1.000000005 * w[1]) + return (x > 0) ? w[0] : -w[0]; + + return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); } /**************************************************************************/ @@ -813,17 +746,15 @@ slow2 (double x) res = do_cos_slow (u, y, del, 0, &cor); if (res == res + cor) return (x > 0) ? res : -res; - else - { - y = fabs (x) - hp0; - y1 = y - hp1; - y2 = (y - y1) - hp1; - __docos (y1, y2, w); - if (w[0] == w[0] + 1.000000005 * w[1]) - return (x > 0) ? w[0] : -w[0]; - else - return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); - } + + y = fabs (x) - hp0; + y1 = y - hp1; + y2 = (y - y1) - hp1; + __docos (y1, y2, w); + if (w[0] == w[0] + 1.000000005 * w[1]) + return (x > 0) ? w[0] : -w[0]; + + return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); } /***************************************************************************/ @@ -836,12 +767,13 @@ slow2 (double x) static double SECTION -sloww (double x, double dx, double orig) +sloww (double x, double dx, double orig, int k) { double y, t, res, cor, w[2], a, da, xn; mynumber v; int4 n; res = TAYLOR_SLOW (x, dx, cor); + if (cor > 0) cor = 1.0005 * cor + fabs (orig) * 3.1e-30; else @@ -849,46 +781,43 @@ sloww (double x, double dx, double orig) if (res == res + cor) return res; + + (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); + if (w[1] > 0) + cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30; else + cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30; + + if (w[0] == w[0] + cor) + return (x > 0) ? w[0] : -w[0]; + + t = (orig * hpinv + toint); + xn = t - toint; + v.x = t; + y = (orig - xn * mp1) - xn * mp2; + n = (v.i[LOW_HALF] + k) & 3; + da = xn * pp3; + t = y - da; + da = (y - t) - da; + y = xn * pp4; + a = t - y; + da = ((t - a) - y) + da; + + if (n == 2 || n == 1) { - (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); - if (w[1] > 0) - cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30; - else - cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30; + a = -a; + da = -da; + } + (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w); + if (w[1] > 0) + cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40; + else + cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40; - if (w[0] == w[0] + cor) - return (x > 0) ? w[0] : -w[0]; - else - { - t = (orig * hpinv + toint); - xn = t - toint; - v.x = t; - y = (orig - xn * mp1) - xn * mp2; - n = v.i[LOW_HALF] & 3; - da = xn * pp3; - t = y - da; - da = (y - t) - da; - y = xn * pp4; - a = t - y; - da = ((t - a) - y) + da; - if (n & 2) - { - a = -a; - da = -da; - } - (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w); - if (w[1] > 0) - cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40; - else - cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40; + if (w[0] == w[0] + cor) + return (a > 0) ? w[0] : -w[0]; - if (w[0] == w[0] + cor) - return (a > 0) ? w[0] : -w[0]; - else - return __mpsin (orig, 0, true); - } - } + return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); } /***************************************************************************/ @@ -900,7 +829,7 @@ sloww (double x, double dx, double orig) static double SECTION -sloww1 (double x, double dx, double orig, int m) +sloww1 (double x, double dx, double orig, int m, int k) { mynumber u; double w[2], y, cor, res; @@ -911,20 +840,18 @@ sloww1 (double x, double dx, double orig, int m) if (res == res + cor) return (m > 0) ? res : -res; + + __dubsin (x, dx, w); + + if (w[1] > 0) + cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); else - { - __dubsin (x, dx, w); + cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); - if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); - else - cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); + if (w[0] == w[0] + cor) + return (m > 0) ? w[0] : -w[0]; - if (w[0] == w[0] + cor) - return (m > 0) ? w[0] : -w[0]; - else - return __mpsin (orig, 0, true); - } + return (k == 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); } /***************************************************************************/ @@ -947,20 +874,18 @@ sloww2 (double x, double dx, double orig, int n) if (res == res + cor) return (n & 2) ? -res : res; + + __docos (x, dx, w); + + if (w[1] > 0) + cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); else - { - __docos (x, dx, w); + cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); - if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); - else - cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); + if (w[0] == w[0] + cor) + return (n & 2) ? -w[0] : w[0]; - if (w[0] == w[0] + cor) - return (n & 2) ? -w[0] : w[0]; - else - return __mpsin (orig, 0, true); - } + return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true); } /***************************************************************************/ @@ -981,18 +906,17 @@ bsloww (double x, double dx, double orig, int n) cor = (cor > 0) ? 1.0005 * cor + 1.1e-24 : 1.0005 * cor - 1.1e-24; if (res == res + cor) return res; + + (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); + if (w[1] > 0) + cor = 1.000000001 * w[1] + 1.1e-24; else - { - (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); - if (w[1] > 0) - cor = 1.000000001 * w[1] + 1.1e-24; - else - cor = 1.000000001 * w[1] - 1.1e-24; - if (w[0] == w[0] + cor) - return (x > 0) ? w[0] : -w[0]; - else - return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); - } + cor = 1.000000001 * w[1] - 1.1e-24; + + if (w[0] == w[0] + cor) + return (x > 0) ? w[0] : -w[0]; + + return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); } /***************************************************************************/ @@ -1016,20 +940,18 @@ bsloww1 (double x, double dx, double orig, int n) res = do_sin_slow (u, y, dx, 1.1e-24, &cor); if (res == res + cor) return (x > 0) ? res : -res; + + __dubsin (fabs (x), dx, w); + + if (w[1] > 0) + cor = 1.000000005 * w[1] + 1.1e-24; else - { - __dubsin (fabs (x), dx, w); + cor = 1.000000005 * w[1] - 1.1e-24; - if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-24; - else - cor = 1.000000005 * w[1] - 1.1e-24; + if (w[0] == w[0] + cor) + return (x > 0) ? w[0] : -w[0]; - if (w[0] == w[0] + cor) - return (x > 0) ? w[0] : -w[0]; - else - return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); - } + return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); } /***************************************************************************/ @@ -1053,20 +975,18 @@ bsloww2 (double x, double dx, double orig, int n) res = do_cos_slow (u, y, dx, 1.1e-24, &cor); if (res == res + cor) return (n & 2) ? -res : res; + + __docos (fabs (x), dx, w); + + if (w[1] > 0) + cor = 1.000000005 * w[1] + 1.1e-24; else - { - __docos (fabs (x), dx, w); + cor = 1.000000005 * w[1] - 1.1e-24; - if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-24; - else - cor = 1.000000005 * w[1] - 1.1e-24; + if (w[0] == w[0] + cor) + return (n & 2) ? -w[0] : w[0]; - if (w[0] == w[0] + cor) - return (n & 2) ? -w[0] : w[0]; - else - return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true); - } + return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true); } /************************************************************************/ @@ -1087,154 +1007,13 @@ cslow2 (double x) res = do_cos_slow (u, y, 0, 0, &cor); if (res == res + cor) return res; - else - { - y = fabs (x); - __docos (y, 0, w); - if (w[0] == w[0] + 1.000000005 * w[1]) - return w[0]; - else - return __mpcos (x, 0, false); - } -} - -/***************************************************************************/ -/* Routine compute cos(x+dx) (Double-Length number) where x is small enough*/ -/* to use Taylor series around zero and (x+dx) .Routine receive also */ -/* (right argument) the original value of x for computing error of */ -/* result.And if result not accurate enough routine calls other routines */ -/***************************************************************************/ - - -static double -SECTION -csloww (double x, double dx, double orig) -{ - double y, t, res, cor, w[2], a, da, xn; - mynumber v; - int4 n; - - /* Taylor series */ - res = TAYLOR_SLOW (x, dx, cor); - - if (cor > 0) - cor = 1.0005 * cor + fabs (orig) * 3.1e-30; - else - cor = 1.0005 * cor - fabs (orig) * 3.1e-30; - - if (res == res + cor) - return res; - else - { - (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); - - if (w[1] > 0) - cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30; - else - cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30; - - if (w[0] == w[0] + cor) - return (x > 0) ? w[0] : -w[0]; - else - { - t = (orig * hpinv + toint); - xn = t - toint; - v.x = t; - y = (orig - xn * mp1) - xn * mp2; - n = v.i[LOW_HALF] & 3; - da = xn * pp3; - t = y - da; - da = (y - t) - da; - y = xn * pp4; - a = t - y; - da = ((t - a) - y) + da; - if (n == 1) - { - a = -a; - da = -da; - } - (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w); - - if (w[1] > 0) - cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40; - else - cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40; - - if (w[0] == w[0] + cor) - return (a > 0) ? w[0] : -w[0]; - else - return __mpcos (orig, 0, true); - } - } -} - -/***************************************************************************/ -/* Routine compute sin(x+dx) (Double-Length number) where x in first or */ -/* third quarter of unit circle.Routine receive also (right argument) the */ -/* original value of x for computing error of result.And if result not */ -/* accurate enough routine calls other routines */ -/***************************************************************************/ -static double -SECTION -csloww1 (double x, double dx, double orig, int m) -{ - mynumber u; - double w[2], y, cor, res; - - u.x = big + x; - y = x - (u.x - big); - res = do_sin_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); - - if (res == res + cor) - return (m > 0) ? res : -res; - else - { - __dubsin (x, dx, w); - if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); - else - cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); - if (w[0] == w[0] + cor) - return (m > 0) ? w[0] : -w[0]; - else - return __mpcos (orig, 0, true); - } -} - - -/***************************************************************************/ -/* Routine compute sin(x+dx) (Double-Length number) where x in second or */ -/* fourth quarter of unit circle.Routine receive also the original value */ -/* and quarter(n= 1or 3)of x for computing error of result.And if result not*/ -/* accurate enough routine calls other routines */ -/***************************************************************************/ - -static double -SECTION -csloww2 (double x, double dx, double orig, int n) -{ - mynumber u; - double w[2], y, cor, res; - - u.x = big + x; - y = x - (u.x - big); - res = do_cos_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); + y = fabs (x); + __docos (y, 0, w); + if (w[0] == w[0] + 1.000000005 * w[1]) + return w[0]; - if (res == res + cor) - return (n) ? -res : res; - else - { - __docos (x, dx, w); - if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); - else - cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); - if (w[0] == w[0] + cor) - return (n) ? -w[0] : w[0]; - else - return __mpcos (orig, 0, true); - } + return __mpcos (x, 0, false); } #ifndef __cos diff --git a/sysdeps/ieee754/dbl-64/s_sincos.c b/sysdeps/ieee754/dbl-64/s_sincos.c index c8a99991cc..c389226b04 100644 --- a/sysdeps/ieee754/dbl-64/s_sincos.c +++ b/sysdeps/ieee754/dbl-64/s_sincos.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. @@ -22,29 +22,89 @@ #include <math_private.h> +#define __sin __sin_local +#define __cos __cos_local +#define IN_SINCOS 1 +#include "s_sin.c" + +/* Consolidated version of reduce_and_compute in s_sin.c that does range + reduction only once and computes sin and cos together. */ +static inline void +__always_inline +reduce_and_compute_sincos (double x, double *sinx, double *cosx) +{ + double a, da; + unsigned int n = __branred (x, &a, &da); + + n = n & 3; + + if (n == 1 || n == 2) + { + a = -a; + da = -da; + } + + if (n & 1) + { + double *temp = cosx; + cosx = sinx; + sinx = temp; + } + + if (a * a < 0.01588) + *sinx = bsloww (a, da, x, n); + else + *sinx = bsloww1 (a, da, x, n); + *cosx = bsloww2 (a, da, x, n); +} void __sincos (double x, double *sinx, double *cosx) { - int32_t ix; + mynumber u; + int k; - /* High word of x. */ - GET_HIGH_WORD (ix, x); + SET_RESTORE_ROUND_53BIT (FE_TONEAREST); - /* |x| ~< pi/4 */ - ix &= 0x7fffffff; - if (ix >= 0x7ff00000) + u.x = x; + k = 0x7fffffff & u.i[HIGH_HALF]; + + if (k < 0x400368fd) { - /* sin(Inf or NaN) is NaN */ - *sinx = *cosx = x - x; - if (__isinf_ns (x)) - __set_errno (EDOM); + *sinx = __sin_local (x); + *cosx = __cos_local (x); + return; } - else + if (k < 0x419921FB) + { + double a, da; + int4 n = reduce_sincos_1 (x, &a, &da); + + *sinx = do_sincos_1 (a, da, x, n, 0); + *cosx = do_sincos_1 (a, da, x, n, 1); + + return; + } + if (k < 0x42F00000) { - *sinx = __sin (x); - *cosx = __cos (x); + double a, da; + int4 n = reduce_sincos_2 (x, &a, &da); + + *sinx = do_sincos_2 (a, da, x, n, 0); + *cosx = do_sincos_2 (a, da, x, n, 1); + + return; + } + if (k < 0x7ff00000) + { + reduce_and_compute_sincos (x, sinx, cosx); + return; } + + if (isinf (x)) + __set_errno (EDOM); + + *sinx = *cosx = x / x; } weak_alias (__sincos, sincos) #ifdef NO_LONG_DOUBLE diff --git a/sysdeps/ieee754/dbl-64/s_tan.c b/sysdeps/ieee754/dbl-64/s_tan.c index dcb4aca2de..ccffdd62d7 100644 --- a/sysdeps/ieee754/dbl-64/s_tan.c +++ b/sysdeps/ieee754/dbl-64/s_tan.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 @@ -34,6 +34,7 @@ /*********************************************************************/ #include <errno.h> +#include <float.h> #include "endian.h" #include <dla.h> #include "mpa.h" @@ -91,6 +92,7 @@ tan (double x) /* (I) The case abs(x) <= 1.259e-8 */ if (w <= g1.d) { + math_check_force_underflow_nonneg (w); retval = x; goto ret; } diff --git a/sysdeps/ieee754/dbl-64/s_tanh.c b/sysdeps/ieee754/dbl-64/s_tanh.c index 23cfcdead5..344a2f0330 100644 --- a/sysdeps/ieee754/dbl-64/s_tanh.c +++ b/sysdeps/ieee754/dbl-64/s_tanh.c @@ -38,6 +38,7 @@ static char rcsid[] = "$NetBSD: s_tanh.c,v 1.7 1995/05/10 20:48:22 jtc Exp $"; * only tanh(0)=0 is exact for finite argument. */ +#include <float.h> #include <math.h> #include <math_private.h> @@ -68,7 +69,10 @@ __tanh (double x) if ((ix | lx) == 0) return x; /* x == +-0 */ if (ix < 0x3c800000) /* |x|<2**-55 */ - return x * (one + x); /* tanh(small) = small */ + { + math_check_force_underflow (x); + return x * (one + x); /* tanh(small) = small */ + } if (ix >= 0x3ff00000) /* |x|>=1 */ { t = __expm1 (two * fabs (x)); diff --git a/sysdeps/ieee754/dbl-64/s_trunc.c b/sysdeps/ieee754/dbl-64/s_trunc.c index ac6c7597b5..f64e097f16 100644 --- a/sysdeps/ieee754/dbl-64/s_trunc.c +++ b/sysdeps/ieee754/dbl-64/s_trunc.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. diff --git a/sysdeps/ieee754/dbl-64/sincos32.c b/sysdeps/ieee754/dbl-64/sincos32.c index d942a1f29f..ad43465f64 100644 --- a/sysdeps/ieee754/dbl-64/sincos32.c +++ b/sysdeps/ieee754/dbl-64/sincos32.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/dbl-64/sincos32.h b/sysdeps/ieee754/dbl-64/sincos32.h index ee70b08743..5a5f2e223a 100644 --- a/sysdeps/ieee754/dbl-64/sincos32.h +++ b/sysdeps/ieee754/dbl-64/sincos32.h @@ -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/dbl-64/sincostab.c b/sysdeps/ieee754/dbl-64/sincostab.c index 44cb525f42..cd88af7331 100644 --- a/sysdeps/ieee754/dbl-64/sincostab.c +++ b/sysdeps/ieee754/dbl-64/sincostab.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/dbl-64/slowexp.c b/sysdeps/ieee754/dbl-64/slowexp.c index 7950111412..b602262783 100644 --- a/sysdeps/ieee754/dbl-64/slowexp.c +++ b/sysdeps/ieee754/dbl-64/slowexp.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/dbl-64/slowpow.c b/sysdeps/ieee754/dbl-64/slowpow.c index 117225ea46..266910c979 100644 --- a/sysdeps/ieee754/dbl-64/slowpow.c +++ b/sysdeps/ieee754/dbl-64/slowpow.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/dbl-64/t_exp.c b/sysdeps/ieee754/dbl-64/t_exp.c index 8d5c7b81a3..3c921f1d3d 100644 --- a/sysdeps/ieee754/dbl-64/t_exp.c +++ b/sysdeps/ieee754/dbl-64/t_exp.c @@ -1,5 +1,5 @@ /* Accurate tables for exp(). - Copyright (C) 1998-2015 Free Software Foundation, Inc. + Copyright (C) 1998-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Geoffrey Keating <geoffk@ozemail.com.au> diff --git a/sysdeps/ieee754/dbl-64/uasncs.h b/sysdeps/ieee754/dbl-64/uasncs.h index c60e93c0a2..e88bd1bf3e 100644 --- a/sysdeps/ieee754/dbl-64/uasncs.h +++ b/sysdeps/ieee754/dbl-64/uasncs.h @@ -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/dbl-64/uatan.tbl b/sysdeps/ieee754/dbl-64/uatan.tbl index f1337fce21..979d65d8b2 100644 --- a/sysdeps/ieee754/dbl-64/uatan.tbl +++ b/sysdeps/ieee754/dbl-64/uatan.tbl @@ -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/dbl-64/uexp.h b/sysdeps/ieee754/dbl-64/uexp.h index 6817eaf074..eb7aa5f5fb 100644 --- a/sysdeps/ieee754/dbl-64/uexp.h +++ b/sysdeps/ieee754/dbl-64/uexp.h @@ -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 @@ -29,7 +29,7 @@ #include "mydefs.h" -const static double one = 1.0, zero = 0.0, hhuge = 1.0e300, tiny = 1.0e-300, +const static double zero = 0.0, hhuge = 1.0e300, tiny = 1.0e-300, err_0 = 1.000014, err_1 = 0.000016; const static int4 bigint = 0x40862002, badint = 0x40876000,smallint = 0x3C8fffff; diff --git a/sysdeps/ieee754/dbl-64/uexp.tbl b/sysdeps/ieee754/dbl-64/uexp.tbl index 859e908c5d..9abfc7ef71 100644 --- a/sysdeps/ieee754/dbl-64/uexp.tbl +++ b/sysdeps/ieee754/dbl-64/uexp.tbl @@ -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/dbl-64/ulog.h b/sysdeps/ieee754/dbl-64/ulog.h index ce4db49024..d507c693d8 100644 --- a/sysdeps/ieee754/dbl-64/ulog.h +++ b/sysdeps/ieee754/dbl-64/ulog.h @@ -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/dbl-64/ulog.tbl b/sysdeps/ieee754/dbl-64/ulog.tbl index e5d2fdd34b..dd92cbb886 100644 --- a/sysdeps/ieee754/dbl-64/ulog.tbl +++ b/sysdeps/ieee754/dbl-64/ulog.tbl @@ -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/dbl-64/upow.h b/sysdeps/ieee754/dbl-64/upow.h index c8569a9456..dd9ba32b71 100644 --- a/sysdeps/ieee754/dbl-64/upow.h +++ b/sysdeps/ieee754/dbl-64/upow.h @@ -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 @@ -34,7 +34,6 @@ /**/ nZERO = {{0x80000000, 0}}, /* -0.0 */ /**/ INF = {{0x7ff00000, 0x00000000}}, /* INF */ /**/ nINF = {{0xfff00000, 0x00000000}}, /* -INF */ -/**/ sqrt_2 = {{0x3ff6a09e, 0x667f3bcc}}, /* sqrt(2) */ /**/ ln2a = {{0x3fe62e42, 0xfefa3800}}, /* ln(2) 43 bits */ /**/ ln2b = {{0x3d2ef357, 0x93c76730}}, /* ln(2)-ln2a */ /**/ bigu = {{0x4297ffff, 0xfffffd2c}}, /* 1.5*2**42 -724*2**-10 */ @@ -48,7 +47,6 @@ /**/ nZERO = {{0, 0x80000000}}, /* -0.0 */ /**/ INF = {{0x00000000, 0x7ff00000}}, /* INF */ /**/ nINF = {{0x00000000, 0xfff00000}}, /* -INF */ -/**/ sqrt_2 = {{0x667f3bcc, 0x3ff6a09e}}, /* sqrt(2) */ /**/ ln2a = {{0xfefa3800, 0x3fe62e42}}, /* ln(2) 43 bits */ /**/ ln2b = {{0x93c76730, 0x3d2ef357}}, /* ln(2)-ln2a */ /**/ bigu = {{0xfffffd2c, 0x4297ffff}}, /* 1.5*2**42 -724*2**-10 */ diff --git a/sysdeps/ieee754/dbl-64/upow.tbl b/sysdeps/ieee754/dbl-64/upow.tbl index 1d506b5ad7..55d878edad 100644 --- a/sysdeps/ieee754/dbl-64/upow.tbl +++ b/sysdeps/ieee754/dbl-64/upow.tbl @@ -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/dbl-64/urem.h b/sysdeps/ieee754/dbl-64/urem.h index f612c135f8..42b26273fd 100644 --- a/sysdeps/ieee754/dbl-64/urem.h +++ b/sysdeps/ieee754/dbl-64/urem.h @@ -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/dbl-64/uroot.h b/sysdeps/ieee754/dbl-64/uroot.h deleted file mode 100644 index 4bbcc3bac7..0000000000 --- a/sysdeps/ieee754/dbl-64/uroot.h +++ /dev/null @@ -1,43 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * Written by International Business Machines Corp. - * Copyright (C) 2001-2015 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 - * the Free Software Foundation; either version 2.1 of the License, or - * (at your option) any later version. - * - * This program 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 this program; if not, see <http://www.gnu.org/licenses/>. - */ - -/******************************************************************/ -/* */ -/* MODULE_NAME:uroot.h */ -/* */ -/* common data and variables prototype and definition */ -/******************************************************************/ - -#ifndef UROOT_H -#define UROOT_H - -#ifdef BIG_ENDI - static const mynumber -/**/ t512 = {{0x5ff00000, 0x00000000 }}, /* 2^512 */ -/**/ tm256 = {{0x2ff00000, 0x00000000 }}; /* 2^-256 */ - -#else -#ifdef LITTLE_ENDI - static const mynumber -/**/ t512 = {{0x00000000, 0x5ff00000 }}, /* 2^512 */ -/**/ tm256 = {{0x00000000, 0x2ff00000 }}; /* 2^-256 */ -#endif -#endif - -#endif diff --git a/sysdeps/ieee754/dbl-64/usncs.h b/sysdeps/ieee754/dbl-64/usncs.h index e3a8193c40..1408695c85 100644 --- a/sysdeps/ieee754/dbl-64/usncs.h +++ b/sysdeps/ieee754/dbl-64/usncs.h @@ -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/dbl-64/utan.h b/sysdeps/ieee754/dbl-64/utan.h index efa13ca5cf..0a99812754 100644 --- a/sysdeps/ieee754/dbl-64/utan.h +++ b/sysdeps/ieee754/dbl-64/utan.h @@ -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/dbl-64/utan.tbl b/sysdeps/ieee754/dbl-64/utan.tbl index c088588229..8a8c2afa53 100644 --- a/sysdeps/ieee754/dbl-64/utan.tbl +++ b/sysdeps/ieee754/dbl-64/utan.tbl @@ -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/dbl-64/w_exp.c b/sysdeps/ieee754/dbl-64/w_exp.c index 61b2dbbd0f..595ed992d2 100644 --- a/sysdeps/ieee754/dbl-64/w_exp.c +++ b/sysdeps/ieee754/dbl-64/w_exp.c @@ -1,4 +1,4 @@ -/* 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 Ulrich Drepper <drepper@gmail.com>, 2011. diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/math_private.h b/sysdeps/ieee754/dbl-64/wordsize-64/math_private.h deleted file mode 100644 index 4f9219934a..0000000000 --- a/sysdeps/ieee754/dbl-64/wordsize-64/math_private.h +++ /dev/null @@ -1,36 +0,0 @@ -#ifndef _MATH_PRIVATE_H_ - -#include_next <math_private.h> -#include <stdint.h> - -#ifndef __isnan -extern __always_inline int -__isnan (double d) -{ - uint64_t di; - EXTRACT_WORDS64 (di, d); - return (di & 0x7fffffffffffffffull) > 0x7ff0000000000000ull; -} -#endif - -#ifndef __isinf_ns -extern __always_inline int -__isinf_ns (double d) -{ - uint64_t di; - EXTRACT_WORDS64 (di, d); - return (di & 0x7fffffffffffffffull) == 0x7ff0000000000000ull; -} -#endif - -#ifndef __finite -extern __always_inline int -__finite (double d) -{ - uint64_t di; - EXTRACT_WORDS64 (di, d); - return (di & 0x7fffffffffffffffull) < 0x7ff0000000000000ull; -} -#endif - -#endif /* _MATH_PRIVATE_H_ */ diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_finite.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_finite.c index fcf2e6d5b6..ef51608f6e 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_finite.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_finite.c @@ -16,6 +16,7 @@ #include <math.h> #include <math_private.h> +#include <shlib-compat.h> #include <stdint.h> #undef __finite @@ -24,11 +25,18 @@ __finite(double x) { int64_t lx; EXTRACT_WORDS64(lx,x); - return (int)((uint64_t)((lx&INT64_C(0x7fffffffffffffff))-INT64_C(0x7ff0000000000000))>>63); + return (int)((uint64_t)((lx&INT64_C(0x7ff0000000000000))-INT64_C(0x7ff0000000000000))>>63); } hidden_def (__finite) weak_alias (__finite, finite) #ifdef NO_LONG_DOUBLE -strong_alias (__finite, __finitel) +# ifdef LDBL_CLASSIFY_COMPAT +# if SHLIB_COMPAT (libc, GLIBC_2_0, GLIBC_2_23) +compat_symbol (libc, __finite, __finitel, GLIBC_2_0); +# endif +# if SHLIB_COMPAT (libm, GLIBC_2_1, GLIBC_2_23) +compat_symbol (libm, __finite, __finitel, GLIBC_2_1); +# endif +# endif weak_alias (__finite, finitel) #endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c index 0dbadb15eb..b7ed14bfa2 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c @@ -1,5 +1,5 @@ /* Round double to integer away from zero. - 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 Ulrich Drepper <drepper@cygnus.com>, 2011. diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c index f1fa3154ee..42068f8187 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c @@ -1,4 +1,4 @@ -/* 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 Ulrich Drepper <drepper@gmail.com>, 2011. diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_isinf.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_isinf.c index 163fc31efc..951fb73239 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_isinf.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_isinf.c @@ -11,6 +11,7 @@ #include <math.h> #include <math_private.h> +#include <shlib-compat.h> int __isinf (double x) @@ -25,6 +26,8 @@ __isinf (double x) hidden_def (__isinf) weak_alias (__isinf, isinf) #ifdef NO_LONG_DOUBLE -strong_alias (__isinf, __isinfl) +# if defined LDBL_CLASSIFY_COMPAT && SHLIB_COMPAT (libc, GLIBC_2_0, GLIBC_2_23) +compat_symbol (libc, __isinf, __isinfl, GLIBC_2_0); +# endif weak_alias (__isinf, isinfl) #endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_isinf_ns.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_isinf_ns.c deleted file mode 100644 index 9d78ed13ae..0000000000 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_isinf_ns.c +++ /dev/null @@ -1,20 +0,0 @@ -/* - * Written by Ulrich Drepper <drepper@gmail.com>. - */ - -/* - * __isinf_ns(x) returns != 0 if x is ±inf, else 0; - * no branching! - */ - -#include <math.h> -#include <math_private.h> - -#undef __isinf_ns -int -__isinf_ns (double x) -{ - int64_t ix; - EXTRACT_WORDS64(ix,x); - return (ix & UINT64_C(0x7fffffffffffffff)) == UINT64_C(0x7ff0000000000000); -} diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c index e80b84ca02..bcff9e3b67 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c @@ -17,6 +17,7 @@ #include <math.h> #include <math_private.h> +#include <shlib-compat.h> #include <stdint.h> #undef __isnan @@ -31,6 +32,8 @@ int __isnan(double x) hidden_def (__isnan) weak_alias (__isnan, isnan) #ifdef NO_LONG_DOUBLE -strong_alias (__isnan, __isnanl) +# if defined LDBL_CLASSIFY_COMPAT && SHLIB_COMPAT (libc, GLIBC_2_0, GLIBC_2_23) +compat_symbol (libc, __isnan, __isnanl, GLIBC_2_0); +# endif weak_alias (__isnan, isnanl) #endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c index 67d77d5732..c22e608c6e 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.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 @@ -25,7 +25,6 @@ __issignaling (double x) u_int64_t xi; EXTRACT_WORDS64 (xi, x); #ifdef HIGH_ORDER_BIT_IS_SET_FOR_SNAN -# error untested /* We only have to care about the high-order bit of x's significand, because having it set (sNaN) already makes the significand different from that used to designate infinity. */ diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c index 8bb52524b7..73ddd8dd9f 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c @@ -1,5 +1,5 @@ /* Round 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. @@ -69,23 +69,10 @@ strong_alias (__llround, __llroundl) weak_alias (__llround, llroundl) #endif -/* long has the same width as long long on LP64 machines, so use an alias. - If building for ILP32 on a machine with 64-bit registers, however, - use a cast if necessary. */ +/* long has the same width as long long on LP64 machines, so use an alias. */ #undef lround #undef __lround -#if !defined (_LP64) && REGISTER_CAST_INT32_TO_INT64 -long int -__lround (double x) -{ - return __llround (x); -} -weak_alias (__lround, lround) -# ifdef NO_LONG_DOUBLE -strong_alias (__lround, __lroundl) -weak_alias (__lround, lroundl) -# endif -#else +#ifdef _LP64 strong_alias (__llround, __lround) weak_alias (__llround, lround) # ifdef NO_LONG_DOUBLE diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c index 1d20e93e6d..c7846054a6 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c @@ -1,5 +1,5 @@ /* Compute radix independent exponent. - 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 Ulrich Drepper <drepper@gmail.com>, 2011. diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c index 1e06fcb16e..9e665b01ec 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c @@ -1 +1,89 @@ -/* The code is the same as llround. Use an alias, see s_llround.c. */ +/* Round double value to long int. + 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/>. */ + +#include <fenv.h> +#include <limits.h> +#include <math.h> + +#include <math_private.h> + +/* For LP64, lround is an alias for llround. */ +#ifndef _LP64 + +long int +__lround (double x) +{ + int32_t j0; + int64_t i0; + long int result; + int sign; + + EXTRACT_WORDS64 (i0, x); + j0 = ((i0 >> 52) & 0x7ff) - 0x3ff; + sign = i0 < 0 ? -1 : 1; + i0 &= UINT64_C(0xfffffffffffff); + i0 |= UINT64_C(0x10000000000000); + + if (j0 < (int32_t) (8 * sizeof (long int)) - 1) + { + if (j0 < 0) + return j0 < -1 ? 0 : sign; + else if (j0 >= 52) + result = i0 << (j0 - 52); + else + { + i0 += UINT64_C(0x8000000000000) >> j0; + + result = i0 >> (52 - 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 + { + /* 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 (sizeof (long int) == 4 + && x <= (double) LONG_MIN - 0.5) + { + /* If truncation produces LONG_MIN, the cast will not raise + the exception, but may raise "inexact". */ + feraiseexcept (FE_INVALID); + return LONG_MIN; + } +#endif + return (long int) x; + } + + return sign * result; +} + +weak_alias (__lround, lround) +# ifdef NO_LONG_DOUBLE +strong_alias (__lround, __lroundl) +weak_alias (__lround, lroundl) +# endif + +#endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c index 304e9d0fb4..716e07e54d 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.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. diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c index 565a4622d4..f1d75bbe68 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c @@ -1,5 +1,5 @@ /* Round 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. diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_scalbn.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_scalbn.c index 3822e1737d..d517a919c8 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_scalbn.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_scalbn.c @@ -55,8 +55,6 @@ __scalbn (double x, int n) INSERT_WORDS64(x,(ix&INT64_C(0x800fffffffffffff))|(k<<52)); return x*twom54; } -weak_alias (__scalbn, scalbn) #ifdef NO_LONG_DOUBLE strong_alias (__scalbn, __scalbnl) -weak_alias (__scalbn, scalbnl) #endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c index 3810f8134d..81ac55e2f6 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.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. diff --git a/sysdeps/ieee754/dbl-64/x2y2m1.c b/sysdeps/ieee754/dbl-64/x2y2m1.c index c96dae59e1..96078888a7 100644 --- a/sysdeps/ieee754/dbl-64/x2y2m1.c +++ b/sysdeps/ieee754/dbl-64/x2y2m1.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. */ double __x2y2m1 (double x, double y) { - double vals[4]; + double vals[5]; SET_RESTORE_ROUND (FE_TONEAREST); mul_split (&vals[1], &vals[0], x, x); mul_split (&vals[3], &vals[2], y, y); - if (x >= 0.75) - vals[1] -= 1.0; - else - { - vals[1] -= 0.5; - vals[3] -= 0.5; - } - qsort (vals, 4, sizeof (double), compare); + vals[4] = -1.0; + qsort (vals, 5, 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 <= 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 (double), compare); + qsort (vals + i + 1, 4 - i, sizeof (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]; } diff --git a/sysdeps/ieee754/dbl-64/x2y2m1f.c b/sysdeps/ieee754/dbl-64/x2y2m1f.c index 43a8acf8a8..480b4c0c67 100644 --- a/sysdeps/ieee754/dbl-64/x2y2m1f.c +++ b/sysdeps/ieee754/dbl-64/x2y2m1f.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 @@ -21,8 +21,8 @@ #include <float.h> /* 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. */ float __x2y2m1f (float x, float y) |