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)

8 files changed, 139 insertions, 182 deletions

diff --git a/src/math/exp.c b/src/math/exp.c
index 5ec8f8a7..2b182acd 100644
--- a/src/math/exp.c
+++ b/src/math/exp.c
@@ -80,7 +80,7 @@ P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
 
 double exp(double x)
 {
-	double hi, lo, c, xx;
+	double_t hi, lo, c, xx, y;
 	int k, sign;
 	uint32_t hx;
 
@@ -89,20 +89,19 @@ double exp(double x)
 	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;
 		}
 	}
 
@@ -128,8 +127,8 @@ double exp(double x)
 	/* 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);
 }
diff --git a/src/math/exp2.c b/src/math/exp2.c
index 8e252280..2e078fb0 100644
--- a/src/math/exp2.c
+++ b/src/math/exp2.c
@@ -305,7 +305,7 @@ static const double tbl[TBLSIZE * 2] = {
  * 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:
@@ -330,41 +330,41 @@ static const double tbl[TBLSIZE * 2] = {
  */
 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] */
diff --git a/src/math/exp2f.c b/src/math/exp2f.c
index ea50db4a..91738f04 100644
--- a/src/math/exp2f.c
+++ b/src/math/exp2f.c
@@ -63,7 +63,7 @@ static const double exp2ft[TBLSIZE] = {
  * 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:
@@ -83,46 +83,42 @@ static const double exp2ft[TBLSIZE] = {
  */
 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;
 }
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];
diff --git a/src/math/expf.c b/src/math/expf.c
index 8aefc917..5572bbf6 100644
--- a/src/math/expf.c
+++ b/src/math/expf.c
@@ -29,7 +29,7 @@ P2 = -2.7667332906e-3f; /* -0xb55215.0p-32 */
 
 float expf(float x)
 {
-	float hi, lo, c, xx;
+	float_t hi, lo, c, xx, y;
 	int k, sign;
 	uint32_t hx;
 
@@ -38,20 +38,17 @@ float expf(float x)
 	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;
 		}
 	}
 
@@ -77,8 +74,8 @@ float expf(float x)
 	/* 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);
 }
diff --git a/src/math/expl.c b/src/math/expl.c
index 50a04297..b62980fa 100644
--- a/src/math/expl.c
+++ b/src/math/expl.c
@@ -100,7 +100,7 @@ long double expl(long double x)
 	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))
diff --git a/src/math/expm1.c b/src/math/expm1.c
index f8f32c46..a7eb2c0b 100644
--- a/src/math/expm1.c
+++ b/src/math/expm1.c
@@ -31,7 +31,7 @@
  *          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,
@@ -107,8 +107,6 @@
 #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 */
@@ -122,39 +120,27 @@ Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
 
 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;
@@ -164,7 +150,7 @@ double expm1(double x)
 				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;
@@ -172,9 +158,9 @@ double expm1(double x)
 		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;
 
@@ -186,9 +172,9 @@ double expm1(double x)
 	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) {
@@ -196,24 +182,20 @@ double expm1(double x)
 			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;
 }
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;
 }