Re: PowerPC floating point little-endian [4 of 15]

[Alan Modra <amodra at gmail dot com>]

On Wed, Jul 10, 2013 at 05:26:49PM +0000, Joseph S. Myers wrote:
> When fixing user-visible bugs, please file them in Bugzilla, and if 
> possible include testcases for them (that fail before, pass after the 
> patch) in libm-test.inc.

OK, I have some libm-test.inc additions.

> My guess is that fmod/remainder tests would be conditional on "#if defined 
> TEST_LDOUBLE && LDBL_MANT_DIG >= 106"; I'm less sure about erf/erfc, but 
> presumably by knowing where the error is, you should be able to make the 
> incorrect negation give errors of up to around 2^52 ulps, and cases with 
> large errors make suitable testcases.  The erf issue with large negative 
> input is certainly easy to test - libm-test.inc is currently lacking any 
> coverage of finite negative values.
> 
> >  	if(__builtin_expect(hx<=hy,0)) {
> > -	    if((hx<hy)||(lx<ly)) return x;	/* |x|<|y| return x */
> > +	    if (hx < hy
> > +		|| (int64_t) (lx ^ sx) < (int64_t) (ly ^ sy))
> > +		return x;			/* |x|<|y| return x */
> >  	    if(lx==ly)
> > -		return Zero[(u_int64_t)sx>>63];	/* |x|=|y| return x*0*/
> > +		return Zero[(uint64_t)sx>>63];	/* |x|=|y| return x*0*/
> 
> It looks to me like the (lx==ly) conditional needs to check (lx^sx) and 
> (ly^sy) instead to be correct, further illustrating the value of a range 
> of new testcases.

What's more, the new testcases revealed that my change wasn't even
correct for the |x|<|y| test.  The following is correct I think, and
may even be optimal.

	    if (hx < hy
		|| (((ly ^ sy) & 0x8000000000000000LL) == 0
		    && (int64_t) (lx ^ sx) < (int64_t) (ly ^ sy))
		|| (((lx ^ sx) & 0x8000000000000000LL) != 0
		    && (int64_t) (lx ^ sx) > (int64_t) (ly ^ sy)))
		return x;			/* |x|<|y| return x */
	    if ((lx ^ sx) == (ly ^ sy))
		return Zero[(uint64_t)sx>>63];	/* |x|=|y| return x*0*/

However, I found that fmodl still didn't give correct results due to
using the exponents before ldbl_extract_mantissa makes an adjustment
for the "borrow from the hidden bit" case.  I also found it was
possible to make fmodl hang due to a bogus mask of a significant bit
before the final nomalization.

	* math/libm-test.inc: Add erf, erfc, fmod and remainder tests.
	* sysdeps/ieee754/ldbl-128ibm/e_fmodl.c (__ieee754_fmodl): Rewrite
	all uses of ieee875 long double macros and unions.  Simplify test
	for 0.0L.  Correct |x|<|y| and |x|=|y| test.  Use
	ldbl_extract_mantissa value for ix,iy exponents.  Properly
	normalize after ldbl_extract_mantissa, and don't add hidden bit
	already handled.  Don't treat low word of ieee854 mantissa like
	low word of IBM long double and mask off bit when testing for
	zero.

diff --git a/math/libm-test.inc b/math/libm-test.inc
index 2e5237b..fb68e1e 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -7771,6 +7771,11 @@ static const struct test_f_f_data erf_test_data[] =
     TEST_f_f (erf, 2.0L, 0.995322265018952734162069256367252929L),
     TEST_f_f (erf, 4.125L, 0.999999994576599200434933994687765914L),
     TEST_f_f (erf, 27.0L, 1.0L),
+    TEST_f_f (erf, -27.0L, -1.0L),
+#if defined TEST_LDOUBLE && LDBL_MANT_DIG >= 54
+    /* The input is not exactly representable as a double.  */
+    TEST_f_f (erf, -0x1.fffffffffffff8p-2L, -0.5204998778130465132916303345518417673509L),
+#endif
   };
 
 static void
@@ -7799,6 +7804,10 @@ static const struct test_f_f_data erfc_test_data[] =
     TEST_f_f (erfc, 0x1.ffa002p+2L, 1.233585992097580296336099501489175967033e-29L),
     TEST_f_f (erfc, 0x1.ffffc8p+2L, 1.122671365033056305522366683719541099329e-29L),
 #ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 54
+    /* The input is not exactly representable as a double.  */
+    TEST_f_f (erfc, -0x1.fffffffffffff8p-2L, 1.52049987781304651329163033455184176735L),
+# endif
     /* The result can only be represented in long double.  */
 # if LDBL_MIN_10_EXP < -319
     TEST_f_f (erfc, 27.0L, 0.523704892378925568501606768284954709e-318L),
@@ -9293,6 +9302,13 @@ static const struct test_ff_f_data fmod_test_data[] =
 #if defined TEST_LDOUBLE && LDBL_MIN_EXP <= -16381
     TEST_ff_f (fmod, 0x0.fffffffffffffffep-16382L, 0x1p-16445L, plus_zero, NO_INEXACT_EXCEPTION),
 #endif
+#if defined TEST_LDOUBLE && LDBL_MANT_DIG >= 56
+    TEST_ff_f (fmod, -0x1.00000000000004p+0L, 0x1.fffffffffffff8p-1L, -0x1p-53L, NO_INEXACT_EXCEPTION),
+    TEST_ff_f (fmod, 0x1.fffffffffffffap-1L, 0x1.fffffffffffff8p-1L, 0x1p-56L, NO_INEXACT_EXCEPTION),
+    TEST_ff_f (fmod, -0x1.fffffffffffffap-1L, 0x1.fffffffffffff8p-1L, -0x1p-56L, NO_INEXACT_EXCEPTION),
+    TEST_ff_f (fmod, 0x1.fffffffffffffap-1L, -0x1.fffffffffffff8p-1L, 0x1p-56L, NO_INEXACT_EXCEPTION),
+    TEST_ff_f (fmod, -0x1.fffffffffffffap-1L, -0x1.fffffffffffff8p-1L, -0x1p-56L, NO_INEXACT_EXCEPTION),
+#endif
   };
 
 static void
@@ -12193,6 +12209,9 @@ static const struct test_ff_f_data remainder_test_data[] =
     TEST_ff_f (remainder, -1.625, -1.0, 0.375, NO_INEXACT_EXCEPTION),
     TEST_ff_f (remainder, 5.0, 2.0, 1.0, NO_INEXACT_EXCEPTION),
     TEST_ff_f (remainder, 3.0, 2.0, -1.0, NO_INEXACT_EXCEPTION),
+#if defined TEST_LDOUBLE && LDBL_MANT_DIG >= 56
+    TEST_ff_f (remainder, -0x1.80000000000002p1L, 2.0, 0x1.fffffffffffff8p-1L, NO_INEXACT_EXCEPTION),
+#endif
   };
 
 static void
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c b/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c
index a60963c..ffb3bec 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c
@@ -27,76 +27,65 @@ static const long double one = 1.0, Zero[] = {0.0, -0.0,};
 long double
 __ieee754_fmodl (long double x, long double y)
 {
-	int64_t n,hx,hy,hz,ix,iy,sx, i;
-	u_int64_t lx,ly,lz;
-	int temp;
+	int64_t hx, hy, hz, sx, sy;
+	uint64_t lx, ly, lz;
+	int n, ix, iy;
+	double xhi, xlo, yhi, ylo;
 
-	GET_LDOUBLE_WORDS64(hx,lx,x);
-	GET_LDOUBLE_WORDS64(hy,ly,y);
+	ldbl_unpack (x, &xhi, &xlo);
+	EXTRACT_WORDS64 (hx, xhi);
+	EXTRACT_WORDS64 (lx, xlo);
+	ldbl_unpack (y, &yhi, &ylo);
+	EXTRACT_WORDS64 (hy, yhi);
+	EXTRACT_WORDS64 (ly, ylo);
 	sx = hx&0x8000000000000000ULL;		/* sign of x */
-	hx ^=sx;				/* |x| */
-	hy &= 0x7fffffffffffffffLL;		/* |y| */
+	hx ^= sx;				/* |x| */
+	sy = hy&0x8000000000000000ULL;		/* sign of y */
+	hy ^= sy;				/* |y| */
 
     /* purge off exception values */
-	if(__builtin_expect((hy|(ly&0x7fffffffffffffff))==0 ||
+	if(__builtin_expect(hy==0 ||
 			    (hx>=0x7ff0000000000000LL)|| /* y=0,or x not finite */
 			    (hy>0x7ff0000000000000LL),0))	/* or y is NaN */
 	    return (x*y)/(x*y);
 	if(__builtin_expect(hx<=hy,0)) {
-	    if((hx<hy)||(lx<ly)) return x;	/* |x|<|y| return x */
-	    if(lx==ly)
-		return Zero[(u_int64_t)sx>>63];	/* |x|=|y| return x*0*/
+	    if (hx < hy
+		|| (((ly ^ sy) & 0x8000000000000000LL) == 0
+		    && (int64_t) (lx ^ sx) < (int64_t) (ly ^ sy))
+		|| (((lx ^ sx) & 0x8000000000000000LL) != 0
+		    && (int64_t) (lx ^ sx) > (int64_t) (ly ^ sy)))
+		return x;			/* |x|<|y| return x */
+	    if ((lx ^ sx) == (ly ^ sy))
+		return Zero[(uint64_t)sx>>63];	/* |x|=|y| return x*0*/
 	}
 
-    /* determine ix = ilogb(x) */
-	if(__builtin_expect(hx<0x0010000000000000LL,0)) {	/* subnormal x */
-	    if(hx==0) {
-		for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
-	    } else {
-		for (ix = -1022, i=(hx<<11); i>0; i<<=1) ix -=1;
-	    }
-	} else ix = (hx>>52)-0x3ff;
-
-    /* determine iy = ilogb(y) */
-	if(__builtin_expect(hy<0x0010000000000000LL,0)) {	/* subnormal y */
-	    if(hy==0) {
-		for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
-	    } else {
-		for (iy = -1022, i=(hy<<11); i>0; i<<=1) iy -=1;
-	    }
-	} else iy = (hy>>52)-0x3ff;
-
     /* Make the IBM extended format 105 bit mantissa look like the ieee854 112
        bit mantissa so the following operations will give the correct
        result.  */
-	ldbl_extract_mantissa(&hx, &lx, &temp, x);
-	ldbl_extract_mantissa(&hy, &ly, &temp, y);
+	ldbl_extract_mantissa(&hx, &lx, &ix, x);
+	ldbl_extract_mantissa(&hy, &ly, &iy, y);
 
-    /* set up {hx,lx}, {hy,ly} and align y to x */
-	if(__builtin_expect(ix >= -1022, 1))
-	    hx = 0x0001000000000000LL|(0x0000ffffffffffffLL&hx);
-	else {		/* subnormal x, shift x to normal */
-	    n = -1022-ix;
-	    if(n<=63) {
-		hx = (hx<<n)|(lx>>(64-n));
-		lx <<= n;
-	    } else {
-		hx = lx<<(n-64);
-		lx = 0;
-	    }
-	}
-	if(__builtin_expect(iy >= -1022, 1))
-	    hy = 0x0001000000000000LL|(0x0000ffffffffffffLL&hy);
-	else {		/* subnormal y, shift y to normal */
-	    n = -1022-iy;
-	    if(n<=63) {
-		hy = (hy<<n)|(ly>>(64-n));
-		ly <<= n;
-	    } else {
-		hy = ly<<(n-64);
-		ly = 0;
-	    }
-	}
+	if (__builtin_expect (ix == -IEEE754_DOUBLE_BIAS, 0))
+	  {
+	    /* subnormal x, shift x to normal.  */
+	    while ((hx & (1LL << 48)) == 0)
+	      {
+		hx = (hx << 1) | (lx >> 63);
+		lx = lx << 1;
+		ix -= 1;
+	      }
+	  }
+
+	if (__builtin_expect (iy == -IEEE754_DOUBLE_BIAS, 0))
+	  {
+	    /* subnormal y, shift y to normal.  */
+	    while ((hy & (1LL << 48)) == 0)
+	      {
+		hy = (hy << 1) | (ly >> 63);
+		ly = ly << 1;
+		iy -= 1;
+	      }
+	  }
 
     /* fix point fmod */
 	n = ix - iy;
@@ -104,7 +93,7 @@ __ieee754_fmodl (long double x, long double y)
 	    hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
 	    if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;}
 	    else {
-		if((hz|(lz&0x7fffffffffffffff))==0)		/* return sign(x)*0 */
+		if((hz|lz)==0)		/* return sign(x)*0 */
 		    return Zero[(u_int64_t)sx>>63];
 		hx = hz+hz+(lz>>63); lx = lz+lz;
 	    }
@@ -113,7 +102,7 @@ __ieee754_fmodl (long double x, long double y)
 	if(hz>=0) {hx=hz;lx=lz;}
 
     /* convert back to floating value and restore the sign */
-	if((hx|(lx&0x7fffffffffffffff))==0)			/* return sign(x)*0 */
+	if((hx|lx)==0)			/* return sign(x)*0 */
 	    return Zero[(u_int64_t)sx>>63];
 	while(hx<0x0001000000000000LL) {	/* normalize x */
 	    hx = hx+hx+(lx>>63); lx = lx+lx;


-- 
Alan Modra
Australia Development Lab, IBM