[committed] libquadmath: Small fmaq and lgamma update

2012-11-01  Tobias Burnus  <burnus@net-b.de>
	    Joseph Myers  <joseph@codesourcery.com>

        * math/fmaq.c (fmaq): Merge from GLIBC. Handle cases
	with small x * y using scaling, not as x * y + z.
	* math/lgammaq.c (lgammaq): Fix to signgam handling.

diff --git a/libquadmath/math/fmaq.c b/libquadmath/math/fmaq.c
index 5616c1a..e5a9d37 100644
--- a/libquadmath/math/fmaq.c
+++ b/libquadmath/math/fmaq.c
@@ -73,6 +73,37 @@  fmaq (__float128 x, __float128 y, __float128 z)
 	  || u.ieee.exponent + v.ieee.exponent
 	     < IEEE854_FLOAT128_BIAS - FLT128_MANT_DIG - 2)
 	return x * y + z;
+      /* If x * y is less than 1/4 of FLT128_DENORM_MIN, neither the
+	 result nor whether there is underflow depends on its exact
+	 value, only on its sign.  */
+      if (u.ieee.exponent + v.ieee.exponent
+	  < IEEE854_FLT128_DOUBLE_BIAS - FLT128_MANT_DIG - 2)
+	{
+	  int neg = u.ieee.negative ^ v.ieee.negative;
+	  __float128 tiny = neg ? -0x1p-16494L : 0x1p-16494L;
+	  if (w.ieee.exponent >= 3)
+	    return tiny + z;
+	  /* Scaling up, adding TINY and scaling down produces the
+	     correct result, because in round-to-nearest mode adding
+	     TINY has no effect and in other modes double rounding is
+	     harmless.  But it may not produce required underflow
+	     exceptions.  */
+	  v.value = z * 0x1p114L + tiny;
+	  if (TININESS_AFTER_ROUNDING
+	      ? v.ieee.exponent < 115
+	      : (w.ieee.exponent == 0
+		 || (w.ieee.exponent == 1
+		     && w.ieee.negative != neg
+		     && w.ieee.mantissa3 == 0
+		     && w.ieee.mantissa2 == 0
+		     && w.ieee.mantissa1 == 0
+		     && w.ieee.mantissa0 == 0)))
+	    {
+	      volatile __float128 force_underflow = x * y;
+	      (void) force_underflow;
+	    }
+	  return v.value * 0x1p-114L;
+	}
       if (u.ieee.exponent + v.ieee.exponent
 	  >= 0x7fff + IEEE854_FLOAT128_BIAS - FLT128_MANT_DIG)
 	{
@@ -228,17 +259,17 @@  fmaq (__float128 x, __float128 y, __float128 z)
 	 for proper rounding.  */
       if (v.ieee.exponent == 226)
 	{
-         /* If the exponent would be in the normal range when
-            rounding to normal precision with unbounded exponent
-            range, the exact result is known and spurious underflows
-            must be avoided on systems detecting tininess after
-            rounding.  */
-         if (TININESS_AFTER_ROUNDING)
-           {
-             w.value = a1 + u.value;
-             if (w.ieee.exponent == 227)
-               return w.value * 0x1p-226L;
-           }
+	  /* If the exponent would be in the normal range when
+	     rounding to normal precision with unbounded exponent
+	     range, the exact result is known and spurious underflows
+	     must be avoided on systems detecting tininess after
+	     rounding.  */
+	  if (TININESS_AFTER_ROUNDING)
+	    {
+	      w.value = a1 + u.value;
+	      if (w.ieee.exponent == 227)
+		return w.value * 0x1p-226L;
+	    }
 	  /* 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. */
diff --git a/libquadmath/math/lgammaq.c b/libquadmath/math/lgammaq.c
index 361f703..485939a 100644
--- a/libquadmath/math/lgammaq.c
+++ b/libquadmath/math/lgammaq.c
@@ -759,7 +759,8 @@  lgammaq (__float128 x)
 {
   __float128 p, q, w, z, nx;
   int i, nn;
-  int sign;
+
+  signgam = 1;
 
   if (! finiteq (x))
     return x * x;
@@ -767,11 +768,9 @@  lgammaq (__float128 x)
   if (x == 0.0Q)
     {
       if (signbitq (x))
-	sign = -1;
+	signgam = -1;
     }
 
-  signgam = sign;
-
   if (x < 0.0Q)
     {
       q = -x;
@@ -780,9 +779,9 @@  lgammaq (__float128 x)
 	return (one / (p - p));
       i = p;
       if ((i & 1) == 0)
-	sign = -1;
+	signgam = -1;
       else
-	sign = 1;
+	signgam = 1;
       z = q - p;
       if (z > 0.5Q)
 	{
@@ -790,10 +789,8 @@  lgammaq (__float128 x)
 	  z = p - q;
 	}
       z = q * sinq (PIQ * z);
-      signgam = sign;
-
       if (z == 0.0Q)
-	return (sign * huge * huge);
+	return (signgam * huge * huge);
       w = lgammaq (q);
       z = logq (PIQ / z) - w;
       return (z);
@@ -1025,7 +1022,7 @@  lgammaq (__float128 x)
     }
 
   if (x > MAXLGM)
-    return (sign * huge * huge);
+    return (signgam * huge * huge);
 
   q = ls2pi - x;
   q = (x - 0.5Q) * logq (x) + q;