diff options
| author | Jakub Jelinek <jakub@redhat.com> | 2017-07-19 15:12:58 +0200 |
|---|---|---|
| committer | Jakub Jelinek <jakub@gcc.gnu.org> | 2017-07-19 15:12:58 +0200 |
| commit | 1eba0867060b1643d71d4185fbc23995248092bf (patch) | |
| tree | b933cf0474996b72f099f09c134fd321c9590600 /libquadmath | |
| parent | 564e405c138dc1891ca92462a585bf89444d7dd1 (diff) | |
re PR libquadmath/65757 (gfortran gives incorrect result for anint with real*16 argument)
PR libquadmath/65757
* quadmath-imp.h (math_opt_barrier, math_force_eval,
math_narrow_eval, math_check_force_underflow,
math_check_force_underflow_nonneg): Define.
* math/ceilq.c: Backport changes from upstream glibc
between 2012-11-01 and 2017-07-13.
* math/remquoq.c: Likewise.
* math/expq.c: Likewise.
* math/llroundq.c: Likewise.
* math/logq.c: Likewise.
* math/atanq.c: Likewise.
* math/nearbyintq.c: Likewise.
* math/scalblnq.c: Likewise.
* math/finiteq.c: Likewise.
* math/atanhq.c: Likewise.
* math/expm1q.c: Likewise.
* math/sinhq.c: Likewise.
* math/log10q.c: Likewise.
* math/rintq.c: Likewise.
* math/roundq.c: Likewise.
* math/fmaq.c: Likewise.
* math/erfq.c: Likewise.
* math/log2q.c: Likewise.
* math/lroundq.c: Likewise.
* math/j1q.c: Likewise.
* math/scalbnq.c: Likewise.
* math/truncq.c: Likewise.
* math/frexpq.c: Likewise.
* math/sincosq.c: Likewise.
* math/tanhq.c: Likewise.
* math/asinq.c: Likewise.
* math/coshq.c: Likewise.
* math/j0q.c: Likewise.
* math/asinhq.c: Likewise.
* math/floorq.c: Likewise.
* math/sinq_kernel.c: Likewise.
* math/powq.c: Likewise.
* math/hypotq.c: Likewise.
* math/sincos_table.c: Likewise.
* math/rem_pio2q.c: Likewise.
* math/nextafterq.c: Likewise.
* math/log1pq.c: Likewise.
* math/sincosq_kernel.c: Likewise.
* math/tanq.c: Likewise.
* math/acosq.c: Likewise.
* math/lrintq.c: Likewise.
* math/llrintq.c: Likewise.
From-SVN: r250343
Diffstat (limited to 'libquadmath')
44 files changed, 789 insertions, 477 deletions
diff --git a/libquadmath/ChangeLog b/libquadmath/ChangeLog index 4bd166b0ff9..73d0e3077c5 100644 --- a/libquadmath/ChangeLog +++ b/libquadmath/ChangeLog @@ -1,3 +1,53 @@ +2017-07-19 Jakub Jelinek <jakub@redhat.com> + + PR libquadmath/65757 + * quadmath-imp.h (math_opt_barrier, math_force_eval, + math_narrow_eval, math_check_force_underflow, + math_check_force_underflow_nonneg): Define. + * math/ceilq.c: Backport changes from upstream glibc + between 2012-11-01 and 2017-07-13. + * math/remquoq.c: Likewise. + * math/expq.c: Likewise. + * math/llroundq.c: Likewise. + * math/logq.c: Likewise. + * math/atanq.c: Likewise. + * math/nearbyintq.c: Likewise. + * math/scalblnq.c: Likewise. + * math/finiteq.c: Likewise. + * math/atanhq.c: Likewise. + * math/expm1q.c: Likewise. + * math/sinhq.c: Likewise. + * math/log10q.c: Likewise. + * math/rintq.c: Likewise. + * math/roundq.c: Likewise. + * math/fmaq.c: Likewise. + * math/erfq.c: Likewise. + * math/log2q.c: Likewise. + * math/lroundq.c: Likewise. + * math/j1q.c: Likewise. + * math/scalbnq.c: Likewise. + * math/truncq.c: Likewise. + * math/frexpq.c: Likewise. + * math/sincosq.c: Likewise. + * math/tanhq.c: Likewise. + * math/asinq.c: Likewise. + * math/coshq.c: Likewise. + * math/j0q.c: Likewise. + * math/asinhq.c: Likewise. + * math/floorq.c: Likewise. + * math/sinq_kernel.c: Likewise. + * math/powq.c: Likewise. + * math/hypotq.c: Likewise. + * math/sincos_table.c: Likewise. + * math/rem_pio2q.c: Likewise. + * math/nextafterq.c: Likewise. + * math/log1pq.c: Likewise. + * math/sincosq_kernel.c: Likewise. + * math/tanq.c: Likewise. + * math/acosq.c: Likewise. + * math/lrintq.c: Likewise. + * math/llrintq.c: Likewise. + 2017-02-09 Gerald Pfeifer <gerald@pfeifer.com> * configure.ac (ACX_BUGURL): Update. diff --git a/libquadmath/math/acosq.c b/libquadmath/math/acosq.c index 7ef79474651..7f2ed2725d1 100644 --- a/libquadmath/math/acosq.c +++ b/libquadmath/math/acosq.c @@ -172,7 +172,7 @@ acosq (__float128 x) } else if (ix < 0x3ffe0000) /* |x| < 0.5 */ { - if (ix < 0x3fc60000) /* |x| < 2**-57 */ + if (ix < 0x3f8e0000) /* |x| < 2**-113 */ return pio2_hi + pio2_lo; if (ix < 0x3ffde000) /* |x| < .4375 */ { diff --git a/libquadmath/math/asinhq.c b/libquadmath/math/asinhq.c index 9be0aa1f053..f3644939b3a 100644 --- a/libquadmath/math/asinhq.c +++ b/libquadmath/math/asinhq.c @@ -46,6 +46,7 @@ asinhq (__float128 x) return x + x; /* x is inf or NaN */ if (ix < 0x3fc70000) { /* |x| < 2^ -56 */ + math_check_force_underflow (x); if (huge + x > one) return x; /* return x inexact except 0 */ } diff --git a/libquadmath/math/asinq.c b/libquadmath/math/asinq.c index 7bd4d768c97..5dff2817694 100644 --- a/libquadmath/math/asinq.c +++ b/libquadmath/math/asinq.c @@ -151,8 +151,10 @@ asinq (__float128 x) { if (ix < 0x3fc60000) /* |x| < 2**-57 */ { - if (huge + x > one) - return x; /* return x with inexact if x!=0 */ + math_check_force_underflow (x); + __float128 force_inexact = huge + x; + math_force_eval (force_inexact); + return x; /* return x with inexact if x!=0 */ } else { diff --git a/libquadmath/math/atanhq.c b/libquadmath/math/atanhq.c index b34036715d7..652138d495c 100644 --- a/libquadmath/math/atanhq.c +++ b/libquadmath/math/atanhq.c @@ -55,7 +55,11 @@ atanhq (__float128 x) else return (x-x)/(x-x); } - if(ix<0x3fc60000 && (huge+x)>zero) return x; /* x < 2^-57 */ + if(ix<0x3fc60000 && (huge+x)>zero) /* x < 2^-57 */ + { + math_check_force_underflow (x); + return x; + } if(ix<0x3ffe0000) { /* x < 0.5 */ t = u.value+u.value; diff --git a/libquadmath/math/atanq.c b/libquadmath/math/atanq.c index 8eccdc3317d..01fd9d69e57 100644 --- a/libquadmath/math/atanq.c +++ b/libquadmath/math/atanq.c @@ -42,7 +42,7 @@ * */ -/* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov> +/* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov> This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -167,7 +167,7 @@ static const __float128 q4 = 2.173623741810414221251136181221172551416E1Q; /* q5 = 1.000000000000000000000000000000000000000E0 */ -static const long double huge = 1.0e4930Q; +static const __float128 huge = 1.0e4930Q; __float128 atanq (__float128 x) @@ -200,6 +200,7 @@ atanq (__float128 x) if (k <= 0x3fc50000) /* |x| < 2**-58 */ { + math_check_force_underflow (x); /* Raise inexact. */ if (huge + x > 0.0) return x; diff --git a/libquadmath/math/ceilq.c b/libquadmath/math/ceilq.c index 0d9bb8b87f3..1adc1e1b9f0 100644 --- a/libquadmath/math/ceilq.c +++ b/libquadmath/math/ceilq.c @@ -15,8 +15,6 @@ #include "quadmath-imp.h" -static const __float128 huge = 1.0e4930Q; - __float128 ceilq (__float128 x) { @@ -25,18 +23,15 @@ ceilq (__float128 x) GET_FLT128_WORDS64(i0,i1,x); j0 = ((i0>>48)&0x7fff)-0x3fff; if(j0<48) { - if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0<0) {i0=0x8000000000000000ULL;i1=0;} - else if((i0|i1)!=0) { i0=0x3fff000000000000ULL;i1=0;} - } + if(j0<0) { + /* return 0*sign(x) if |x|<1 */ + if(i0<0) {i0=0x8000000000000000ULL;i1=0;} + else if((i0|i1)!=0) { i0=0x3fff000000000000ULL;i1=0;} } else { i = (0x0000ffffffffffffULL)>>j0; if(((i0&i)|i1)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0>0) i0 += (0x0001000000000000LL)>>j0; - i0 &= (~i); i1=0; - } + if(i0>0) i0 += (0x0001000000000000LL)>>j0; + i0 &= (~i); i1=0; } } else if (j0>111) { if(j0==0x4000) return x+x; /* inf or NaN */ @@ -44,17 +39,15 @@ ceilq (__float128 x) } else { i = -1ULL>>(j0-48); if((i1&i)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0>0) { - if(j0==48) i0+=1; - else { - j = i1+(1LL<<(112-j0)); - if(j<i1) i0 +=1 ; /* got a carry */ - i1=j; - } + if(i0>0) { + if(j0==48) i0+=1; + else { + j = i1+(1LL<<(112-j0)); + if(j<i1) i0 +=1 ; /* got a carry */ + i1=j; } - i1 &= (~i); } + i1 &= (~i); } SET_FLT128_WORDS64(x,i0,i1); return x; diff --git a/libquadmath/math/coshq.c b/libquadmath/math/coshq.c index 77f4b98338b..61397509527 100644 --- a/libquadmath/math/coshq.c +++ b/libquadmath/math/coshq.c @@ -76,10 +76,10 @@ coshq (__float128 x) /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expq(|x|)) */ if (ex < 0x3ffd62e4) /* 0.3465728759765625 */ { + if (ex < 0x3fb80000) /* |x| < 2^-116 */ + return one; /* cosh(tiny) = 1 */ t = expm1q (u.value); w = one + t; - if (ex < 0x3fb80000) /* |x| < 2^-116 */ - return w; /* cosh(tiny) = 1 */ return one + (t * t) / (w + w); } diff --git a/libquadmath/math/erfq.c b/libquadmath/math/erfq.c index 8d383e9ca70..45a8c014e5c 100644 --- a/libquadmath/math/erfq.c +++ b/libquadmath/math/erfq.c @@ -11,9 +11,9 @@ /* Modifications and expansions for 128-bit long double are Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov> - and are incorporated herein by permission of the author. The author + and are incorporated herein by permission of the author. The author reserves the right to distribute this material elsewhere under different - copying permissions. These modifications are distributed here under + copying permissions. These modifications are distributed here under the following terms: This library is free software; you can redistribute it and/or @@ -96,6 +96,7 @@ * erfc/erf(NaN) is NaN */ +#include <errno.h> #include "quadmath-imp.h" @@ -142,13 +143,10 @@ deval (__float128 x, const __float128 *p, int n) static const __float128 tiny = 1e-4931Q, - half = 0.5Q, one = 1.0Q, two = 2.0Q, /* 2/sqrt(pi) - 1 */ - efx = 1.2837916709551257389615890312154517168810E-1Q, - /* 8 * (2/sqrt(pi) - 1) */ - efx8 = 1.0270333367641005911692712249723613735048E0Q; + efx = 1.2837916709551257389615890312154517168810E-1Q; /* erf(x) = x + x R(x^2) @@ -773,6 +771,8 @@ erfq (__float128 x) if (ix >= 0x3fff0000) /* |x| >= 1.0 */ { + if (ix >= 0x40030000 && sign > 0) + return one; /* x >= 16, avoid spurious underflow from erfc. */ y = erfcq (x); return (one - y); /* return (one - erfcq (x)); */ @@ -785,7 +785,12 @@ erfq (__float128 x) if (ix < 0x3fc60000) /* |x|<2**-57 */ { if (ix < 0x00080000) - return 0.125 * (8.0 * x + efx8 * x); /*avoid underflow */ + { + /* Avoid spurious underflow. */ + __float128 ret = 0.0625 * (16.0 * x + (16.0 * efx) * x); + math_check_force_underflow (ret); + return ret; + } return x + efx * x; } y = a + a * neval (z, TN1, NTN1) / deval (z, TD1, NTD1); @@ -867,7 +872,7 @@ erfcq (__float128 x) y = C19b + z * neval (z, RNr19, NRNr19) / deval (z, RDr19, NRDr19); y += C19a; break; - case 9: + default: /* i == 9. */ z = x - 1.125Q; y = C20b + z * neval (z, RNr20, NRNr20) / deval (z, RDr20, NRDr20); y += C20a; @@ -921,14 +926,22 @@ erfcq (__float128 x) z = u.value; r = expq (-z * z - 0.5625) * expq ((z - x) * (z + x) + p); if ((sign & 0x80000000) == 0) - return r / x; + { + __float128 ret = r / x; + if (ret == 0) + errno = ERANGE; + return ret; + } else return two - r / x; } else { if ((sign & 0x80000000) == 0) - return tiny * tiny; + { + errno = ERANGE; + return tiny * tiny; + } else return two - tiny; } diff --git a/libquadmath/math/expm1q.c b/libquadmath/math/expm1q.c index 8cfdd8eec94..9060d480858 100644 --- a/libquadmath/math/expm1q.c +++ b/libquadmath/math/expm1q.c @@ -35,7 +35,7 @@ * */ -/* Copyright 2001 by Stephen L. Moshier +/* Copyright 2001 by Stephen L. Moshier This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -82,8 +82,6 @@ static const __float128 C1 = 6.93145751953125E-1Q, C2 = 1.428606820309417232121458176568075500134E-6Q, -/* ln (2^16384 * (1 - 2^-113)) */ - maxlog = 1.1356523406294143949491931077970764891253E4Q, /* ln 2^-114 */ minarg = -7.9018778583833765273564461846232128760607E1Q; @@ -108,33 +106,30 @@ expm1q (__float128 x) } if (ix >= 0x7fff0000) { - /* Infinity. */ + /* Infinity (which must be negative infinity). */ if (((ix & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3) == 0) - { - if (sign) - return -1.0Q; - else - return x; - } - /* NaN. No invalid exception. */ - return x; + return -1.0Q; + /* NaN. Invalid exception if signaling. */ + return x + x; } /* expm1(+- 0) = +- 0. */ if ((ix == 0) && (u.words32.w1 | u.words32.w2 | u.words32.w3) == 0) return x; - /* Overflow. */ - if (x > maxlog) - { - errno = ERANGE; - return (HUGE_VALQ * HUGE_VALQ); - } - /* Minimum value. */ if (x < minarg) return (4.0/HUGE_VALQ - 1.0Q); + /* Avoid internal underflow when result does not underflow, while + ensuring underflow (without returning a zero of the wrong sign) + when the result does underflow. */ + if (fabsq (x) < 0x1p-113Q) + { + math_check_force_underflow (x); + return x; + } + /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */ xx = C1 + C2; /* ln 2. */ px = floorq (0.5 + x / xx); diff --git a/libquadmath/math/expq.c b/libquadmath/math/expq.c index 70c638d7a08..5df6cd8e192 100644 --- a/libquadmath/math/expq.c +++ b/libquadmath/math/expq.c @@ -1,5 +1,5 @@ /* Quad-precision floating point e^x. - Copyright (C) 1999 Free Software Foundation, Inc. + Copyright (C) 1999-2017 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jj@ultra.linux.cz> Partly based on double-precision code @@ -1075,7 +1075,7 @@ static const __float128 C[] = { #define TWO15 C[11] 32768.0Q, -/* Chebyshev polynom coeficients for (exp(x)-1)/x */ +/* Chebyshev polynom coefficients for (exp(x)-1)/x */ #define P1 C[12] #define P2 C[13] #define P3 C[14] @@ -1142,7 +1142,7 @@ expq (__float128 x) * __expq_table[T_EXPL_RES2 + tval2]; n_i = (int)n; /* 'unsafe' is 1 iff n_1 != 0. */ - unsafe = abs(n_i) >= -FLT128_MIN_EXP - 1; + unsafe = abs(n_i) >= 15000; ex2_u.ieee.exponent += n_i >> unsafe; /* Compute scale = 2^n_1. */ @@ -1179,7 +1179,7 @@ expq (__float128 x) ex3_u.d = (result - ex2_u.d) - x22 * ex2_u.d; ex2_u.d = result; ex3_u.ieee.exponent += LDBL_MANT_DIG + 15 + IEEE854_LONG_DOUBLE_BIAS - - ex2_u.ieee.exponent; + - ex2_u.ieee.exponent; n_i = abs (ex3_u.d); n_i = (n_i + 1) / 2; #ifdef USE_FENV_H @@ -1196,7 +1196,11 @@ expq (__float128 x) if (!unsafe) return result; else - return result * scale_u.value; + { + result *= scale_u.value; + math_check_force_underflow_nonneg (result); + return result; + } } /* Exceptional cases: */ else if (__builtin_isless (x, himark)) diff --git a/libquadmath/math/finiteq.c b/libquadmath/math/finiteq.c index c5326d2194b..e6703fb2261 100644 --- a/libquadmath/math/finiteq.c +++ b/libquadmath/math/finiteq.c @@ -25,6 +25,6 @@ finiteq (const __float128 x) { int64_t hx; GET_FLT128_MSW64(hx,x); - return (int)((uint64_t)((hx&0x7fffffffffffffffLL) + return (int)((uint64_t)((hx&0x7fff000000000000LL) -0x7fff000000000000LL)>>63); } diff --git a/libquadmath/math/floorq.c b/libquadmath/math/floorq.c index 0d74f94e27d..41b993fa7a0 100644 --- a/libquadmath/math/floorq.c +++ b/libquadmath/math/floorq.c @@ -15,8 +15,6 @@ #include "quadmath-imp.h" -static const __float128 huge = 1.0e4930Q; - __float128 floorq (__float128 x) { @@ -25,19 +23,16 @@ floorq (__float128 x) GET_FLT128_WORDS64(i0,i1,x); j0 = ((i0>>48)&0x7fff)-0x3fff; if(j0<48) { - if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0>=0) {i0=i1=0;} - else if(((i0&0x7fffffffffffffffLL)|i1)!=0) - { i0=0xbfff000000000000ULL;i1=0;} - } + if(j0<0) { + /* return 0*sign(x) if |x|<1 */ + if(i0>=0) {i0=i1=0;} + else if(((i0&0x7fffffffffffffffLL)|i1)!=0) + { i0=0xbfff000000000000ULL;i1=0;} } else { i = (0x0000ffffffffffffULL)>>j0; if(((i0&i)|i1)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0<0) i0 += (0x0001000000000000LL)>>j0; - i0 &= (~i); i1=0; - } + if(i0<0) i0 += (0x0001000000000000LL)>>j0; + i0 &= (~i); i1=0; } } else if (j0>111) { if(j0==0x4000) return x+x; /* inf or NaN */ @@ -45,17 +40,15 @@ floorq (__float128 x) } else { i = -1ULL>>(j0-48); if((i1&i)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0<0) { - if(j0==48) i0+=1; - else { - j = i1+(1LL<<(112-j0)); - if(j<i1) i0 +=1 ; /* got a carry */ - i1=j; - } + if(i0<0) { + if(j0==48) i0+=1; + else { + j = i1+(1LL<<(112-j0)); + if(j<i1) i0 +=1 ; /* got a carry */ + i1=j; } - i1 &= (~i); } + i1 &= (~i); } SET_FLT128_WORDS64(x,i0,i1); return x; diff --git a/libquadmath/math/fmaq.c b/libquadmath/math/fmaq.c index d3c5fb3901a..68a63cf8fd6 100644 --- a/libquadmath/math/fmaq.c +++ b/libquadmath/math/fmaq.c @@ -1,5 +1,5 @@ /* Compute x * y + z as ternary operation. - Copyright (C) 2010-2012 Free Software Foundation, Inc. + Copyright (C) 2010-2017 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jakub@redhat.com>, 2010. @@ -97,8 +97,8 @@ fmaq (__float128 x, __float128 y, __float128 z) && w.ieee.mant_low == 0 && w.ieee.mant_high == 0))) { - volatile __float128 force_underflow = x * y; - (void) force_underflow; + __float128 force_underflow = x * y; + math_force_eval (force_underflow); } return v.value * 0x1p-114Q; } @@ -161,15 +161,15 @@ fmaq (__float128 x, __float128 y, __float128 z) <= IEEE854_FLOAT128_BIAS + FLT128_MANT_DIG) */ { if (u.ieee.exponent > v.ieee.exponent) - u.ieee.exponent += 2 * FLT128_MANT_DIG; + u.ieee.exponent += 2 * FLT128_MANT_DIG + 2; else - v.ieee.exponent += 2 * FLT128_MANT_DIG; - if (w.ieee.exponent <= 4 * FLT128_MANT_DIG + 4) + v.ieee.exponent += 2 * FLT128_MANT_DIG + 2; + if (w.ieee.exponent <= 4 * FLT128_MANT_DIG + 6) { if (w.ieee.exponent) - w.ieee.exponent += 2 * FLT128_MANT_DIG; + w.ieee.exponent += 2 * FLT128_MANT_DIG + 2; else - w.value *= 0x1p226Q; + w.value *= 0x1p228Q; adjust = -1; } /* Otherwise x * y should just affect inexact @@ -182,7 +182,10 @@ fmaq (__float128 x, __float128 y, __float128 z) /* Ensure correct sign of exact 0 + 0. */ if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0)) - return x * y + z; + { + x = math_opt_barrier (x); + return x * y + z; + } #ifdef USE_FENV_H fenv_t env; @@ -208,24 +211,24 @@ fmaq (__float128 x, __float128 y, __float128 z) t1 = m1 - t1; t2 = z - t2; __float128 a2 = t1 + t2; + /* Ensure the arithmetic is not scheduled after feclearexcept call. */ + math_force_eval (m2); + math_force_eval (a2); #ifdef USE_FENV_H feclearexcept (FE_INEXACT); #endif - /* If the result is an exact zero, ensure it has the correct - sign. */ + /* If the result is an exact zero, ensure it has the correct sign. */ if (a1 == 0 && m2 == 0) { #ifdef USE_FENV_H feupdateenv (&env); #endif - /* Ensure that round-to-nearest value of z + m1 is not - reused. */ - asm volatile ("" : "=m" (z) : "m" (z)); + /* Ensure that round-to-nearest value of z + m1 is not reused. */ + z = math_opt_barrier (z); return z + m1; } - #ifdef USE_FENV_H fesetround (FE_TOWARDZERO); #endif @@ -273,19 +276,19 @@ fmaq (__float128 x, __float128 y, __float128 z) /* If a1 + u.value is exact, the only rounding happens during scaling down. */ if (j == 0) - return v.value * 0x1p-226Q; + return v.value * 0x1p-228Q; /* If result rounded to zero is not subnormal, no double rounding will occur. */ - if (v.ieee.exponent > 226) - return (a1 + u.value) * 0x1p-226Q; - /* If v.value * 0x1p-226Q with round to zero is a subnormal above - or equal to FLT128_MIN / 2, then v.value * 0x1p-226Q shifts mantissa + if (v.ieee.exponent > 228) + return (a1 + u.value) * 0x1p-228Q; + /* If v.value * 0x1p-228Q with round to zero is a subnormal above + or equal to FLT128_MIN / 2, then v.value * 0x1p-228Q shifts mantissa down just by 1 bit, which means v.ieee.mant_low |= j would change the round bit, not sticky or guard bit. - v.value * 0x1p-226Q never normalizes by shifting up, + v.value * 0x1p-228Q never normalizes by shifting up, so round bit plus sticky bit should be already enough for proper rounding. */ - if (v.ieee.exponent == 226) + if (v.ieee.exponent == 228) { /* If the exponent would be in the normal range when rounding to normal precision with unbounded exponent @@ -295,8 +298,8 @@ fmaq (__float128 x, __float128 y, __float128 z) if (TININESS_AFTER_ROUNDING) { w.value = a1 + u.value; - if (w.ieee.exponent == 227) - return w.value * 0x1p-226Q; + if (w.ieee.exponent == 229) + return w.value * 0x1p-228Q; } /* v.ieee.mant_low & 2 is LSB bit of the result before rounding, v.ieee.mant_low & 1 is the round bit and j is our sticky @@ -305,11 +308,11 @@ fmaq (__float128 x, __float128 y, __float128 z) w.ieee.mant_low = ((v.ieee.mant_low & 3) << 1) | j; w.ieee.negative = v.ieee.negative; v.ieee.mant_low &= ~3U; - v.value *= 0x1p-226Q; + v.value *= 0x1p-228Q; w.value *= 0x1p-2Q; return v.value + w.value; } v.ieee.mant_low |= j; - return v.value * 0x1p-226Q; + return v.value * 0x1p-228Q; } } diff --git a/libquadmath/math/frexpq.c b/libquadmath/math/frexpq.c index b0b305af574..2bd77829bf2 100644 --- a/libquadmath/math/frexpq.c +++ b/libquadmath/math/frexpq.c @@ -35,7 +35,7 @@ frexpq (__float128 x, int *eptr) GET_FLT128_WORDS64(hx,lx,x); ix = 0x7fffffffffffffffULL&hx; *eptr = 0; - if(ix>=0x7fff000000000000ULL||((ix|lx)==0)) return x; /* 0,inf,nan */ + if(ix>=0x7fff000000000000ULL||((ix|lx)==0)) return x + x;/* 0,inf,nan */ if (ix<0x0001000000000000ULL) { /* subnormal */ x *= two114; GET_FLT128_MSW64(hx,x); diff --git a/libquadmath/math/hypotq.c b/libquadmath/math/hypotq.c index 2df317f3681..057901073dc 100644 --- a/libquadmath/math/hypotq.c +++ b/libquadmath/math/hypotq.c @@ -89,6 +89,17 @@ hypotq (__float128 x, __float128 y) b *= t1; a *= t1; k -= 16382; + GET_FLT128_MSW64 (ha, a); + GET_FLT128_MSW64 (hb, b); + if (hb > ha) + { + t1 = a; + a = b; + b = t1; + j = ha; + ha = hb; + hb = j; + } } else { /* scale a and b by 2^9600 */ ha += 0x2580000000000000LL; /* a *= 2^9600 */ hb += 0x2580000000000000LL; /* b *= 2^9600 */ @@ -119,6 +130,8 @@ hypotq (__float128 x, __float128 y) t1 = 1.0Q; GET_FLT128_MSW64(high,t1); SET_FLT128_MSW64(t1,high+(k<<48)); - return t1*w; + w *= t1; + math_check_force_underflow_nonneg (w); + return w; } else return w; } diff --git a/libquadmath/math/j0q.c b/libquadmath/math/j0q.c index 8c6f811125e..c6e482b1c51 100644 --- a/libquadmath/math/j0q.c +++ b/libquadmath/math/j0q.c @@ -681,7 +681,7 @@ j0q (__float128 x) if (! finiteq (x)) { if (x != x) - return x; + return x + x; else return 0.0Q; } @@ -691,6 +691,8 @@ j0q (__float128 x) xx = fabsq (x); if (xx <= 2.0Q) { + if (xx < 0x1p-57Q) + return 1.0Q; /* 0 <= x <= 2 */ z = xx * xx; p = z * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D); @@ -699,6 +701,28 @@ j0q (__float128 x) return p; } + /* X = x - pi/4 + cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4) + = 1/sqrt(2) * (cos(x) + sin(x)) + sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4) + = 1/sqrt(2) * (sin(x) - cos(x)) + sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) + cf. Fdlibm. */ + sincosq (xx, &s, &c); + ss = s - c; + cc = s + c; + if (xx <= FLT128_MAX / 2.0Q) + { + z = -cosq (xx + xx); + if ((s * c) < 0) + cc = z / ss; + else + ss = z / cc; + } + + if (xx > 0x1p256Q) + return ONEOSQPI * cc / sqrtq (xx); + xinv = 1.0Q / xx; z = xinv * xinv; if (xinv <= 0.25) @@ -760,21 +784,6 @@ j0q (__float128 x) p = 1.0Q + z * p; q = z * xinv * q; q = q - 0.125Q * xinv; - /* X = x - pi/4 - cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4) - = 1/sqrt(2) * (cos(x) + sin(x)) - sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4) - = 1/sqrt(2) * (sin(x) - cos(x)) - sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - cf. Fdlibm. */ - sincosq (xx, &s, &c); - ss = s - c; - cc = s + c; - z = - cosq (xx + xx); - if ((s * c) < 0) - cc = z / ss; - else - ss = z / cc; z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx); return z; } @@ -817,17 +826,12 @@ y0q (__float128 x) __float128 xx, xinv, z, p, q, c, s, cc, ss; if (! finiteq (x)) - { - if (x != x) - return x; - else - return 0.0Q; - } + return 1 / (x + x * x); if (x <= 0.0Q) { if (x < 0.0Q) return (zero / (zero * x)); - return -HUGE_VALQ + x; + return -1 / zero; /* -inf and divide by zero exception. */ } xx = fabsq (x); if (xx <= 0x1p-57) @@ -841,6 +845,28 @@ y0q (__float128 x) return p; } + /* X = x - pi/4 + cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4) + = 1/sqrt(2) * (cos(x) + sin(x)) + sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4) + = 1/sqrt(2) * (sin(x) - cos(x)) + sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) + cf. Fdlibm. */ + sincosq (x, &s, &c); + ss = s - c; + cc = s + c; + if (xx <= FLT128_MAX / 2.0Q) + { + z = -cosq (x + x); + if ((s * c) < 0) + cc = z / ss; + else + ss = z / cc; + } + + if (xx > 0x1p256Q) + return ONEOSQPI * ss / sqrtq (x); + xinv = 1.0Q / xx; z = xinv * xinv; if (xinv <= 0.25) @@ -902,21 +928,6 @@ y0q (__float128 x) p = 1.0Q + z * p; q = z * xinv * q; q = q - 0.125Q * xinv; - /* X = x - pi/4 - cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4) - = 1/sqrt(2) * (cos(x) + sin(x)) - sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4) - = 1/sqrt(2) * (sin(x) - cos(x)) - sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - cf. Fdlibm. */ - sincosq (x, &s, &c); - ss = s - c; - cc = s + c; - z = - cosq (x + x); - if ((s * c) < 0) - cc = z / ss; - else - ss = z / cc; z = ONEOSQPI * (p * ss + q * cc) / sqrtq (x); return z; } diff --git a/libquadmath/math/j1q.c b/libquadmath/math/j1q.c index eb599c949a9..5eb705084e2 100644 --- a/libquadmath/math/j1q.c +++ b/libquadmath/math/j1q.c @@ -95,6 +95,7 @@ License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ +#include <errno.h> #include "quadmath-imp.h" /* 1 / sqrt(pi) */ @@ -687,13 +688,21 @@ j1q (__float128 x) if (! finiteq (x)) { if (x != x) - return x; + return x + x; else return 0.0Q; } if (x == 0.0Q) return x; xx = fabsq (x); + if (xx <= 0x1p-58Q) + { + __float128 ret = x * 0.5Q; + math_check_force_underflow (ret); + if (ret == 0) + errno = ERANGE; + return ret; + } if (xx <= 2.0Q) { /* 0 <= x <= 2 */ @@ -705,6 +714,32 @@ j1q (__float128 x) return p; } + /* X = x - 3 pi/4 + cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4) + = 1/sqrt(2) * (-cos(x) + sin(x)) + sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4) + = -1/sqrt(2) * (sin(x) + cos(x)) + cf. Fdlibm. */ + sincosq (xx, &s, &c); + ss = -s - c; + cc = s - c; + if (xx <= FLT128_MAX / 2.0Q) + { + z = cosq (xx + xx); + if ((s * c) > 0) + cc = z / ss; + else + ss = z / cc; + } + + if (xx > 0x1p256Q) + { + z = ONEOSQPI * cc / sqrtq (xx); + if (x < 0) + z = -z; + return z; + } + xinv = 1.0Q / xx; z = xinv * xinv; if (xinv <= 0.25) @@ -766,20 +801,6 @@ j1q (__float128 x) p = 1.0Q + z * p; q = z * q; q = q * xinv + 0.375Q * xinv; - /* X = x - 3 pi/4 - cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4) - = 1/sqrt(2) * (-cos(x) + sin(x)) - sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4) - = -1/sqrt(2) * (sin(x) + cos(x)) - cf. Fdlibm. */ - sincosq (xx, &s, &c); - ss = -s - c; - cc = s - c; - z = cosq (xx + xx); - if ((s * c) > 0) - cc = z / ss; - else - ss = z / cc; z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx); if (x < 0) z = -z; @@ -823,24 +844,25 @@ y1q (__float128 x) __float128 xx, xinv, z, p, q, c, s, cc, ss; if (! finiteq (x)) - { - if (x != x) - return x; - else - return 0.0Q; - } + return 1 / (x + x * x); if (x <= 0.0Q) { if (x < 0.0Q) return (zero / (zero * x)); - return -HUGE_VALQ + x; + return -1 / zero; /* -inf and divide by zero exception. */ } xx = fabsq (x); if (xx <= 0x1p-114) - return -TWOOPI / x; + { + z = -TWOOPI / x; + if (isinfq (z)) + errno = ERANGE; + return z; + } if (xx <= 2.0Q) { /* 0 <= x <= 2 */ + /* FIXME: SET_RESTORE_ROUNDL (FE_TONEAREST); */ z = xx * xx; p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D); p = -TWOOPI / xx + p; @@ -848,6 +870,27 @@ y1q (__float128 x) return p; } + /* X = x - 3 pi/4 + cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4) + = 1/sqrt(2) * (-cos(x) + sin(x)) + sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4) + = -1/sqrt(2) * (sin(x) + cos(x)) + cf. Fdlibm. */ + sincosq (xx, &s, &c); + ss = -s - c; + cc = s - c; + if (xx <= FLT128_MAX / 2.0Q) + { + z = cosq (xx + xx); + if ((s * c) > 0) + cc = z / ss; + else + ss = z / cc; + } + + if (xx > 0x1p256Q) + return ONEOSQPI * ss / sqrtq (xx); + xinv = 1.0Q / xx; z = xinv * xinv; if (xinv <= 0.25) @@ -909,20 +952,6 @@ y1q (__float128 x) p = 1.0Q + z * p; q = z * q; q = q * xinv + 0.375Q * xinv; - /* X = x - 3 pi/4 - cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4) - = 1/sqrt(2) * (-cos(x) + sin(x)) - sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4) - = -1/sqrt(2) * (sin(x) + cos(x)) - cf. Fdlibm. */ - sincosq (xx, &s, &c); - ss = -s - c; - cc = s - c; - z = cosq (xx + xx); - if ((s * c) > 0) - cc = z / ss; - else - ss = z / cc; z = ONEOSQPI * (p * ss + q * cc) / sqrtq (xx); return z; } diff --git a/libquadmath/math/llrintq.c b/libquadmath/math/llrintq.c index eef31d823b6..a6a0ae64bd3 100644 --- a/libquadmath/math/llrintq.c +++ b/libquadmath/math/llrintq.c @@ -1,9 +1,9 @@ /* Round argument to nearest integral value according to current rounding direction. - Copyright (C) 1997, 1999, 2006 Free Software Foundation, Inc. + Copyright (C) 1997-2017 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and - Jakub Jelinek <jj@ultra.linux.cz>, 1999. + Jakub Jelinek <jj@ultra.linux.cz>, 1999. The GNU C Library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -33,7 +33,7 @@ llrintq (__float128 x) { int32_t j0; uint64_t i0,i1; - volatile __float128 w; + __float128 w; __float128 t; long long int result; int sx; @@ -46,8 +46,23 @@ llrintq (__float128 x) if (j0 < (int32_t) (8 * sizeof (long long int)) - 1) { - w = two112[sx] + x; - t = w - two112[sx]; +#if defined FE_INVALID || defined FE_INEXACT + /* X < LLONG_MAX + 1 implied by J0 < 63. */ + if (x > (__float128) LLONG_MAX) + { + /* In the event of overflow we must raise the "invalid" + exception, but not "inexact". */ + t = nearbyintq (x); +#ifdef USE_FENV_H + feraiseexcept (t == LLONG_MAX ? FE_INEXACT : FE_INVALID); +#endif + } + else +#endif + { + w = two112[sx] + x; + t = w - two112[sx]; + } GET_FLT128_WORDS64 (i0, i1, t); j0 = ((i0 >> 48) & 0x7fff) - 0x3fff; i0 &= 0x0000ffffffffffffLL; @@ -62,6 +77,24 @@ llrintq (__float128 x) } else { + /* The number is too large. Unless it rounds to LLONG_MIN, + FE_INVALID must be raised and the return value is + unspecified. */ +#if defined FE_INVALID || defined FE_INEXACT + if (x < (__float128) LLONG_MIN + && x > (__float128) LLONG_MIN - 1.0Q) + { + /* If truncation produces LLONG_MIN, the cast will not raise + the exception, but may raise "inexact". */ + t = nearbyintq (x); +#ifdef USE_FENV_H + feraiseexcept (t == LLONG_MIN ? FE_INEXACT : FE_INVALID); +#endif + return LLONG_MIN; + } + +#endif + /* The number is too large. It is left implementation defined what happens. */ return (long long int) x; diff --git a/libquadmath/math/llroundq.c b/libquadmath/math/llroundq.c index d22180d6bba..098fb9ef72b 100644 --- a/libquadmath/math/llroundq.c +++ b/libquadmath/math/llroundq.c @@ -1,8 +1,8 @@ /* Round __float128 value to long long int. - Copyright (C) 1997, 1999, 2004 Free Software Foundation, Inc. + Copyright (C) 1997-2017 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and - Jakub Jelinek <jj@ultra.linux.cz>, 1999. + Jakub Jelinek <jj@ultra.linux.cz>, 1999. The GNU C Library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -59,13 +59,32 @@ llroundq (__float128 x) if (j0 == 48) result = (long long int) i0; else - result = ((long long int) i0 << (j0 - 48)) | (j >> (112 - j0)); + { + result = ((long long int) i0 << (j0 - 48)) | (j >> (112 - j0)); +#if defined FE_INVALID && defined USE_FENV_H + if (sign == 1 && result == LLONG_MIN) + /* Rounding brought the value out of range. */ + feraiseexcept (FE_INVALID); +#endif + } } } else { - /* The number is too large. It is left implementation defined - what happens. */ + /* The number is too large. Unless it rounds to LLONG_MIN, + FE_INVALID must be raised and the return value is + unspecified. */ +#ifdef FE_INVALID + if (x <= (__float128) LLONG_MIN - 0.5Q) + { + /* If truncation produces LLONG_MIN, the cast will not raise + the exception, but may raise "inexact". */ +#ifdef USE_FENV_H + feraiseexcept (FE_INVALID); +#endif + return LLONG_MIN; + } +#endif return (long long int) x; } diff --git a/libquadmath/math/log10q.c b/libquadmath/math/log10q.c index 9eeb9ae3fc4..3afb1121267 100644 --- a/libquadmath/math/log10q.c +++ b/libquadmath/math/log10q.c @@ -188,12 +188,15 @@ log10q (__float128 x) /* Test for domain */ GET_FLT128_WORDS64 (hx, lx, x); if (((hx & 0x7fffffffffffffffLL) | lx) == 0) - return (-1.0Q / (x - x)); + return (-1.0Q / fabsq (x)); /* log10l(+-0)=-inf */ if (hx < 0) return (x - x) / (x - x); if (hx >= 0x7fff000000000000LL) return (x + x); + if (x == 1.0Q) + return 0.0Q; + /* separate mantissa from exponent */ /* Note, frexp is used so that denormal numbers diff --git a/libquadmath/math/log1pq.c b/libquadmath/math/log1pq.c index d8bff405dff..c59ceef5c9e 100644 --- a/libquadmath/math/log1pq.c +++ b/libquadmath/math/log1pq.c @@ -36,7 +36,7 @@ * IEEE -1, 8 100000 1.9e-34 4.3e-35 */ -/* Copyright 2001 by Stephen L. Moshier +/* Copyright 2001 by Stephen L. Moshier This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -128,8 +128,8 @@ log1pq (__float128 xm1) /* Test for NaN or infinity input. */ u.value = xm1; hx = u.words32.w0; - if (hx >= 0x7fff0000) - return xm1; + if ((hx & 0x7fffffff) >= 0x7fff0000) + return xm1 + fabsq (xm1); /* log1p(+- 0) = +- 0. */ if (((hx & 0x7fffffff) == 0) @@ -138,17 +138,21 @@ log1pq (__float128 xm1) if ((hx & 0x7fffffff) < 0x3f8e0000) { + math_check_force_underflow (xm1); if ((int) xm1 == 0) return xm1; } - x = xm1 + 1.0Q; + if (xm1 >= 0x1p113Q) + x = xm1; + else + x = xm1 + 1.0Q; /* log1p(-1) = -inf */ if (x <= 0.0Q) { if (x == 0.0Q) - return (-1.0Q / (x - x)); + return (-1.0Q / zero); /* log1p(-1) = -inf */ else return (zero / (x - x)); } diff --git a/libquadmath/math/log2q.c b/libquadmath/math/log2q.c index f8275369b37..865f341f03d 100644 --- a/libquadmath/math/log2q.c +++ b/libquadmath/math/log2q.c @@ -181,12 +181,15 @@ log2q (__float128 x) /* Test for domain */ GET_FLT128_WORDS64 (hx, lx, x); if (((hx & 0x7fffffffffffffffLL) | lx) == 0) - return (-1.0Q / (x - x)); + return (-1.0Q / fabsq (x)); /* log2l(+-0)=-inf */ if (hx < 0) return (x - x) / (x - x); if (hx >= 0x7fff000000000000LL) return (x + x); + if (x == 1.0Q) + return 0.0Q; + /* separate mantissa from exponent */ /* Note, frexp is used so that denormal numbers diff --git a/libquadmath/math/logq.c b/libquadmath/math/logq.c index 7aae9b101ad..43249bb498a 100644 --- a/libquadmath/math/logq.c +++ b/libquadmath/math/logq.c @@ -212,9 +212,8 @@ logq (__float128 x) } /* Extract exponent and reduce domain to 0.703125 <= u < 1.40625 */ - e = (int) (m >> 16) - (int) 0x3ffe; - m &= 0xffff; - u.words32.w0 = m | 0x3ffe0000; + u.value = frexpq (x, &e); + m = u.words32.w0 & 0xffff; m |= 0x10000; /* Find lookup table index k from high order bits of the significand. */ if (m < 0x16800) @@ -241,6 +240,8 @@ logq (__float128 x) /* On this interval the table is not used due to cancellation error. */ if ((x <= 1.0078125Q) && (x >= 0.9921875Q)) { + if (x == 1.0Q) + return 0.0Q; z = x - 1.0Q; k = 64; t.value = 1.0Q; diff --git a/libquadmath/math/lrintq.c b/libquadmath/math/lrintq.c index d1497ae3882..50a300554b9 100644 --- a/libquadmath/math/lrintq.c +++ b/libquadmath/math/lrintq.c @@ -1,9 +1,9 @@ /* Round argument to nearest integral value according to current rounding direction. - Copyright (C) 1997, 1999, 2004, 2006 Free Software Foundation, Inc. + Copyright (C) 1997-2017 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and - Jakub Jelinek <jj@ultra.linux.cz>, 1999. + Jakub Jelinek <jj@ultra.linux.cz>, 1999. The GNU C Library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -33,7 +33,7 @@ lrintq (__float128 x) { int32_t j0; uint64_t i0,i1; - volatile __float128 w; + __float128 w; __float128 t; long int result; int sx; @@ -44,25 +44,57 @@ lrintq (__float128 x) i0 &= 0x0000ffffffffffffLL; i0 |= 0x0001000000000000LL; - if (j0 < 48) + if (j0 < (int32_t) (8 * sizeof (long int)) - 1) { - w = two112[sx] + x; - t = w - two112[sx]; - GET_FLT128_WORDS64 (i0, i1, t); - j0 = ((i0 >> 48) & 0x7fff) - 0x3fff; - i0 &= 0x0000ffffffffffffLL; - i0 |= 0x0001000000000000LL; + if (j0 < 48) + { +#if defined FE_INVALID || defined FE_INEXACT + /* X < LONG_MAX + 1 implied by J0 < 31. */ + if (sizeof (long int) == 4 + && x > (__float128) LONG_MAX) + { + /* In the event of overflow we must raise the "invalid" + exception, but not "inexact". */ + t = nearbyintq (x); +#ifdef USE_FENV_H + feraiseexcept (t == LONG_MAX ? FE_INEXACT : FE_INVALID); +#endif + } + else +#endif + { + w = two112[sx] + x; + t = w - two112[sx]; + } + GET_FLT128_WORDS64 (i0, i1, t); + j0 = ((i0 >> 48) & 0x7fff) - 0x3fff; + i0 &= 0x0000ffffffffffffLL; + i0 |= 0x0001000000000000LL; - result = (j0 < 0 ? 0 : i0 >> (48 - j0)); - } - else if (j0 < (int32_t) (8 * sizeof (long int)) - 1) - { - if (j0 >= 112) + result = (j0 < 0 ? 0 : i0 >> (48 - j0)); + } + else if (j0 >= 112) result = ((long int) i0 << (j0 - 48)) | (i1 << (j0 - 112)); else { - w = two112[sx] + x; - t = w - two112[sx]; +#if defined FE_INVALID || defined FE_INEXACT + /* X < LONG_MAX + 1 implied by J0 < 63. */ + if (sizeof (long int) == 8 + && x > (__float128) LONG_MAX) + { + /* In the event of overflow we must raise the "invalid" + exception, but not "inexact". */ + t = nearbyintq (x); +#ifdef USE_FENV_H + feraiseexcept (t == LONG_MAX ? FE_INEXACT : FE_INVALID); +#endif + } + else +#endif + { + w = two112[sx] + x; + t = w - two112[sx]; + } GET_FLT128_WORDS64 (i0, i1, t); j0 = ((i0 >> 48) & 0x7fff) - 0x3fff; i0 &= 0x0000ffffffffffffLL; @@ -76,8 +108,22 @@ lrintq (__float128 x) } else { - /* The number is too large. It is left implementation defined - what happens. */ + /* The number is too large. Unless it rounds to LONG_MIN, + FE_INVALID must be raised and the return value is + unspecified. */ +#if defined FE_INVALID || defined FE_INEXACT + if (x < (__float128) LONG_MIN + && x > (__float128) LONG_MIN - 1.0Q) + { + /* If truncation produces LONG_MIN, the cast will not raise + the exception, but may raise "inexact". */ + t = nearbyintq (x); +#ifdef USE_FENV_H + feraiseexcept (t == LONG_MIN ? FE_INEXACT : FE_INVALID); +#endif + return LONG_MIN; + } +#endif return (long int) x; } diff --git a/libquadmath/math/lroundq.c b/libquadmath/math/lroundq.c index 59c883a1464..55285034cec 100644 --- a/libquadmath/math/lroundq.c +++ b/libquadmath/math/lroundq.c @@ -1,8 +1,8 @@ /* Round __float128 value to long int. - Copyright (C) 1997, 1999, 2004 Free Software Foundation, Inc. + Copyright (C) 1997-2017 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and - Jakub Jelinek <jj@ultra.linux.cz>, 1999. + Jakub Jelinek <jj@ultra.linux.cz>, 1999. The GNU C Library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -36,19 +36,26 @@ lroundq (__float128 x) i0 &= 0x0000ffffffffffffLL; i0 |= 0x0001000000000000LL; - if (j0 < 48) + if (j0 < (int32_t) (8 * sizeof (long int)) - 1) { - if (j0 < 0) - return j0 < -1 ? 0 : sign; - else + if (j0 < 48) { - i0 += 0x0000800000000000LL >> j0; - result = i0 >> (48 - j0); + if (j0 < 0) + return j0 < -1 ? 0 : sign; + else + { + i0 += 0x0000800000000000LL >> j0; + result = i0 >> (48 - j0); +#if defined FE_INVALID && defined USE_FENV_H + if (sizeof (long int) == 4 + && sign == 1 + && result == LONG_MIN) + /* Rounding brought the value out of range. */ + feraiseexcept (FE_INVALID); +#endif + } } - } - else if (j0 < (int32_t) (8 * sizeof (long int)) - 1) - { - if (j0 >= 112) + else if (j0 >= 112) result = ((long int) i0 << (j0 - 48)) | (i1 << (j0 - 112)); else { @@ -59,13 +66,34 @@ lroundq (__float128 x) if (j0 == 48) result = (long int) i0; else - result = ((long int) i0 << (j0 - 48)) | (j >> (112 - j0)); + { + result = ((long int) i0 << (j0 - 48)) | (j >> (112 - j0)); +#if defined FE_INVALID && defined USE_FENV_H + if (sizeof (long int) == 8 + && sign == 1 + && result == LONG_MIN) + /* Rounding brought the value out of range. */ + feraiseexcept (FE_INVALID); +#endif + } } } else { - /* The number is too large. 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 (x <= (__float128) LONG_MIN - 0.5Q) + { + /* If truncation produces LONG_MIN, the cast will not raise + the exception, but may raise "inexact". */ +#ifdef USE_FENV_H + feraiseexcept (FE_INVALID); +#endif + return LONG_MIN; + } +#endif return (long int) x; } diff --git a/libquadmath/math/nearbyintq.c b/libquadmath/math/nearbyintq.c index 207124808e7..b250927ea2a 100644 --- a/libquadmath/math/nearbyintq.c +++ b/libquadmath/math/nearbyintq.c @@ -44,7 +44,7 @@ nearbyintq(__float128 x) fenv_t env; #endif int64_t i0,j0,sx; - uint64_t i1; + uint64_t i1 __attribute__ ((unused)); __float128 w,t; GET_FLT128_WORDS64(i0,i1,x); sx = (((uint64_t)i0)>>63); @@ -56,6 +56,7 @@ nearbyintq(__float128 x) #endif w = TWO112[sx]+x; t = w-TWO112[sx]; + math_force_eval (t); #ifdef USE_FENV_H fesetenv (&env); #endif @@ -72,6 +73,7 @@ nearbyintq(__float128 x) #endif w = TWO112[sx]+x; t = w-TWO112[sx]; + math_force_eval (t); #ifdef USE_FENV_H fesetenv (&env); #endif diff --git a/libquadmath/math/nextafterq.c b/libquadmath/math/nextafterq.c index 04d63deb8e2..a030e9c6444 100644 --- a/libquadmath/math/nextafterq.c +++ b/libquadmath/math/nextafterq.c @@ -13,6 +13,7 @@ * ==================================================== */ +#include <errno.h> #include "quadmath-imp.h" __float128 @@ -54,9 +55,15 @@ nextafterq (__float128 x, __float128 y) } } hy = hx&0x7fff000000000000LL; - if(hy==0x7fff000000000000LL) return x+x;/* overflow */ + if(hy==0x7fff000000000000LL) { + __float128 u = x + x; /* overflow */ + math_force_eval (u); + errno = ERANGE; + } if(hy==0) { - /* here we should raise an underflow flag */ + __float128 u = x*x; /* underflow */ + math_force_eval (u); /* raise underflow flag */ + errno = ERANGE; } SET_FLT128_WORDS64(x,hx,lx); return x; diff --git a/libquadmath/math/powq.c b/libquadmath/math/powq.c index dd44b7c175a..acad2057516 100644 --- a/libquadmath/math/powq.c +++ b/libquadmath/math/powq.c @@ -147,7 +147,7 @@ __float128 powq (__float128 x, __float128 y) { __float128 z, ax, z_h, z_l, p_h, p_l; - __float128 y1, t1, t2, r, s, t, u, v, w; + __float128 y1, t1, t2, r, s, sgn, t, u, v, w; __float128 s2, s_h, s_l, t_h, t_l, ay; int32_t i, j, k, yisint, n; uint32_t ix, iy; @@ -261,6 +261,11 @@ powq (__float128 x, __float128 y) if (((((uint32_t) hx >> 31) - 1) | yisint) == 0) return (x - x) / (x - x); + /* sgn (sign of result -ve**odd) = -1 else = 1 */ + sgn = one; + if (((((u_int32_t) hx >> 31) - 1) | (yisint - 1)) == 0) + sgn = -one; /* (-ve)**(odd int) */ + /* |y| is huge. 2^-16495 = 1/2 of smallest representable value. If (1 - 1/131072)^y underflows, y > 1.4986e9 */ @@ -276,9 +281,9 @@ powq (__float128 x, __float128 y) } /* over/underflow if x is not close to one */ if (ix < 0x3ffeffff) - return (hy < 0) ? huge * huge : tiny * tiny; + return (hy < 0) ? sgn * huge * huge : sgn * tiny * tiny; if (ix > 0x3fff0000) - return (hy > 0) ? huge * huge : tiny * tiny; + return (hy > 0) ? sgn * huge * huge : sgn * tiny * tiny; } ay = y > 0 ? y : -y; @@ -365,11 +370,6 @@ powq (__float128 x, __float128 y) t1 = o.value; t2 = z_l - (((t1 - t) - dp_h[k]) - z_h); - /* s (sign of result -ve**odd) = -1 else = 1 */ - s = one; - if (((((uint32_t) hx >> 31) - 1) | (yisint - 1)) == 0) - s = -one; /* (-ve)**(odd int) */ - /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ y1 = y; o.value = y1; @@ -385,11 +385,11 @@ powq (__float128 x, __float128 y) { /* if z > 16384 */ if (((j - 0x400d0000) | o.words32.w1 | o.words32.w2 | o.words32.w3) != 0) - return s * huge * huge; /* overflow */ + return sgn * huge * huge; /* overflow */ else { if (p_l + ovt > z - p_h) - return s * huge * huge; /* overflow */ + return sgn * huge * huge; /* overflow */ } } else if ((j & 0x7fffffff) >= 0x400d01b9) /* z <= -16495 */ @@ -397,11 +397,11 @@ powq (__float128 x, __float128 y) /* z < -16495 */ if (((j - 0xc00d01bc) | o.words32.w1 | o.words32.w2 | o.words32.w3) != 0) - return s * tiny * tiny; /* underflow */ + return sgn * tiny * tiny; /* underflow */ else { if (p_l <= z - p_h) - return s * tiny * tiny; /* underflow */ + return sgn * tiny * tiny; /* underflow */ } } /* compute 2**(p_h+p_l) */ @@ -434,11 +434,15 @@ powq (__float128 x, __float128 y) j = o.words32.w0; j += (n << 16); if ((j >> 16) <= 0) - z = scalbnq (z, n); /* subnormal output */ + { + z = scalbnq (z, n); /* subnormal output */ + __float128 force_underflow = z * z; + math_force_eval (force_underflow); + } else { o.words32.w0 = j; z = o.value; } - return s * z; + return sgn * z; } diff --git a/libquadmath/math/rem_pio2q.c b/libquadmath/math/rem_pio2q.c index 60653c8d1d3..3308b218473 100644 --- a/libquadmath/math/rem_pio2q.c +++ b/libquadmath/math/rem_pio2q.c @@ -312,7 +312,7 @@ recompute: /* Quad-precision floating point argument reduction. - Copyright (C) 1999 Free Software Foundation, Inc. + Copyright (C) 1999-2017 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jj@ultra.linux.cz> @@ -332,176 +332,176 @@ recompute: 02111-1307 USA. */ /* - * Table of constants for 2/pi, 5628 hexadecimal digits of 2/pi + * Table of constants for 2/pi, 5628 hexadecimal digits of 2/pi */ static const int32_t two_over_pi[] = { -0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62, -0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a, -0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129, -0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41, -0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8, -0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf, -0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5, -0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08, -0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3, -0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880, -0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b, -0x47c419, 0xc367cd, 0xdce809, 0x2a8359, 0xc4768b, 0x961ca6, -0xddaf44, 0xd15719, 0x053ea5, 0xff0705, 0x3f7e33, 0xe832c2, -0xde4f98, 0x327dbb, 0xc33d26, 0xef6b1e, 0x5ef89f, 0x3a1f35, -0xcaf27f, 0x1d87f1, 0x21907c, 0x7c246a, 0xfa6ed5, 0x772d30, -0x433b15, 0xc614b5, 0x9d19c3, 0xc2c4ad, 0x414d2c, 0x5d000c, -0x467d86, 0x2d71e3, 0x9ac69b, 0x006233, 0x7cd2b4, 0x97a7b4, -0xd55537, 0xf63ed7, 0x1810a3, 0xfc764d, 0x2a9d64, 0xabd770, -0xf87c63, 0x57b07a, 0xe71517, 0x5649c0, 0xd9d63b, 0x3884a7, -0xcb2324, 0x778ad6, 0x23545a, 0xb91f00, 0x1b0af1, 0xdfce19, -0xff319f, 0x6a1e66, 0x615799, 0x47fbac, 0xd87f7e, 0xb76522, -0x89e832, 0x60bfe6, 0xcdc4ef, 0x09366c, 0xd43f5d, 0xd7de16, -0xde3b58, 0x929bde, 0x2822d2, 0xe88628, 0x4d58e2, 0x32cac6, -0x16e308, 0xcb7de0, 0x50c017, 0xa71df3, 0x5be018, 0x34132e, -0x621283, 0x014883, 0x5b8ef5, 0x7fb0ad, 0xf2e91e, 0x434a48, -0xd36710, 0xd8ddaa, 0x425fae, 0xce616a, 0xa4280a, 0xb499d3, -0xf2a606, 0x7f775c, 0x83c2a3, 0x883c61, 0x78738a, 0x5a8caf, -0xbdd76f, 0x63a62d, 0xcbbff4, 0xef818d, 0x67c126, 0x45ca55, -0x36d9ca, 0xd2a828, 0x8d61c2, 0x77c912, 0x142604, 0x9b4612, -0xc459c4, 0x44c5c8, 0x91b24d, 0xf31700, 0xad43d4, 0xe54929, -0x10d5fd, 0xfcbe00, 0xcc941e, 0xeece70, 0xf53e13, 0x80f1ec, -0xc3e7b3, 0x28f8c7, 0x940593, 0x3e71c1, 0xb3092e, 0xf3450b, -0x9c1288, 0x7b20ab, 0x9fb52e, 0xc29247, 0x2f327b, 0x6d550c, -0x90a772, 0x1fe76b, 0x96cb31, 0x4a1679, 0xe27941, 0x89dff4, -0x9794e8, 0x84e6e2, 0x973199, 0x6bed88, 0x365f5f, 0x0efdbb, -0xb49a48, 0x6ca467, 0x427271, 0x325d8d, 0xb8159f, 0x09e5bc, -0x25318d, 0x3974f7, 0x1c0530, 0x010c0d, 0x68084b, 0x58ee2c, -0x90aa47, 0x02e774, 0x24d6bd, 0xa67df7, 0x72486e, 0xef169f, -0xa6948e, 0xf691b4, 0x5153d1, 0xf20acf, 0x339820, 0x7e4bf5, -0x6863b2, 0x5f3edd, 0x035d40, 0x7f8985, 0x295255, 0xc06437, -0x10d86d, 0x324832, 0x754c5b, 0xd4714e, 0x6e5445, 0xc1090b, -0x69f52a, 0xd56614, 0x9d0727, 0x50045d, 0xdb3bb4, 0xc576ea, -0x17f987, 0x7d6b49, 0xba271d, 0x296996, 0xacccc6, 0x5414ad, -0x6ae290, 0x89d988, 0x50722c, 0xbea404, 0x940777, 0x7030f3, -0x27fc00, 0xa871ea, 0x49c266, 0x3de064, 0x83dd97, 0x973fa3, -0xfd9443, 0x8c860d, 0xde4131, 0x9d3992, 0x8c70dd, 0xe7b717, -0x3bdf08, 0x2b3715, 0xa0805c, 0x93805a, 0x921110, 0xd8e80f, -0xaf806c, 0x4bffdb, 0x0f9038, 0x761859, 0x15a562, 0xbbcb61, -0xb989c7, 0xbd4010, 0x04f2d2, 0x277549, 0xf6b6eb, 0xbb22db, -0xaa140a, 0x2f2689, 0x768364, 0x333b09, 0x1a940e, 0xaa3a51, -0xc2a31d, 0xaeedaf, 0x12265c, 0x4dc26d, 0x9c7a2d, 0x9756c0, -0x833f03, 0xf6f009, 0x8c402b, 0x99316d, 0x07b439, 0x15200c, -0x5bc3d8, 0xc492f5, 0x4badc6, 0xa5ca4e, 0xcd37a7, 0x36a9e6, -0x9492ab, 0x6842dd, 0xde6319, 0xef8c76, 0x528b68, 0x37dbfc, -0xaba1ae, 0x3115df, 0xa1ae00, 0xdafb0c, 0x664d64, 0xb705ed, -0x306529, 0xbf5657, 0x3aff47, 0xb9f96a, 0xf3be75, 0xdf9328, -0x3080ab, 0xf68c66, 0x15cb04, 0x0622fa, 0x1de4d9, 0xa4b33d, -0x8f1b57, 0x09cd36, 0xe9424e, 0xa4be13, 0xb52333, 0x1aaaf0, -0xa8654f, 0xa5c1d2, 0x0f3f0b, 0xcd785b, 0x76f923, 0x048b7b, -0x721789, 0x53a6c6, 0xe26e6f, 0x00ebef, 0x584a9b, 0xb7dac4, -0xba66aa, 0xcfcf76, 0x1d02d1, 0x2df1b1, 0xc1998c, 0x77adc3, -0xda4886, 0xa05df7, 0xf480c6, 0x2ff0ac, 0x9aecdd, 0xbc5c3f, -0x6dded0, 0x1fc790, 0xb6db2a, 0x3a25a3, 0x9aaf00, 0x9353ad, -0x0457b6, 0xb42d29, 0x7e804b, 0xa707da, 0x0eaa76, 0xa1597b, -0x2a1216, 0x2db7dc, 0xfde5fa, 0xfedb89, 0xfdbe89, 0x6c76e4, -0xfca906, 0x70803e, 0x156e85, 0xff87fd, 0x073e28, 0x336761, -0x86182a, 0xeabd4d, 0xafe7b3, 0x6e6d8f, 0x396795, 0x5bbf31, -0x48d784, 0x16df30, 0x432dc7, 0x356125, 0xce70c9, 0xb8cb30, -0xfd6cbf, 0xa200a4, 0xe46c05, 0xa0dd5a, 0x476f21, 0xd21262, -0x845cb9, 0x496170, 0xe0566b, 0x015299, 0x375550, 0xb7d51e, -0xc4f133, 0x5f6e13, 0xe4305d, 0xa92e85, 0xc3b21d, 0x3632a1, -0xa4b708, 0xd4b1ea, 0x21f716, 0xe4698f, 0x77ff27, 0x80030c, -0x2d408d, 0xa0cd4f, 0x99a520, 0xd3a2b3, 0x0a5d2f, 0x42f9b4, -0xcbda11, 0xd0be7d, 0xc1db9b, 0xbd17ab, 0x81a2ca, 0x5c6a08, -0x17552e, 0x550027, 0xf0147f, 0x8607e1, 0x640b14, 0x8d4196, -0xdebe87, 0x2afdda, 0xb6256b, 0x34897b, 0xfef305, 0x9ebfb9, -0x4f6a68, 0xa82a4a, 0x5ac44f, 0xbcf82d, 0x985ad7, 0x95c7f4, -0x8d4d0d, 0xa63a20, 0x5f57a4, 0xb13f14, 0x953880, 0x0120cc, -0x86dd71, 0xb6dec9, 0xf560bf, 0x11654d, 0x6b0701, 0xacb08c, -0xd0c0b2, 0x485551, 0x0efb1e, 0xc37295, 0x3b06a3, 0x3540c0, -0x7bdc06, 0xcc45e0, 0xfa294e, 0xc8cad6, 0x41f3e8, 0xde647c, -0xd8649b, 0x31bed9, 0xc397a4, 0xd45877, 0xc5e369, 0x13daf0, -0x3c3aba, 0x461846, 0x5f7555, 0xf5bdd2, 0xc6926e, 0x5d2eac, -0xed440e, 0x423e1c, 0x87c461, 0xe9fd29, 0xf3d6e7, 0xca7c22, -0x35916f, 0xc5e008, 0x8dd7ff, 0xe26a6e, 0xc6fdb0, 0xc10893, -0x745d7c, 0xb2ad6b, 0x9d6ecd, 0x7b723e, 0x6a11c6, 0xa9cff7, -0xdf7329, 0xbac9b5, 0x5100b7, 0x0db2e2, 0x24ba74, 0x607de5, -0x8ad874, 0x2c150d, 0x0c1881, 0x94667e, 0x162901, 0x767a9f, -0xbefdfd, 0xef4556, 0x367ed9, 0x13d9ec, 0xb9ba8b, 0xfc97c4, -0x27a831, 0xc36ef1, 0x36c594, 0x56a8d8, 0xb5a8b4, 0x0ecccf, -0x2d8912, 0x34576f, 0x89562c, 0xe3ce99, 0xb920d6, 0xaa5e6b, -0x9c2a3e, 0xcc5f11, 0x4a0bfd, 0xfbf4e1, 0x6d3b8e, 0x2c86e2, -0x84d4e9, 0xa9b4fc, 0xd1eeef, 0xc9352e, 0x61392f, 0x442138, -0xc8d91b, 0x0afc81, 0x6a4afb, 0xd81c2f, 0x84b453, 0x8c994e, -0xcc2254, 0xdc552a, 0xd6c6c0, 0x96190b, 0xb8701a, 0x649569, -0x605a26, 0xee523f, 0x0f117f, 0x11b5f4, 0xf5cbfc, 0x2dbc34, -0xeebc34, 0xcc5de8, 0x605edd, 0x9b8e67, 0xef3392, 0xb817c9, -0x9b5861, 0xbc57e1, 0xc68351, 0x103ed8, 0x4871dd, 0xdd1c2d, -0xa118af, 0x462c21, 0xd7f359, 0x987ad9, 0xc0549e, 0xfa864f, -0xfc0656, 0xae79e5, 0x362289, 0x22ad38, 0xdc9367, 0xaae855, -0x382682, 0x9be7ca, 0xa40d51, 0xb13399, 0x0ed7a9, 0x480569, -0xf0b265, 0xa7887f, 0x974c88, 0x36d1f9, 0xb39221, 0x4a827b, -0x21cf98, 0xdc9f40, 0x5547dc, 0x3a74e1, 0x42eb67, 0xdf9dfe, -0x5fd45e, 0xa4677b, 0x7aacba, 0xa2f655, 0x23882b, 0x55ba41, -0x086e59, 0x862a21, 0x834739, 0xe6e389, 0xd49ee5, 0x40fb49, -0xe956ff, 0xca0f1c, 0x8a59c5, 0x2bfa94, 0xc5c1d3, 0xcfc50f, -0xae5adb, 0x86c547, 0x624385, 0x3b8621, 0x94792c, 0x876110, -0x7b4c2a, 0x1a2c80, 0x12bf43, 0x902688, 0x893c78, 0xe4c4a8, -0x7bdbe5, 0xc23ac4, 0xeaf426, 0x8a67f7, 0xbf920d, 0x2ba365, -0xb1933d, 0x0b7cbd, 0xdc51a4, 0x63dd27, 0xdde169, 0x19949a, -0x9529a8, 0x28ce68, 0xb4ed09, 0x209f44, 0xca984e, 0x638270, -0x237c7e, 0x32b90f, 0x8ef5a7, 0xe75614, 0x08f121, 0x2a9db5, -0x4d7e6f, 0x5119a5, 0xabf9b5, 0xd6df82, 0x61dd96, 0x023616, -0x9f3ac4, 0xa1a283, 0x6ded72, 0x7a8d39, 0xa9b882, 0x5c326b, -0x5b2746, 0xed3400, 0x7700d2, 0x55f4fc, 0x4d5901, 0x8071e0, -0xe13f89, 0xb295f3, 0x64a8f1, 0xaea74b, 0x38fc4c, 0xeab2bb, -0x47270b, 0xabc3a7, 0x34ba60, 0x52dd34, 0xf8563a, 0xeb7e8a, -0x31bb36, 0x5895b7, 0x47f7a9, 0x94c3aa, 0xd39225, 0x1e7f3e, -0xd8974e, 0xbba94f, 0xd8ae01, 0xe661b4, 0x393d8e, 0xa523aa, -0x33068e, 0x1633b5, 0x3bb188, 0x1d3a9d, 0x4013d0, 0xcc1be5, -0xf862e7, 0x3bf28f, 0x39b5bf, 0x0bc235, 0x22747e, 0xa247c0, -0xd52d1f, 0x19add3, 0x9094df, 0x9311d0, 0xb42b25, 0x496db2, -0xe264b2, 0x5ef135, 0x3bc6a4, 0x1a4ad0, 0xaac92e, 0x64e886, -0x573091, 0x982cfb, 0x311b1a, 0x08728b, 0xbdcee1, 0x60e142, -0xeb641d, 0xd0bba3, 0xe559d4, 0x597b8c, 0x2a4483, 0xf332ba, -0xf84867, 0x2c8d1b, 0x2fa9b0, 0x50f3dd, 0xf9f573, 0xdb61b4, -0xfe233e, 0x6c41a6, 0xeea318, 0x775a26, 0xbc5e5c, 0xcea708, -0x94dc57, 0xe20196, 0xf1e839, 0xbe4851, 0x5d2d2f, 0x4e9555, -0xd96ec2, 0xe7d755, 0x6304e0, 0xc02e0e, 0xfc40a0, 0xbbf9b3, -0x7125a7, 0x222dfb, 0xf619d8, 0x838c1c, 0x6619e6, 0xb20d55, -0xbb5137, 0x79e809, 0xaf9149, 0x0d73de, 0x0b0da5, 0xce7f58, -0xac1934, 0x724667, 0x7a1a13, 0x9e26bc, 0x4555e7, 0x585cb5, -0x711d14, 0x486991, 0x480d60, 0x56adab, 0xd62f64, 0x96ee0c, -0x212ff3, 0x5d6d88, 0xa67684, 0x95651e, 0xab9e0a, 0x4ddefe, -0x571010, 0x836a39, 0xf8ea31, 0x9e381d, 0xeac8b1, 0xcac96b, -0x37f21e, 0xd505e9, 0x984743, 0x9fc56c, 0x0331b7, 0x3b8bf8, -0x86e56a, 0x8dc343, 0x6230e7, 0x93cfd5, 0x6a8f2d, 0x733005, -0x1af021, 0xa09fcb, 0x7415a1, 0xd56b23, 0x6ff725, 0x2f4bc7, -0xb8a591, 0x7fac59, 0x5c55de, 0x212c38, 0xb13296, 0x5cff50, -0x366262, 0xfa7b16, 0xf4d9a6, 0x2acfe7, 0xf07403, 0xd4d604, -0x6fd916, 0x31b1bf, 0xcbb450, 0x5bd7c8, 0x0ce194, 0x6bd643, -0x4fd91c, 0xdf4543, 0x5f3453, 0xe2b5aa, 0xc9aec8, 0x131485, -0xf9d2bf, 0xbadb9e, 0x76f5b9, 0xaf15cf, 0xca3182, 0x14b56d, -0xe9fe4d, 0x50fc35, 0xf5aed5, 0xa2d0c1, 0xc96057, 0x192eb6, -0xe91d92, 0x07d144, 0xaea3c6, 0x343566, 0x26d5b4, 0x3161e2, -0x37f1a2, 0x209eff, 0x958e23, 0x493798, 0x35f4a6, 0x4bdc02, -0xc2be13, 0xbe80a0, 0x0b72a3, 0x115c5f, 0x1e1bd1, 0x0db4d3, -0x869e85, 0x96976b, 0x2ac91f, 0x8a26c2, 0x3070f0, 0x041412, -0xfc9fa5, 0xf72a38, 0x9c6878, 0xe2aa76, 0x50cfe1, 0x559274, -0x934e38, 0x0a92f7, 0x5533f0, 0xa63db4, 0x399971, 0xe2b755, -0xa98a7c, 0x008f19, 0xac54d2, 0x2ea0b4, 0xf5f3e0, 0x60c849, -0xffd269, 0xae52ce, 0x7a5fdd, 0xe9ce06, 0xfb0ae8, 0xa50cce, -0xea9d3e, 0x3766dd, 0xb834f5, 0x0da090, 0x846f88, 0x4ae3d5, -0x099a03, 0x2eae2d, 0xfcb40a, 0xfb9b33, 0xe281dd, 0x1b16ba, -0xd8c0af, 0xd96b97, 0xb52dc9, 0x9c277f, 0x5951d5, 0x21ccd6, -0xb6496b, 0x584562, 0xb3baf2, 0xa1a5c4, 0x7ca2cf, 0xa9b93d, -0x7b7b89, 0x483d38, +0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62, +0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a, +0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129, +0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41, +0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8, +0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf, +0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5, +0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08, +0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3, +0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880, +0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b, +0x47c419, 0xc367cd, 0xdce809, 0x2a8359, 0xc4768b, 0x961ca6, +0xddaf44, 0xd15719, 0x053ea5, 0xff0705, 0x3f7e33, 0xe832c2, +0xde4f98, 0x327dbb, 0xc33d26, 0xef6b1e, 0x5ef89f, 0x3a1f35, +0xcaf27f, 0x1d87f1, 0x21907c, 0x7c246a, 0xfa6ed5, 0x772d30, +0x433b15, 0xc614b5, 0x9d19c3, 0xc2c4ad, 0x414d2c, 0x5d000c, +0x467d86, 0x2d71e3, 0x9ac69b, 0x006233, 0x7cd2b4, 0x97a7b4, +0xd55537, 0xf63ed7, 0x1810a3, 0xfc764d, 0x2a9d64, 0xabd770, +0xf87c63, 0x57b07a, 0xe71517, 0x5649c0, 0xd9d63b, 0x3884a7, +0xcb2324, 0x778ad6, 0x23545a, 0xb91f00, 0x1b0af1, 0xdfce19, +0xff319f, 0x6a1e66, 0x615799, 0x47fbac, 0xd87f7e, 0xb76522, +0x89e832, 0x60bfe6, 0xcdc4ef, 0x09366c, 0xd43f5d, 0xd7de16, +0xde3b58, 0x929bde, 0x2822d2, 0xe88628, 0x4d58e2, 0x32cac6, +0x16e308, 0xcb7de0, 0x50c017, 0xa71df3, 0x5be018, 0x34132e, +0x621283, 0x014883, 0x5b8ef5, 0x7fb0ad, 0xf2e91e, 0x434a48, +0xd36710, 0xd8ddaa, 0x425fae, 0xce616a, 0xa4280a, 0xb499d3, +0xf2a606, 0x7f775c, 0x83c2a3, 0x883c61, 0x78738a, 0x5a8caf, +0xbdd76f, 0x63a62d, 0xcbbff4, 0xef818d, 0x67c126, 0x45ca55, +0x36d9ca, 0xd2a828, 0x8d61c2, 0x77c912, 0x142604, 0x9b4612, +0xc459c4, 0x44c5c8, 0x91b24d, 0xf31700, 0xad43d4, 0xe54929, +0x10d5fd, 0xfcbe00, 0xcc941e, 0xeece70, 0xf53e13, 0x80f1ec, +0xc3e7b3, 0x28f8c7, 0x940593, 0x3e71c1, 0xb3092e, 0xf3450b, +0x9c1288, 0x7b20ab, 0x9fb52e, 0xc29247, 0x2f327b, 0x6d550c, +0x90a772, 0x1fe76b, 0x96cb31, 0x4a1679, 0xe27941, 0x89dff4, +0x9794e8, 0x84e6e2, 0x973199, 0x6bed88, 0x365f5f, 0x0efdbb, +0xb49a48, 0x6ca467, 0x427271, 0x325d8d, 0xb8159f, 0x09e5bc, +0x25318d, 0x3974f7, 0x1c0530, 0x010c0d, 0x68084b, 0x58ee2c, +0x90aa47, 0x02e774, 0x24d6bd, 0xa67df7, 0x72486e, 0xef169f, +0xa6948e, 0xf691b4, 0x5153d1, 0xf20acf, 0x339820, 0x7e4bf5, +0x6863b2, 0x5f3edd, 0x035d40, 0x7f8985, 0x295255, 0xc06437, +0x10d86d, 0x324832, 0x754c5b, 0xd4714e, 0x6e5445, 0xc1090b, +0x69f52a, 0xd56614, 0x9d0727, 0x50045d, 0xdb3bb4, 0xc576ea, +0x17f987, 0x7d6b49, 0xba271d, 0x296996, 0xacccc6, 0x5414ad, +0x6ae290, 0x89d988, 0x50722c, 0xbea404, 0x940777, 0x7030f3, +0x27fc00, 0xa871ea, 0x49c266, 0x3de064, 0x83dd97, 0x973fa3, +0xfd9443, 0x8c860d, 0xde4131, 0x9d3992, 0x8c70dd, 0xe7b717, +0x3bdf08, 0x2b3715, 0xa0805c, 0x93805a, 0x921110, 0xd8e80f, +0xaf806c, 0x4bffdb, 0x0f9038, 0x761859, 0x15a562, 0xbbcb61, +0xb989c7, 0xbd4010, 0x04f2d2, 0x277549, 0xf6b6eb, 0xbb22db, +0xaa140a, 0x2f2689, 0x768364, 0x333b09, 0x1a940e, 0xaa3a51, +0xc2a31d, 0xaeedaf, 0x12265c, 0x4dc26d, 0x9c7a2d, 0x9756c0, +0x833f03, 0xf6f009, 0x8c402b, 0x99316d, 0x07b439, 0x15200c, +0x5bc3d8, 0xc492f5, 0x4badc6, 0xa5ca4e, 0xcd37a7, 0x36a9e6, +0x9492ab, 0x6842dd, 0xde6319, 0xef8c76, 0x528b68, 0x37dbfc, +0xaba1ae, 0x3115df, 0xa1ae00, 0xdafb0c, 0x664d64, 0xb705ed, +0x306529, 0xbf5657, 0x3aff47, 0xb9f96a, 0xf3be75, 0xdf9328, +0x3080ab, 0xf68c66, 0x15cb04, 0x0622fa, 0x1de4d9, 0xa4b33d, +0x8f1b57, 0x09cd36, 0xe9424e, 0xa4be13, 0xb52333, 0x1aaaf0, +0xa8654f, 0xa5c1d2, 0x0f3f0b, 0xcd785b, 0x76f923, 0x048b7b, +0x721789, 0x53a6c6, 0xe26e6f, 0x00ebef, 0x584a9b, 0xb7dac4, +0xba66aa, 0xcfcf76, 0x1d02d1, 0x2df1b1, 0xc1998c, 0x77adc3, +0xda4886, 0xa05df7, 0xf480c6, 0x2ff0ac, 0x9aecdd, 0xbc5c3f, +0x6dded0, 0x1fc790, 0xb6db2a, 0x3a25a3, 0x9aaf00, 0x9353ad, +0x0457b6, 0xb42d29, 0x7e804b, 0xa707da, 0x0eaa76, 0xa1597b, +0x2a1216, 0x2db7dc, 0xfde5fa, 0xfedb89, 0xfdbe89, 0x6c76e4, +0xfca906, 0x70803e, 0x156e85, 0xff87fd, 0x073e28, 0x336761, +0x86182a, 0xeabd4d, 0xafe7b3, 0x6e6d8f, 0x396795, 0x5bbf31, +0x48d784, 0x16df30, 0x432dc7, 0x356125, 0xce70c9, 0xb8cb30, +0xfd6cbf, 0xa200a4, 0xe46c05, 0xa0dd5a, 0x476f21, 0xd21262, +0x845cb9, 0x496170, 0xe0566b, 0x015299, 0x375550, 0xb7d51e, +0xc4f133, 0x5f6e13, 0xe4305d, 0xa92e85, 0xc3b21d, 0x3632a1, +0xa4b708, 0xd4b1ea, 0x21f716, 0xe4698f, 0x77ff27, 0x80030c, +0x2d408d, 0xa0cd4f, 0x99a520, 0xd3a2b3, 0x0a5d2f, 0x42f9b4, +0xcbda11, 0xd0be7d, 0xc1db9b, 0xbd17ab, 0x81a2ca, 0x5c6a08, +0x17552e, 0x550027, 0xf0147f, 0x8607e1, 0x640b14, 0x8d4196, +0xdebe87, 0x2afdda, 0xb6256b, 0x34897b, 0xfef305, 0x9ebfb9, +0x4f6a68, 0xa82a4a, 0x5ac44f, 0xbcf82d, 0x985ad7, 0x95c7f4, +0x8d4d0d, 0xa63a20, 0x5f57a4, 0xb13f14, 0x953880, 0x0120cc, +0x86dd71, 0xb6dec9, 0xf560bf, 0x11654d, 0x6b0701, 0xacb08c, +0xd0c0b2, 0x485551, 0x0efb1e, 0xc37295, 0x3b06a3, 0x3540c0, +0x7bdc06, 0xcc45e0, 0xfa294e, 0xc8cad6, 0x41f3e8, 0xde647c, +0xd8649b, 0x31bed9, 0xc397a4, 0xd45877, 0xc5e369, 0x13daf0, +0x3c3aba, 0x461846, 0x5f7555, 0xf5bdd2, 0xc6926e, 0x5d2eac, +0xed440e, 0x423e1c, 0x87c461, 0xe9fd29, 0xf3d6e7, 0xca7c22, +0x35916f, 0xc5e008, 0x8dd7ff, 0xe26a6e, 0xc6fdb0, 0xc10893, +0x745d7c, 0xb2ad6b, 0x9d6ecd, 0x7b723e, 0x6a11c6, 0xa9cff7, +0xdf7329, 0xbac9b5, 0x5100b7, 0x0db2e2, 0x24ba74, 0x607de5, +0x8ad874, 0x2c150d, 0x0c1881, 0x94667e, 0x162901, 0x767a9f, +0xbefdfd, 0xef4556, 0x367ed9, 0x13d9ec, 0xb9ba8b, 0xfc97c4, +0x27a831, 0xc36ef1, 0x36c594, 0x56a8d8, 0xb5a8b4, 0x0ecccf, +0x2d8912, 0x34576f, 0x89562c, 0xe3ce99, 0xb920d6, 0xaa5e6b, +0x9c2a3e, 0xcc5f11, 0x4a0bfd, 0xfbf4e1, 0x6d3b8e, 0x2c86e2, +0x84d4e9, 0xa9b4fc, 0xd1eeef, 0xc9352e, 0x61392f, 0x442138, +0xc8d91b, 0x0afc81, 0x6a4afb, 0xd81c2f, 0x84b453, 0x8c994e, +0xcc2254, 0xdc552a, 0xd6c6c0, 0x96190b, 0xb8701a, 0x649569, +0x605a26, 0xee523f, 0x0f117f, 0x11b5f4, 0xf5cbfc, 0x2dbc34, +0xeebc34, 0xcc5de8, 0x605edd, 0x9b8e67, 0xef3392, 0xb817c9, +0x9b5861, 0xbc57e1, 0xc68351, 0x103ed8, 0x4871dd, 0xdd1c2d, +0xa118af, 0x462c21, 0xd7f359, 0x987ad9, 0xc0549e, 0xfa864f, +0xfc0656, 0xae79e5, 0x362289, 0x22ad38, 0xdc9367, 0xaae855, +0x382682, 0x9be7ca, 0xa40d51, 0xb13399, 0x0ed7a9, 0x480569, +0xf0b265, 0xa7887f, 0x974c88, 0x36d1f9, 0xb39221, 0x4a827b, +0x21cf98, 0xdc9f40, 0x5547dc, 0x3a74e1, 0x42eb67, 0xdf9dfe, +0x5fd45e, 0xa4677b, 0x7aacba, 0xa2f655, 0x23882b, 0x55ba41, +0x086e59, 0x862a21, 0x834739, 0xe6e389, 0xd49ee5, 0x40fb49, +0xe956ff, 0xca0f1c, 0x8a59c5, 0x2bfa94, 0xc5c1d3, 0xcfc50f, +0xae5adb, 0x86c547, 0x624385, 0x3b8621, 0x94792c, 0x876110, +0x7b4c2a, 0x1a2c80, 0x12bf43, 0x902688, 0x893c78, 0xe4c4a8, +0x7bdbe5, 0xc23ac4, 0xeaf426, 0x8a67f7, 0xbf920d, 0x2ba365, +0xb1933d, 0x0b7cbd, 0xdc51a4, 0x63dd27, 0xdde169, 0x19949a, +0x9529a8, 0x28ce68, 0xb4ed09, 0x209f44, 0xca984e, 0x638270, +0x237c7e, 0x32b90f, 0x8ef5a7, 0xe75614, 0x08f121, 0x2a9db5, +0x4d7e6f, 0x5119a5, 0xabf9b5, 0xd6df82, 0x61dd96, 0x023616, +0x9f3ac4, 0xa1a283, 0x6ded72, 0x7a8d39, 0xa9b882, 0x5c326b, +0x5b2746, 0xed3400, 0x7700d2, 0x55f4fc, 0x4d5901, 0x8071e0, +0xe13f89, 0xb295f3, 0x64a8f1, 0xaea74b, 0x38fc4c, 0xeab2bb, +0x47270b, 0xabc3a7, 0x34ba60, 0x52dd34, 0xf8563a, 0xeb7e8a, +0x31bb36, 0x5895b7, 0x47f7a9, 0x94c3aa, 0xd39225, 0x1e7f3e, +0xd8974e, 0xbba94f, 0xd8ae01, 0xe661b4, 0x393d8e, 0xa523aa, +0x33068e, 0x1633b5, 0x3bb188, 0x1d3a9d, 0x4013d0, 0xcc1be5, +0xf862e7, 0x3bf28f, 0x39b5bf, 0x0bc235, 0x22747e, 0xa247c0, +0xd52d1f, 0x19add3, 0x9094df, 0x9311d0, 0xb42b25, 0x496db2, +0xe264b2, 0x5ef135, 0x3bc6a4, 0x1a4ad0, 0xaac92e, 0x64e886, +0x573091, 0x982cfb, 0x311b1a, 0x08728b, 0xbdcee1, 0x60e142, +0xeb641d, 0xd0bba3, 0xe559d4, 0x597b8c, 0x2a4483, 0xf332ba, +0xf84867, 0x2c8d1b, 0x2fa9b0, 0x50f3dd, 0xf9f573, 0xdb61b4, +0xfe233e, 0x6c41a6, 0xeea318, 0x775a26, 0xbc5e5c, 0xcea708, +0x94dc57, 0xe20196, 0xf1e839, 0xbe4851, 0x5d2d2f, 0x4e9555, +0xd96ec2, 0xe7d755, 0x6304e0, 0xc02e0e, 0xfc40a0, 0xbbf9b3, +0x7125a7, 0x222dfb, 0xf619d8, 0x838c1c, 0x6619e6, 0xb20d55, +0xbb5137, 0x79e809, 0xaf9149, 0x0d73de, 0x0b0da5, 0xce7f58, +0xac1934, 0x724667, 0x7a1a13, 0x9e26bc, 0x4555e7, 0x585cb5, +0x711d14, 0x486991, 0x480d60, 0x56adab, 0xd62f64, 0x96ee0c, +0x212ff3, 0x5d6d88, 0xa67684, 0x95651e, 0xab9e0a, 0x4ddefe, +0x571010, 0x836a39, 0xf8ea31, 0x9e381d, 0xeac8b1, 0xcac96b, +0x37f21e, 0xd505e9, 0x984743, 0x9fc56c, 0x0331b7, 0x3b8bf8, +0x86e56a, 0x8dc343, 0x6230e7, 0x93cfd5, 0x6a8f2d, 0x733005, +0x1af021, 0xa09fcb, 0x7415a1, 0xd56b23, 0x6ff725, 0x2f4bc7, +0xb8a591, 0x7fac59, 0x5c55de, 0x212c38, 0xb13296, 0x5cff50, +0x366262, 0xfa7b16, 0xf4d9a6, 0x2acfe7, 0xf07403, 0xd4d604, +0x6fd916, 0x31b1bf, 0xcbb450, 0x5bd7c8, 0x0ce194, 0x6bd643, +0x4fd91c, 0xdf4543, 0x5f3453, 0xe2b5aa, 0xc9aec8, 0x131485, +0xf9d2bf, 0xbadb9e, 0x76f5b9, 0xaf15cf, 0xca3182, 0x14b56d, +0xe9fe4d, 0x50fc35, 0xf5aed5, 0xa2d0c1, 0xc96057, 0x192eb6, +0xe91d92, 0x07d144, 0xaea3c6, 0x343566, 0x26d5b4, 0x3161e2, +0x37f1a2, 0x209eff, 0x958e23, 0x493798, 0x35f4a6, 0x4bdc02, +0xc2be13, 0xbe80a0, 0x0b72a3, 0x115c5f, 0x1e1bd1, 0x0db4d3, +0x869e85, 0x96976b, 0x2ac91f, 0x8a26c2, 0x3070f0, 0x041412, +0xfc9fa5, 0xf72a38, 0x9c6878, 0xe2aa76, 0x50cfe1, 0x559274, +0x934e38, 0x0a92f7, 0x5533f0, 0xa63db4, 0x399971, 0xe2b755, +0xa98a7c, 0x008f19, 0xac54d2, 0x2ea0b4, 0xf5f3e0, 0x60c849, +0xffd269, 0xae52ce, 0x7a5fdd, 0xe9ce06, 0xfb0ae8, 0xa50cce, +0xea9d3e, 0x3766dd, 0xb834f5, 0x0da090, 0x846f88, 0x4ae3d5, +0x099a03, 0x2eae2d, 0xfcb40a, 0xfb9b33, 0xe281dd, 0x1b16ba, +0xd8c0af, 0xd96b97, 0xb52dc9, 0x9c277f, 0x5951d5, 0x21ccd6, +0xb6496b, 0x584562, 0xb3baf2, 0xa1a5c4, 0x7ca2cf, 0xa9b93d, +0x7b7b89, 0x483d38, }; static const __float128 c[] = { -/* 93 bits of pi/2 */ +/* 113 bits of pi/2 */ #define PI_2_1 c[0] - 1.57079632679489661923132169155131424e+00Q, /* 3fff921fb54442d18469898cc5100000 */ + 0x1.921fb54442d18469898cc51701b8p+0Q, /* pi/2 - PI_2_1 */ #define PI_2_1t c[1] - 8.84372056613570112025531863263659260e-29Q, /* 3fa1c06e0e68948127044533e63a0106 */ + 0x3.9a252049c1114cf98e804177d4c8p-116Q, }; @@ -525,8 +525,8 @@ __quadmath_rem_pio2q (__float128 x, __float128 *y) if (ix < 0x40002d97c7f3321dLL) /* |x| in <pi/4, 3pi/4) */ { if (hx > 0) - { - /* 113 + 93 bit PI is ok */ + { + /* 113 + 113 bit PI is ok */ z = x - PI_2_1; y[0] = z - PI_2_1t; y[1] = (z - y[0]) - PI_2_1t; @@ -534,7 +534,7 @@ __quadmath_rem_pio2q (__float128 x, __float128 *y) } else { - /* 113 + 93 bit PI is ok */ + /* 113 + 113 bit PI is ok */ z = x + PI_2_1; y[0] = z + PI_2_1t; y[1] = (z - y[0]) + PI_2_1t; diff --git a/libquadmath/math/remquoq.c b/libquadmath/math/remquoq.c index f7001afc3e5..fa85c5b2d50 100644 --- a/libquadmath/math/remquoq.c +++ b/libquadmath/math/remquoq.c @@ -1,5 +1,5 @@ /* Compute remainder and a congruent to the quotient. - Copyright (C) 1997, 1999, 2002 Free Software Foundation, Inc. + Copyright (C) 1997-2017 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and Jakub Jelinek <jj@ultra.linux.cz>, 1999. @@ -49,7 +49,7 @@ remquoq (__float128 x, __float128 y, int *quo) if (hy <= 0x7ffbffffffffffffLL) x = fmodq (x, 8 * y); /* now x < 8y */ - + if (((hx - hy) | (lx - ly)) == 0) { *quo = qs ? -1 : 1; @@ -60,12 +60,12 @@ remquoq (__float128 x, __float128 y, int *quo) y = fabsq (y); cquo = 0; - if (x >= 4 * y) + if (hy <= 0x7ffcffffffffffffLL && x >= 4 * y) { x -= 4 * y; cquo += 4; } - if (x >= 2 * y) + if (hy <= 0x7ffdffffffffffffLL && x >= 2 * y) { x -= 2 * y; cquo += 2; @@ -101,6 +101,9 @@ remquoq (__float128 x, __float128 y, int *quo) *quo = qs ? -cquo : cquo; + /* Ensure correct sign of zero result in round-downward mode. */ + if (x == 0.0Q) + x = 0.0Q; if (sx) x = -x; return x; diff --git a/libquadmath/math/rintq.c b/libquadmath/math/rintq.c index 15d4c781ba4..7bc9684e5ea 100644 --- a/libquadmath/math/rintq.c +++ b/libquadmath/math/rintq.c @@ -35,7 +35,7 @@ __float128 rintq (__float128 x) { int64_t i0,j0,sx; - uint64_t i1; + uint64_t i1 __attribute__ ((unused)); __float128 w,t; GET_FLT128_WORDS64(i0,i1,x); sx = (((uint64_t)i0)>>63); diff --git a/libquadmath/math/roundq.c b/libquadmath/math/roundq.c index 7c9d640e933..b48366315f4 100644 --- a/libquadmath/math/roundq.c +++ b/libquadmath/math/roundq.c @@ -1,5 +1,5 @@ /* Round __float128 to integer away from zero. - Copyright (C) 1997, 1999 Free Software Foundation, Inc. + Copyright (C) 1997-2017 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and Jakub Jelinek <jj@ultra.linux.cz>, 1999. @@ -21,9 +21,6 @@ #include "quadmath-imp.h" -static const __float128 huge = 1.0E4930Q; - - __float128 roundq (__float128 x) { @@ -32,17 +29,14 @@ roundq (__float128 x) GET_FLT128_WORDS64 (i0, i1, x); j0 = ((i0 >> 48) & 0x7fff) - 0x3fff; - if (j0 < 31) + if (j0 < 48) { if (j0 < 0) { - if (huge + x > 0.0) - { - i0 &= 0x8000000000000000ULL; - if (j0 == -1) - i0 |= 0x3fff000000000000LL; - i1 = 0; - } + i0 &= 0x8000000000000000ULL; + if (j0 == -1) + i0 |= 0x3fff000000000000LL; + i1 = 0; } else { @@ -50,13 +44,9 @@ roundq (__float128 x) if (((i0 & i) | i1) == 0) /* X is integral. */ return x; - if (huge + x > 0.0) - { - /* Raise inexact if x != 0. */ - i0 += 0x0000800000000000LL >> j0; - i0 &= ~i; - i1 = 0; - } + i0 += 0x0000800000000000LL >> j0; + i0 &= ~i; + i1 = 0; } } else if (j0 > 111) @@ -74,14 +64,10 @@ roundq (__float128 x) /* X is integral. */ return x; - if (huge + x > 0.0) - { - /* Raise inexact if x != 0. */ - uint64_t j = i1 + (1LL << (111 - j0)); - if (j < i1) - i0 += 1; - i1 = j; - } + uint64_t j = i1 + (1LL << (111 - j0)); + if (j < i1) + i0 += 1; + i1 = j; i1 &= ~i; } diff --git a/libquadmath/math/scalblnq.c b/libquadmath/math/scalblnq.c index a414045b51f..6aac21125f8 100644 --- a/libquadmath/math/scalblnq.c +++ b/libquadmath/math/scalblnq.c @@ -1,7 +1,7 @@ /* scalblnq.c -- __float128 version of s_scalbn.c. * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz. */ - + /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. diff --git a/libquadmath/math/scalbnq.c b/libquadmath/math/scalbnq.c index 9975a47154c..9ed5fe61799 100644 --- a/libquadmath/math/scalbnq.c +++ b/libquadmath/math/scalbnq.c @@ -1,7 +1,7 @@ /* scalbnq.c -- __float128 version of s_scalbn.c. * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz. */ - + /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. diff --git a/libquadmath/math/sincos_table.c b/libquadmath/math/sincos_table.c index b7e7c750314..362f517a6ef 100644 --- a/libquadmath/math/sincos_table.c +++ b/libquadmath/math/sincos_table.c @@ -1,5 +1,5 @@ /* Quad-precision floating point sine and cosine tables. - Copyright (C) 1999 Free Software Foundation, Inc. + Copyright (C) 1999-2017 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jj@ultra.linux.cz> @@ -37,7 +37,7 @@ const __float128 __sincosq_table[] = { /* sin(x) = 0.25dc50bc95711d0d9787d108fd438cf5959ee0bfb7a1e36e8b1a112968f356657420e9cc9ea */ 1.47892995873409608580026675734609314e-01Q, /* 3ffc2ee285e4ab88e86cbc3e8847ea1c */ 9.74950446464233268291647449768590886e-36Q, /* 3f8a9eb2b3dc17f6f43c6dd16342252d */ - + /* x = 1.56250000000000000000000000000000000e-01 3ffc4000000000000000000000000000 */ /* cos(x) = 0.fce1a053e621438b6d60c76e8c45bf0a9dc71aa16f922acc10e95144ec796a249813c9cb649 */ 9.87817783816471944100503034363211317e-01Q, /* 3ffef9c340a7cc428716dac18edd188b */ diff --git a/libquadmath/math/sincosq.c b/libquadmath/math/sincosq.c index b7c221486d7..d83b1a6b757 100644 --- a/libquadmath/math/sincosq.c +++ b/libquadmath/math/sincosq.c @@ -1,5 +1,5 @@ /* Compute sine and cosine of argument. - Copyright (C) 1997, 1999 Free Software Foundation, Inc. + Copyright (C) 1997-2017 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and Jakub Jelinek <jj@ultra.linux.cz>. @@ -19,6 +19,7 @@ Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. */ +#include <errno.h> #include "quadmath-imp.h" void @@ -37,6 +38,8 @@ sincosq (__float128 x, __float128 *sinx, __float128 *cosx) { /* sin(Inf or NaN) is NaN */ *sinx = *cosx = x - x; + if (isinfq (x)) + errno = EDOM; } else { diff --git a/libquadmath/math/sincosq_kernel.c b/libquadmath/math/sincosq_kernel.c index f6341a4d948..171ed6fc117 100644 --- a/libquadmath/math/sincosq_kernel.c +++ b/libquadmath/math/sincosq_kernel.c @@ -1,5 +1,5 @@ /* Quad-precision floating point sine and cosine on <-pi/4,pi/4>. - Copyright (C) 1999 Free Software Foundation, Inc. + Copyright (C) 1999-2017 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jj@ultra.linux.cz> @@ -110,12 +110,15 @@ __quadmath_kernel_sincosq(__float128 x, __float128 y, __float128 *sinx, /* Argument is small enough to approximate it by a Chebyshev polynomial of degree 16(17). */ if (tix < 0x3fc60000) /* |x| < 2^-57 */ - if (!((int)x)) /* generate inexact */ - { - *sinx = x; - *cosx = ONE; - return; - } + { + math_check_force_underflow (x); + if (!((int)x)) /* generate inexact */ + { + *sinx = x; + *cosx = ONE; + return; + } + } z = x * x; *sinx = x + (x * (z*(SIN1+z*(SIN2+z*(SIN3+z*(SIN4+ z*(SIN5+z*(SIN6+z*(SIN7+z*SIN8))))))))); diff --git a/libquadmath/math/sinhq.c b/libquadmath/math/sinhq.c index eff8149b515..a4191c0fd3d 100644 --- a/libquadmath/math/sinhq.c +++ b/libquadmath/math/sinhq.c @@ -85,8 +85,11 @@ sinhq (__float128 x) if (ix <= 0x40044000) { if (ix < 0x3fc60000) /* |x| < 2^-57 */ - if (shuge + x > one) - return x; /* sinh(tiny) = tiny with inexact */ + { + math_check_force_underflow (x); + if (shuge + x > one) + return x; /* sinh(tiny) = tiny with inexact */ + } t = expm1q (u.value); if (ix < 0x3fff0000) return h * (2.0Q * t - t * t / (t + one)); diff --git a/libquadmath/math/sinq_kernel.c b/libquadmath/math/sinq_kernel.c index 86034551d43..873341a7617 100644 --- a/libquadmath/math/sinq_kernel.c +++ b/libquadmath/math/sinq_kernel.c @@ -90,7 +90,10 @@ __quadmath_kernel_sinq (__float128 x, __float128 y, int iy) /* Argument is small enough to approximate it by a Chebyshev polynomial of degree 17. */ if (tix < 0x3fc60000) /* |x| < 2^-57 */ - if (!((int)x)) return x; /* generate inexact */ + { + math_check_force_underflow (x); + if (!((int)x)) return x; /* generate inexact */ + } z = x * x; return x + (x * (z*(SIN1+z*(SIN2+z*(SIN3+z*(SIN4+ z*(SIN5+z*(SIN6+z*(SIN7+z*SIN8))))))))); diff --git a/libquadmath/math/tanhq.c b/libquadmath/math/tanhq.c index 4ef4fd02143..6b60dbf1b2f 100644 --- a/libquadmath/math/tanhq.c +++ b/libquadmath/math/tanhq.c @@ -72,7 +72,10 @@ tanhq (__float128 x) if (u.value == 0) return x; /* x == +- 0 */ if (ix < 0x3fc60000) /* |x| < 2^-57 */ - return x * (one + tiny); /* tanh(small) = small */ + { + math_check_force_underflow (x); + return x * (one + tiny); /* tanh(small) = small */ + } u.words32.w0 = ix; /* Absolute value of x. */ if (ix >= 0x3fff0000) { /* |x| >= 1 */ diff --git a/libquadmath/math/tanq.c b/libquadmath/math/tanq.c index 690d94b782c..d2864f61c7d 100644 --- a/libquadmath/math/tanq.c +++ b/libquadmath/math/tanq.c @@ -12,9 +12,9 @@ /* Long double expansions are Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov> - and are incorporated herein by permission of the author. The author + and are incorporated herein by permission of the author. The author reserves the right to distribute this material elsewhere under different - copying permissions. These modifications are distributed here under + copying permissions. These modifications are distributed here under the following terms: This library is free software; you can redistribute it and/or @@ -99,8 +99,13 @@ __quadmath_kernel_tanq (__float128 x, __float128 y, int iy) if ((ix | u.words32.w1 | u.words32.w2 | u.words32.w3 | (iy + 1)) == 0) return one / fabsq (x); + else if (iy == 1) + { + math_check_force_underflow (x); + return x; + } else - return (iy == 1) ? x : -one / x; + return -one / x; } } if (ix >= 0x3ffe5942) /* |x| >= 0.6743316650390625 */ @@ -163,7 +168,7 @@ __quadmath_kernel_tanq (__float128 x, __float128 y, int iy) /* tanq.c -- __float128 version of s_tan.c. * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz. */ - + /* @(#)s_tan.c 5.1 93/09/24 */ /* * ==================================================== diff --git a/libquadmath/math/truncq.c b/libquadmath/math/truncq.c index 608871785d1..607c24985e3 100644 --- a/libquadmath/math/truncq.c +++ b/libquadmath/math/truncq.c @@ -1,8 +1,8 @@ /* Truncate argument to nearest integral value not larger than the argument. - Copyright (C) 1997, 1999 Free Software Foundation, Inc. + Copyright (C) 1997-2017 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and - Jakub Jelinek <jj@ultra.linux.cz>, 1999. + Jakub Jelinek <jj@ultra.linux.cz>, 1999. The GNU C Library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/libquadmath/quadmath-imp.h b/libquadmath/quadmath-imp.h index 40b346b6ff2..ed186144e75 100644 --- a/libquadmath/quadmath-imp.h +++ b/libquadmath/quadmath-imp.h @@ -186,4 +186,45 @@ do { \ __builtin_fpclassify (QUADFP_NAN, QUADFP_INFINITE, QUADFP_NORMAL, \ QUADFP_SUBNORMAL, QUADFP_ZERO, x) +#ifndef math_opt_barrier +# define math_opt_barrier(x) \ +({ __typeof (x) __x = (x); __asm ("" : "+m" (__x)); __x; }) +# define math_force_eval(x) \ +({ __typeof (x) __x = (x); __asm __volatile__ ("" : : "m" (__x)); }) +#endif + +/* math_narrow_eval reduces its floating-point argument to the range + and precision of its semantic type. (The original evaluation may + still occur with excess range and precision, so the result may be + affected by double rounding.) */ +#define math_narrow_eval(x) (x) + +/* If X (which is not a NaN) is subnormal, force an underflow + exception. */ +#define math_check_force_underflow(x) \ + do \ + { \ + __float128 force_underflow_tmp = (x); \ + if (fabsq (force_underflow_tmp) < FLT128_MIN) \ + { \ + __float128 force_underflow_tmp2 \ + = force_underflow_tmp * force_underflow_tmp; \ + math_force_eval (force_underflow_tmp2); \ + } \ + } \ + while (0) +/* Likewise, but X is also known to be nonnegative. */ +#define math_check_force_underflow_nonneg(x) \ + do \ + { \ + __float128 force_underflow_tmp = (x); \ + if (force_underflow_tmp < FLT128_MIN) \ + { \ + __float128 force_underflow_tmp2 \ + = force_underflow_tmp * force_underflow_tmp; \ + math_force_eval (force_underflow_tmp2); \ + } \ + } \ + while (0) + #endif |