[musl] [PATCH 2/2] math: fix ld80 powl(x,huge) and powl(LDBL_MAX,small)

[Szabolcs Nagy <nsz@port70.net>]

powl used >= LDBL_MAX as infinity check, but LDBL_MAX is finite, so
this can cause wrong results e.g. powl(LDBL_MAX, 0.5) returned inf
or powl(2, LDBL_MAX) returned inf without raising overflow.

huge y values (close to LDBL_MAX) could cause intermediate results to
overflow (computing y * log2(x) with more than long double precision)
and e.g. powl(0.5, 0x1p16380L) or powl(10, 0x1p16380L) returned nan.
this is fixed by handling huge y early since that always overflows or
underflows.

reported by Paul Zimmermann against expl10 (which uses powl).
---
 src/math/powl.c | 34 +++++++++++++++++++++-------------
 1 file changed, 21 insertions(+), 13 deletions(-)

diff --git a/src/math/powl.c b/src/math/powl.c
index 5b6da07b..6f64ea71 100644
--- a/src/math/powl.c
+++ b/src/math/powl.c
@@ -212,25 +212,33 @@ long double powl(long double x, long double y)
 	}
 	if (x == 1.0)
 		return 1.0; /* 1**y = 1, even if y is nan */
-	if (x == -1.0 && !isfinite(y))
-		return 1.0; /* -1**inf = 1 */
 	if (y == 0.0)
 		return 1.0; /* x**0 = 1, even if x is nan */
 	if (y == 1.0)
 		return x;
-	if (y >= LDBL_MAX) {
-		if (x > 1.0 || x < -1.0)
-			return INFINITY;
-		if (x != 0.0)
-			return 0.0;
-	}
-	if (y <= -LDBL_MAX) {
-		if (x > 1.0 || x < -1.0)
+	/* if y*log2(x) < log2(LDBL_TRUE_MIN)-1 then x^y uflows to 0
+	   if y*log2(x) > -log2(LDBL_TRUE_MIN)+1 > LDBL_MAX_EXP then x^y oflows
+	   if |x|!=1 then |log2(x)| > |log(x)| > LDBL_EPSILON/2 so
+	   x^y oflows/uflows if |y|*LDBL_EPSILON/2 > -log2(LDBL_TRUE_MIN)+1 */
+	if (fabsl(y) > 2*(-LDBL_MIN_EXP+LDBL_MANT_DIG+1)/LDBL_EPSILON) {
+		/* y is not an odd int */
+		if (x == -1.0)
+			return 1.0;
+		if (y == INFINITY) {
+			if (x > 1.0 || x < -1.0)
+				return INFINITY;
 			return 0.0;
-		if (x != 0.0 || y == -INFINITY)
+		}
+		if (y == -INFINITY) {
+			if (x > 1.0 || x < -1.0)
+				return 0.0;
 			return INFINITY;
+		}
+		if ((x > 1.0 || x < -1.0) == (y > 0))
+			return huge * huge;
+		return twom10000 * twom10000;
 	}
-	if (x >= LDBL_MAX) {
+	if (x == INFINITY) {
 		if (y > 0.0)
 			return INFINITY;
 		return 0.0;
@@ -253,7 +261,7 @@ long double powl(long double x, long double y)
 			yoddint = 1;
 	}
 
-	if (x <= -LDBL_MAX) {
+	if (x == -INFINITY) {
 		if (y > 0.0) {
 			if (yoddint)
 				return -INFINITY;
-- 
2.41.0