summaryrefslogtreecommitdiff
path: root/src/math/exp2l.c
diff options
context:
space:
mode:
authorSzabolcs Nagy <nsz@port70.net>2013-09-04 07:51:11 +0000
committerSzabolcs Nagy <nsz@port70.net>2013-09-05 11:30:08 +0000
commit39c910fb061114e6aa5c3bf2c94b1d7262d62221 (patch)
tree5ca9b82746e8ac224f93dcda1b8d19f95b68eda2 /src/math/exp2l.c
parentea9bb95a5b36c0a3d2ed8fb03808745b406c2633 (diff)
downloadmusl-39c910fb061114e6aa5c3bf2c94b1d7262d62221.tar.gz
musl-39c910fb061114e6aa5c3bf2c94b1d7262d62221.tar.bz2
musl-39c910fb061114e6aa5c3bf2c94b1d7262d62221.tar.xz
musl-39c910fb061114e6aa5c3bf2c94b1d7262d62221.zip
math: fix underflow in exp*.c and long double handling in exp2l
* don't care about inexact flag * use double_t and float_t (faster, smaller, more precise on x86) * exp: underflow when result is zero or subnormal and not -inf * exp2: underflow when result is zero or subnormal and not exact * expm1: underflow when result is zero or subnormal * expl: don't underflow on -inf * exp2: fix incorrect comment * expm1: simplify special case handling and overflow properly * expm1: cleanup final scaling and fix negative left shift ub (twopk)
Diffstat (limited to 'src/math/exp2l.c')
-rw-r--r--src/math/exp2l.c44
1 files changed, 20 insertions, 24 deletions
diff --git a/src/math/exp2l.c b/src/math/exp2l.c
index 145c20fe..8fc4037c 100644
--- a/src/math/exp2l.c
+++ b/src/math/exp2l.c
@@ -33,13 +33,9 @@ long double exp2l(long double x)
return exp2(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
-
#define TBLBITS 7
#define TBLSIZE (1 << TBLBITS)
-#define BIAS (LDBL_MAX_EXP - 1)
-#define EXPMASK (BIAS + LDBL_MAX_EXP)
-
static const double
redux = 0x1.8p63 / TBLSIZE,
P1 = 0x1.62e42fefa39efp-1,
@@ -203,29 +199,29 @@ static const double tbl[TBLSIZE * 2] = {
*/
long double exp2l(long double x)
{
- union IEEEl2bits u, v;
+ union ldshape u = {x};
+ int e = u.i.se & 0x7fff;
long double r, z;
- uint32_t hx, ix, i0;
+ uint32_t i0;
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 == EXPMASK) {
- if (u.xbits.man == 1ULL << 63 && hx == 0xffff) /* -inf */
+ if (e >= 0x3fff + 13) { /* |x| >= 8192 or x is NaN */
+ if (u.i.se >= 0x3fff + 14 && u.i.se >> 15 == 0)
+ /* overflow */
+ return x * 0x1p16383L;
+ if (e == 0x7fff) /* -inf or -nan */
+ return -1/x;
+ if (x < -16382) {
+ if (x <= -16446 || x - 0x1p63 + 0x1p63 != x)
+ /* underflow */
+ FORCE_EVAL((float)(-0x1p-149/x));
+ if (x <= -16446)
return 0;
- return x;
- }
- if (x >= 16384) {
- x *= 0x1p16383L;
- return x;
}
- if (x <= -16446)
- return 0x1p-10000L*0x1p-10000L;
- } else if (ix < BIAS - 64) /* |x| < 0x1p-64 */
+ } else if (e < 0x3fff - 64) {
return 1 + x;
+ }
/*
* Reduce x, computing z, i0, and k. The low bits of x + redux
@@ -238,13 +234,13 @@ long double exp2l(long double x)
* We split this into k = 0xabc and i0 = 0x12 (adjusted to
* index into the table), then we compute z = 0x0.003456p0.
*/
- u.e = x + redux;
- i0 = u.bits.manl + TBLSIZE / 2;
+ u.f = x + redux;
+ i0 = u.i.m + TBLSIZE / 2;
k.u = i0 / TBLSIZE * TBLSIZE;
k.i /= TBLSIZE;
i0 %= TBLSIZE;
- u.e -= redux;
- z = x - u.e;
+ u.f -= redux;
+ z = x - u.f;
/* Compute r = exp2l(y) = exp2lt[i0] * p(z). */
long double t_hi = tbl[2*i0];