[libm] Fix undefined behavior of a left shifted of a signed integer

The patch fixes a few instances of left shifts on
signed integer entities.  A 'static inline' helper function
'subnormal_ilogb()' has been added to math_private.h.  This
function is then used e_fmod.c, s_ilogb(), and s_remquo.c.

The change in s_remquo.c has only been compile tested.

The change to e_fmod.c has been test on over 3 billion pairs
of subnormal numbers where testing included x > y and x < y
pairs.  The test compared the output from fmod() with the
output from mpfr_fmod() from MPFR.  There were no difference.

The change to s_ilogb() has had limited testing where its
output was compared against frexp().  In this testing, no
differences in output were detected.

PR:	288778
MFC after:	1 week

lib/msun/src/e_fmod.c

@@ -26,14 +26,14 @@ static const double one = 1.0, Zero[] = {0.0, -0.0,};
 double
 fmod(double x, double y)
 {
-	int32_t n,hx,hy,hz,ix,iy,sx,i;
-	u_int32_t lx,ly,lz;
+	int32_t hx, hy, hz, ix, iy, n, sx;
+	u_int32_t lx, ly, lz;
 
 	EXTRACT_WORDS(hx,lx,x);
 	EXTRACT_WORDS(hy,ly,y);
 	sx = hx&0x80000000;		/* sign of x */
-	hx ^=sx;		/* |x| */
-	hy &= 0x7fffffff;	/* |y| */
+	hx ^= sx;			/* |x| */
+	hy &= 0x7fffffff;		/* |y| */
 
     /* purge off exception values */
 	if((hy|ly)==0||(hx>=0x7ff00000)||	/* y=0,or x not finite */
@@ -46,22 +46,16 @@ fmod(double x, double y)
 	}
 
     /* determine ix = ilogb(x) */
-	if(hx<0x00100000) {	/* 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>>20)-1023;
+	if(hx<0x00100000)
+	    ix = subnormal_ilogb(hx, lx);
+	else
+	    ix = (hx>>20)-1023;
 
     /* determine iy = ilogb(y) */
-	if(hy<0x00100000) {	/* 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>>20)-1023;
+	if(hy<0x00100000)
+	    iy = subnormal_ilogb(hy, ly);
+	else
+	    iy = (hy>>20)-1023;
 
     /* set up {hx,lx}, {hy,ly} and align y to x */
 	if(ix >= -1022) 

lib/msun/src/math_private.h

@@ -739,6 +739,27 @@ irintl(long double x)
 	(ar) = (x) - (ai);				\
 } while (0)
 
+/*
+ * For a double entity split into high and low parts, compute ilogb.
+ */
+static inline int32_t
+subnormal_ilogb(int32_t hi, int32_t lo)
+{
+	int32_t j;
+	uint32_t i;
+
+	j = -1022;
+	if (hi == 0) {
+	    j -= 21;
+	    i = (uint32_t)lo;
+	} else
+	    i = (uint32_t)hi << 11;
+
+	for (; i < 0x7fffffff; i <<= 1) j -= 1;
+
+	return (j);
+}
+
 #ifdef DEBUG
 #if defined(__amd64__) || defined(__i386__)
 #define	breakpoint()	asm("int $3")

lib/msun/src/s_ilogb.c

@@ -21,21 +21,18 @@
 #include "math.h"
 #include "math_private.h"
 
-	int ilogb(double x)
+int
+ilogb(double x)
 {
-	int32_t hx,lx,ix;
+	int32_t hx, ix, lx;
 
 	EXTRACT_WORDS(hx,lx,x);
 	hx &= 0x7fffffff;
 	if(hx<0x00100000) {
 	    if((hx|lx)==0)
 		return 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;
-		}
+	    else
+		ix = subnormal_ilogb(hx, lx);
 	    return ix;
 	}
 	else if (hx<0x7ff00000) return (hx>>20)-1023;

lib/msun/src/s_remquo.c

@@ -4,7 +4,7 @@
  *
  * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -27,7 +27,7 @@ static const double Zero[] = {0.0, -0.0,};
 double
 remquo(double x, double y, int *quo)
 {
-	int32_t n,hx,hy,hz,ix,iy,sx,i;
+	int32_t hx,hy,hz,ix,iy,n,sx;
 	u_int32_t lx,ly,lz,q,sxy;
 
 	EXTRACT_WORDS(hx,lx,x);
@@ -53,25 +53,19 @@ remquo(double x, double y, int *quo)
 	}
 
     /* determine ix = ilogb(x) */
-	if(hx<0x00100000) {	/* 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>>20)-1023;
+	if(hx<0x00100000)
+	    ix = subnormal_ilogb(hx, lx);
+	else
+	    ix = (hx>>20)-1023;
 
     /* determine iy = ilogb(y) */
-	if(hy<0x00100000) {	/* 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>>20)-1023;
+	if(hy<0x00100000)
+	    iy = subnormal_ilogb(hy, ly);
+	else
+	    iy = (hy>>20)-1023;
 
     /* set up {hx,lx}, {hy,ly} and align y to x */
-	if(ix >= -1022) 
+	if(ix >= -1022)
 	    hx = 0x00100000|(0x000fffff&hx);
 	else {		/* subnormal x, shift x to normal */
 	    n = -1022-ix;
@@ -83,7 +77,7 @@ remquo(double x, double y, int *quo)
 		lx = 0;
 	    }
 	}
-	if(iy >= -1022) 
+	if(iy >= -1022)
 	    hy = 0x00100000|(0x000fffff&hy);
 	else {		/* subnormal y, shift y to normal */
 	    n = -1022-iy;