diff options
| -rw-r--r-- | src/math/__rem_pio2.c | 15 | ||||
| -rw-r--r-- | src/math/__rem_pio2f.c | 13 | ||||
| -rw-r--r-- | src/math/__rem_pio2l.c | 16 |
3 files changed, 41 insertions, 3 deletions
diff --git a/src/math/__rem_pio2.c b/src/math/__rem_pio2.c index d403f81c..dcf672fb 100644 --- a/src/math/__rem_pio2.c +++ b/src/math/__rem_pio2.c @@ -36,6 +36,7 @@ */ static const double toint = 1.5/EPS, +pio4 = 0x1.921fb54442d18p-1, invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */ @@ -117,11 +118,23 @@ int __rem_pio2(double x, double *y) } if (ix < 0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */ medium: - /* rint(x/(pi/2)), Assume round-to-nearest. */ + /* rint(x/(pi/2)) */ fn = (double_t)x*invpio2 + toint - toint; n = (int32_t)fn; r = x - fn*pio2_1; w = fn*pio2_1t; /* 1st round, good to 85 bits */ + /* Matters with directed rounding. */ + if (predict_false(r - w < -pio4)) { + n--; + fn--; + r = x - fn*pio2_1; + w = fn*pio2_1t; + } else if (predict_false(r - w > pio4)) { + n++; + fn++; + r = x - fn*pio2_1; + w = fn*pio2_1t; + } y[0] = r - w; u.f = y[0]; ey = u.i>>52 & 0x7ff; diff --git a/src/math/__rem_pio2f.c b/src/math/__rem_pio2f.c index 4473c1c4..e6765643 100644 --- a/src/math/__rem_pio2f.c +++ b/src/math/__rem_pio2f.c @@ -35,6 +35,7 @@ */ static const double toint = 1.5/EPS, +pio4 = 0x1.921fb6p-1, invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ pio2_1 = 1.57079631090164184570e+00, /* 0x3FF921FB, 0x50000000 */ pio2_1t = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */ @@ -50,10 +51,20 @@ int __rem_pio2f(float x, double *y) ix = u.i & 0x7fffffff; /* 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. */ + /* Use a specialized rint() to get fn. */ fn = (double_t)x*invpio2 + toint - toint; n = (int32_t)fn; *y = x - fn*pio2_1 - fn*pio2_1t; + /* Matters with directed rounding. */ + if (predict_false(*y < -pio4)) { + n--; + fn--; + *y = x - fn*pio2_1 - fn*pio2_1t; + } else if (predict_false(*y > pio4)) { + n++; + fn++; + *y = x - fn*pio2_1 - fn*pio2_1t; + } return n; } if(ix>=0x7f800000) { /* x is inf or NaN */ diff --git a/src/math/__rem_pio2l.c b/src/math/__rem_pio2l.c index 77255bd8..236b2def 100644 --- a/src/math/__rem_pio2l.c +++ b/src/math/__rem_pio2l.c @@ -44,6 +44,7 @@ pio2_1 = 1.57079632679597125389e+00, /* 0x3FF921FB, 0x54444000 */ pio2_2 = -1.07463465549783099519e-12, /* -0x12e7b967674000.0p-92 */ pio2_3 = 6.36831716351370313614e-25; /* 0x18a2e037074000.0p-133 */ static const long double +pio4 = 0x1.921fb54442d1846ap-1L, invpio2 = 6.36619772367581343076e-01L, /* 0xa2f9836e4e44152a.0p-64 */ pio2_1t = -1.07463465549719416346e-12L, /* -0x973dcb3b399d747f.0p-103 */ pio2_2t = 6.36831716351095013979e-25L, /* 0xc51701b839a25205.0p-144 */ @@ -57,6 +58,7 @@ pio2_3t = -2.75299651904407171810e-37L; /* -0xbb5bf6c7ddd660ce.0p-185 */ #define NX 5 #define NY 3 static const long double +pio4 = 0x1.921fb54442d18469898cc51701b8p-1L, invpio2 = 6.3661977236758134307553505349005747e-01L, /* 0x145f306dc9c882a53f84eafa3ea6a.0p-113 */ pio2_1 = 1.5707963267948966192292994253909555e+00L, /* 0x1921fb54442d18469800000000000.0p-112 */ pio2_1t = 2.0222662487959507323996846200947577e-21L, /* 0x13198a2e03707344a4093822299f3.0p-181 */ @@ -76,11 +78,23 @@ int __rem_pio2l(long double x, long double *y) u.f = x; ex = u.i.se & 0x7fff; if (SMALL(u)) { - /* rint(x/(pi/2)), Assume round-to-nearest. */ + /* rint(x/(pi/2)) */ fn = x*invpio2 + toint - toint; n = QUOBITS(fn); r = x-fn*pio2_1; w = fn*pio2_1t; /* 1st round good to 102/180 bits (ld80/ld128) */ + /* Matters with directed rounding. */ + if (predict_false(r - w < -pio4)) { + n--; + fn--; + r = x - fn*pio2_1; + w = fn*pio2_1t; + } else if (predict_false(r - w > pio4)) { + n++; + fn++; + r = x - fn*pio2_1; + w = fn*pio2_1t; + } y[0] = r-w; u.f = y[0]; ey = u.i.se & 0x7fff; |