#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,
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,
*/
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;
}
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);
}
#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,
{
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;
}
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;
#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
#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,
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
* 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
projects / oweals / musl.git / commitdiff