This is the mail archive of the
gcc-patches@gcc.gnu.org
mailing list for the GCC project.
[PATCH] Misc libquadmath backports from upstram glibc (PR libquadmath/65757)
- From: Jakub Jelinek <jakub at redhat dot com>
- To: Tobias Burnus <tobias dot burnus at physik dot fu-berlin dot de>, "Joseph S. Myers" <joseph at codesourcery dot com>, "Carlos O'Donell" <carlos at redhat dot com>, Florian Weimer <fweimer at redhat dot com>
- Cc: gcc-patches at gcc dot gnu dot org
- Date: Thu, 13 Jul 2017 22:47:06 +0200
- Subject: [PATCH] Misc libquadmath backports from upstram glibc (PR libquadmath/65757)
- Authentication-results: sourceware.org; auth=none
- Authentication-results: ext-mx04.extmail.prod.ext.phx2.redhat.com; dmarc=none (p=none dis=none) header.from=redhat.com
- Authentication-results: ext-mx04.extmail.prod.ext.phx2.redhat.com; spf=pass smtp.mailfrom=jakub at redhat dot com
- Dkim-filter: OpenDKIM Filter v2.11.0 mx1.redhat.com AFAF98B13D
- Dmarc-filter: OpenDMARC Filter v1.3.2 mx1.redhat.com AFAF98B13D
- Reply-to: Jakub Jelinek <jakub at redhat dot com>
Hi!
This patch is a manual backport of the 2012-2017 sysdeps/ieee754/ldbl-128/
glibc changes into libquadmath. As mentioned in the PR, which has
detailed git diff commands, I've left *jnl.c, *lgamma* and x2y2m1l.c
changes so far, those were too large, and the long double -> _Float128
and 123.456L -> L(123.456) changes.
Bootstrapped/regtested on x86_64-linux and i686-linux, no further testing
done though. I'll commit this next week unless I hear objections.
Of course, further testing e.g. using glibc testsuite would be greatly
appreciated.
2017-07-13 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.
--- libquadmath/quadmath-imp.h.jj 2013-02-07 08:59:55.000000000 +0100
+++ libquadmath/quadmath-imp.h 2017-07-13 13:04:21.928641461 +0200
@@ -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
--- libquadmath/math/ceilq.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/ceilq.c 2017-07-13 13:04:21.932641408 +0200
@@ -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;
--- libquadmath/math/remquoq.c.jj 2013-02-07 08:59:54.000000000 +0100
+++ libquadmath/math/remquoq.c 2017-07-13 13:04:21.936641356 +0200
@@ -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
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
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 = 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;
--- libquadmath/math/expq.c.jj 2013-02-07 08:59:54.000000000 +0100
+++ libquadmath/math/expq.c 2017-07-13 13:04:21.939641316 +0200
@@ -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))
--- libquadmath/math/llroundq.c.jj 2013-02-07 08:59:54.000000000 +0100
+++ libquadmath/math/llroundq.c 2017-07-13 13:04:21.944641251 +0200
@@ -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;
}
--- libquadmath/math/logq.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/logq.c 2017-07-13 13:04:21.948641198 +0200
@@ -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;
--- libquadmath/math/atanq.c.jj 2012-11-22 11:04:24.000000000 +0100
+++ libquadmath/math/atanq.c 2017-07-13 13:04:21.952641145 +0200
@@ -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;
--- libquadmath/math/nearbyintq.c.jj 2012-11-01 09:33:30.000000000 +0100
+++ libquadmath/math/nearbyintq.c 2017-07-13 13:06:03.533307283 +0200
@@ -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
--- libquadmath/math/scalblnq.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/scalblnq.c 2017-07-13 13:04:21.960641040 +0200
@@ -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.
--- libquadmath/math/finiteq.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/finiteq.c 2017-07-13 13:04:21.964640988 +0200
@@ -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);
}
--- libquadmath/math/atanhq.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/atanhq.c 2017-07-13 13:04:21.968640935 +0200
@@ -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;
--- libquadmath/math/expm1q.c.jj 2012-11-01 09:33:30.000000000 +0100
+++ libquadmath/math/expm1q.c 2017-07-13 13:04:21.972640883 +0200
@@ -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);
--- libquadmath/math/sinhq.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/sinhq.c 2017-07-13 13:04:21.976640830 +0200
@@ -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));
--- libquadmath/math/log10q.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/log10q.c 2017-07-13 13:55:14.279634475 +0200
@@ -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
--- libquadmath/math/rintq.c.jj 2012-12-13 13:33:09.000000000 +0100
+++ libquadmath/math/rintq.c 2017-07-13 13:04:21.984640725 +0200
@@ -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);
--- libquadmath/math/roundq.c.jj 2013-02-07 08:59:54.000000000 +0100
+++ libquadmath/math/roundq.c 2017-07-13 13:04:21.988640673 +0200
@@ -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;
}
--- libquadmath/math/fmaq.c.jj 2013-02-07 08:59:54.000000000 +0100
+++ libquadmath/math/fmaq.c 2017-07-13 13:04:21.991640633 +0200
@@ -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, __floa
&& 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, __floa
<= 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, __floa
/* 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, __floa
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, __floa
/* 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, __floa
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, __floa
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;
}
}
--- libquadmath/math/erfq.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/erfq.c 2017-07-13 13:04:21.995640581 +0200
@@ -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
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;
}
--- libquadmath/math/log2q.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/log2q.c 2017-07-13 13:55:56.478098300 +0200
@@ -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
--- libquadmath/math/lroundq.c.jj 2013-02-07 08:59:54.000000000 +0100
+++ libquadmath/math/lroundq.c 2017-07-13 13:04:22.003640476 +0200
@@ -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;
}
--- libquadmath/math/j1q.c.jj 2012-11-22 11:04:24.000000000 +0100
+++ libquadmath/math/j1q.c 2017-07-13 13:53:36.232880261 +0200
@@ -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;
}
--- libquadmath/math/scalbnq.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/scalbnq.c 2017-07-13 13:04:22.010640384 +0200
@@ -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.
--- libquadmath/math/truncq.c.jj 2013-02-07 08:59:54.000000000 +0100
+++ libquadmath/math/truncq.c 2017-07-13 13:04:22.014640331 +0200
@@ -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
--- libquadmath/math/frexpq.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/frexpq.c 2017-07-13 13:04:22.017640292 +0200
@@ -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);
--- libquadmath/math/sincosq.c.jj 2013-02-07 08:59:54.000000000 +0100
+++ libquadmath/math/sincosq.c 2017-07-13 13:04:22.021640239 +0200
@@ -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,
{
/* sin(Inf or NaN) is NaN */
*sinx = *cosx = x - x;
+ if (isinfq (x))
+ errno = EDOM;
}
else
{
--- libquadmath/math/tanhq.c.jj 2010-11-19 20:56:37.000000000 +0100
+++ libquadmath/math/tanhq.c 2017-07-13 13:04:22.025640187 +0200
@@ -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 */
--- libquadmath/math/asinq.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/asinq.c 2017-07-13 13:33:01.676595863 +0200
@@ -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
{
--- libquadmath/math/coshq.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/coshq.c 2017-07-13 13:04:22.032640095 +0200
@@ -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);
}
--- libquadmath/math/j0q.c.jj 2012-11-22 11:04:24.000000000 +0100
+++ libquadmath/math/j0q.c 2017-07-13 13:53:27.655989239 +0200
@@ -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;
}
--- libquadmath/math/asinhq.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/asinhq.c 2017-07-13 13:04:22.039640003 +0200
@@ -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 */
}
--- libquadmath/math/floorq.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/floorq.c 2017-07-13 13:04:22.043639951 +0200
@@ -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;
--- libquadmath/math/sinq_kernel.c.jj 2013-02-07 08:59:54.000000000 +0100
+++ libquadmath/math/sinq_kernel.c 2017-07-13 13:04:22.047639898 +0200
@@ -90,7 +90,10 @@ __quadmath_kernel_sinq (__float128 x, __
/* 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)))))))));
--- libquadmath/math/powq.c.jj 2012-11-22 11:04:24.000000000 +0100
+++ libquadmath/math/powq.c 2017-07-13 13:04:22.051639845 +0200
@@ -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;
}
--- libquadmath/math/hypotq.c.jj 2011-01-18 12:18:53.000000000 +0100
+++ libquadmath/math/hypotq.c 2017-07-13 13:04:22.055639793 +0200
@@ -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;
}
--- libquadmath/math/sincos_table.c.jj 2013-02-07 08:59:54.000000000 +0100
+++ libquadmath/math/sincos_table.c 2017-07-13 13:04:22.058639754 +0200
@@ -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 */
--- libquadmath/math/rem_pio2q.c.jj 2013-02-07 08:59:54.000000000 +0100
+++ libquadmath/math/rem_pio2q.c 2017-07-13 13:56:40.754535721 +0200
@@ -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, __fl
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, __fl
}
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;
--- libquadmath/math/nextafterq.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/nextafterq.c 2017-07-13 13:04:22.066639649 +0200
@@ -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;
--- libquadmath/math/log1pq.c.jj 2012-11-22 11:04:24.000000000 +0100
+++ libquadmath/math/log1pq.c 2017-07-13 13:04:22.069639609 +0200
@@ -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));
}
--- libquadmath/math/sincosq_kernel.c.jj 2013-02-07 08:59:54.000000000 +0100
+++ libquadmath/math/sincosq_kernel.c 2017-07-13 13:04:22.072639570 +0200
@@ -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,
/* 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)))))))));
--- libquadmath/math/tanq.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/tanq.c 2017-07-13 13:04:22.076639517 +0200
@@ -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, __
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, __
/* 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 */
/*
* ====================================================
--- libquadmath/math/acosq.c.jj 2012-11-02 09:01:48.000000000 +0100
+++ libquadmath/math/acosq.c 2017-07-13 13:04:22.080639465 +0200
@@ -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 */
{
--- libquadmath/math/lrintq.c.jj 2013-02-07 08:59:54.000000000 +0100
+++ libquadmath/math/lrintq.c 2017-07-13 13:04:22.083639425 +0200
@@ -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;
}
--- libquadmath/math/llrintq.c.jj 2013-02-07 08:59:54.000000000 +0100
+++ libquadmath/math/llrintq.c 2017-07-13 13:34:20.858586606 +0200
@@ -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;
Jakub