fmaq.c (fmaq): Merge from GLIBC.

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 signgam handling.

Co-Authored-By: Joseph Myers <joseph@codesourcery.com>
From-SVN: r193099

diff --git a/libquadmath/ChangeLog b/libquadmath/ChangeLog
index efb95d179d24fa974681a11fa3162eb59c6388d1..f1b4ce9b816ef1bfe556525e0efa103682bd81bd 100644 (file)
--- a/libquadmath/ChangeLog
+++ b/libquadmath/ChangeLog
@@ -1,3 +1,10 @@
+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 signgam handling.
+
 2012-11-01  Tobias Burnus  <burnus@net-b.de>
 
        * Makefile.am (libquadmath_la_SOURCES): Add new math/* files.

diff --git a/libquadmath/math/fmaq.c b/libquadmath/math/fmaq.c
index 5616c1a278176aa43580cb197c9748b8aa67457b..e5a9d37627d6ef4cf87dbbd100a5bfcbeaae5a4f 100644 (file)
--- 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 361f7037bc3cf69e9a7cfcfce08ea8462f61725d..485939af1b43f972631b62603f4dd62214c49f58 100644 (file)
--- 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;