hare

commit f3015828b05cee7631259bfa0c8cfa55fe2f6dbe
parent 6ec62fc6657f2010540da02574697c257b14e102
Author: Lassi Pulkkinen <lassi@pulk.fi>
Date:   Thu,  6 Apr 2023 20:02:20 +0300

math::powf64: Check for overflowing exponent

Add two special cases from Go's implementation (already cited as a reference)
that were previously missing. The (p_frac != 0f64 && x < 0f64) one was
already in the test suite and seemingly working, but I'd rather have an
explicit check.

Add tests for overflowing exponents.

Signed-off-by: Lassi Pulkkinen <lassi@pulk.fi>

Diffstat:
M
math/+test/math.ha
 | 
6
++++++
M
math/math.ha
 | 
21
++++++++++++++-------

2 files changed, 20 insertions(+), 7 deletions(-)
diff --git a/math/+test/math.ha b/math/+test/math.ha
@@ -110,12 +110,18 @@
 	// Very large positive integer
 	assert(!isinf(powf64(2f64, 1020f64)));
 	assert(isinf(powf64(2f64, 1050f64)));
+	// Very very large positive integer
+	assert(isinf(powf64(2f64, F64_MAX_NORMAL)));
+	assert(powf64(0.5f64, F64_MAX_NORMAL) == 0f64);
 	// Negative integer
 	assert(powf64(2f64, -1f64) == 0.5f64);
 	assert(powf64(2f64, -2f64) == 0.25f64);
 	// Very small negative integer
 	assert(powf64(2f64, -1020f64) > 0f64);
 	assert(powf64(2f64, -1080f64) == 0f64);
+	// Very very small negative integer
+	assert(powf64(2f64, -F64_MAX_NORMAL) == 0f64);
+	assert(isinf(powf64(0.5f64, -F64_MAX_NORMAL)));
 	// Positive fractional powers
 	assert(eqwithin(powf64(2f64, 1.5f64),
 		2.8284271247461900976033774f64,
diff --git a/math/math.ha b/math/math.ha
@@ -731,13 +731,19 @@ export fn powf64(x: f64, p: f64) f64 = {
 		return 1f64 / sqrtf64(x);
 	};
 
-	const x_parts = frexp(x);
-	let x_mantissa = x_parts.0;
-	let x_exp = x_parts.1;
-
-	const p_int_frac = modfracf64(absf64(p));
-	let p_int = p_int_frac.0;
-	let p_frac = p_int_frac.1;
+	let (p_int, p_frac) = modfracf64(absf64(p));
+	if (p_frac != 0f64 && x < 0f64) {
+		return NAN;
+	};
+	if (p_int > types::I64_MAX: f64) {
+		if (x == -1f64) {
+			return 1f64;
+		} else if ((absf64(x) < 1f64) == (p > 0f64)) {
+			return 0f64;
+		} else {
+			return INF;
+		};
+	};
 
 	let res_mantissa = 1f64;
 	let res_exp = 0i64;
@@ -756,6 +762,7 @@ export fn powf64(x: f64, p: f64) f64 = {
 	// Repeatedly square our number x, for each bit in our power p.
 	// If the current bit is 1 in p, add the respective power of x to our
 	// result.
+	let (x_mantissa, x_exp) = frexp(x);
 	for (let i = p_int: i64; i != 0; i >>= 1) {
 		// Check for over/underflow.
 		if (x_exp <= -1i64 << (F64_EXPONENT_BITS: i64)) {