mirror of
git://git.musl-libc.org/musl
synced 2025-02-16 02:46:53 +00:00
math: cleanup exp2.c exp2f.c and exp2l.c
* old code relied on sign extension on right shift * exp2l ld64 wrapper was wrong * use scalbn instead of bithacks
This commit is contained in:
Szabolcs Nagy
parent
bbbf045ce9
commit
159c7655d0
@ -27,11 +27,9 @@
|
||||
|
||||
#include "libm.h"
|
||||
|
||||
#define TBLBITS 8
|
||||
#define TBLSIZE (1 << TBLBITS)
|
||||
#define TBLSIZE 256
|
||||
|
||||
static const double
|
||||
huge = 0x1p1000,
|
||||
redux = 0x1.8p52 / TBLSIZE,
|
||||
P1 = 0x1.62e42fefa39efp-1,
|
||||
P2 = 0x1.ebfbdff82c575p-3,
|
||||
@ -39,8 +37,6 @@ P3 = 0x1.c6b08d704a0a6p-5,
|
||||
P4 = 0x1.3b2ab88f70400p-7,
|
||||
P5 = 0x1.5d88003875c74p-10;
|
||||
|
||||
static const volatile double twom1000 = 0x1p-1000;
|
||||
|
||||
static const double tbl[TBLSIZE * 2] = {
|
||||
/* exp2(z + eps) eps */
|
||||
0x1.6a09e667f3d5dp-1, 0x1.9880p-44,
|
||||
@ -334,25 +330,28 @@ static const double tbl[TBLSIZE * 2] = {
|
||||
*/
|
||||
double exp2(double x)
|
||||
{
|
||||
double r, t, twopk, twopkp1000, z;
|
||||
uint32_t hx, ix, lx, i0;
|
||||
int k;
|
||||
double r, t, z;
|
||||
uint32_t hx, ix, i0;
|
||||
union {uint32_t u; int32_t i;} k;
|
||||
|
||||
/* Filter out exceptional cases. */
|
||||
GET_HIGH_WORD(hx, x);
|
||||
ix = hx & 0x7fffffff;
|
||||
if (ix >= 0x40900000) { /* |x| >= 1024 */
|
||||
if (ix >= 0x7ff00000) {
|
||||
GET_LOW_WORD(lx, x);
|
||||
if (((ix & 0xfffff) | lx) != 0 || (hx & 0x80000000) == 0)
|
||||
return x + x; /* x is NaN or +Inf */
|
||||
else
|
||||
return 0.0; /* x is -Inf */
|
||||
GET_LOW_WORD(ix, x);
|
||||
if (hx == 0xfff00000 && ix == 0) /* -inf */
|
||||
return 0;
|
||||
return x;
|
||||
}
|
||||
if (x >= 1024) {
|
||||
STRICT_ASSIGN(double, x, x * 0x1p1023);
|
||||
return x;
|
||||
}
|
||||
if (x <= -1075) {
|
||||
STRICT_ASSIGN(double, x, 0x1p-1000*0x1p-1000);
|
||||
return x;
|
||||
}
|
||||
if (x >= 0x1.0p10)
|
||||
return huge * huge; /* overflow */
|
||||
if (x <= -0x1.0ccp10)
|
||||
return twom1000 * twom1000; /* underflow */
|
||||
} else if (ix < 0x3c900000) { /* |x| < 0x1p-54 */
|
||||
return 1.0 + x;
|
||||
}
|
||||
@ -361,24 +360,16 @@ double exp2(double x)
|
||||
STRICT_ASSIGN(double, t, x + redux);
|
||||
GET_LOW_WORD(i0, t);
|
||||
i0 += TBLSIZE / 2;
|
||||
k = (i0 >> TBLBITS) << 20;
|
||||
i0 = (i0 & (TBLSIZE - 1)) << 1;
|
||||
k.u = i0 / TBLSIZE * TBLSIZE;
|
||||
k.i /= TBLSIZE;
|
||||
i0 %= TBLSIZE;
|
||||
t -= redux;
|
||||
z = x - t;
|
||||
|
||||
/* Compute r = exp2(y) = exp2t[i0] * p(z - eps[i]). */
|
||||
t = tbl[i0]; /* exp2t[i0] */
|
||||
z -= tbl[i0 + 1]; /* eps[i0] */
|
||||
if (k >= -1021 << 20)
|
||||
INSERT_WORDS(twopk, 0x3ff00000 + k, 0);
|
||||
else
|
||||
INSERT_WORDS(twopkp1000, 0x3ff00000 + k + (1000 << 20), 0);
|
||||
t = tbl[2*i0]; /* exp2t[i0] */
|
||||
z -= tbl[2*i0 + 1]; /* eps[i0] */
|
||||
r = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * P5))));
|
||||
|
||||
/* Scale by 2**(k>>20). */
|
||||
if (k < -1021 << 20)
|
||||
return r * twopkp1000 * twom1000;
|
||||
if (k == 1024 << 20)
|
||||
return r * 2.0 * 0x1p1023;
|
||||
return r * twopk;
|
||||
return scalbn(r, k.i);
|
||||
}
|
||||
|
||||
@ -27,19 +27,15 @@
|
||||
|
||||
#include "libm.h"
|
||||
|
||||
#define TBLBITS 4
|
||||
#define TBLSIZE (1 << TBLBITS)
|
||||
#define TBLSIZE 16
|
||||
|
||||
static const float
|
||||
huge = 0x1p100f,
|
||||
redux = 0x1.8p23f / TBLSIZE,
|
||||
P1 = 0x1.62e430p-1f,
|
||||
P2 = 0x1.ebfbe0p-3f,
|
||||
P3 = 0x1.c6b348p-5f,
|
||||
P4 = 0x1.3b2c9cp-7f;
|
||||
|
||||
static const volatile float twom100 = 0x1p-100f;
|
||||
|
||||
static const double exp2ft[TBLSIZE] = {
|
||||
0x1.6a09e667f3bcdp-1,
|
||||
0x1.7a11473eb0187p-1,
|
||||
@ -89,23 +85,25 @@ float exp2f(float x)
|
||||
{
|
||||
double tv, twopk, u, z;
|
||||
float t;
|
||||
uint32_t hx, ix, i0;
|
||||
int32_t k;
|
||||
uint32_t hx, ix, i0, k;
|
||||
|
||||
/* Filter out exceptional cases. */
|
||||
GET_FLOAT_WORD(hx, x);
|
||||
ix = hx & 0x7fffffff;
|
||||
if (ix >= 0x43000000) { /* |x| >= 128 */
|
||||
if (ix >= 0x7f800000) {
|
||||
if ((ix & 0x7fffff) != 0 || (hx & 0x80000000) == 0)
|
||||
return x + x; /* x is NaN or +Inf */
|
||||
else
|
||||
return 0.0; /* x is -Inf */
|
||||
if (hx == 0xff800000) /* -inf */
|
||||
return 0;
|
||||
return x;
|
||||
}
|
||||
if (x >= 128) {
|
||||
STRICT_ASSIGN(float, x, x * 0x1p127);
|
||||
return x;
|
||||
}
|
||||
if (x <= -150) {
|
||||
STRICT_ASSIGN(float, x, 0x1p-100*0x1p-100);
|
||||
return x;
|
||||
}
|
||||
if (x >= 0x1.0p7f)
|
||||
return huge * huge; /* overflow */
|
||||
if (x <= -0x1.2cp7f)
|
||||
return twom100 * twom100; /* underflow */
|
||||
} else if (ix <= 0x33000000) { /* |x| <= 0x1p-25 */
|
||||
return 1.0f + x;
|
||||
}
|
||||
@ -114,7 +112,7 @@ float exp2f(float x)
|
||||
STRICT_ASSIGN(float, t, x + redux);
|
||||
GET_FLOAT_WORD(i0, t);
|
||||
i0 += TBLSIZE / 2;
|
||||
k = (i0 >> TBLBITS) << 20;
|
||||
k = (i0 / TBLSIZE) << 20;
|
||||
i0 &= TBLSIZE - 1;
|
||||
t -= redux;
|
||||
z = x - t;
|
||||
|
||||
@ -30,7 +30,7 @@
|
||||
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
|
||||
long double exp2l(long double x)
|
||||
{
|
||||
return exp2l(x);
|
||||
return exp2(x);
|
||||
}
|
||||
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
|
||||
|
||||
@ -40,10 +40,6 @@ long double exp2l(long double x)
|
||||
#define BIAS (LDBL_MAX_EXP - 1)
|
||||
#define EXPMASK (BIAS + LDBL_MAX_EXP)
|
||||
|
||||
static const long double huge = 0x1p10000L;
|
||||
/* XXX Prevent gcc from erroneously constant folding this. */
|
||||
static const volatile long double twom10000 = 0x1p-10000L;
|
||||
|
||||
static const double
|
||||
redux = 0x1.8p63 / TBLSIZE,
|
||||
P1 = 0x1.62e42fefa39efp-1,
|
||||
@ -208,27 +204,28 @@ static const double tbl[TBLSIZE * 2] = {
|
||||
long double exp2l(long double x)
|
||||
{
|
||||
union IEEEl2bits u, v;
|
||||
long double r, twopk, twopkp10000, z;
|
||||
long double r, z;
|
||||
uint32_t hx, ix, i0;
|
||||
int k;
|
||||
union {uint32_t u; int32_t i;} k;
|
||||
|
||||
/* Filter out exceptional cases. */
|
||||
u.e = x;
|
||||
hx = u.xbits.expsign;
|
||||
ix = hx & EXPMASK;
|
||||
if (ix >= BIAS + 14) { /* |x| >= 16384 or x is NaN */
|
||||
if (ix == BIAS + LDBL_MAX_EXP) {
|
||||
if (u.xbits.man != 1ULL << 63 || (hx & 0x8000) == 0)
|
||||
return x + x; /* x is +Inf or NaN */
|
||||
return 0.0; /* x is -Inf */
|
||||
if (ix == EXPMASK) {
|
||||
if (u.xbits.man == 1ULL << 63 && hx == 0xffff) /* -inf */
|
||||
return 0;
|
||||
return x;
|
||||
}
|
||||
if (x >= 16384) {
|
||||
x *= 0x1p16383L;
|
||||
return x;
|
||||
}
|
||||
if (x >= 16384)
|
||||
return huge * huge; /* overflow */
|
||||
if (x <= -16446)
|
||||
return twom10000 * twom10000; /* underflow */
|
||||
} else if (ix <= BIAS - 66) { /* |x| < 0x1p-66 */
|
||||
return 1.0 + x;
|
||||
}
|
||||
return 0x1p-10000L*0x1p-10000L;
|
||||
} else if (ix < BIAS - 64) /* |x| < 0x1p-64 */
|
||||
return 1 + x;
|
||||
|
||||
/*
|
||||
* Reduce x, computing z, i0, and k. The low bits of x + redux
|
||||
@ -240,38 +237,22 @@ long double exp2l(long double x)
|
||||
* Then the low-order word of x + redux is 0x000abc12,
|
||||
* We split this into k = 0xabc and i0 = 0x12 (adjusted to
|
||||
* index into the table), then we compute z = 0x0.003456p0.
|
||||
*
|
||||
* XXX If the exponent is negative, the computation of k depends on
|
||||
* '>>' doing sign extension.
|
||||
*/
|
||||
u.e = x + redux;
|
||||
i0 = u.bits.manl + TBLSIZE / 2;
|
||||
k = (int)i0 >> TBLBITS;
|
||||
i0 = (i0 & (TBLSIZE - 1)) << 1;
|
||||
k.u = i0 / TBLSIZE * TBLSIZE;
|
||||
k.i /= TBLSIZE;
|
||||
i0 %= TBLSIZE;
|
||||
u.e -= redux;
|
||||
z = x - u.e;
|
||||
v.xbits.man = 1ULL << 63;
|
||||
if (k >= LDBL_MIN_EXP) {
|
||||
v.xbits.expsign = LDBL_MAX_EXP - 1 + k;
|
||||
twopk = v.e;
|
||||
} else {
|
||||
v.xbits.expsign = LDBL_MAX_EXP - 1 + k + 10000;
|
||||
twopkp10000 = v.e;
|
||||
}
|
||||
|
||||
/* Compute r = exp2l(y) = exp2lt[i0] * p(z). */
|
||||
long double t_hi = tbl[i0];
|
||||
long double t_lo = tbl[i0 + 1];
|
||||
long double t_hi = tbl[2*i0];
|
||||
long double t_lo = tbl[2*i0 + 1];
|
||||
/* XXX This gives > 1 ulp errors outside of FE_TONEAREST mode */
|
||||
r = t_lo + (t_hi + t_lo) * z * (P1 + z * (P2 + z * (P3 + z * (P4
|
||||
+ z * (P5 + z * P6))))) + t_hi;
|
||||
|
||||
/* Scale by 2**k. */
|
||||
if (k >= LDBL_MIN_EXP) {
|
||||
if (k == LDBL_MAX_EXP)
|
||||
return r * 2.0 * 0x1p16383L;
|
||||
return r * twopk;
|
||||
}
|
||||
return r * twopkp10000 * twom10000;
|
||||
return scalbnl(r, k.i);
|
||||
}
|
||||
#endif
|
||||
|
||||
Loading…
Reference in New Issue
Block a user