summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--libm/s_rint.c54
-rw-r--r--test/math/rint.c33
-rw-r--r--test/math/signgam.c27
3 files changed, 80 insertions, 34 deletions
diff --git a/libm/s_rint.c b/libm/s_rint.c
index 358ce76..06432c6 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 04c1953..c7bfab9 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 c60375a..d79d6af 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;
}