projects
/
musl
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
add more arch-specific MAP_ macros to bits/mman.h
[musl]
/
src
/
math
/
exp2.c
diff --git
a/src/math/exp2.c
b/src/math/exp2.c
index
08c21a6
..
8e25228
100644
(file)
--- a/
src/math/exp2.c
+++ b/
src/math/exp2.c
@@
-27,11
+27,9
@@
#include "libm.h"
#include "libm.h"
-#define TBLBITS 8
-#define TBLSIZE (1 << TBLBITS)
+#define TBLSIZE 256
static const double
static const double
-huge = 0x1p1000,
redux = 0x1.8p52 / TBLSIZE,
P1 = 0x1.62e42fefa39efp-1,
P2 = 0x1.ebfbdff82c575p-3,
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;
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,
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 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) {
/* 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;
}
} 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;
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 -= 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))));
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);
}
}