Patchwork [committed] libquadmath: Small fmaq and lgamma update

login
register
mail settings
Submitter Tobias Burnus
Date Nov. 2, 2012, 4:59 p.m.
Message ID <5093FC03.7000306@net-b.de>
Download mbox | patch
Permalink /patch/196604/
State New
Headers show

Comments

Tobias Burnus - Nov. 2, 2012, 4:59 p.m.
The attached patch ports Joseph's GLIBC fmal patch from today; 
additionally, I made a small adjustment to lgamma's signgam handling.

Build and tested on x86-64-gnu-linux and applied as Rev. 193099


Otherwise, it still applies what I wrote at 
http://gcc.gnu.org/ml/gcc-patches/2012-11/msg00042.html:

Besides those isses [i.e. the test failures], one should - as mentioned 
in the yesterday's thread - update strtod and printf; the GLIBC test 
should be modified in such a way that they work properly. I would also 
like to see the test automatically be run (with "make test") and 
possibly to have some additional DejaGNU support for 
libquadmath-specific tests. (One could also add support for __complex__ 
__float128 as Andreas suggested.)

Tobias
Paolo Carlini - Nov. 2, 2012, 6:01 p.m.
On 11/02/2012 05:59 PM, Tobias Burnus wrote:
> The attached patch ports Joseph's GLIBC fmal patch from today; 
> additionally, I made a small adjustment to lgamma's signgam handling.
>
> Build and tested on x86-64-gnu-linux and applied as Rev. 193099
My builds are now failing in libquadmath with IEEE854_FLT128_DOUBLE_BIAS 
undeclared. I'm not doing anything special...

Paolo.
Tobias Burnus - Nov. 2, 2012, 6:10 p.m.
Paolo Carlini wrote:
> On 11/02/2012 05:59 PM, Tobias Burnus wrote:
>> The attached patch ports Joseph's GLIBC fmal patch from today; 
>> additionally, I made a small adjustment to lgamma's signgam handling.
>>
>> Build and tested on x86-64-gnu-linux and applied as Rev. 193099
> My builds are now failing in libquadmath with 
> IEEE854_FLT128_DOUBLE_BIAS undeclared. I'm not doing anything special...

Sorry, I seemingly applied an old premature version of the patch instead 
of the one I had build and tested. Fixed by the attached patch (Rev. 
193100).

Tobias
Paolo Carlini - Nov. 2, 2012, 6:24 p.m.
On 11/02/2012 07:10 PM, Tobias Burnus wrote:
> Paolo Carlini wrote:
>> On 11/02/2012 05:59 PM, Tobias Burnus wrote:
>>> The attached patch ports Joseph's GLIBC fmal patch from today; 
>>> additionally, I made a small adjustment to lgamma's signgam handling.
>>>
>>> Build and tested on x86-64-gnu-linux and applied as Rev. 193099
>> My builds are now failing in libquadmath with 
>> IEEE854_FLT128_DOUBLE_BIAS undeclared. I'm not doing anything special...
> Sorry, I seemingly applied an old premature version of the patch 
> instead of the one I had build and tested. Fixed by the attached patch 
> (Rev. 193100).
Thanks!

Paolo.

Patch

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;