2012-10-31 Tobias Burnus <burnus@net-b.de>

            Joseph Myers <joseph@codesourcery.com>
            David S. Miller <davem@davemloft.net>
            Ulrich Drepper <drepper@redhat.com>
            Marek Polacek <polacek@redhat.com>:
            Petr Baudis <pasky@suse.cz>

        * math/complex.c (csqrtq): NaN and INF fixes.
        * math/sqrtq.c (sqrt): NaN, INF and < 0 fixes.
        * math/expm1q.c (expm1q): Changes from GLIBC. Use expq for
        large parameters. Fix errno for boundary conditions.
        * math/finiteq.c (finiteq): Add comment.
        * math/fmaq.c (fmaq): Changes from GLIBC. Fix missing underflows
        and bad results for some subnormal results. Fix sign of inexact
        zero return. Fix sign of exact zero return.
        Ensure additions are not scheduled after fetestexcept.
        * math/jnq.c (jnq): Changes from GLIBC. Set up errno properly
        for ynq. Fix jnq precision.
        * math/nearbyintq.c (nearbyintq): Changes from GLIBC. Do not
        manipulate bits before adding and subtracting TWO112[sx].
        * math/rintq.c (rintq): Ditto.
        * math/scalbnq.c (scalbnq): Changes from GLIBC. Fix integer
        overflow.



git-svn-id: svn+ssh://gcc.gnu.org/svn/gcc/trunk@193037 138bc75d-0d04-0410-961f-82ee72b054a4

11 files changed, 131 insertions, 74 deletions

diff --git a/libquadmath/ChangeLog b/libquadmath/ChangeLog
index dd37cd348a4..37cc5a25028 100644
--- a/libquadmath/ChangeLog
+++ b/libquadmath/ChangeLog
@@ -1,3 +1,27 @@
+2012-10-31  Tobias Burnus  <burnus@net-b.de>
+	    Joseph Myers <joseph@codesourcery.com>
+	    David S. Miller <davem@davemloft.net>
+	    Ulrich Drepper <drepper@redhat.com>
+	    Marek Polacek <polacek@redhat.com>:
+	    Petr Baudis <pasky@suse.cz>
+
+	* math/complex.c (csqrtq): NaN and INF fixes. 
+	* math/sqrtq.c (sqrt): NaN, INF and < 0 fixes.
+	* math/expm1q.c (expm1q): Changes from GLIBC. Use expq for
+	large parameters. Fix errno for boundary conditions.
+	* math/finiteq.c (finiteq): Add comment.
+	* math/fmaq.c (fmaq): Changes from GLIBC. Fix missing underflows
+	and bad results for some subnormal results. Fix sign of inexact
+	zero return. Fix sign of exact zero return. 
+	Ensure additions are not scheduled after fetestexcept.
+	* math/jnq.c (jnq): Changes from GLIBC. Set up errno properly
+	for ynq. Fix jnq precision.
+	* math/nearbyintq.c (nearbyintq): Changes from GLIBC. Do not
+	manipulate bits before adding and subtracting TWO112[sx].
+	* math/rintq.c (rintq): Ditto.
+	* math/scalbnq.c (scalbnq): Changes from GLIBC. Fix integer
+	overflow.
+
 2012-09-14  David Edelsohn  <dje.gcc@gmail.com>
 
 	* configure: Regenerated.
diff --git a/libquadmath/math/complex.c b/libquadmath/math/complex.c
index f67448a2c12..9953f52ca6b 100644
--- a/libquadmath/math/complex.c
+++ b/libquadmath/math/complex.c
@@ -177,7 +177,11 @@ csqrtq (__complex128 z)
 
   if (im == 0)
   {
-    if (re < 0)
+    if (isnanq (re))
+    {
+      COMPLEX_ASSIGN (v, -re, -re);
+    }
+    else if (re < 0)
     {
       COMPLEX_ASSIGN (v, 0, copysignq (sqrtq (-re), im));
     }
@@ -186,6 +190,10 @@ csqrtq (__complex128 z)
       COMPLEX_ASSIGN (v, fabsq (sqrtq (re)), copysignq (0, im));
     }
   }
+  else if (isinfq (im))
+  {
+    COMPLEX_ASSIGN (v, fabsq (im), im);
+  }
   else if (re == 0)
   {
     __float128 r = sqrtq (0.5 * fabsq (im));
diff --git a/libquadmath/math/expm1q.c b/libquadmath/math/expm1q.c
index 510c98fe493..8cfdd8eec94 100644
--- a/libquadmath/math/expm1q.c
+++ b/libquadmath/math/expm1q.c
@@ -53,6 +53,7 @@
 
 
 
+#include <errno.h>
 #include "quadmath-imp.h"
 
 /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x)
@@ -100,6 +101,11 @@ expm1q (__float128 x)
   ix = u.words32.w0;
   sign = ix & 0x80000000;
   ix &= 0x7fffffff;
+  if (!sign && ix >= 0x40060000)
+    {
+      /* If num is positive and exp >= 6 use plain exp.  */
+      return expq (x);
+    }
   if (ix >= 0x7fff0000)
     {
       /* Infinity. */
@@ -120,7 +126,10 @@ expm1q (__float128 x)
 
   /* Overflow.  */
   if (x > maxlog)
-    return (HUGE_VALQ * HUGE_VALQ);
+    {
+      errno = ERANGE;
+      return (HUGE_VALQ * HUGE_VALQ);
+    }
 
   /* Minimum value.  */
   if (x < minarg)
diff --git a/libquadmath/math/finiteq.c b/libquadmath/math/finiteq.c
index f22e9d7f20e..663c9289263 100644
--- a/libquadmath/math/finiteq.c
+++ b/libquadmath/math/finiteq.c
@@ -15,6 +15,11 @@
 
 #include "quadmath-imp.h"
 
+/*
+ * finiteq(x) returns 1 is x is finite, else 0;
+ * no branching!
+ */
+
 int
 finiteq (const __float128 x)
 {
diff --git a/libquadmath/math/fmaq.c b/libquadmath/math/fmaq.c
index 126b0a2d26b..23e3188669e 100644
--- a/libquadmath/math/fmaq.c
+++ b/libquadmath/math/fmaq.c
@@ -1,5 +1,5 @@
 /* Compute x * y + z as ternary operation.
-   Copyright (C) 2010 Free Software Foundation, Inc.
+   Copyright (C) 2010-2012 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Jakub Jelinek <jakub@redhat.com>, 2010.
 
@@ -57,6 +57,11 @@ fmaq (__float128 x, __float128 y, __float128 z)
 	  && u.ieee.exponent != 0x7fff
           && v.ieee.exponent != 0x7fff)
 	return (z + x) + y;
+      /* If z is zero and x are y are nonzero, compute the result
+	 as x * y to avoid the wrong sign of a zero result if x * y
+	 underflows to 0.  */
+      if (z == 0 && x != 0 && y != 0)
+	return x * y;
       /* If x or y or z is Inf/NaN, or if fma will certainly overflow,
 	 or if x * y is less than half of FLT128_DENORM_MIN,
 	 compute as x * y + z.  */
@@ -136,6 +141,11 @@ fmaq (__float128 x, __float128 y, __float128 z)
       y = v.value;
       z = w.value;
     }
+
+  /* Ensure correct sign of exact 0 + 0.  */
+  if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0))
+    return x * y + z;
+
   /* Multiplication m1 + m2 = x * y using Dekker's algorithm.  */
 #define C ((1LL << (FLT128_MANT_DIG + 1) / 2) + 1)
   __float128 x1 = x * C;
@@ -191,7 +201,7 @@ fmaq (__float128 x, __float128 y, __float128 z)
 #endif
       v.value = a1 + u.value;
       /* Ensure the addition is not scheduled after fetestexcept call.  */
-      asm volatile ("" : : "m" (v));
+      asm volatile ("" : : "m" (v.value));
 #ifdef USE_FENV_H
       int j = fetestexcept (FE_INEXACT) != 0;
       feupdateenv (&env);
@@ -220,20 +230,14 @@ fmaq (__float128 x, __float128 y, __float128 z)
 	{
 	  /* 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
-	     bit.  In round-to-nearest 001 rounds down like 00,
-	     011 rounds up, even though 01 rounds down (thus we need
-	     to adjust), 101 rounds down like 10 and 111 rounds up
-	     like 11.  */
-	  if ((v.ieee.mant_low & 3) == 1)
-	    {
-	      v.value *= 0x1p-226Q;
-	      if (v.ieee.negative)
-		return v.value - 0x1p-16494Q /* __FLT128_DENORM_MIN__ */;
-	      else
-		return v.value + 0x1p-16494Q /* __FLT128_DENORM_MIN__ */;
-	    }
-	  else
-	    return v.value * 0x1p-226Q;
+	     bit. */
+	  w.value = 0.0Q;
+	  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-226L;
+	  w.value *= 0x1p-2L;
+	  return v.value + w.value;
 	}
       v.ieee.mant_low |= j;
       return v.value * 0x1p-226Q;
diff --git a/libquadmath/math/jnq.c b/libquadmath/math/jnq.c
index d82947a3cca..56a183604c1 100644
--- a/libquadmath/math/jnq.c
+++ b/libquadmath/math/jnq.c
@@ -11,9 +11,9 @@
 
 /* Modifications 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
@@ -56,6 +56,7 @@
  *
  */
 
+#include <errno.h>
 #include "quadmath-imp.h"
 
 static const __float128
@@ -273,7 +274,16 @@ jnq (int n, __float128 x)
 		    }
 		}
 	    }
-	  b = (t * j0q (x) / b);
+	  /* j0() and j1() suffer enormous loss of precision at and
+	   * near zero; however, we know that their zero points never
+	   * coincide, so just choose the one further away from zero.
+	   */
+	  z = j0q (x);
+	  w = j1q (x);
+	  if (fabsq (z) >= fabsq (w))
+	    b = (t * z / b);
+	  else
+	    b = (t * w / a);
 	}
     }
   if (sgn == 1)
@@ -374,6 +384,9 @@ ynq (int n, __float128 x)
 	  a = temp;
 	}
     }
+  /* If B is +-Inf, set up errno accordingly.  */
+  if (! finiteq (b))
+    errno = ERANGE;
   if (sign > 0)
     return b;
   else
diff --git a/libquadmath/math/nearbyintq.c b/libquadmath/math/nearbyintq.c
index 8e92c5afd4b..207124808e7 100644
--- a/libquadmath/math/nearbyintq.c
+++ b/libquadmath/math/nearbyintq.c
@@ -44,18 +44,13 @@ nearbyintq(__float128 x)
 	fenv_t env;
 #endif
 	int64_t i0,j0,sx;
-	uint64_t i,i1;
+	uint64_t i1;
 	__float128 w,t;
 	GET_FLT128_WORDS64(i0,i1,x);
 	sx = (((uint64_t)i0)>>63);
 	j0 = ((i0>>48)&0x7fff)-0x3fff;
-	if(j0<48) {
+	if(j0<112) {
 	    if(j0<0) {
-		if(((i0&0x7fffffffffffffffLL)|i1)==0) return x;
-		i1 |= (i0&0x0000ffffffffffffLL);
-		i0 &= 0xffffe00000000000ULL;
-		i0 |= ((i1|-i1)>>16)&0x0000800000000000LL;
-		SET_FLT128_MSW64(x,i0);
 #ifdef USE_FENV_H
 		feholdexcept (&env);
 #endif
@@ -67,25 +62,11 @@ nearbyintq(__float128 x)
 		GET_FLT128_MSW64(i0,t);
 		SET_FLT128_MSW64(t,(i0&0x7fffffffffffffffLL)|(sx<<63));
 	        return t;
-	    } else {
-		i = (0x0000ffffffffffffLL)>>j0;
-		if(((i0&i)|i1)==0) return x; /* x is integral */
-		i>>=1;
-		if(((i0&i)|i1)!=0) {
-		    if(j0==47) i1 = 0x4000000000000000ULL; else
-		    i0 = (i0&(~i))|((0x0000200000000000LL)>>j0);
-		}
 	    }
-	} else if (j0>111) {
+	} else {
 	    if(j0==0x4000) return x+x;	/* inf or NaN */
 	    else return x;		/* x is integral */
-	} else {
-	    i = -1ULL>>(j0-48);
-	    if((i1&i)==0) return x;	/* x is integral */
-	    i>>=1;
-	    if((i1&i)!=0) i1 = (i1&(~i))|((0x4000000000000000LL)>>(j0-48));
 	}
-	SET_FLT128_WORDS64(x,i0,i1);
 #ifdef USE_FENV_H
 	feholdexcept (&env);
 #endif
diff --git a/libquadmath/math/rintq.c b/libquadmath/math/rintq.c
index fe7a59188b4..8a93fdbfa78 100644
--- a/libquadmath/math/rintq.c
+++ b/libquadmath/math/rintq.c
@@ -13,6 +13,16 @@
  * ====================================================
  */
 
+/*
+ * rintq(x)
+ * Return x rounded to integral value according to the prevailing
+ * rounding mode.
+ * Method:
+ *	Using floating addition.
+ * Exception:
+ *	Inexact flag raised if x not equal to rintq(x).
+ */
+
 #include "quadmath-imp.h"
 
 static const __float128
@@ -25,42 +35,23 @@ __float128
 rintq (__float128 x)
 {
 	int64_t i0,j0,sx;
-	uint64_t i,i1;
+	uint64_t i1;
 	__float128 w,t;
 	GET_FLT128_WORDS64(i0,i1,x);
 	sx = (((uint64_t)i0)>>63);
 	j0 = ((i0>>48)&0x7fff)-0x3fff;
-	if(j0<48) {
+	if(j0<112) {
 	    if(j0<0) {
-		if(((i0&0x7fffffffffffffffLL)|i1)==0) return x;
-		i1 |= (i0&0x0000ffffffffffffLL);
-		i0 &= 0xffffe00000000000ULL;
-		i0 |= ((i1|-i1)>>16)&0x0000800000000000LL;
-		SET_FLT128_MSW64(x,i0);
 	        w = TWO112[sx]+x;
 	        t = w-TWO112[sx];
 		GET_FLT128_MSW64(i0,t);
 		SET_FLT128_MSW64(t,(i0&0x7fffffffffffffffLL)|(sx<<63));
 	        return t;
-	    } else {
-		i = (0x0000ffffffffffffLL)>>j0;
-		if(((i0&i)|i1)==0) return x; /* x is integral */
-		i>>=1;
-		if(((i0&i)|i1)!=0) {
-		    if(j0==47) i1 = 0x4000000000000000ULL; else
-		    i0 = (i0&(~i))|((0x0000200000000000LL)>>j0);
-		}
 	    }
-	} else if (j0>111) {
+	} else {
 	    if(j0==0x4000) return x+x;	/* inf or NaN */
 	    else return x;		/* x is integral */
-	} else {
-	    i = -1ULL>>(j0-48);
-	    if((i1&i)==0) return x;	/* x is integral */
-	    i>>=1;
-	    if((i1&i)!=0) i1 = (i1&(~i))|((0x4000000000000000LL)>>(j0-48));
 	}
-	SET_FLT128_WORDS64(x,i0,i1);
 	w = TWO112[sx]+x;
 	return w-TWO112[sx];
 }
diff --git a/libquadmath/math/scalblnq.c b/libquadmath/math/scalblnq.c
index 75997f688c4..b2332250c61 100644
--- a/libquadmath/math/scalblnq.c
+++ b/libquadmath/math/scalblnq.c
@@ -13,6 +13,13 @@
  * ====================================================
  */
 
+/*
+ * scalblnq (_float128 x, long int n)
+ * scalblnq(x,n) returns x* 2**n  computed by  exponent
+ * manipulation rather than by actually performing an
+ * exponentiation or a multiplication.
+ */
+
 #include "quadmath-imp.h"
 
 static const __float128
@@ -34,10 +41,12 @@ scalblnq (__float128 x, long int n)
 	    k = ((hx>>48)&0x7fff) - 114;
 	}
         if (k==0x7fff) return x+x;		/* NaN or Inf */
-        k = k+n;
-        if (n> 50000 || k > 0x7ffe)
-	  return huge*copysignq(huge,x); /* overflow  */
 	if (n< -50000) return tiny*copysignq(tiny,x); /*underflow*/
+        if (n> 50000 || k+n > 0x7ffe)
+	  return huge*copysignq(huge,x); /* overflow  */
+	/* Now k and n are bounded we know that k = k+n does not
+	   overflow.  */
+        k = k+n;
         if (k > 0) 				/* normal result */
 	    {SET_FLT128_MSW64(x,(hx&0x8000ffffffffffffULL)|(k<<48)); return x;}
         if (k <= -114)
diff --git a/libquadmath/math/scalbnq.c b/libquadmath/math/scalbnq.c
index b7049a70e45..f0852ee038d 100644
--- a/libquadmath/math/scalbnq.c
+++ b/libquadmath/math/scalbnq.c
@@ -13,6 +13,14 @@
  * ====================================================
  */
 
+
+/*
+ * scalbnq (__float128 x, int n)
+ * scalbnq(x,n) returns x* 2**n  computed by  exponent
+ * manipulation rather than by actually performing an
+ * exponentiation or a multiplication.
+ */
+
 #include "quadmath-imp.h"
 
 static const __float128
@@ -34,10 +42,12 @@ scalbnq (__float128 x, int n)
 	    k = ((hx>>48)&0x7fff) - 114;
 	}
         if (k==0x7fff) return x+x;		/* NaN or Inf */
-        k = k+n;
-        if (n> 50000 || k > 0x7ffe)
-	  return huge*copysignq(huge,x); /* overflow  */
 	if (n< -50000) return tiny*copysignq(tiny,x); /*underflow*/
+        if (n> 50000 || k+n > 0x7ffe)
+	  return huge*copysignq(huge,x); /* overflow  */
+	/* Now k and n are bounded we know that k = k+n does not
+	   overflow.  */
+        k = k+n;
         if (k > 0) 				/* normal result */
 	    {SET_FLT128_MSW64(x,(hx&0x8000ffffffffffffULL)|(k<<48)); return x;}
         if (k <= -114)
diff --git a/libquadmath/math/sqrtq.c b/libquadmath/math/sqrtq.c
index 6ed4605ed5c..f63c0d1f6d2 100644
--- a/libquadmath/math/sqrtq.c
+++ b/libquadmath/math/sqrtq.c
@@ -8,14 +8,17 @@ sqrtq (const __float128 x)
   __float128 y;
   int exp;
 
-  if (x == 0)
+  if (isnanq (x) || (isinfq (x) && x > 0))
     return x;
 
-  if (isnanq (x))
+  if (x == 0)
     return x;
 
   if (x < 0)
-    return nanq ("");
+    {
+      /* Return NaN with invalid signal.  */
+      return (x - x) / (x - x);
+    }
 
   if (x <= DBL_MAX && x >= DBL_MIN)
   {