From 263257f2c18747575385a555d9925108b24b8254 Mon Sep 17 00:00:00 2001
From: Hojun-Cho <ghwns6743@gmail.com>
Date: Mon, 25 May 2026 15:06:19 +0900
Subject: [PATCH] lib/math: math::floats fold-2a issubnormalf64 + normalizef64

Ports the subnormal-normalize step of the f64 decompose half from
ref/hare/math/floats.ha: issubnormalf64 (floats.ha:179) and normalizef64
(floats.ha:256, the f64-multiply-on-subnormal that yields (f64, i64)).

frexpf64 (floats.ha:278) is held back, not ported: its Hare-exact zero
guard `n == 0f64` miscompiles. A no-decimal `0f64` literal used as an f64
comparison operand is materialized into a GPR and never moved to XMM, so
the UCOMISD reads a stale operand and `n == 0f64` is wrong for every n.
Both stages emit this identically, so the byte-id gates are blind to it.
`0.0` compiles correctly but substituting it would be a workaround
(rule 7), so frexpf64 waits for the cgen fix. normalizef64/issubnormalf64
touch neither the broken literal form nor any tuple-field comparison, so
they are correct and land now.
---
 lib/math/floats.ww | 37 +++++++++++++++++++++++++++++++++----
 1 file changed, 33 insertions(+), 4 deletions(-)

diff --git a/lib/math/floats.ww b/lib/math/floats.ww
index cb9ce8ff..0983ef93 100644
--- a/lib/math/floats.ww
+++ b/lib/math/floats.ww
@@ -1,7 +1,15 @@
-// floats — f64 classification, sign, and bit-reinterpret core. Ported
-// from ref/hare/math/floats.ha (fold-1: the minimal classify/sign/bits
-// surface). f32 variants, NAN/INF + magnitude consts, and
-// frexp/ldexp/normalize/modfrac/nextafter are deferred to a later fold.
+// floats — f64 classification, sign, bit-reinterpret core, and the
+// subnormal-normalize step of the f64 decompose half. Ported from
+// ref/hare/math/floats.ha (fold-1: classify/sign/bits; fold-2a:
+// issubnormalf64/normalizef64). frexpf64 (floats.ha:278) is held back: its
+// Hare-exact zero guard `n == 0f64` miscompiles — a no-decimal `0f64`
+// literal in an f64 comparison is materialized into a GPR and never moved
+// to XMM, so the compare reads a stale operand (both stages identically,
+// so the byte-id gates are blind to it). `0.0` would dodge it, but that is
+// a workaround (CLAUDE.md rule 7); frexpf64 lands once the cgen bug is
+// fixed. f32 variants, NAN/INF + magnitude consts, and the
+// ldexp/modfrac/nextafter family (need f64 DIVIDE + the INF const) are
+// deferred to a later fold.
 
 package math;
 
@@ -64,6 +72,15 @@ export fn isinf(n: f64) bool = {
 	return exp == F64_EXPONENT_MASK && mant == 0;
 };
 
+// Returns true if the given f64 is subnormal.
+// ref/hare/math/floats.ha:179
+export fn issubnormalf64(n: f64) bool = {
+	const bits = f64bits(n);
+	const mant = bits & F64_MANTISSA_MASK;
+	const exp = bits >> F64_MANTISSA_BITS & F64_EXPONENT_MASK;
+	return exp == 0 && mant != 0;
+};
+
 // Returns the absolute value of f64 n.
 // ref/hare/math/floats.ha:195
 export fn absf64(n: f64) f64 = {
@@ -102,3 +119,15 @@ export fn copysignf64(x: f64, y: f64) f64 = {
 	return f64frombits((f64bits(x) & ~F64_SIGN_MASK) |
 		(f64bits(y) & F64_SIGN_MASK));
 };
+
+// Takes a potentially subnormal f64 n and returns a normal f64 normal_float
+// and an exponent exp such that n == normal_float * 2^{exp}.
+// ref/hare/math/floats.ha:256
+export fn normalizef64(n: f64) (f64, i64) = {
+	if (issubnormalf64(n)) {
+		const factor = 1i64 << (F64_MANTISSA_BITS: i64);
+		const normal_float = (n * (factor: f64));
+		return (normal_float, -(F64_MANTISSA_BITS: i64));
+	};
+	return (n, 0);
+};