diff --git a/src/math/__rem_pio2.c b/src/math/__rem_pio2.c index a40db9f..5fafc4a 100644 --- a/src/math/__rem_pio2.c +++ b/src/math/__rem_pio2.c @@ -19,12 +19,6 @@ #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 33 bit of pi/2 @@ -35,7 +29,6 @@ * 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 */ @@ -48,8 +41,8 @@ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ 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]; + double_t z,w,t,r; + double tx[3],ty[2],fn; uint32_t ix; int sign, n, ex, ey, i; @@ -118,7 +111,8 @@ 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. */ - fn = x*invpio2 + toint - toint; + fn = x*invpio2 + 0x1.8p52; + fn = fn - 0x1.8p52; n = (int32_t)fn; r = x - fn*pio2_1; w = fn*pio2_1t; /* 1st round, good to 85 bits */ diff --git a/src/math/__rem_pio2f.c b/src/math/__rem_pio2f.c index f516385..838e1fc 100644 --- a/src/math/__rem_pio2f.c +++ b/src/math/__rem_pio2f.c @@ -22,19 +22,12 @@ #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 */ @@ -42,8 +35,7 @@ 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]; - double_t fn; + double tx[1],ty[1],fn; uint32_t ix; int n, sign, e0; @@ -51,7 +43,8 @@ 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 + toint - toint; + fn = x*invpio2 + 0x1.8p52; + fn = fn - 0x1.8p52; n = (int32_t)fn; *y = x - fn*pio2_1 - fn*pio2_1t; return n; diff --git a/src/math/__rem_pio2l.c b/src/math/__rem_pio2l.c index 77255bd..8b15b7b 100644 --- a/src/math/__rem_pio2l.c +++ b/src/math/__rem_pio2l.c @@ -20,11 +20,10 @@ * use __rem_pio2_large() for large x */ -static const long double toint = 1.5/LDBL_EPSILON; - #if LDBL_MANT_DIG == 64 /* u ~< 0x1p25*pi/2 */ #define SMALL(u) (((u.i.se & 0x7fffU)<<16 | u.i.m>>48) < ((0x3fff + 25)<<16 | 0x921f>>1 | 0x8000)) +#define TOINT 0x1.8p63 #define QUOBITS(x) ((uint32_t)(int32_t)x & 0x7fffffff) #define ROUND1 22 #define ROUND2 61 @@ -51,6 +50,7 @@ pio2_3t = -2.75299651904407171810e-37L; /* -0xbb5bf6c7ddd660ce.0p-185 */ #elif LDBL_MANT_DIG == 113 /* u ~< 0x1p45*pi/2 */ #define SMALL(u) (((u.i.se & 0x7fffU)<<16 | u.i.top) < ((0x3fff + 45)<<16 | 0x921f)) +#define TOINT 0x1.8p112 #define QUOBITS(x) ((uint32_t)(int64_t)x & 0x7fffffff) #define ROUND1 51 #define ROUND2 119 @@ -77,7 +77,7 @@ int __rem_pio2l(long double x, long double *y) ex = u.i.se & 0x7fff; if (SMALL(u)) { /* rint(x/(pi/2)), Assume round-to-nearest. */ - fn = x*invpio2 + toint - toint; + 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) */ diff --git a/src/math/ceil.c b/src/math/ceil.c index b13e6f2..22dc224 100644 --- a/src/math/ceil.c +++ b/src/math/ceil.c @@ -1,12 +1,5 @@ #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 -static const double_t toint = 1/EPS; - double ceil(double x) { union {double f; uint64_t i;} u = {x}; @@ -17,9 +10,9 @@ double ceil(double x) return x; /* y = int(x) - x, where int(x) is an integer neighbor of x */ if (u.i >> 63) - y = x - toint + toint - x; + y = (double)(x - 0x1p52) + 0x1p52 - x; else - y = x + toint - toint - x; + y = (double)(x + 0x1p52) - 0x1p52 - x; /* special case because of non-nearest rounding modes */ if (e <= 0x3ff-1) { FORCE_EVAL(y); diff --git a/src/math/ceill.c b/src/math/ceill.c index 60a8302..a2cb0a7 100644 --- a/src/math/ceill.c +++ b/src/math/ceill.c @@ -6,9 +6,11 @@ long double ceill(long double x) return ceil(x); } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 - -static const long double toint = 1/LDBL_EPSILON; - +#if LDBL_MANT_DIG == 64 +#define TOINT 0x1p63 +#elif LDBL_MANT_DIG == 113 +#define TOINT 0x1p112 +#endif long double ceill(long double x) { union ldshape u = {x}; @@ -19,9 +21,9 @@ long double ceill(long double x) return x; /* y = int(x) - x, where int(x) is an integer neighbor of x */ if (u.i.se >> 15) - y = x - toint + toint - x; + y = x - TOINT + TOINT - x; else - y = x + toint - toint - x; + y = x + TOINT - TOINT - x; /* special case because of non-nearest rounding modes */ if (e <= 0x3fff-1) { FORCE_EVAL(y); diff --git a/src/math/floor.c b/src/math/floor.c index 14a31cd..ebc9fab 100644 --- a/src/math/floor.c +++ b/src/math/floor.c @@ -1,12 +1,5 @@ #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 -static const double_t toint = 1/EPS; - double floor(double x) { union {double f; uint64_t i;} u = {x}; @@ -17,9 +10,9 @@ double floor(double x) return x; /* y = int(x) - x, where int(x) is an integer neighbor of x */ if (u.i >> 63) - y = x - toint + toint - x; + y = (double)(x - 0x1p52) + 0x1p52 - x; else - y = x + toint - toint - x; + y = (double)(x + 0x1p52) - 0x1p52 - x; /* special case because of non-nearest rounding modes */ if (e <= 0x3ff-1) { FORCE_EVAL(y); diff --git a/src/math/floorl.c b/src/math/floorl.c index 16aaec4..961f9e8 100644 --- a/src/math/floorl.c +++ b/src/math/floorl.c @@ -6,9 +6,11 @@ long double floorl(long double x) return floor(x); } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 - -static const long double toint = 1/LDBL_EPSILON; - +#if LDBL_MANT_DIG == 64 +#define TOINT 0x1p63 +#elif LDBL_MANT_DIG == 113 +#define TOINT 0x1p112 +#endif long double floorl(long double x) { union ldshape u = {x}; @@ -19,9 +21,9 @@ long double floorl(long double x) return x; /* y = int(x) - x, where int(x) is an integer neighbor of x */ if (u.i.se >> 15) - y = x - toint + toint - x; + y = x - TOINT + TOINT - x; else - y = x + toint - toint - x; + y = x + TOINT - TOINT - x; /* special case because of non-nearest rounding modes */ if (e <= 0x3fff-1) { FORCE_EVAL(y); diff --git a/src/math/modfl.c b/src/math/modfl.c index a47b192..4b03a4b 100644 --- a/src/math/modfl.c +++ b/src/math/modfl.c @@ -11,9 +11,11 @@ long double modfl(long double x, long double *iptr) return r; } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 - -static const long double toint = 1/LDBL_EPSILON; - +#if LDBL_MANT_DIG == 64 +#define TOINT 0x1p63 +#elif LDBL_MANT_DIG == 113 +#define TOINT 0x1p112 +#endif long double modfl(long double x, long double *iptr) { union ldshape u = {x}; @@ -38,7 +40,7 @@ long double modfl(long double x, long double *iptr) /* raises spurious inexact */ absx = s ? -x : x; - y = absx + toint - toint - absx; + y = absx + TOINT - TOINT - absx; if (y == 0) { *iptr = x; return s ? -0.0 : 0.0; diff --git a/src/math/rintl.c b/src/math/rintl.c index 374327d..2672507 100644 --- a/src/math/rintl.c +++ b/src/math/rintl.c @@ -6,9 +6,11 @@ long double rintl(long double x) return rint(x); } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 - -static const long double toint = 1/LDBL_EPSILON; - +#if LDBL_MANT_DIG == 64 +#define TOINT 0x1p63 +#elif LDBL_MANT_DIG == 113 +#define TOINT 0x1p112 +#endif long double rintl(long double x) { union ldshape u = {x}; @@ -19,9 +21,9 @@ long double rintl(long double x) if (e >= 0x3fff+LDBL_MANT_DIG-1) return x; if (s) - y = x - toint + toint; + y = x - TOINT + TOINT; else - y = x + toint - toint; + y = x + TOINT - TOINT; if (y == 0) return 0*x; return y; diff --git a/src/math/round.c b/src/math/round.c index 130d58d..4b38d1f 100644 --- a/src/math/round.c +++ b/src/math/round.c @@ -1,12 +1,5 @@ #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 -static const double_t toint = 1/EPS; - double round(double x) { union {double f; uint64_t i;} u = {x}; @@ -19,10 +12,10 @@ double round(double x) x = -x; if (e < 0x3ff-1) { /* raise inexact if x!=0 */ - FORCE_EVAL(x + toint); + FORCE_EVAL(x + 0x1p52); return 0*u.f; } - y = x + toint - toint - x; + y = (double)(x + 0x1p52) - 0x1p52 - x; if (y > 0.5) y = y + x - 1; else if (y <= -0.5) diff --git a/src/math/roundf.c b/src/math/roundf.c index e8210af..c6b2779 100644 --- a/src/math/roundf.c +++ b/src/math/roundf.c @@ -1,14 +1,5 @@ #include "libm.h" -#if FLT_EVAL_METHOD==0 -#define EPS FLT_EPSILON -#elif FLT_EVAL_METHOD==1 -#define EPS DBL_EPSILON -#elif FLT_EVAL_METHOD==2 -#define EPS LDBL_EPSILON -#endif -static const float_t toint = 1/EPS; - float roundf(float x) { union {float f; uint32_t i;} u = {x}; @@ -20,10 +11,10 @@ float roundf(float x) if (u.i >> 31) x = -x; if (e < 0x7f-1) { - FORCE_EVAL(x + toint); + FORCE_EVAL(x + 0x1p23f); return 0*u.f; } - y = x + toint - toint - x; + y = (float)(x + 0x1p23f) - 0x1p23f - x; if (y > 0.5f) y = y + x - 1; else if (y <= -0.5f) diff --git a/src/math/roundl.c b/src/math/roundl.c index f4ff682..8f3f28b 100644 --- a/src/math/roundl.c +++ b/src/math/roundl.c @@ -6,9 +6,11 @@ long double roundl(long double x) return round(x); } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 - -static const long double toint = 1/LDBL_EPSILON; - +#if LDBL_MANT_DIG == 64 +#define TOINT 0x1p63 +#elif LDBL_MANT_DIG == 113 +#define TOINT 0x1p112 +#endif long double roundl(long double x) { union ldshape u = {x}; @@ -20,10 +22,10 @@ long double roundl(long double x) if (u.i.se >> 15) x = -x; if (e < 0x3fff-1) { - FORCE_EVAL(x + toint); + FORCE_EVAL(x + TOINT); return 0*u.f; } - y = x + toint - toint - x; + y = x + TOINT - TOINT - x; if (y > 0.5) y = y + x - 1; else if (y <= -0.5) diff --git a/src/math/truncl.c b/src/math/truncl.c index f07b193..3eedb08 100644 --- a/src/math/truncl.c +++ b/src/math/truncl.c @@ -6,9 +6,11 @@ long double truncl(long double x) return trunc(x); } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 - -static const long double toint = 1/LDBL_EPSILON; - +#if LDBL_MANT_DIG == 64 +#define TOINT 0x1p63 +#elif LDBL_MANT_DIG == 113 +#define TOINT 0x1p112 +#endif long double truncl(long double x) { union ldshape u = {x}; @@ -25,7 +27,7 @@ long double truncl(long double x) /* y = int(|x|) - |x|, where int(|x|) is an integer neighbor of |x| */ if (s) x = -x; - y = x + toint - toint - x; + y = x + TOINT - TOINT - x; if (y > 0) y -= 1; x += y;