diff options
Diffstat (limited to 'src/math/__rem_pio2f.c')
| -rw-r--r-- | src/math/__rem_pio2f.c | 13 |
1 files changed, 10 insertions, 3 deletions
diff --git a/src/math/__rem_pio2f.c b/src/math/__rem_pio2f.c index 838e1fca..f5163856 100644 --- a/src/math/__rem_pio2f.c +++ b/src/math/__rem_pio2f.c @@ -22,12 +22,19 @@ #include "libm.h" +#if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1 +#define EPS DBL_EPSILON +#elif FLT_EVAL_METHOD==2 +#define EPS LDBL_EPSILON +#endif + /* * invpio2: 53 bits of 2/pi * pio2_1: first 25 bits of pi/2 * pio2_1t: pi/2 - pio2_1 */ static const double +toint = 1.5/EPS, invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ pio2_1 = 1.57079631090164184570e+00, /* 0x3FF921FB, 0x50000000 */ pio2_1t = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */ @@ -35,7 +42,8 @@ pio2_1t = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */ int __rem_pio2f(float x, double *y) { union {float f; uint32_t i;} u = {x}; - double tx[1],ty[1],fn; + double tx[1],ty[1]; + double_t fn; uint32_t ix; int n, sign, e0; @@ -43,8 +51,7 @@ int __rem_pio2f(float x, double *y) /* 25+53 bit pi is good enough for medium size */ if (ix < 0x4dc90fdb) { /* |x| ~< 2^28*(pi/2), medium size */ /* Use a specialized rint() to get fn. Assume round-to-nearest. */ - fn = x*invpio2 + 0x1.8p52; - fn = fn - 0x1.8p52; + fn = x*invpio2 + toint - toint; n = (int32_t)fn; *y = x - fn*pio2_1 - fn*pio2_1t; return n; |