diff options
Diffstat (limited to 'src/math/expm1f.c')
| -rw-r--r-- | src/math/expm1f.c | 65 |
1 files changed, 26 insertions, 39 deletions
diff --git a/src/math/expm1f.c b/src/math/expm1f.c index d9568466..698ab45f 100644 --- a/src/math/expm1f.c +++ b/src/math/expm1f.c @@ -16,8 +16,6 @@ #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 */ @@ -32,35 +30,27 @@ Q2 = 1.5807170421e-3; /* 0xcf3010.0p-33 */ 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; @@ -70,7 +60,7 @@ float expm1f(float x) 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; @@ -78,8 +68,9 @@ float expm1f(float x) 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; @@ -91,9 +82,9 @@ float expm1f(float x) 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) { @@ -101,24 +92,20 @@ float expm1f(float x) 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; } |