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

libquadmath/ChangeLog

+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.

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 */
 	{

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 */
     }

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
 	{

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;

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;

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;

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);
     }

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;
     }

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);

libquadmath/math/expq.c

 /* 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/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);
 }

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;

libquadmath/math/fmaq.c

 /* 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;
     }
 }

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);

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;
 }

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;
 }

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;
 }

libquadmath/math/llrintq.c

 /* 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;

libquadmath/math/llroundq.c

 /* 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/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

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));
     }

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

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;

libquadmath/math/lrintq.c

 /* 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/lroundq.c

 /* 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/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

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;

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;
 }

libquadmath/math/remquoq.c

 /* 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;

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);

libquadmath/math/roundq.c

 /* 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/scalblnq.c

 /* 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/scalbnq.c

 /* 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/sincos_table.c

 /* 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/sincosq.c

 /* 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
     {

libquadmath/math/sincosq_kernel.c

 /* 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)))))))));

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));

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)))))))));

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  */

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 */
 /*
  * ====================================================

libquadmath/math/truncq.c

 /* 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/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