| Index: fusl/src/math/__rem_pio2.c
|
| diff --git a/fusl/src/math/__rem_pio2.c b/fusl/src/math/__rem_pio2.c
|
| index d403f81c792217a195d9680d9356eec17f734f41..025ce4e3799be149bfc917c15aa9973aef83c2f7 100644
|
| --- a/fusl/src/math/__rem_pio2.c
|
| +++ b/fusl/src/math/__rem_pio2.c
|
| @@ -19,9 +19,9 @@
|
|
|
| #include "libm.h"
|
|
|
| -#if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1
|
| +#if FLT_EVAL_METHOD == 0 || FLT_EVAL_METHOD == 1
|
| #define EPS DBL_EPSILON
|
| -#elif FLT_EVAL_METHOD==2
|
| +#elif FLT_EVAL_METHOD == 2
|
| #define EPS LDBL_EPSILON
|
| #endif
|
|
|
| @@ -34,144 +34,146 @@
|
| * pio2_3: third 33 bit of pi/2
|
| * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
|
| */
|
| -static const double
|
| -toint = 1.5/EPS,
|
| -invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
|
| -pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
|
| -pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
|
| -pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
|
| -pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
|
| -pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
|
| -pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
|
| +static const double toint = 1.5 / EPS,
|
| + invpio2 =
|
| + 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
|
| + pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
|
| + pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
|
| + pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
|
| + pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
|
| + pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
|
| + pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
|
|
|
| /* caller must handle the case when reduction is not needed: |x| ~<= pi/4 */
|
| -int __rem_pio2(double x, double *y)
|
| -{
|
| - union {double f; uint64_t i;} u = {x};
|
| - double_t z,w,t,r,fn;
|
| - double tx[3],ty[2];
|
| - uint32_t ix;
|
| - int sign, n, ex, ey, i;
|
| +int __rem_pio2(double x, double* y) {
|
| + union {
|
| + double f;
|
| + uint64_t i;
|
| + } u = {x};
|
| + double_t z, w, t, r, fn;
|
| + double tx[3], ty[2];
|
| + uint32_t ix;
|
| + int sign, n, ex, ey, i;
|
|
|
| - sign = u.i>>63;
|
| - ix = u.i>>32 & 0x7fffffff;
|
| - if (ix <= 0x400f6a7a) { /* |x| ~<= 5pi/4 */
|
| - if ((ix & 0xfffff) == 0x921fb) /* |x| ~= pi/2 or 2pi/2 */
|
| - goto medium; /* cancellation -- use medium case */
|
| - if (ix <= 0x4002d97c) { /* |x| ~<= 3pi/4 */
|
| - if (!sign) {
|
| - z = x - pio2_1; /* one round good to 85 bits */
|
| - y[0] = z - pio2_1t;
|
| - y[1] = (z-y[0]) - pio2_1t;
|
| - return 1;
|
| - } else {
|
| - z = x + pio2_1;
|
| - y[0] = z + pio2_1t;
|
| - y[1] = (z-y[0]) + pio2_1t;
|
| - return -1;
|
| - }
|
| - } else {
|
| - if (!sign) {
|
| - z = x - 2*pio2_1;
|
| - y[0] = z - 2*pio2_1t;
|
| - y[1] = (z-y[0]) - 2*pio2_1t;
|
| - return 2;
|
| - } else {
|
| - z = x + 2*pio2_1;
|
| - y[0] = z + 2*pio2_1t;
|
| - y[1] = (z-y[0]) + 2*pio2_1t;
|
| - return -2;
|
| - }
|
| - }
|
| - }
|
| - if (ix <= 0x401c463b) { /* |x| ~<= 9pi/4 */
|
| - if (ix <= 0x4015fdbc) { /* |x| ~<= 7pi/4 */
|
| - if (ix == 0x4012d97c) /* |x| ~= 3pi/2 */
|
| - goto medium;
|
| - if (!sign) {
|
| - z = x - 3*pio2_1;
|
| - y[0] = z - 3*pio2_1t;
|
| - y[1] = (z-y[0]) - 3*pio2_1t;
|
| - return 3;
|
| - } else {
|
| - z = x + 3*pio2_1;
|
| - y[0] = z + 3*pio2_1t;
|
| - y[1] = (z-y[0]) + 3*pio2_1t;
|
| - return -3;
|
| - }
|
| - } else {
|
| - if (ix == 0x401921fb) /* |x| ~= 4pi/2 */
|
| - goto medium;
|
| - if (!sign) {
|
| - z = x - 4*pio2_1;
|
| - y[0] = z - 4*pio2_1t;
|
| - y[1] = (z-y[0]) - 4*pio2_1t;
|
| - return 4;
|
| - } else {
|
| - z = x + 4*pio2_1;
|
| - y[0] = z + 4*pio2_1t;
|
| - y[1] = (z-y[0]) + 4*pio2_1t;
|
| - return -4;
|
| - }
|
| - }
|
| - }
|
| - if (ix < 0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */
|
| -medium:
|
| - /* rint(x/(pi/2)), Assume round-to-nearest. */
|
| - 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 */
|
| - y[0] = r - w;
|
| - u.f = y[0];
|
| - ey = u.i>>52 & 0x7ff;
|
| - ex = ix>>20;
|
| - if (ex - ey > 16) { /* 2nd round, good to 118 bits */
|
| - t = r;
|
| - w = fn*pio2_2;
|
| - r = t - w;
|
| - w = fn*pio2_2t - ((t-r)-w);
|
| - y[0] = r - w;
|
| - u.f = y[0];
|
| - ey = u.i>>52 & 0x7ff;
|
| - if (ex - ey > 49) { /* 3rd round, good to 151 bits, covers all cases */
|
| - t = r;
|
| - w = fn*pio2_3;
|
| - r = t - w;
|
| - w = fn*pio2_3t - ((t-r)-w);
|
| - y[0] = r - w;
|
| - }
|
| - }
|
| - y[1] = (r - y[0]) - w;
|
| - return n;
|
| - }
|
| - /*
|
| - * all other (large) arguments
|
| - */
|
| - if (ix >= 0x7ff00000) { /* x is inf or NaN */
|
| - y[0] = y[1] = x - x;
|
| - return 0;
|
| - }
|
| - /* set z = scalbn(|x|,-ilogb(x)+23) */
|
| - u.f = x;
|
| - u.i &= (uint64_t)-1>>12;
|
| - u.i |= (uint64_t)(0x3ff + 23)<<52;
|
| - z = u.f;
|
| - for (i=0; i < 2; i++) {
|
| - tx[i] = (double)(int32_t)z;
|
| - z = (z-tx[i])*0x1p24;
|
| - }
|
| - tx[i] = z;
|
| - /* skip zero terms, first term is non-zero */
|
| - while (tx[i] == 0.0)
|
| - i--;
|
| - n = __rem_pio2_large(tx,ty,(int)(ix>>20)-(0x3ff+23),i+1,1);
|
| - if (sign) {
|
| - y[0] = -ty[0];
|
| - y[1] = -ty[1];
|
| - return -n;
|
| - }
|
| - y[0] = ty[0];
|
| - y[1] = ty[1];
|
| - return n;
|
| + sign = u.i >> 63;
|
| + ix = u.i >> 32 & 0x7fffffff;
|
| + if (ix <= 0x400f6a7a) { /* |x| ~<= 5pi/4 */
|
| + if ((ix & 0xfffff) == 0x921fb) /* |x| ~= pi/2 or 2pi/2 */
|
| + goto medium; /* cancellation -- use medium case */
|
| + if (ix <= 0x4002d97c) { /* |x| ~<= 3pi/4 */
|
| + if (!sign) {
|
| + z = x - pio2_1; /* one round good to 85 bits */
|
| + y[0] = z - pio2_1t;
|
| + y[1] = (z - y[0]) - pio2_1t;
|
| + return 1;
|
| + } else {
|
| + z = x + pio2_1;
|
| + y[0] = z + pio2_1t;
|
| + y[1] = (z - y[0]) + pio2_1t;
|
| + return -1;
|
| + }
|
| + } else {
|
| + if (!sign) {
|
| + z = x - 2 * pio2_1;
|
| + y[0] = z - 2 * pio2_1t;
|
| + y[1] = (z - y[0]) - 2 * pio2_1t;
|
| + return 2;
|
| + } else {
|
| + z = x + 2 * pio2_1;
|
| + y[0] = z + 2 * pio2_1t;
|
| + y[1] = (z - y[0]) + 2 * pio2_1t;
|
| + return -2;
|
| + }
|
| + }
|
| + }
|
| + if (ix <= 0x401c463b) { /* |x| ~<= 9pi/4 */
|
| + if (ix <= 0x4015fdbc) { /* |x| ~<= 7pi/4 */
|
| + if (ix == 0x4012d97c) /* |x| ~= 3pi/2 */
|
| + goto medium;
|
| + if (!sign) {
|
| + z = x - 3 * pio2_1;
|
| + y[0] = z - 3 * pio2_1t;
|
| + y[1] = (z - y[0]) - 3 * pio2_1t;
|
| + return 3;
|
| + } else {
|
| + z = x + 3 * pio2_1;
|
| + y[0] = z + 3 * pio2_1t;
|
| + y[1] = (z - y[0]) + 3 * pio2_1t;
|
| + return -3;
|
| + }
|
| + } else {
|
| + if (ix == 0x401921fb) /* |x| ~= 4pi/2 */
|
| + goto medium;
|
| + if (!sign) {
|
| + z = x - 4 * pio2_1;
|
| + y[0] = z - 4 * pio2_1t;
|
| + y[1] = (z - y[0]) - 4 * pio2_1t;
|
| + return 4;
|
| + } else {
|
| + z = x + 4 * pio2_1;
|
| + y[0] = z + 4 * pio2_1t;
|
| + y[1] = (z - y[0]) + 4 * pio2_1t;
|
| + return -4;
|
| + }
|
| + }
|
| + }
|
| + if (ix < 0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */
|
| + medium:
|
| + /* rint(x/(pi/2)), Assume round-to-nearest. */
|
| + 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 */
|
| + y[0] = r - w;
|
| + u.f = y[0];
|
| + ey = u.i >> 52 & 0x7ff;
|
| + ex = ix >> 20;
|
| + if (ex - ey > 16) { /* 2nd round, good to 118 bits */
|
| + t = r;
|
| + w = fn * pio2_2;
|
| + r = t - w;
|
| + w = fn * pio2_2t - ((t - r) - w);
|
| + y[0] = r - w;
|
| + u.f = y[0];
|
| + ey = u.i >> 52 & 0x7ff;
|
| + if (ex - ey > 49) { /* 3rd round, good to 151 bits, covers all cases */
|
| + t = r;
|
| + w = fn * pio2_3;
|
| + r = t - w;
|
| + w = fn * pio2_3t - ((t - r) - w);
|
| + y[0] = r - w;
|
| + }
|
| + }
|
| + y[1] = (r - y[0]) - w;
|
| + return n;
|
| + }
|
| + /*
|
| + * all other (large) arguments
|
| + */
|
| + if (ix >= 0x7ff00000) { /* x is inf or NaN */
|
| + y[0] = y[1] = x - x;
|
| + return 0;
|
| + }
|
| + /* set z = scalbn(|x|,-ilogb(x)+23) */
|
| + u.f = x;
|
| + u.i &= (uint64_t)-1 >> 12;
|
| + u.i |= (uint64_t)(0x3ff + 23) << 52;
|
| + z = u.f;
|
| + for (i = 0; i < 2; i++) {
|
| + tx[i] = (double)(int32_t)z;
|
| + z = (z - tx[i]) * 0x1p24;
|
| + }
|
| + tx[i] = z;
|
| + /* skip zero terms, first term is non-zero */
|
| + while (tx[i] == 0.0)
|
| + i--;
|
| + n = __rem_pio2_large(tx, ty, (int)(ix >> 20) - (0x3ff + 23), i + 1, 1);
|
| + if (sign) {
|
| + y[0] = -ty[0];
|
| + y[1] = -ty[1];
|
| + return -n;
|
| + }
|
| + y[0] = ty[0];
|
| + y[1] = ty[1];
|
| + return n;
|
| }
|
|
|