perf: reduce checked determinant query overhead

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"

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<const D: usize> Ldlt<D> {
     pub const fn det(&self) -> Result<f64, LaError> {
         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.

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<const D: usize> Lu<D> {
     pub const fn det(&self) -> Result<f64, LaError> {
         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.

src/scaled_product.rs

@@ -10,39 +10,45 @@ 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;
     let product_exponent = (product.to_bits() >> FRACTION_BITS) & EXPONENT_MASK;
     // 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());
+        }
     }
 }