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>

diff --git a/libm/s_rint.c b/libm/s_rint.c
index 358ce76..06432c6 100644 (file)
--- 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;

diff --git a/test/math/rint.c b/test/math/rint.c
index 04c1953..c7bfab9 100644 (file)
--- a/test/math/rint.c
+++ b/test/math/rint.c
@@ -1,11 +1,32 @@
 #include <math.h>
 #include <stdlib.h>
+#include <stdint.h>
 #include <stdio.h>
 
-int main(void) {
-    double d1, d2;
-    d1 = 0.6;  d2 = rint(d1);
-    printf("d1 = %f, d2 = %f\n", d1, d2);
-    return 0;
-}
+#define check_d1(func, param, expected) \
+do { \
+       int err; hex_union ur; hex_union up; \
+       up.f = param; ur.f = result = func(param); \
+       errors += (err = (result != expected)); \
+       err \
+       ? printf("FAIL: %s(%g/"HEXFMT")=%g/"HEXFMT" (expected %g)\n", \
+               #func, (param), (long long)up.hex, result, (long long)ur.hex, expected) \
+       : printf("PASS: %s(%g)=%g\n", #func, (param), result); \
+} while (0)
+
+#define HEXFMT "%08llx"
+typedef union {
+       double f;
+       uint64_t hex;
+} hex_union;
+double result;
+
+int errors = 0;
 
+int main(void)
+{
+       check_d1(rint, 0.6, 1.0);
+
+        printf("Errors: %d\n", errors);
+        return errors;
+}

diff --git a/test/math/signgam.c b/test/math/signgam.c
index c60375a..d79d6af 100644 (file)
--- a/test/math/signgam.c
+++ b/test/math/signgam.c
@@ -5,14 +5,23 @@
 double zero = 0.0;
 double mzero;
 
-int
-main (void)
+int main(void)
 {
-  double d;
-  mzero = copysign (zero, -1.0);
-  d = lgamma (zero);
-  printf ("%g %d\n", d, signgam);
-  d = lgamma (mzero);
-  printf ("%g %d\n", d, signgam);
-  return 0;
+       double d;
+       int errors = 0;
+
+       mzero = copysign(zero, -1.0);
+
+       d = lgamma(zero);
+       printf("%g %d\n", d, signgam);
+       errors += !(d == HUGE_VAL);
+       errors += !(signgam == 1);
+
+       d = lgamma(mzero);
+       printf("%g %d\n", d, signgam);
+       errors += !(d == HUGE_VAL);
+       errors += !(signgam == -1);
+
+       printf("Errors: %d\n", errors);
+       return errors;
 }