Fix some special cases in log-overkill related to subnormals.

2 files changed, 31 insertions, 12 deletions

diff --git a/lib/math/log-overkill.myr b/lib/math/log-overkill.myr
index ecac77d..6c74486 100644
--- a/lib/math/log-overkill.myr
+++ b/lib/math/log-overkill.myr
@@ -206,17 +206,30 @@ const logoverkill64 = {x : flt64
 	/* Special cases */
 	if std.isnan(x)
 		-> (std.flt64frombits(0x7ff8000000000000), std.flt64frombits(0x7ff8000000000000))
-	elif xe == -1024 && xs == 0
+	elif xe <= -1023 && xs == 0
 		/* log (+/- 0) is -infinity */
 		-> (std.flt64frombits(0xfff0000000000000), std.flt64frombits(0xfff0000000000000))
-	elif xe == 1023
-		-> (std.flt64frombits(0xfff8000000000000), std.flt64frombits(0xfff8000000000000))
 	elif xn
 		/* log(-anything) is NaN */
 		-> (std.flt64frombits(0x7ff8000000000000), std.flt64frombits(0x7ff8000000000000))
 	;;
 
+	if xe <= -1023
+		/*
+		   We depend on being able to pick bits out of xs as if it were normal, so
+		   normalize any subnormals.
+		 */
+		xe++
+		var check = 1 << 52
+		while xs & check == 0
+			xs <<= 1
+			xe--
+		;;
+		xsf = std.flt64assem(false, 0, xs)
+	;;
+
 	var shift = 0
+	var non_trivial = 0
 	var then_invert = false
 	var j1, F1, f1, logF1_hi, logF1_lo, F1_inv_hi, F1_inv_lo, fF1_hi, fF1_lo
 
@@ -237,8 +250,9 @@ const logoverkill64 = {x : flt64
 		var xinv_hi = 1.0 / x
 		var xinv_lo = fma64(-1.0 * xinv_hi, x, 1.0) / x
 		(tn, te, ts) = std.flt64explode(xinv_hi)
-		shift = ((47 >= te) : uint64) * ((47 - te) : uint64) + ((47 < te) : uint64) * 64
-		j1 = (ts >> shift) & 0x1f
+		non_trivial = ((47 >= te) : uint64) * ((47 - te < 64) : uint64)
+		shift = non_trivial * ((47 - te) : uint64)
+		j1 = non_trivial * ((ts >> shift) & 0x1f)
 		var F1m1 = scale2((j1 : flt64), -5)
 		F1 = 1.0 + F1m1
 		var f1_hi, f1_lo
@@ -266,8 +280,9 @@ const logoverkill64 = {x : flt64
 
 	/* F2 */
 	(tn, te, ts) = std.flt64explode(fF1_hi)
-	shift = ((42 >= te) : uint64) * ((42 - te) : uint64) + ((42 < te) : uint64) * 64
-	var j2 = (ts >> shift) & 0x1f
+	non_trivial = ((42 >= te) : uint64) * ((42 - te < 64) : uint64)
+	shift = non_trivial * ((42 - te) : uint64)
+	var j2 = non_trivial * ((ts >> shift) & 0x1f)
 	var F2m1 = scale2((j2 : flt64), -10)
 	var F2 = 1.0 + F2m1
 	var f2_hi, f2_lo
@@ -284,8 +299,9 @@ const logoverkill64 = {x : flt64
 
 	/* F3 (just like F2) */
 	(tn, te, ts) = std.flt64explode(fF2_hi)
-	shift = ((37 >= te) : uint64) * ((37 - te) : uint64) + ((37 < te) : uint64) * 64
-	var j3 = (ts >> shift) & 0x1f
+	non_trivial = ((37 >= te) : uint64) * ((37 - te < 64) : uint64)
+	shift = non_trivial * ((37 - te) : uint64)
+	var j3 = non_trivial * ((ts >> shift) & 0x1f)
 	var F3m1 = scale2((j3 : flt64), -15)
 	var F3 = 1.0 + F3m1
 	var f3_hi, f3_lo
@@ -302,8 +318,9 @@ const logoverkill64 = {x : flt64
 
 	/* F4 (just like F2) */
 	(tn, te, ts) = std.flt64explode(fF3_hi)
-	shift = ((32 >= te) : uint64) * ((32 - te) : uint64) + ((32 < te) : uint64) * 64
-	var j4 = (ts >> shift) & 0x1f
+	non_trivial = ((32 >= te) : uint64) * ((32 - te < 64) : uint64)
+	shift = non_trivial * ((32 - te) : uint64)
+	var j4 = non_trivial * ((ts >> shift) & 0x1f)
 	var F4m1 = scale2((j4 : flt64), -20)
 	var F4 = 1.0 + F4m1
 	var f4_hi, f4_lo
@@ -376,7 +393,6 @@ const logoverkill64 = {x : flt64
 	var xel2_hi, xel2_lo, lx_hi, lx_lo
 	(xel2_hi, xel2_lo) = hl_mult((xe : flt64), 0.0, std.flt64frombits(0x3fe62e42fefa39ef), std.flt64frombits(0x3c7abc9e3b39803f))
 
-
 	(t1, t2) = slow2sum(xel2_lo, lsig_lo)
 	(lx_lo, t1) = slow2sum(lsig_hi, t1)
 	(lx_hi, lx_lo) = slow2sum(xel2_hi, lx_lo)
diff --git a/lib/math/test/log-overkill.myr b/lib/math/test/log-overkill.myr
index b407a8d..03e67ef 100644
--- a/lib/math/test/log-overkill.myr
+++ b/lib/math/test/log-overkill.myr
@@ -227,6 +227,9 @@ const log02 = {c
 		(0x4014a5ce06d7df05, 0x3ffa42cc8df38d10, 0xbc5e54e0ca2ed44c),
 		(0x3f68d1447d5a29e8, 0xc017328d1195bac4, 0x3cb4fd5db024d1da),
 		(0x404963f9a80919eb, 0x400f6b9160b05ec4, 0x3c4526374db12c53),
+		(0x3fe78000b3cf1a39, 0xbfd3c2508d81ebf9, 0x3c7d6df43454d213),
+		(0x7fe6c53d8cef3d27, 0x40862b8a1ec909c8, 0x3cf9a8752da53a7e),
+		(0x000342cdeeb18fc9, 0xc0862fe5598ee7e6, 0xbd2bf7df7d1e9517),
 	][:]
 
 	for (x, y1, y2) : inputs