Re: [PATCH] Fix ldbl128ibm fmodl for subnormals

[Andreas Schwab <schwab at linux-m68k dot org>]

Adhemerval Zanella <azanella@linux.vnet.ibm.com> writes:

> Current subnormal tests for fmod (TEST_ff_f (fmod, 0x0.fffffffffffffp-1022L, 0x1p-1074L, plus_zero)) 
> fails for ldbl-128ibm. This patch fixes it and it is related to commit 'Fix fmod for subnormals' 
> (c5bfe3d5ba29d36563f1e4bd4f8d7336093ee6fc). 
>
> Tested on ppc32 and ppc64.

How did you test it?  I get an infinite loop here:

	while(hx<0x0001000000000000LL) {	/* normalize x */
	    hx = hx+hx+(lx>>63); lx = lx+lx;
	    iy -= 1;
	}

where hx == 0xfffe000000000000 (ie. negative).

The following doesn't look correct:

	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;
	    }

First, lx is the image of the low double, not the low half of a 128 bit
IEEE long double.  Second, since i is unsigned, the loops shift out all
bits, not just the leading zero bits.  In ldbl-128/e_fmodl.c, i is
signed, relying on undefined behaviour.  The same bug is present in
ilogbl for both ldbl-128 and ldbl-128ibm.

David, please test the ldbl-128 fmodl/ilogbl, I may have added some
off-by-one.

Andreas.

	* sysdeps/ieee754/ldbl-128/e_fmodl.c (__ieee754_fmodl): Use
	__builtin_clzll to calculate exponent of subnormal number.
	* sysdeps/ieee754/ldbl-128ibm/e_fmodl.c (__ieee754_fmodl):
	Likewise.
	* sysdeps/ieee754/ldbl-128/e_ilogbl.c (__ieee754_ilogbl):
	Likewise.
	* sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c (__ieee754_ilogbl):
	Likewise.
---
 sysdeps/ieee754/ldbl-128/e_fmodl.c     |   12 +++++++-----
 sysdeps/ieee754/ldbl-128/e_ilogbl.c    |    7 +++++--
 sysdeps/ieee754/ldbl-128ibm/e_fmodl.c  |   14 +++-----------
 sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c |    6 +-----
 4 files changed, 16 insertions(+), 23 deletions(-)

diff --git a/sysdeps/ieee754/ldbl-128/e_fmodl.c b/sysdeps/ieee754/ldbl-128/e_fmodl.c
index 3328003..9b27a4f 100644
--- a/sysdeps/ieee754/ldbl-128/e_fmodl.c
+++ b/sysdeps/ieee754/ldbl-128/e_fmodl.c
@@ -26,7 +26,7 @@ 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;
+	int64_t n,hx,hy,hz,ix,iy,sx;
 	u_int64_t lx,ly,lz;
 
 	GET_LDOUBLE_WORDS64(hx,lx,x);
@@ -47,19 +47,21 @@ __ieee754_fmodl (long double x, long double y)
 
     /* determine ix = ilogb(x) */
 	if(hx<0x0001000000000000LL) {	/* subnormal x */
+	    ix = -16382;
 	    if(hx==0) {
-		for (ix = -16431, i=lx; i>0; i<<=1) ix -=1;
+		ix -= __builtin_clzll (lx) + 49;
 	    } else {
-		for (ix = -16382, i=hx<<15; i>0; i<<=1) ix -=1;
+		ix -= __builtin_clzll (hx << 15);
 	    }
 	} else ix = (hx>>48)-0x3fff;
 
     /* determine iy = ilogb(y) */
 	if(hy<0x0001000000000000LL) {	/* subnormal y */
+	    iy = -16382;
 	    if(hy==0) {
-		for (iy = -16431, i=ly; i>0; i<<=1) iy -=1;
+		iy -= __builtin_clzll (ly) + 49;
 	    } else {
-		for (iy = -16382, i=hy<<15; i>0; i<<=1) iy -=1;
+		iy -= __builtin_clzll (hy << 15);
 	    }
 	} else iy = (hy>>48)-0x3fff;
 
diff --git a/sysdeps/ieee754/ldbl-128/e_ilogbl.c b/sysdeps/ieee754/ldbl-128/e_ilogbl.c
index 0a47649..17e4d12 100644
--- a/sysdeps/ieee754/ldbl-128/e_ilogbl.c
+++ b/sysdeps/ieee754/ldbl-128/e_ilogbl.c
@@ -39,11 +39,14 @@ int __ieee754_ilogbl (long double x)
 	    if((hx|lx)==0)
 		return FP_ILOGB0;	/* ilogbl(0) = FP_ILOGB0 */
 	    else			/* subnormal x */
+	      {
+		ix = -16382;
 		if(hx==0) {
-		    for (ix = -16431; lx>0; lx<<=1) ix -=1;
+		    ix -= __builtin_clzll (lx) + 49;
 		} else {
-		    for (ix = -16382, hx<<=15; hx>0; hx<<=1) ix -=1;
+		    ix -= __builtin_clzll (hx << 15);
 		}
+	      }
 	    return ix;
 	}
 	else if (hx<0x7fff000000000000LL) return (hx>>48)-0x3fff;
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c b/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c
index 31e709b..7e075df 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c
@@ -28,7 +28,7 @@ long double
 __ieee754_fmodl (long double x, long double y)
 {
 	int64_t n,hx,hy,hz,ix,iy,sx;
-	u_int64_t lx,ly,lz, i;
+	u_int64_t lx,ly,lz;
 	int temp;
 
 	GET_LDOUBLE_WORDS64(hx,lx,x);
@@ -50,20 +50,12 @@ __ieee754_fmodl (long double x, long double y)
 
     /* 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;
-	    }
+	  ix = -1022 - __builtin_clzll (hx << 11);
 	} 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;
-	    }
+	  iy = -1022 - __builtin_clzll (hy << 11);
 	} else iy = (hy>>52)-0x3ff;
 
     /* Make the IBM extended format 105 bit mantissa look like the ieee854 112
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c b/sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c
index 55f87ed..d483c7a 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c
@@ -40,11 +40,7 @@ int __ieee754_ilogbl(long double x)
 	    if((hx|(lx&0x7fffffffffffffffLL))==0)
 		return FP_ILOGB0;	/* ilogbl(0) = FP_ILOGB0 */
 	    else			/* subnormal x */
-		if(hx==0) {
-		    for (ix = -1043; lx>0; lx<<=1) ix -=1;
-		} else {
-		    for (ix = -1022, hx<<=11; hx>0; hx<<=1) ix -=1;
-		}
+		ix = -1022 - __builtin_clzll (hx << 11);
 	    return ix;
 	}
 	else if (hx<0x7ff0000000000000LL) return (hx>>52)-0x3ff;
-- 
1.7.10.4


-- 
Andreas Schwab, schwab@linux-m68k.org
GPG Key fingerprint = 58CA 54C7 6D53 942B 1756  01D3 44D5 214B 8276 4ED5
"And now for something completely different."