diff --git a/justfile b/justfile index 1f55b67..30c63ae 100644 --- a/justfile +++ b/justfile @@ -24,7 +24,7 @@ clippy_sarif_version := "0.8.0" dprint_version := "0.55.2" git_cliff_version := "2.13.1" just_version := "1.58.0" -rumdl_version := "0.2.52" +rumdl_version := "0.2.53" sarif_fmt_version := "0.8.0" taplo_version := "0.10.0" typos_version := "1.49.0" diff --git a/src/ldlt.rs b/src/ldlt.rs index ae4afa1..0508b44 100644 --- a/src/ldlt.rs +++ b/src/ldlt.rs @@ -20,7 +20,7 @@ use core::hint::cold_path; use crate::matrix::SymmetricMatrix; -use crate::scaled_product::{RangeCheckedProduct, ScaledProduct, range_checked_product}; +use crate::scaled_product::{ScaledProduct, range_checked_product}; use crate::vector::Vector; use crate::{ArithmeticOperation, FactorizationKind, LaError, Tolerance}; @@ -276,18 +276,36 @@ impl Ldlt { pub const fn det(&self) -> Result { let mut det = 1.0; let mut i = 0; - while i < D { - let factor = self.factors.diag(i); - match range_checked_product(det, factor) { - RangeCheckedProduct::Safe(next) => det = next, - RangeCheckedProduct::NeedsScaling => { + + if D <= 4 { + // Tiny determinants compose better with factorization when range + // loss exits immediately instead of carrying an aggregate proof. + while i < D { + let step = range_checked_product(det, self.factors.diag(i)); + if !step.range_preserved() { cold_path(); return self.scaled_det(); } + det = step.product(); + i += 1; } + return Ok(det); + } + + let mut range_preserved = true; + while i < D { + let factor = self.factors.diag(i); + let step = range_checked_product(det, factor); + det = step.product(); + range_preserved &= step.range_preserved(); i += 1; } - Ok(det) + if range_preserved { + Ok(det) + } else { + cold_path(); + self.scaled_det() + } } /// Recompute the determinant with normalized mantissa/exponent scaling. diff --git a/src/lu.rs b/src/lu.rs index 8e799a1..665431d 100644 --- a/src/lu.rs +++ b/src/lu.rs @@ -10,7 +10,7 @@ use core::hint::cold_path; use crate::matrix::Matrix; -use crate::scaled_product::{RangeCheckedProduct, ScaledProduct, range_checked_product}; +use crate::scaled_product::{ScaledProduct, range_checked_product}; use crate::vector::Vector; use crate::{ArithmeticOperation, FactorizationKind, LaError, Tolerance}; @@ -364,18 +364,36 @@ impl Lu { pub const fn det(&self) -> Result { let mut det = if self.permutation.is_odd() { -1.0 } else { 1.0 }; let mut i = 0; - while i < D { - let factor = self.factors.diag(i); - match range_checked_product(det, factor) { - RangeCheckedProduct::Safe(next) => det = next, - RangeCheckedProduct::NeedsScaling => { + + if D <= 4 { + // Tiny determinants compose better with factorization when range + // loss exits immediately instead of carrying an aggregate proof. + while i < D { + let step = range_checked_product(det, self.factors.diag(i)); + if !step.range_preserved() { cold_path(); return self.scaled_det(); } + det = step.product(); + i += 1; } + return Ok(det); + } + + let mut range_preserved = true; + while i < D { + let factor = self.factors.diag(i); + let step = range_checked_product(det, factor); + det = step.product(); + range_preserved &= step.range_preserved(); i += 1; } - Ok(det) + if range_preserved { + Ok(det) + } else { + cold_path(); + self.scaled_det() + } } /// Recompute the determinant with normalized mantissa/exponent scaling. diff --git a/src/scaled_product.rs b/src/scaled_product.rs index 65e2b76..0a2397d 100644 --- a/src/scaled_product.rs +++ b/src/scaled_product.rs @@ -10,24 +10,35 @@ const EXPONENT_BIAS: i128 = 1023; const MIN_NORMAL_EXPONENT: i128 = -1022; const MIN_SUBNORMAL_EXPONENT: i128 = -1074; -/// Result of multiplying one direct product step with its range proof attached. +/// One direct product step with its range proof attached. #[derive(Clone, Copy, Debug, PartialEq)] -pub(crate) enum RangeCheckedProduct { - /// The value is finite and normal, or is an exact zero caused by a zero - /// operand, so direct accumulation may continue. - Safe(f64), - /// Direct multiplication overflowed or lost range through gradual - /// underflow, so all factors must be recomputed with scaling. - NeedsScaling, +pub(crate) struct RangeCheckedProduct { + product: f64, + range_preserved: bool, } -/// Multiply one direct product step and retain only values proven safe for -/// sequential accumulation. +impl RangeCheckedProduct { + /// Return the directly accumulated product. + #[inline] + pub(crate) const fn product(self) -> f64 { + self.product + } + + /// Return whether direct accumulation preserved the non-zero product's range. + #[inline] + pub(crate) const fn range_preserved(self) -> bool { + self.range_preserved + } +} + +/// Multiply one non-zero direct product step and attach its range proof. /// -/// The common normal case needs only the result's binary exponent field. Zero -/// operands are checked only when that field is zero, preserving signed-zero -/// multiplication without treating underflow from two non-zero operands as an -/// exact zero. +/// Successful LU and LDLT construction proves every diagonal factor is finite +/// and non-zero. Starting from `±1.0`, a normal result therefore proves that +/// direct accumulation has not overflowed or lost range through gradual +/// underflow. Callers combine every step's proof and replay the complete product +/// with scaling if any step fails, keeping the success path branch-free between +/// factors. #[inline] pub(crate) const fn range_checked_product(accumulator: f64, factor: f64) -> RangeCheckedProduct { let product = accumulator * factor; @@ -35,14 +46,9 @@ pub(crate) const fn range_checked_product(accumulator: f64, factor: f64) -> Rang // Subtracting one maps the valid normal fields 1..=0x7fe to // 0..=0x7fd. Zero wraps high and 0x7ff maps to the exclusive upper // bound, so the common normal path needs one unsigned comparison. - if product_exponent.wrapping_sub(1) < EXPONENT_MASK - 1 { - return RangeCheckedProduct::Safe(product); - } - - if product_exponent == 0 && (accumulator == 0.0 || factor == 0.0) { - RangeCheckedProduct::Safe(product) - } else { - RangeCheckedProduct::NeedsScaling + RangeCheckedProduct { + product, + range_preserved: product_exponent.wrapping_sub(1) < EXPONENT_MASK - 1, } } @@ -216,7 +222,7 @@ impl ScaledProduct { #[cfg(test)] mod tests { - use super::{RangeCheckedProduct, SIGN_MASK, ScaledProduct, range_checked_product}; + use super::{SIGN_MASK, ScaledProduct, range_checked_product}; const TWO_NEG_800: f64 = f64::from_bits(223_u64 << 52); const TWO_POS_800: f64 = f64::from_bits(1823_u64 << 52); @@ -334,10 +340,7 @@ mod tests { product.multiply(factor); } - assert_eq!( - range_checked_product(TWO_NEG_800, TWO_NEG_800), - RangeCheckedProduct::NeedsScaling - ); + assert!(!range_checked_product(TWO_NEG_800, TWO_NEG_800).range_preserved()); assert_eq!(product.finish().map(f64::to_bits), Some(1)); } @@ -356,26 +359,17 @@ mod tests { } #[test] - fn direct_product_range_check_distinguishes_exact_zero_from_range_loss() { - assert_eq!( - range_checked_product(1.5, 2.0), - RangeCheckedProduct::Safe(3.0) - ); - assert_eq!( - range_checked_product(-0.0, -2.0), - RangeCheckedProduct::Safe(0.0) - ); - assert_eq!( + fn direct_product_range_proof_distinguishes_ordinary_values_from_range_loss() { + let ordinary = range_checked_product(1.5, 2.0); + assert_eq!(ordinary.product().to_bits(), 3.0_f64.to_bits()); + assert!(ordinary.range_preserved()); + + for step in [ range_checked_product(TWO_NEG_800, TWO_NEG_800), - RangeCheckedProduct::NeedsScaling - ); - assert_eq!( range_checked_product(f64::MIN_POSITIVE, 0.5), - RangeCheckedProduct::NeedsScaling - ); - assert_eq!( range_checked_product(TWO_POS_800, TWO_POS_800), - RangeCheckedProduct::NeedsScaling - ); + ] { + assert!(!step.range_preserved()); + } } }