From a164c4c7527d2b74809a0e948ca280beb532cc9d Mon Sep 17 00:00:00 2001 From: Cursor Agent Date: Tue, 25 Aug 2026 22:10:04 +0000 Subject: [PATCH 1/2] Scale RMS residuals by 2^exp without underflowing the power first. 2.0.powi(-1075) is zero even when (nf/df)*2^-1075 is the smallest subnormal. Apply MIN_POSITIVE in 1022-bit steps so a scaled sample just above 2^-1075 still publishes rms_lsb = 2^-1074. This is not MAX_ERR_LSB. Co-authored-by: Steven Giacomelli --- CHANGELOG.md | 5 ++-- crates/surfaces-bake/src/bound.rs | 41 ++++++++++++++++++++++++---- crates/surfaces-bake/src/quantize.rs | 24 ++++++++++++++++ 3 files changed, 63 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d5b12b8..d8a81e7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -34,8 +34,9 @@ 1_048_576 cells, and `i32` overflow are closed `BakeError`s. On-knot lookup is binary search on the ordered knot lists. The operator RMS statistic scales before squaring so a representable tiny residual does - not underflow to 0, and converting the residual to `f64` keeps - subnormals through `2^-1074` finite. The bound is an + not underflow to 0, and converting the residual to `f64` applies the + binary exponent without first underflowing `2^exp`, so a significand + just above `2^-1075` still becomes the smallest subnormal. The bound is an upper bound on deviation from the supplied samples, not a device or accuracy claim. This is not a runtime API change. - Reviewed host crates on `ph-surfaces-bake` (`num-bigint`, `num-rational`, diff --git a/crates/surfaces-bake/src/bound.rs b/crates/surfaces-bake/src/bound.rs index 15c74ee..96324d8 100644 --- a/crates/surfaces-bake/src/bound.rs +++ b/crates/surfaces-bake/src/bound.rs @@ -45,10 +45,10 @@ impl Ratio { /// Approximate the ratio as `f64` for the RMS statistic only. /// - /// Take a 53-bit window of each limb and restore `2^{n_bits - d_bits}` so - /// `1 / 2^1074` (`f64::from_bits(1)`) stays the smallest subnormal instead - /// of becoming infinity. Independent limb conversion still overflows - /// around `2^1024` (`1 + 1e-300`). This is not the bound. + /// Take a 53-bit window of each limb and restore `2^{n_shift - d_shift}` + /// without first materializing `2^exp` when `exp` is below `-1022`: + /// `2.0.powi(-1075)` underflows even when `(nf / df) * 2^exp` is the + /// smallest subnormal. This is not the bound. pub(crate) fn to_f64(&self) -> f64 { if self.is_zero() { return 0.0; @@ -71,7 +71,7 @@ impl Ratio { _ if n_shift > d_shift => return f64::INFINITY, _ => return 0.0, }; - let x = (nf / df) * 2.0_f64.powi(exp); + let x = scale_by_pow2(nf / df, exp); if self.0.is_negative() { -x } else { x } } @@ -91,6 +91,28 @@ impl Ratio { } } +/// `x * 2^exp` without letting `2^exp` underflow before `x` is applied. +fn scale_by_pow2(mut x: f64, mut exp: i32) -> f64 { + if x == 0.0 || !x.is_finite() || exp == 0 { + return x; + } + while exp <= -1022 { + x *= f64::MIN_POSITIVE; + exp += 1022; + if x == 0.0 { + return x; + } + } + while exp >= 1024 { + x *= 2.0_f64.powi(1023); + exp -= 1023; + if !x.is_finite() { + return x; + } + } + x * 2.0_f64.powi(exp) +} + pub(crate) fn ratio_from_f64(x: f64) -> Option { Some(Ratio(dyadic(x)?)) } @@ -190,4 +212,13 @@ mod tests { assert!(x.is_finite()); assert_eq!(x, tiny); } + + #[test] + fn a_near_half_min_subnormal_times_scale_is_the_smallest_subnormal() { + let value = f64::from_bits(0x0c70_0000_0000_0001); + let scale = f64::from_bits(0x3040_0000_0000_0001); + let x = scaled_sample(value, scale).unwrap().to_f64(); + assert!(x.is_finite()); + assert_eq!(x, f64::from_bits(1)); + } } diff --git a/crates/surfaces-bake/src/quantize.rs b/crates/surfaces-bake/src/quantize.rs index ab0f572..2e7f7aa 100644 --- a/crates/surfaces-bake/src/quantize.rs +++ b/crates/surfaces-bake/src/quantize.rs @@ -871,6 +871,30 @@ mod tests { assert_eq!(table.rms_lsb, tiny); } + #[test] + fn a_scaled_near_half_min_subnormal_rms_is_not_zero() { + let value = f64::from_bits(0x0c70_0000_0000_0001); + let scale = f64::from_bits(0x3040_0000_0000_0001); + let table = BakeInput::new( + vec![ + Sample::new(0.0, 0.0, value), + Sample::new(1.0, 0.0, value), + Sample::new(0.0, 1.0, value), + Sample::new(1.0, 1.0, value), + Sample::new(0.5, 0.5, value), + ], + Axis::knots(vec![0, 1]), + Axis::knots(vec![0, 1]), + scale, + ) + .unwrap() + .quantize() + .unwrap(); + assert_eq!(table.values, vec![vec![0, 0], vec![0, 0]]); + assert_eq!(table.max_err_lsb, 1); + assert_eq!(table.rms_lsb, f64::from_bits(1)); + } + #[test] fn a_finite_residual_too_wide_for_i32_is_bound_overflow() { let err = BakeInput::new( From 798d19bf9f5ea52a21c4a707ade4750f6daa0e12 Mon Sep 17 00:00:00 2001 From: Cursor Agent Date: Tue, 25 Aug 2026 22:23:10 +0000 Subject: [PATCH 2/2] Apply the RMS subnormal remainder before the MIN_POSITIVE chunk. Multiplying by 2^-1022 first can make a still-normal significand subnormal and round again on the leftover exponent, turning a three-unit residual into four and the five-sample RMS into two minimum-subnormal units. Co-authored-by: Steven Giacomelli --- CHANGELOG.md | 6 ++-- crates/surfaces-bake/src/bound.rs | 41 +++++++++++++++++++++++----- crates/surfaces-bake/src/quantize.rs | 26 ++++++++++++++++++ 3 files changed, 64 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d8a81e7..7f65eaa 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -35,8 +35,10 @@ lookup is binary search on the ordered knot lists. The operator RMS statistic scales before squaring so a representable tiny residual does not underflow to 0, and converting the residual to `f64` applies the - binary exponent without first underflowing `2^exp`, so a significand - just above `2^-1075` still becomes the smallest subnormal. The bound is an + binary exponent without first underflowing `2^exp` and applies any + remaining negative exponent while the significand is still normal, so a + significand just above `2^-1075` still becomes the smallest subnormal + with only one subnormal rounding. The bound is an upper bound on deviation from the supplied samples, not a device or accuracy claim. This is not a runtime API change. - Reviewed host crates on `ph-surfaces-bake` (`num-bigint`, `num-rational`, diff --git a/crates/surfaces-bake/src/bound.rs b/crates/surfaces-bake/src/bound.rs index 96324d8..1bcde54 100644 --- a/crates/surfaces-bake/src/bound.rs +++ b/crates/surfaces-bake/src/bound.rs @@ -48,7 +48,9 @@ impl Ratio { /// Take a 53-bit window of each limb and restore `2^{n_shift - d_shift}` /// without first materializing `2^exp` when `exp` is below `-1022`: /// `2.0.powi(-1075)` underflows even when `(nf / df) * 2^exp` is the - /// smallest subnormal. This is not the bound. + /// smallest subnormal. The remainder is applied while the significand is + /// still normal so there is only one subnormal rounding. This is not the + /// bound. pub(crate) fn to_f64(&self) -> f64 { if self.is_zero() { return 0.0; @@ -92,16 +94,20 @@ impl Ratio { } /// `x * 2^exp` without letting `2^exp` underflow before `x` is applied. +/// +/// When `exp <= -1022`, apply `2^{exp + 1022}` first so the significand stays +/// normal, then multiply by `2^-1022` once. Multiplying by `MIN_POSITIVE` +/// first can make `x` subnormal and round a second time on the remainder. fn scale_by_pow2(mut x: f64, mut exp: i32) -> f64 { if x == 0.0 || !x.is_finite() || exp == 0 { return x; } - while exp <= -1022 { - x *= f64::MIN_POSITIVE; - exp += 1022; - if x == 0.0 { - return x; + if exp <= -1022 { + let rem = exp + 1022; + if rem != 0 { + x *= 2.0_f64.powi(rem); } + return x * f64::MIN_POSITIVE; } while exp >= 1024 { x *= 2.0_f64.powi(1023); @@ -161,7 +167,7 @@ fn dyadic(x: f64) -> Option { #[cfg(test)] mod tests { - use super::{Ratio, scaled_sample}; + use super::{Ratio, lerp_ratio, ratio_from_f64, scaled_sample}; #[test] fn ceil_of_half_is_one() { @@ -221,4 +227,25 @@ mod tests { assert!(x.is_finite()); assert_eq!(x, f64::from_bits(1)); } + + #[test] + fn a_subnormal_lerp_residual_rounds_once_to_three_units() { + // Knots [0, 17] × [0, 1], rows [0, 1] / [0, 1]. The off-knot sample + // at x = 0x0088_0000_0000_0000, y = 0.5 has an exact residual of three + // minimum-subnormal units. Applying 2^-1022 first rounds the + // intermediate to seven units and the remaining ×0.5 to four. + let x = ratio_from_f64(f64::from_bits(0x0088_0000_0000_0000)).unwrap(); + let y = ratio_from_f64(0.5).unwrap(); + let sample = scaled_sample(f64::from_bits(0x0046_9696_9696_9697), 1.0).unwrap(); + let x0 = Ratio::from_i128(0); + let x1 = Ratio::from_i128(17); + let y0 = Ratio::from_i128(0); + let y1 = Ratio::from_i128(1); + let v0 = Ratio::from_i128(0); + let v1 = Ratio::from_i128(1); + let row = lerp_ratio(&x, &x0, &x1, &v0, &v1).unwrap(); + let reconstructed = lerp_ratio(&y, &y0, &y1, &row, &row).unwrap(); + let residual = sample.sub(&reconstructed).unwrap(); + assert_eq!(residual.to_f64().abs(), f64::from_bits(3)); + } } diff --git a/crates/surfaces-bake/src/quantize.rs b/crates/surfaces-bake/src/quantize.rs index 2e7f7aa..74a256d 100644 --- a/crates/surfaces-bake/src/quantize.rs +++ b/crates/surfaces-bake/src/quantize.rs @@ -895,6 +895,32 @@ mod tests { assert_eq!(table.rms_lsb, f64::from_bits(1)); } + #[test] + fn a_subnormal_lerp_rms_does_not_double_round() { + // Five-sample RMS of a three-unit residual is one minimum-subnormal + // unit. Double-rounding the residual to four units would report two. + let x = f64::from_bits(0x0088_0000_0000_0000); + let value = f64::from_bits(0x0046_9696_9696_9697); + let table = BakeInput::new( + vec![ + Sample::new(0.0, 0.0, 0.0), + Sample::new(17.0, 0.0, 1.0), + Sample::new(0.0, 1.0, 0.0), + Sample::new(17.0, 1.0, 1.0), + Sample::new(x, 0.5, value), + ], + Axis::knots(vec![0, 17]), + Axis::knots(vec![0, 1]), + 1.0, + ) + .unwrap() + .quantize() + .unwrap(); + assert_eq!(table.values, vec![vec![0, 1], vec![0, 1]]); + assert_eq!(table.max_err_lsb, 1); + assert_eq!(table.rms_lsb, f64::from_bits(1)); + } + #[test] fn a_finite_residual_too_wide_for_i32_is_bound_overflow() { let err = BakeInput::new(