Patchwork Update libquadmath from GLIBC

login
register
mail settings
Submitter Tobias Burnus
Date Oct. 31, 2012, 3:52 p.m.
Message ID <50914933.8070300@net-b.de>
Download mbox | patch
Permalink /patch/195905/
State New
Headers show

Comments

Tobias Burnus - Oct. 31, 2012, 3:52 p.m.
Tobias Burnus wrote:
> libquadmath's math functions are based on (but not identical to) 
> GLIBC's sysdeps/ieee754/ldbl-128 functions. In the attached patch, I 
> have ported the bug fixes from GLIBC over to libquadmath. Hopefully, 
> the port is complete and correct.

Slightly updated version, committed as Rev. 193037.

Jakub suggested to use a modified version of GLIBC's test-ldouble.c, 
which I did. Before the patch, the result was (using the too strict 
sysdeps/x86_64/fpu/libm-test-ulps):

   5644 test cases plus 4809 tests for exception flags executed.
   815 errors occurred.

With the GLIBC merge, it went down to 776 errors (and reduced ulp for 
some tests) - and with the sqrtq/csqrtq patch, the failure number went 
down to 708.

Tobias
Jakub Jelinek - Oct. 31, 2012, 3:56 p.m.
On Wed, Oct 31, 2012 at 04:52:19PM +0100, Tobias Burnus wrote:
> Tobias Burnus wrote:
> >libquadmath's math functions are based on (but not identical to)
> >GLIBC's sysdeps/ieee754/ldbl-128 functions. In the attached patch,
> >I have ported the bug fixes from GLIBC over to libquadmath.
> >Hopefully, the port is complete and correct.
> 
> Slightly updated version, committed as Rev. 193037.

Thanks.

> Jakub suggested to use a modified version of GLIBC's test-ldouble.c,
> which I did. Before the patch, the result was (using the too strict
> sysdeps/x86_64/fpu/libm-test-ulps):
> 
>   5644 test cases plus 4809 tests for exception flags executed.
>   815 errors occurred.
> 
> With the GLIBC merge, it went down to 776 errors (and reduced ulp
> for some tests) - and with the sqrtq/csqrtq patch, the failure
> number went down to 708.

I think it would be nice if you also posted the changes you did to
test-ldouble.c and libm-test.inc, so that next time we could more easily
test it again.

	Jakub
Joseph S. Myers - Oct. 31, 2012, 4:58 p.m.
On Wed, 31 Oct 2012, Tobias Burnus wrote:

> Tobias Burnus wrote:
> > libquadmath's math functions are based on (but not identical to) GLIBC's
> > sysdeps/ieee754/ldbl-128 functions. In the attached patch, I have ported the
> > bug fixes from GLIBC over to libquadmath. Hopefully, the port is complete
> > and correct.
> 
> Slightly updated version, committed as Rev. 193037.

At a glance, I don't see any sign of

commit c0df8e693f34b535bd6ee1b691bc4ca6bc3b4579
Author: Joseph Myers <joseph@codesourcery.com>
Date:   Thu Mar 22 12:52:50 2012 +0000

    Fix low-part sign handling in sin/cos for ldbl-128 and ldbl-128ibm.

being included (although the effects of that change were only small 
reductions, on average, in ulps errors); I don't know if other relevant 
commits to the ldbl-128 code (since it was last merged from glibc) are 
missing.  (And unsurprisingly it doesn't have the fma change I committed 
today relating to spurious exceptions on after-rounding architectures - 
FWIW, all the architectures currently using libquadmath are 
after-rounding.)

There have been lots of improvements to the complex functions in glibc 
lately (those are generally in math/ not sysdeps/ieee754/), although 
libquadmath seems to be following glibc's versions of those less closely, 
and fixing those in glibc is still very much a work in progress (inverse 
trig and hyperbolic functions, and cpow, still need major work to be 
reasonably accurate for all inputs and avoid spurious exceptions / produce 
correct exceptions).

Separately, soft-fp in GCC could do with being updated from glibc - that 
might eliminate a few spurious underflow exceptions if you get those in 
your test results.  (Note that current soft-fp should, ideally, have 
FP_TRAPPING_EXCEPTIONS definitions added to sfp-machine.h files, where 
hardware exceptions are supported at all by those files.  The relevant 
architecture maintainers - aarch64, i386, ia64 - could reasonably be asked 
to add those.)

printf and strtod code in libquadmath looks like it could do with being 
updated from glibc as well - strtod in particular has had major changes in 
glibc since it was last updated.

Patch

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.

diff --git a/libquadmath/math/complex.c b/libquadmath/math/complex.c
index f67448a..9953f52 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 510c98f..8cfdd8e 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 f22e9d7..663c928 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 126b0a2..23e3188 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 d82947a..56a1836 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 8e92c5a..2071248 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 fe7a591..8a93fdb 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 75997f6..b233225 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 b7049a7..f0852ee 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 6ed4605..f63c0d1 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)
   {