diff options
| -rw-r--r-- | libm/s_rint.c | 54 | ||||
| -rw-r--r-- | test/math/rint.c | 33 | ||||
| -rw-r--r-- | test/math/signgam.c | 27 |
3 files changed, 80 insertions, 34 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; diff --git a/test/math/rint.c b/test/math/rint.c index 04c195385..c7bfab920 100644 --- 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 c60375aec..d79d6afb2 100644 --- 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; } |