double exp(double x)
{
- double hi, lo, c, xx;
+ double_t hi, lo, c, xx, y;
int k, sign;
uint32_t hx;
hx &= 0x7fffffff; /* high word of |x| */
/* special cases */
- if (hx >= 0x40862e42) { /* if |x| >= 709.78... */
+ if (hx >= 0x4086232b) { /* if |x| >= 708.39... */
if (isnan(x))
return x;
- if (hx == 0x7ff00000 && sign) /* -inf */
- return 0;
if (x > 709.782712893383973096) {
/* overflow if x!=inf */
STRICT_ASSIGN(double, x, 0x1p1023 * x);
return x;
}
- if (x < -745.13321910194110842) {
- /* underflow */
- STRICT_ASSIGN(double, x, 0x1p-1000 * 0x1p-1000);
- return x;
+ if (x < -708.39641853226410622) {
+ /* underflow if x!=-inf */
+ FORCE_EVAL((float)(-0x1p-149/x));
+ if (x < -745.13321910194110842)
+ return 0;
}
}
/* x is now in primary range */
xx = x*x;
c = x - xx*(P1+xx*(P2+xx*(P3+xx*(P4+xx*P5))));
- x = 1 + (x*c/(2-c) - lo + hi);
+ y = 1 + (x*c/(2-c) - lo + hi);
if (k == 0)
- return x;
- return scalbn(x, k);
+ return y;
+ return scalbn(y, k);
}
* Method: (accurate tables)
*
* Reduce x:
- * x = 2**k + y, for integer k and |y| <= 1/2.
+ * x = k + y, for integer k and |y| <= 1/2.
* Thus we have exp2(x) = 2**k * exp2(y).
*
* Reduce y:
*/
double exp2(double x)
{
- double r, t, z;
- uint32_t hx, ix, i0;
+ double_t r, t, z;
+ uint32_t ix, i0;
+ union {double f; uint64_t i;} u = {x};
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(ix, x);
- if (hx == 0xfff00000 && ix == 0) /* -inf */
- return 0;
- return x;
- }
- if (x >= 1024) {
+ ix = u.i>>32 & 0x7fffffff;
+ if (ix >= 0x408ff000) { /* |x| >= 1022 or nan */
+ if (ix >= 0x40900000 && u.i>>63 == 0) { /* x >= 1024 or nan */
+ /* overflow */
STRICT_ASSIGN(double, x, x * 0x1p1023);
return x;
}
- if (x <= -1075) {
- STRICT_ASSIGN(double, x, 0x1p-1000*0x1p-1000);
- return x;
+ if (ix >= 0x7ff00000) /* -inf or -nan */
+ return -1/x;
+ if (u.i>>63) { /* x <= -1022 */
+ /* underflow */
+ if (x <= -1075 || x - 0x1p52 + 0x1p52 != x)
+ FORCE_EVAL((float)(-0x1p-149/x));
+ if (x <= -1075)
+ return 0;
}
} else if (ix < 0x3c900000) { /* |x| < 0x1p-54 */
return 1.0 + x;
}
/* Reduce x, computing z, i0, and k. */
- STRICT_ASSIGN(double, t, x + redux);
- GET_LOW_WORD(i0, t);
+ u.f = x + redux;
+ i0 = u.i;
i0 += TBLSIZE / 2;
k.u = i0 / TBLSIZE * TBLSIZE;
k.i /= TBLSIZE;
i0 %= TBLSIZE;
- t -= redux;
- z = x - t;
+ u.f -= redux;
+ z = x - u.f;
/* Compute r = exp2(y) = exp2t[i0] * p(z - eps[i]). */
t = tbl[2*i0]; /* exp2t[i0] */
* Method: (equally-spaced tables)
*
* Reduce x:
- * x = 2**k + y, for integer k and |y| <= 1/2.
+ * x = k + y, for integer k and |y| <= 1/2.
* Thus we have exp2f(x) = 2**k * exp2(y).
*
* Reduce y:
*/
float exp2f(float x)
{
- double tv, twopk, u, z;
- float t;
- uint32_t hx, ix, i0, k;
+ double_t t, r, z;
+ union {float f; uint32_t i;} u = {x};
+ union {double f; uint64_t i;} uk;
+ uint32_t ix, i0, k;
/* Filter out exceptional cases. */
- GET_FLOAT_WORD(hx, x);
- ix = hx & 0x7fffffff;
- if (ix >= 0x43000000) { /* |x| >= 128 */
- if (ix >= 0x7f800000) {
- if (hx == 0xff800000) /* -inf */
- return 0;
- return x;
- }
- if (x >= 128) {
+ ix = u.i & 0x7fffffff;
+ if (ix > 0x42fc0000) { /* |x| > 126 */
+ if (u.i >= 0x43000000 && u.i < 0x80000000) { /* x >= 128 */
STRICT_ASSIGN(float, x, x * 0x1p127f);
return x;
}
- if (x <= -150) {
- STRICT_ASSIGN(float, x, 0x1p-100f*0x1p-100f);
- return x;
+ if (u.i >= 0x80000000) { /* x < -126 */
+ if (u.i >= 0xc3160000 || (u.i & 0x0000ffff))
+ FORCE_EVAL(-0x1p-149f/x);
+ if (u.i >= 0xc3160000) /* x <= -150 */
+ return 0;
}
} else if (ix <= 0x33000000) { /* |x| <= 0x1p-25 */
return 1.0f + x;
}
/* Reduce x, computing z, i0, and k. */
- STRICT_ASSIGN(float, t, x + redux);
- GET_FLOAT_WORD(i0, t);
+ u.f = x + redux;
+ i0 = u.i;
i0 += TBLSIZE / 2;
- k = (i0 / TBLSIZE) << 20;
+ k = i0 / TBLSIZE;
+ uk.i = (uint64_t)(0x3ff + k)<<52;
i0 &= TBLSIZE - 1;
- t -= redux;
- z = x - t;
- INSERT_WORDS(twopk, 0x3ff00000 + k, 0);
-
+ u.f -= redux;
+ z = x - u.f;
/* Compute r = exp2(y) = exp2ft[i0] * p(z). */
- tv = exp2ft[i0];
- u = tv * z;
- tv = tv + u * (P1 + z * P2) + u * (z * z) * (P3 + z * P4);
+ r = exp2ft[i0];
+ t = r * z;
+ r = r + t * (P1 + z * P2) + t * (z * z) * (P3 + z * P4);
- /* Scale by 2**(k>>20). */
- return tv * twopk;
+ /* Scale by 2**k */
+ return r * uk.f;
}
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,
*/
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
* 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];
float expf(float x)
{
- float hi, lo, c, xx;
+ float_t hi, lo, c, xx, y;
int k, sign;
uint32_t hx;
hx &= 0x7fffffff; /* high word of |x| */
/* special cases */
- if (hx >= 0x42b17218) { /* if |x| >= 88.722839f or NaN */
- if (hx > 0x7f800000) /* NaN */
- return x;
- if (!sign) {
- /* overflow if x!=inf */
+ if (hx >= 0x42aeac50) { /* if |x| >= -87.33655f or NaN */
+ if (hx >= 0x42b17218 && !sign) { /* x >= 88.722839f */
+ /* overflow */
STRICT_ASSIGN(float, x, x * 0x1p127f);
return x;
}
- if (hx == 0x7f800000) /* -inf */
- return 0;
- if (hx >= 0x42cff1b5) { /* x <= -103.972084f */
+ if (sign) {
/* underflow */
- STRICT_ASSIGN(float, x, 0x1p-100f*0x1p-100f);
- return x;
+ FORCE_EVAL(-0x1p-149f/x);
+ if (hx >= 0x42cff1b5) /* x <= -103.972084f */
+ return 0;
}
}
/* x is now in primary range */
xx = x*x;
c = x - xx*(P1+xx*P2);
- x = 1 + (x*c/(2-c) - lo + hi);
+ y = 1 + (x*c/(2-c) - lo + hi);
if (k == 0)
- return x;
- return scalbnf(x, k);
+ return y;
+ return scalbnf(y, k);
}
if (x > 11356.5234062941439488L) /* x > ln(2^16384 - 0.5) */
return x * 0x1p16383L;
if (x < -11399.4985314888605581L) /* x < ln(2^-16446) */
- return 0x1p-10000L * 0x1p-10000L;
+ return -0x1p-16445L/x;
/* Express e**x = e**f 2**k
* = e**(f + k ln(2))
* R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
* = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r))
* = 1 - r^2/60 + r^4/2520 - r^6/100800 + ...
- * We use a special Reme algorithm on [0,0.347] to generate
+ * We use a special Remez algorithm on [0,0.347] to generate
* a polynomial of degree 5 in r*r to approximate R1. The
* maximum error of this polynomial approximation is bounded
* by 2**-61. In other words,
#include "libm.h"
static const double
-huge = 1.0e+300,
-tiny = 1.0e-300,
o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
ln2_hi = 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
ln2_lo = 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
double expm1(double x)
{
- double y,hi,lo,c,t,e,hxs,hfx,r1,twopk;
- int32_t k,xsb;
- uint32_t hx;
-
- GET_HIGH_WORD(hx, x);
- xsb = hx&0x80000000; /* sign bit of x */
- hx &= 0x7fffffff; /* high word of |x| */
+ double_t y,hi,lo,c,t,e,hxs,hfx,r1,twopk;
+ union {double f; uint64_t i;} u = {x};
+ uint32_t hx = u.i>>32 & 0x7fffffff;
+ int k, sign = u.i>>63;
/* filter out huge and non-finite argument */
if (hx >= 0x4043687A) { /* if |x|>=56*ln2 */
- if (hx >= 0x40862E42) { /* if |x|>=709.78... */
- if (hx >= 0x7ff00000) {
- uint32_t low;
-
- GET_LOW_WORD(low, x);
- if (((hx&0xfffff)|low) != 0) /* NaN */
- return x+x;
- return xsb==0 ? x : -1.0; /* exp(+-inf)={inf,-1} */
- }
- if(x > o_threshold)
- return huge*huge; /* overflow */
- }
- if (xsb != 0) { /* x < -56*ln2, return -1.0 with inexact */
- /* raise inexact */
- if(x+tiny<0.0)
- return tiny-1.0; /* return -1 */
+ if (isnan(x))
+ return x;
+ if (sign)
+ return -1;
+ if (x > o_threshold) {
+ x *= 0x1p1023;
+ return x;
}
}
/* argument reduction */
if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
- if (xsb == 0) {
+ if (!sign) {
hi = x - ln2_hi;
lo = ln2_lo;
k = 1;
k = -1;
}
} else {
- k = invln2*x + (xsb==0 ? 0.5 : -0.5);
+ k = invln2*x + (sign ? -0.5 : 0.5);
t = k;
hi = x - t*ln2_hi; /* t*ln2_hi is exact here */
lo = t*ln2_lo;
STRICT_ASSIGN(double, x, hi - lo);
c = (hi-x)-lo;
} else if (hx < 0x3c900000) { /* |x| < 2**-54, return x */
- /* raise inexact flags when x != 0 */
- t = huge+x;
- return x - (t-(huge+x));
+ if (hx < 0x00100000)
+ FORCE_EVAL((float)x);
+ return x;
} else
k = 0;
e = hxs*((r1-t)/(6.0 - x*t));
if (k == 0) /* c is 0 */
return x - (x*e-hxs);
- INSERT_WORDS(twopk, 0x3ff00000+(k<<20), 0); /* 2^k */
e = x*(e-c) - c;
e -= hxs;
+ /* exp(x) ~ 2^k (x_reduced - e + 1) */
if (k == -1)
return 0.5*(x-e) - 0.5;
if (k == 1) {
return -2.0*(e-(x+0.5));
return 1.0+2.0*(x-e);
}
- if (k <= -2 || k > 56) { /* suffice to return exp(x)-1 */
- y = 1.0 - (e-x);
+ u.i = (uint64_t)(0x3ff + k)<<52; /* 2^k */
+ twopk = u.f;
+ if (k < 0 || k > 56) { /* suffice to return exp(x)-1 */
+ y = x - e + 1.0;
if (k == 1024)
y = y*2.0*0x1p1023;
else
y = y*twopk;
return y - 1.0;
}
- t = 1.0;
- if (k < 20) {
- SET_HIGH_WORD(t, 0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */
- y = t-(e-x);
- y = y*twopk;
- } else {
- SET_HIGH_WORD(t, ((0x3ff-k)<<20)); /* 2^-k */
- y = x-(e+t);
- y += 1.0;
- y = y*twopk;
- }
+ u.i = (uint64_t)(0x3ff - k)<<52; /* 2^-k */
+ if (k < 20)
+ y = (x-e+(1-u.f))*twopk;
+ else
+ y = (x-(e+u.f)+1)*twopk;
return y;
}
#include "libm.h"
static const float
-huge = 1.0e+30,
-tiny = 1.0e-30,
o_threshold = 8.8721679688e+01, /* 0x42b17180 */
ln2_hi = 6.9313812256e-01, /* 0x3f317180 */
ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */
float expm1f(float x)
{
- float y,hi,lo,c,t,e,hxs,hfx,r1,twopk;
- int32_t k,xsb;
- uint32_t hx;
-
- GET_FLOAT_WORD(hx, x);
- xsb = hx&0x80000000; /* sign bit of x */
- hx &= 0x7fffffff; /* high word of |x| */
+ float_t y,hi,lo,c,t,e,hxs,hfx,r1,twopk;
+ union {float f; uint32_t i;} u = {x};
+ uint32_t hx = u.i & 0x7fffffff;
+ int k, sign = u.i >> 31;
/* filter out huge and non-finite argument */
if (hx >= 0x4195b844) { /* if |x|>=27*ln2 */
- if (hx >= 0x42b17218) { /* if |x|>=88.721... */
- if (hx > 0x7f800000) /* NaN */
- return x+x;
- if (hx == 0x7f800000) /* exp(+-inf)={inf,-1} */
- return xsb==0 ? x : -1.0;
- if (x > o_threshold)
- return huge*huge; /* overflow */
- }
- if (xsb != 0) { /* x < -27*ln2 */
- /* raise inexact */
- if (x+tiny < 0.0f)
- return tiny-1.0f; /* return -1 */
+ if (hx > 0x7f800000) /* NaN */
+ return x;
+ if (sign)
+ return -1;
+ if (x > o_threshold) {
+ x *= 0x1p127f;
+ return x;
}
}
/* argument reduction */
if (hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */
if (hx < 0x3F851592) { /* and |x| < 1.5 ln2 */
- if (xsb == 0) {
+ if (!sign) {
hi = x - ln2_hi;
lo = ln2_lo;
k = 1;
k = -1;
}
} else {
- k = invln2*x + (xsb==0 ? 0.5f : -0.5f);
+ k = invln2*x + (sign ? -0.5f : 0.5f);
t = k;
hi = x - t*ln2_hi; /* t*ln2_hi is exact here */
lo = t*ln2_lo;
STRICT_ASSIGN(float, x, hi - lo);
c = (hi-x)-lo;
} else if (hx < 0x33000000) { /* when |x|<2**-25, return x */
- t = huge+x; /* return x with inexact flags when x!=0 */
- return x - (t-(huge+x));
+ if (hx < 0x00800000)
+ FORCE_EVAL(x*x);
+ return x;
} else
k = 0;
e = hxs*((r1-t)/(6.0f - x*t));
if (k == 0) /* c is 0 */
return x - (x*e-hxs);
- SET_FLOAT_WORD(twopk, 0x3f800000+(k<<23)); /* 2^k */
e = x*(e-c) - c;
e -= hxs;
+ /* exp(x) ~ 2^k (x_reduced - e + 1) */
if (k == -1)
return 0.5f*(x-e) - 0.5f;
if (k == 1) {
return -2.0f*(e-(x+0.5f));
return 1.0f + 2.0f*(x-e);
}
- if (k <= -2 || k > 56) { /* suffice to return exp(x)-1 */
- y = 1.0f - (e - x);
+ u.i = (0x7f+k)<<23; /* 2^k */
+ twopk = u.f;
+ if (k < 0 || k > 56) { /* suffice to return exp(x)-1 */
+ y = x - e + 1.0f;
if (k == 128)
y = y*2.0f*0x1p127f;
else
y = y*twopk;
return y - 1.0f;
}
- t = 1.0f;
- if (k < 23) {
- SET_FLOAT_WORD(t, 0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */
- y = t - (e - x);
- y = y*twopk;
- } else {
- SET_FLOAT_WORD(t, (0x7f-k)<<23); /* 2^-k */
- y = x - (e + t);
- y += 1.0f;
- y = y*twopk;
- }
+ u.i = (0x7f-k)<<23; /* 2^-k */
+ if (k < 23)
+ y = (x-e+(1-u.f))*twopk;
+ else
+ y = (x-(e+u.f)+1)*twopk;
return y;
}
projects / oweals / musl.git / commitdiff