diff options
| author | Szabolcs Nagy <nsz@port70.net> | 2014-10-29 00:34:37 +0100 |
|---|---|---|
| committer | Rich Felker <dalias@aerifal.cx> | 2014-10-31 11:35:40 -0400 |
| commit | 0ce946cf808274c2d6e5419b139e130c8ad4bd30 (patch) | |
| tree | e6614e756dde0afbcd48e38916d0208fed93ece1 /src/math/__rem_pio2.c | |
| parent | 79ca86094d70f43252b683c3a3ccb572d462cf28 (diff) | |
| download | musl-0ce946cf808274c2d6e5419b139e130c8ad4bd30.tar.gz | |
math: use the rounding idiom consistently
the idiomatic rounding of x is
n = x + toint - toint;
where toint is either 1/EPSILON (x is non-negative) or 1.5/EPSILON
(x may be negative and nearest rounding mode is assumed) and EPSILON is
according to the evaluation precision (the type of toint is not very
important, because single precision float can represent the 1/EPSILON of
ieee binary128).
in case of FLT_EVAL_METHOD!=0 this avoids a useless store to double or
float precision, and the long double code became cleaner with
1/LDBL_EPSILON instead of ifdefs for toint.
__rem_pio2f and __rem_pio2 functions slightly changed semantics:
on i386 a double-rounding is avoided so close to half-way cases may
get evaluated differently eg. as sin(pi/4-eps) instead of cos(pi/4+eps)
Diffstat (limited to 'src/math/__rem_pio2.c')
| -rw-r--r-- | src/math/__rem_pio2.c | 14 |
1 files changed, 10 insertions, 4 deletions
diff --git a/src/math/__rem_pio2.c b/src/math/__rem_pio2.c index 5fafc4a8..a40db9fc 100644 --- a/src/math/__rem_pio2.c +++ b/src/math/__rem_pio2.c @@ -19,6 +19,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 33 bit of pi/2 @@ -29,6 +35,7 @@ * 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 */ @@ -41,8 +48,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; - double tx[3],ty[2],fn; + double_t z,w,t,r,fn; + double tx[3],ty[2]; uint32_t ix; int sign, n, ex, ey, i; @@ -111,8 +118,7 @@ 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 + 0x1.8p52; - fn = fn - 0x1.8p52; + fn = x*invpio2 + toint - toint; n = (int32_t)fn; r = x - fn*pio2_1; w = fn*pio2_1t; /* 1st round, good to 85 bits */ |