diff --git a/lib/math/floats.ww b/lib/math/floats.ww index 0983ef93..19523a83 100644 --- a/lib/math/floats.ww +++ b/lib/math/floats.ww @@ -1,15 +1,11 @@ -// floats — f64 classification, sign, bit-reinterpret core, and the -// subnormal-normalize step of the f64 decompose half. Ported from +// floats — f64 classification, sign, bit-reinterpret core, and the f64 +// decompose half (subnormal-normalize + frexp). 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. +// issubnormalf64/normalizef64/frexpf64). frexpf64's zero guard `n == 0f64` +// rides the #103 fix (no-decimal f64 literal now materialized into XMM) +// and its (f64, i64) tuple return rides the #105 fix (tuple f64-word read). +// f32 variants stay deferred (#104 blocks the f32 call-arg narrowing), as +// do the ldexp/modfrac/nextafter family (need f64 DIVIDE + the INF const). package math; @@ -56,6 +52,18 @@ export def F64_EXPONENT_MASK: u64 = (1 << F64_EXPONENT_BITS) - 1; // ref/hare/math/floats.ha:75 def F64_SIGN_MASK: u64 = 1u64 << 63; +// Mask that clears an f64's exponent field, keeping sign + mantissa. +// ref/hare/math/floats.ha:77. Hare hardcodes the 0x800FFFFFFFFFFFFF binary +// literal because its lexer can't const-fold the expression; ww's #88 +// def-const-fold can, so the readable form is kept. floats.ha:79's NOTE +// expression has an `0u64 &` upstream typo (it would yield 0); the value it +// documents is exactly ~(F64_EXPONENT_MASK << F64_MANTISSA_BITS). +def F64_EXP_REMOVAL_MASK: u64 = ~(F64_EXPONENT_MASK << F64_MANTISSA_BITS); + +// The f64 bit pattern whose exponent field evaluates to zero (0.5 scale). +// ref/hare/math/floats.ha:84 +def F64_EXP_ZERO: u64 = (F64_EXPONENT_BIAS - 1) << F64_MANTISSA_BITS; + // Returns true if the given floating-point number is NaN. // ref/hare/math/floats.ha:144 (Hare's expression body inlined into a // block: ww has no expression-bodied fn form, only brace blocks). @@ -131,3 +139,22 @@ export fn normalizef64(n: f64) (f64, i64) = { }; return (n, 0); }; + +// Breaks a f64 down into its mantissa and exponent. The mantissa will be +// between 0.5 and 1. +// ref/hare/math/floats.ha:278 +export fn frexpf64(n: f64) (f64, i64) = { + if (isnan(n) || isinf(n) || n == 0f64) { + return (n, 0); + }; + const normalized = normalizef64(n); + const normal_float = normalized.0; + const normalization_exp = normalized.1; + const bits = f64bits(normal_float); + const raw_exp: u64 = (bits >> F64_MANTISSA_BITS) & F64_EXPONENT_MASK; + const exp: i64 = normalization_exp + + (raw_exp: i64) - (F64_EXPONENT_BIAS: i64) + 1; + const mantissa: f64 = + f64frombits((bits & F64_EXP_REMOVAL_MASK) | F64_EXP_ZERO); + return (mantissa, exp); +};