Apply RMS exponent without underflowing 2^exp first

CHANGELOG.md

@@ -34,8 +34,11 @@
   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` 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`,

crates/surfaces-bake/src/bound.rs

@@ -45,10 +45,12 @@ 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. 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;
@@ -71,7 +73,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 +93,32 @@ 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;
+    }
+    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);
+        exp -= 1023;
+        if !x.is_finite() {
+            return x;
+        }
+    }
+    x * 2.0_f64.powi(exp)
+}
+
 pub(crate) fn ratio_from_f64(x: f64) -> Option<Ratio> {
     Some(Ratio(dyadic(x)?))
 }
@@ -139,7 +167,7 @@ fn dyadic(x: f64) -> Option<BigRational> {
 
 #[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() {
@@ -190,4 +218,34 @@ 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));
+    }
+
+    #[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));
+    }
 }

crates/surfaces-bake/src/quantize.rs

@@ -871,6 +871,56 @@ 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_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(