libm: fix rint/scalb testcase failures

These failures no longer happen:

Failure: Test: scalb (2.0, 0.5) == NaN plus invalid exception
Failure: Test: scalb (3.0, -2.5) == NaN plus invalid exception
Failure: Test: rint (0.5) == 0.0
Failure: Test: rint (1.5) == 2.0
Failure: Test: rint (2.5) == 2.0
Failure: Test: rint (3.5) == 4.0
Failure: Test: rint (4.5) == 4.0
Failure: Test: rint (-0.5) == -0.0
Failure: Test: rint (-1.5) == -2.0
Failure: Test: rint (-2.5) == -2.0
Failure: Test: rint (-3.5) == -4.0
Failure: Test: rint (-4.5) == -4.0

Signed-off-by: Denys Vlasenko <vda.linux@googlemail.com>

1 files changed, 35 insertions, 19 deletions

diff --git a/libm/s_rint.c b/libm/s_rint.c
index 358ce76b4..06432c622 100644
--- a/libm/s_rint.c
+++ b/libm/s_rint.c
@@ -30,41 +30,57 @@ TWO52[2]={
 
 double rint(double x)
 {
-	int32_t i0,j0,sx;
+	int32_t i0, j0, sx;
 	u_int32_t i,i1;
-	double w,t;
+	double t;
+	/* We use w = x + 2^52; t = w - 2^52; trick to round x to integer.
+	 * This trick requires that compiler does not optimize it
+	 * by keeping intermediate result w in a register wider than double.
+	 * Declaring w volatile assures that value gets truncated to double
+	 * (unfortunately, it also forces store+load):
+	 */
+	volatile double w;
+
 	EXTRACT_WORDS(i0,i1,x);
-	sx = (i0>>31)&1;
-	j0 = ((i0>>20)&0x7ff)-0x3ff;
-	if(j0<20) {
-	    if(j0<0) {
-		if(((i0&0x7fffffff)|i1)==0) return x;
+	/* Unbiased exponent */
+	j0 = ((((u_int32_t)i0) >> 20)&0x7ff)-0x3ff;
+
+	if (j0 > 51) {
+		//Why bother? Just returning x works too
+		//if (j0 == 0x400)  /* inf or NaN */
+		//	return x+x;
+		return x;  /* x is integral */
+	}
+
+	/* Sign */
+	sx = ((u_int32_t)i0) >> 31;
+
+	if (j0<20) {
+	    if (j0<0) { /* |x| < 1 */
+		if (((i0&0x7fffffff)|i1)==0) return x;
 		i1 |= (i0&0x0fffff);
 		i0 &= 0xfffe0000;
 		i0 |= ((i1|-i1)>>12)&0x80000;
 		SET_HIGH_WORD(x,i0);
-	        w = TWO52[sx]+x;
-	        t =  w-TWO52[sx];
+		w = TWO52[sx]+x;
+		t = w-TWO52[sx];
 		GET_HIGH_WORD(i0,t);
 		SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
-	        return t;
+		return t;
 	    } else {
 		i = (0x000fffff)>>j0;
-		if(((i0&i)|i1)==0) return x; /* x is integral */
+		if (((i0&i)|i1)==0) return x; /* x is integral */
 		i>>=1;
-		if(((i0&i)|i1)!=0) {
-		    if(j0==19) i1 = 0x40000000; else
-		    i0 = (i0&(~i))|((0x20000)>>j0);
+		if (((i0&i)|i1)!=0) {
+		    if (j0==19) i1 = 0x40000000;
+		    else i0 = (i0&(~i))|((0x20000)>>j0);
 		}
 	    }
-	} else if (j0>51) {
-	    if(j0==0x400) return x+x;	/* inf or NaN */
-	    else return x;		/* x is integral */
 	} else {
 	    i = ((u_int32_t)(0xffffffff))>>(j0-20);
-	    if((i1&i)==0) return x;	/* x is integral */
+	    if ((i1&i)==0) return x;	/* x is integral */
 	    i>>=1;
-	    if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
+	    if ((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
 	}
 	INSERT_WORDS(x,i0,i1);
 	w = TWO52[sx]+x;