diff --git a/CHANGELOG.md b/CHANGELOG.md index d5b12b8..7f65eaa 100644 --- a/CHANGELOG.md +++ b/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`, diff --git a/crates/surfaces-bake/src/bound.rs b/crates/surfaces-bake/src/bound.rs index 15c74ee..1bcde54 100644 --- a/crates/surfaces-bake/src/bound.rs +++ b/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 { Some(Ratio(dyadic(x)?)) } @@ -139,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() { @@ -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)); + } } diff --git a/crates/surfaces-bake/src/quantize.rs b/crates/surfaces-bake/src/quantize.rs index ab0f572..74a256d 100644 --- a/crates/surfaces-bake/src/quantize.rs +++ b/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(