* the interval [0,0.34658]:
* Write
* R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
- * We use a special Remes algorithm on [0,0.34658] to generate
+ * We use a special Remez algorithm on [0,0.34658] to generate
* a polynomial of degree 5 to approximate R. The maximum error
* of this polynomial approximation is bounded by 2**-59. In
* other words,
* | 2.0+P1*z+...+P5*z - R(z) | <= 2
* | |
* The computation of exp(r) thus becomes
- * 2*r
- * exp(r) = 1 + -------
- * R - r
- * r*R1(r)
+ * 2*r
+ * exp(r) = 1 + ----------
+ * R(r) - r
+ * r*c(r)
* = 1 + r + ----------- (for better accuracy)
- * 2 - R1(r)
+ * 2 - c(r)
* where
- * 2 4 10
- * R1(r) = r - (P1*r + P2*r + ... + P5*r ).
+ * 2 4 10
+ * c(r) = r - (P1*r + P2*r + ... + P5*r ).
*
* 3. Scale back to obtain exp(x):
* From step 1, we have
*
* Misc. info.
* For IEEE double
- * if x > 7.09782712893383973096e+02 then exp(x) overflow
- * if x < -7.45133219101941108420e+02 then exp(x) underflow
- *
- * Constants:
- * The hexadecimal values are the intended ones for the following
- * constants. The decimal values may be used, provided that the
- * compiler will convert from decimal to binary accurately enough
- * to produce the hexadecimal values shown.
+ * if x > 709.782712893383973096 then exp(x) overflows
+ * if x < -745.133219101941108420 then exp(x) underflows
*/
#include "libm.h"
static const double
-halF[2] = {0.5,-0.5,},
-huge = 1.0e+300,
-o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
-u_threshold = -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */
-ln2HI[2] = { 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
- -6.93147180369123816490e-01},/* 0xbfe62e42, 0xfee00000 */
-ln2LO[2] = { 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
- -1.90821492927058770002e-10},/* 0xbdea39ef, 0x35793c76 */
+half[2] = {0.5,-0.5},
+ln2hi = 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
+ln2lo = 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
-static const volatile double
-twom1000 = 9.33263618503218878990e-302; /* 2**-1000=0x01700000,0 */
-
double exp(double x)
{
- double y,hi=0.0,lo=0.0,c,t,twopk;
- int32_t k=0,xsb;
+ double hi, lo, c, xx;
+ int k, sign;
uint32_t hx;
GET_HIGH_WORD(hx, x);
- xsb = (hx>>31)&1; /* sign bit of x */
+ sign = hx>>31;
hx &= 0x7fffffff; /* high word of |x| */
- /* filter out non-finite argument */
- if (hx >= 0x40862E42) { /* if |x| >= 709.78... */
- if (hx >= 0x7ff00000) {
- uint32_t lx;
-
- GET_LOW_WORD(lx,x);
- if (((hx&0xfffff)|lx) != 0) /* NaN */
- return x+x;
- return xsb==0 ? x : 0.0; /* exp(+-inf)={inf,0} */
+ /* special cases */
+ if (hx >= 0x40862e42) { /* if |x| >= 709.78... */
+ 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 > o_threshold)
- return huge*huge; /* overflow */
- if (x < u_threshold)
- return twom1000*twom1000; /* underflow */
}
/* argument reduction */
- if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
- if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
- hi = x-ln2HI[xsb];
- lo = ln2LO[xsb];
- k = 1 - xsb - xsb;
- } else {
- k = (int)(invln2*x+halF[xsb]);
- t = k;
- hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */
- lo = t*ln2LO[0];
- }
+ if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
+ if (hx >= 0x3ff0a2b2) /* if |x| >= 1.5 ln2 */
+ k = (int)(invln2*x + half[sign]);
+ else
+ k = 1 - sign - sign;
+ hi = x - k*ln2hi; /* k*ln2hi is exact here */
+ lo = k*ln2lo;
STRICT_ASSIGN(double, x, hi - lo);
- } else if(hx < 0x3e300000) { /* |x| < 2**-28 */
- /* raise inexact */
- if (huge+x > 1.0)
- return 1.0+x;
- } else
+ } else if (hx > 0x3e300000) { /* if |x| > 2**-28 */
k = 0;
+ hi = x;
+ lo = 0;
+ } else {
+ /* inexact if x!=0 */
+ FORCE_EVAL(0x1p1023 + x);
+ return 1 + x;
+ }
/* x is now in primary range */
- t = x*x;
- if (k >= -1021)
- INSERT_WORDS(twopk, 0x3ff00000+(k<<20), 0);
- else
- INSERT_WORDS(twopk, 0x3ff00000+((k+1000)<<20), 0);
- c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
+ xx = x*x;
+ c = x - xx*(P1+xx*(P2+xx*(P3+xx*(P4+xx*P5))));
+ x = 1 + (x*c/(2-c) - lo + hi);
if (k == 0)
- return 1.0 - ((x*c)/(c-2.0) - x);
- y = 1.0-((lo-(x*c)/(2.0-c))-hi);
- if (k < -1021)
- return y*twopk*twom1000;
- if (k == 1024)
- return y*2.0*0x1p1023;
- return y*twopk;
+ return x;
+ return scalbn(x, k);
}
float exp10f(float x)
{
static const float p10[] = {
- 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1,
+ 1e-7f, 1e-6f, 1e-5f, 1e-4f, 1e-3f, 1e-2f, 1e-1f,
1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7
};
float n, y = modff(x, &n);
#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
#include "libm.h"
static const float
-halF[2] = {0.5,-0.5,},
-huge = 1.0e+30,
-o_threshold = 8.8721679688e+01, /* 0x42b17180 */
-u_threshold = -1.0397208405e+02, /* 0xc2cff1b5 */
-ln2HI[2] = { 6.9314575195e-01, /* 0x3f317200 */
- -6.9314575195e-01,},/* 0xbf317200 */
-ln2LO[2] = { 1.4286067653e-06, /* 0x35bfbe8e */
- -1.4286067653e-06,},/* 0xb5bfbe8e */
-invln2 = 1.4426950216e+00, /* 0x3fb8aa3b */
+half[2] = {0.5,-0.5},
+ln2hi = 6.9314575195e-1f, /* 0x3f317200 */
+ln2lo = 1.4286067653e-6f, /* 0x35bfbe8e */
+invln2 = 1.4426950216e+0f, /* 0x3fb8aa3b */
/*
* Domain [-0.34568, 0.34568], range ~[-4.278e-9, 4.447e-9]:
* |x*(exp(x)+1)/(exp(x)-1) - p(x)| < 2**-27.74
*/
-P1 = 1.6666625440e-1, /* 0xaaaa8f.0p-26 */
-P2 = -2.7667332906e-3; /* -0xb55215.0p-32 */
-
-static const volatile float twom100 = 7.8886090522e-31; /* 2**-100=0x0d800000 */
+P1 = 1.6666625440e-1f, /* 0xaaaa8f.0p-26 */
+P2 = -2.7667332906e-3f; /* -0xb55215.0p-32 */
float expf(float x)
{
- float y,hi=0.0,lo=0.0,c,t,twopk;
- int32_t k=0,xsb;
+ float hi, lo, c, xx;
+ int k, sign;
uint32_t hx;
GET_FLOAT_WORD(hx, x);
- xsb = (hx>>31)&1; /* sign bit of x */
+ sign = hx >> 31; /* sign bit of x */
hx &= 0x7fffffff; /* high word of |x| */
- /* filter out non-finite argument */
- if (hx >= 0x42b17218) { /* if |x|>=88.721... */
+ /* special cases */
+ if (hx >= 0x42b17218) { /* if |x| >= 88.722839f or NaN */
if (hx > 0x7f800000) /* NaN */
- return x+x;
- if (hx == 0x7f800000) /* exp(+-inf)={inf,0} */
- return xsb==0 ? x : 0.0;
- if (x > o_threshold)
- return huge*huge; /* overflow */
- if (x < u_threshold)
- return twom100*twom100; /* underflow */
+ return x;
+ if (!sign) {
+ /* overflow if x!=inf */
+ STRICT_ASSIGN(float, x, x * 0x1p127f);
+ return x;
+ }
+ if (hx == 0x7f800000) /* -inf */
+ return 0;
+ if (hx >= 0x42cff1b5) { /* x <= -103.972084f */
+ /* underflow */
+ STRICT_ASSIGN(float, x, 0x1p-100f*0x1p-100f);
+ return x;
+ }
}
/* argument reduction */
- if (hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */
- if (hx < 0x3F851592) { /* and |x| < 1.5 ln2 */
- hi = x-ln2HI[xsb];
- lo = ln2LO[xsb];
- k = 1 - xsb - xsb;
- } else {
- k = invln2*x + halF[xsb];
- t = k;
- hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */
- lo = t*ln2LO[0];
- }
+ if (hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */
+ if (hx > 0x3f851592) /* if |x| > 1.5 ln2 */
+ k = invln2*x + half[sign];
+ else
+ k = 1 - sign - sign;
+ hi = x - k*ln2hi; /* k*ln2hi is exact here */
+ lo = k*ln2lo;
STRICT_ASSIGN(float, x, hi - lo);
- } else if(hx < 0x39000000) { /* |x|<2**-14 */
- /* raise inexact */
- if (huge+x > 1.0f)
- return 1.0f + x;
- } else
+ } else if (hx > 0x39000000) { /* |x| > 2**-14 */
k = 0;
+ hi = x;
+ lo = 0;
+ } else {
+ /* raise inexact */
+ FORCE_EVAL(0x1p127f + x);
+ return 1 + x;
+ }
/* x is now in primary range */
- t = x*x;
- if (k >= -125)
- SET_FLOAT_WORD(twopk, 0x3f800000+(k<<23));
- else
- SET_FLOAT_WORD(twopk, 0x3f800000+((k+100)<<23));
- c = x - t*(P1+t*P2);
+ xx = x*x;
+ c = x - xx*(P1+xx*P2);
+ x = 1 + (x*c/(2-c) - lo + hi);
if (k == 0)
- return 1.0f - ((x*c)/(c - 2.0f) - x);
- y = 1.0f - ((lo - (x*c)/(2.0f - c)) - hi);
- if (k < -125)
- return y*twopk*twom100;
- if (k == 128)
- return y*2.0f*0x1p127f;
- return y*twopk;
+ return x;
+ return scalbnf(x, k);
}
* x k f
* e = 2 e.
*
- * A Pade' form of degree 2/3 is used to approximate exp(f) - 1
+ * A Pade' form of degree 5/6 is used to approximate exp(f) - 1
* in the basic range [-0.5 ln 2, 0.5 ln 2].
*
*
2.0000000000000000000897E0L,
};
static const long double
-C1 = 6.9314575195312500000000E-1L,
-C2 = 1.4286068203094172321215E-6L,
-MAXLOGL = 1.1356523406294143949492E4L,
-MINLOGL = -1.13994985314888605586758E4L,
-LOG2EL = 1.4426950408889634073599E0L;
+LN2HI = 6.9314575195312500000000E-1L,
+LN2LO = 1.4286068203094172321215E-6L,
+LOG2E = 1.4426950408889634073599E0L;
long double expl(long double x)
{
long double px, xx;
- int n;
+ int k;
if (isnan(x))
return x;
- if (x > MAXLOGL)
- return INFINITY;
- if (x < MINLOGL)
- return 0.0;
+ 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;
- /* Express e**x = e**g 2**n
- * = e**g e**(n loge(2))
- * = e**(g + n loge(2))
+ /* Express e**x = e**f 2**k
+ * = e**(f + k ln(2))
*/
- px = floorl(LOG2EL * x + 0.5); /* floor() truncates toward -infinity. */
- n = px;
- x -= px * C1;
- x -= px * C2;
+ px = floorl(LOG2E * x + 0.5);
+ k = px;
+ x -= px * LN2HI;
+ x -= px * LN2LO;
- /* rational approximation for exponential
- * of the fractional part:
- * e**x = 1 + 2x P(x**2)/(Q(x**2) - P(x**2))
+ /* rational approximation of the fractional part:
+ * e**x = 1 + 2x P(x**2)/(Q(x**2) - x P(x**2))
*/
xx = x * x;
px = x * __polevll(xx, P, 2);
- x = px/(__polevll(xx, Q, 3) - px);
+ x = px/(__polevll(xx, Q, 3) - px);
x = 1.0 + 2.0 * x;
- x = scalbnl(x, n);
- return x;
+ return scalbnl(x, k);
}
#endif
projects / musl / commitdiff