From 39c910fb061114e6aa5c3bf2c94b1d7262d62221 Mon Sep 17 00:00:00 2001
From: Szabolcs Nagy <nsz@port70.net>
Date: Wed, 4 Sep 2013 07:51:11 +0000
Subject: 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)
---
 src/math/expm1.c | 72 +++++++++++++++++++++-----------------------------------
 1 file changed, 27 insertions(+), 45 deletions(-)

(limited to 'src/math/expm1.c')

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;
 }
-- 
cgit v1.2.1