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)

1 files changed, 10 insertions, 11 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);
 }