From b8bfa9fa0c908ea034d294ceb123c8191eff3395 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Sat, 11 Jul 2026 15:35:34 -0700 Subject: [PATCH] perf(exact): restore small-matrix det_sign_exact throughput - Reuse proof-bearing shared minors for D4 determinant and permanent evaluation. - Restore the dense D3 filter while preserving sparse, overflow, and underflow fallbacks. - Require headline benchmarks to exercise the intended filter and document historical harness overhead. Closes #157 --- benches/common/exact.rs | 4 + docs/BENCHMARKING.md | 10 ++ src/matrix.rs | 305 ++++++++++++++++++++++++++---------- tests/exact_bench_config.rs | 23 +++ 4 files changed, 256 insertions(+), 86 deletions(-) diff --git a/benches/common/exact.rs b/benches/common/exact.rs index a21b3f3..328bddb 100644 --- a/benches/common/exact.rs +++ b/benches/common/exact.rs @@ -529,6 +529,10 @@ pub fn validate_f64_determinant_benchmarks(input: &ValidatedExac panic!("the baseline fixture must have a certified D={D} determinant bound"); }; assert_eq!(estimate.determinant().to_bits(), direct.to_bits()); + assert!( + estimate.determinant().abs() > estimate.absolute_error_bound(), + "the headline D={D} det_sign_exact benchmark must exercise the fast filter", + ); let observed_error = (rational_from_f64(direct) - &exact).abs(); let certified_bound = rational_from_f64(estimate.absolute_error_bound()); assert!( diff --git a/docs/BENCHMARKING.md b/docs/BENCHMARKING.md index d61ee40..d1db7d5 100644 --- a/docs/BENCHMARKING.md +++ b/docs/BENCHMARKING.md @@ -352,6 +352,16 @@ claim that every timed sample is revalidated. Criterion closures remain free of oracle work so their measurements cover only the named operation; the operation is deterministic for the already-validated input. +The original v0.4.2-to-v0.4.3 report is not library-only evidence for its +headline `det_sign_exact` rows. The v0.4.3 harness routed those calls through an +operation enum and constructed a complete matrix/RHS input inside each timed +iteration, while the saved v0.4.2 samples used direct closures. The current +shared harness borrows a prevalidated input for both revisions. It also verifies +that the D=2–4 headline fixtures resolve through the floating-point filter, so +those rows continue to measure the intended common path. Use a current +shared-harness comparison before attributing the historical D=2 or D=3 changes +to the library implementation. + For exact-arithmetic comparisons against v0.4.2 or older baselines, rows such as `det_exact_rounded_f64 (vs det_exact_f64)` mean the current rounded API is being compared to the historical lossy `*_exact_f64` benchmark. Rows such as diff --git a/src/matrix.rs b/src/matrix.rs index b43e9a7..2c209c1 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -155,6 +155,52 @@ impl FilterArithmetic { } } +/// A finite D=4 matrix proven safe for shared-minor determinant and permanent +/// evaluation. +/// +/// Construction proves both the fixed dimension and that every coefficient in +/// the first two rows is non-zero. The latter makes every shared 2×2 minor part +/// of an active Leibniz term, so the dense kernel cannot evaluate an overflowing +/// minor solely for a mathematically absent term. +#[repr(transparent)] +#[derive(Clone, Copy)] +struct Det4SharedMinorInput<'a, const D: usize> { + matrix: &'a Matrix, +} + +impl<'a, const D: usize> Det4SharedMinorInput<'a, D> { + /// Parse a matrix into the shared-minor D=4 domain. + /// + /// `None` selects the guarded determinant path; it does not represent an + /// invalid public matrix. + #[expect( + clippy::inline_always, + reason = "the D=4 determinant hot path must eliminate its proof wrapper" + )] + #[inline(always)] + const fn try_new(matrix: &'a Matrix) -> Option { + if D != 4 { + return None; + } + + let r = &matrix.rows; + let shared_minors_are_active = (r[0][0] != 0.0) + && (r[0][1] != 0.0) + && (r[0][2] != 0.0) + && (r[0][3] != 0.0) + && (r[1][0] != 0.0) + && (r[1][1] != 0.0) + && (r[1][2] != 0.0) + && (r[1][3] != 0.0); + + if shared_minors_are_active { + Some(Self { matrix }) + } else { + None + } + } +} + impl SymmetricMatrix { /// Construct a symmetric matrix proof without checking the invariant. /// @@ -872,28 +918,14 @@ impl Matrix { [self.rows[2][0], self.rows[2][1], self.rows[2][2]], )), 4 => { - let r = &self.rows; - - // The dense hot path shares the six unique 2×2 minors across - // all four cofactors. Requiring non-zero coefficients in the - // first two rows keeps this equivalent to the guarded path - // below: no mathematically absent term can force evaluation of - // an overflowing minor. - let dense = (r[0][0] != 0.0) - && (r[0][1] != 0.0) - && (r[0][2] != 0.0) - && (r[0][3] != 0.0) - && (r[1][0] != 0.0) - && (r[1][1] != 0.0) - && (r[1][2] != 0.0) - && (r[1][3] != 0.0); - if !TRACK_UNDERFLOW && dense { + if !TRACK_UNDERFLOW && let Some(input) = Det4SharedMinorInput::try_new(self) { return Some(FilterArithmetic { - value: Self::det4_dense_elements(r), + value: Self::det4_dense_elements(input), underflow_safe: true, }); } + let r = &self.rows; let mut det = if r[0][3] == 0.0 { FilterArithmetic { value: 0.0, @@ -952,12 +984,7 @@ impl Matrix { } } - /// Evaluate the dense 4×4 cofactor expansion with shared 2×2 minors. - /// - /// Callers must establish that the first two rows contain no zero - /// coefficients. That proof makes every minor part of an active Leibniz - /// term, so evaluating all six shared minors cannot introduce a non-finite - /// result through a mathematically absent term. + /// Evaluate the proof-bearing 4×4 cofactor expansion with shared 2×2 minors. /// /// When no intermediate undergoes gradual underflow, the rounding error is /// bounded by `ERR_COEFF_4 · p(|A|)`, where `p(|A|)` is the absolute Leibniz @@ -969,7 +996,8 @@ impl Matrix { reason = "the D=4 determinant hot path must inline its shared-minor expansion" )] #[inline(always)] - const fn det4_dense_elements(r: &[[f64; D]; D]) -> f64 { + const fn det4_dense_elements(input: Det4SharedMinorInput<'_, D>) -> f64 { + let r = &input.matrix.rows; let s23 = r[2][2].mul_add(r[3][3], -(r[2][3] * r[3][2])); let s13 = r[2][1].mul_add(r[3][3], -(r[2][3] * r[3][1])); let s12 = r[2][1].mul_add(r[3][2], -(r[2][2] * r[3][1])); @@ -988,6 +1016,47 @@ impl Matrix { ) } + /// Evaluate the dense 4×4 absolute permanent with shared 2×2 minors. + /// + /// The proof carried by `input` makes every shared minor part of an active + /// Leibniz term. The caller separately establishes a wide exponent margin, + /// so this branch-free kernel cannot hide gradual underflow or evaluate an + /// overflowing minor for a mathematically absent term. + #[expect( + clippy::inline_always, + reason = "the D=4 determinant filter must inline its shared-minor permanent" + )] + #[inline(always)] + const fn det4_dense_abs_permanent_elements(input: Det4SharedMinorInput<'_, D>) -> f64 { + let r = &input.matrix.rows; + let sp23 = (r[2][2] * r[3][3]).abs() + (r[2][3] * r[3][2]).abs(); + let sp13 = (r[2][1] * r[3][3]).abs() + (r[2][3] * r[3][1]).abs(); + let sp12 = (r[2][1] * r[3][2]).abs() + (r[2][2] * r[3][1]).abs(); + let sp03 = (r[2][0] * r[3][3]).abs() + (r[2][3] * r[3][0]).abs(); + let sp02 = (r[2][0] * r[3][2]).abs() + (r[2][2] * r[3][0]).abs(); + let sp01 = (r[2][0] * r[3][1]).abs() + (r[2][1] * r[3][0]).abs(); + + let pc0 = r[1][3] + .abs() + .mul_add(sp12, r[1][2].abs().mul_add(sp13, r[1][1].abs() * sp23)); + let pc1 = r[1][3] + .abs() + .mul_add(sp02, r[1][2].abs().mul_add(sp03, r[1][0].abs() * sp23)); + let pc2 = r[1][3] + .abs() + .mul_add(sp01, r[1][1].abs().mul_add(sp03, r[1][0].abs() * sp13)); + let pc3 = r[1][2] + .abs() + .mul_add(sp01, r[1][1].abs().mul_add(sp02, r[1][0].abs() * sp12)); + + r[0][3].abs().mul_add( + pc3, + r[0][2] + .abs() + .mul_add(pc2, r[0][1].abs().mul_add(pc1, r[0][0].abs() * pc0)), + ) + } + /// Floating-point determinant, using closed-form formulas for D ≤ 4 and /// LU decomposition for D ≥ 5. /// @@ -1265,68 +1334,7 @@ impl Matrix { ); Self::certified_error_bound(ERR_COEFF_3, permanent) } - 4 => { - let r = &self.rows; - let mut permanent = if r[0][3] == 0.0 { - FilterArithmetic { - value: 0.0, - underflow_safe: true, - } - } else { - let pc3 = Self::det3_abs_permanent_elements::( - [r[1][0], r[1][1], r[1][2]], - [r[2][0], r[2][1], r[2][2]], - [r[3][0], r[3][1], r[3][2]], - ); - let mut term = - FilterArithmetic::::multiply(r[0][3].abs(), pc3.value); - term.underflow_safe &= pc3.underflow_safe; - term - }; - if r[0][2] != 0.0 { - let pc2 = Self::det3_abs_permanent_elements::( - [r[1][0], r[1][1], r[1][3]], - [r[2][0], r[2][1], r[2][3]], - [r[3][0], r[3][1], r[3][3]], - ); - let prior_safe = permanent.underflow_safe && pc2.underflow_safe; - permanent = FilterArithmetic::::mul_add( - r[0][2].abs(), - pc2.value, - permanent.value, - ); - permanent.underflow_safe &= prior_safe; - } - if r[0][1] != 0.0 { - let pc1 = Self::det3_abs_permanent_elements::( - [r[1][0], r[1][2], r[1][3]], - [r[2][0], r[2][2], r[2][3]], - [r[3][0], r[3][2], r[3][3]], - ); - let prior_safe = permanent.underflow_safe && pc1.underflow_safe; - permanent = FilterArithmetic::::mul_add( - r[0][1].abs(), - pc1.value, - permanent.value, - ); - permanent.underflow_safe &= prior_safe; - } - if r[0][0] != 0.0 { - let pc0 = Self::det3_abs_permanent_elements::( - [r[1][1], r[1][2], r[1][3]], - [r[2][1], r[2][2], r[2][3]], - [r[3][1], r[3][2], r[3][3]], - ); - let prior_safe = permanent.underflow_safe && pc0.underflow_safe; - permanent = FilterArithmetic::::mul_add( - r[0][0].abs(), - pc0.value, - permanent.value, - ); - permanent.underflow_safe &= prior_safe; - } - Self::certified_error_bound(ERR_COEFF_4, permanent) - } + 4 => self.det4_errbound::(), _ => { cold_path(); Ok(None) @@ -1334,6 +1342,79 @@ impl Matrix { } } + /// Compute the D=4 determinant error bound after the dimension dispatch. + const fn det4_errbound(&self) -> Result, LaError> { + if !TRACK_UNDERFLOW && let Some(input) = Det4SharedMinorInput::try_new(self) { + return Self::certified_error_bound( + ERR_COEFF_4, + FilterArithmetic:: { + value: Self::det4_dense_abs_permanent_elements(input), + underflow_safe: true, + }, + ); + } + + let r = &self.rows; + let mut permanent = if r[0][3] == 0.0 { + FilterArithmetic { + value: 0.0, + underflow_safe: true, + } + } else { + let pc3 = Self::det3_abs_permanent_elements::( + [r[1][0], r[1][1], r[1][2]], + [r[2][0], r[2][1], r[2][2]], + [r[3][0], r[3][1], r[3][2]], + ); + let mut term = FilterArithmetic::::multiply(r[0][3].abs(), pc3.value); + term.underflow_safe &= pc3.underflow_safe; + term + }; + if r[0][2] != 0.0 { + let pc2 = Self::det3_abs_permanent_elements::( + [r[1][0], r[1][1], r[1][3]], + [r[2][0], r[2][1], r[2][3]], + [r[3][0], r[3][1], r[3][3]], + ); + let prior_safe = permanent.underflow_safe && pc2.underflow_safe; + permanent = FilterArithmetic::::mul_add( + r[0][2].abs(), + pc2.value, + permanent.value, + ); + permanent.underflow_safe &= prior_safe; + } + if r[0][1] != 0.0 { + let pc1 = Self::det3_abs_permanent_elements::( + [r[1][0], r[1][2], r[1][3]], + [r[2][0], r[2][2], r[2][3]], + [r[3][0], r[3][2], r[3][3]], + ); + let prior_safe = permanent.underflow_safe && pc1.underflow_safe; + permanent = FilterArithmetic::::mul_add( + r[0][1].abs(), + pc1.value, + permanent.value, + ); + permanent.underflow_safe &= prior_safe; + } + if r[0][0] != 0.0 { + let pc0 = Self::det3_abs_permanent_elements::( + [r[1][1], r[1][2], r[1][3]], + [r[2][1], r[2][2], r[2][3]], + [r[3][1], r[3][2], r[3][3]], + ); + let prior_safe = permanent.underflow_safe && pc0.underflow_safe; + permanent = FilterArithmetic::::mul_add( + r[0][0].abs(), + pc0.value, + permanent.value, + ); + permanent.underflow_safe &= prior_safe; + } + Self::certified_error_bound(ERR_COEFF_4, permanent) + } + /// Evaluate a 3×3 determinant expansion with a guarded sparse fallback. /// /// When all three first-row coefficients are non-zero, one branch-free @@ -1413,6 +1494,19 @@ impl Matrix { r1: [f64; 3], r2: [f64; 3], ) -> FilterArithmetic { + let dense = (r0[0] != 0.0) && (r0[1] != 0.0) && (r0[2] != 0.0); + if !TRACK_UNDERFLOW && dense { + let pm00 = (r1[1] * r2[2]).abs() + (r1[2] * r2[1]).abs(); + let pm01 = (r1[0] * r2[2]).abs() + (r1[2] * r2[0]).abs(); + let pm02 = (r1[0] * r2[1]).abs() + (r1[1] * r2[0]).abs(); + return FilterArithmetic { + value: r0[2] + .abs() + .mul_add(pm02, r0[1].abs().mul_add(pm01, r0[0].abs() * pm00)), + underflow_safe: true, + }; + } + let mut permanent = if r0[2] == 0.0 { FilterArithmetic { value: 0.0, @@ -1769,6 +1863,25 @@ mod tests { assert_eq!(paired.determinant().to_bits(), direct.to_bits()); } + #[test] + fn det_direct_d3_dense_reports_legitimate_overflow() { + // The unscaled matrix has determinant 54, so scaling every entry by + // 1.6e102 gives a determinant of approximately 2.21e308. + let scale = 1.6e102; + let m = black_box( + Matrix::<3>::try_from_rows([ + [4.0 * scale, scale, scale], + [scale, 4.0 * scale, scale], + [scale, scale, 4.0 * scale], + ]) + .unwrap(), + ); + let expected = LaError::non_finite_computation_scalar(ArithmeticOperation::Determinant); + + assert_eq!(m.det_direct(), Err(expected)); + assert_eq!(m.det(), Err(expected)); + } + #[test] fn det_direct_d3_nonsingular() { // [[2,1,0],[0,3,1],[1,0,2]] → det = 2*(6-0) - 1*(0-1) + 0 = 13 @@ -1822,6 +1935,26 @@ mod tests { assert_eq!(paired.determinant().to_bits(), direct.to_bits()); } + #[test] + fn det_direct_d4_dense_reports_legitimate_overflow() { + // The unscaled matrix has determinant 189, so scaling every entry by + // 3.2e76 gives a determinant of approximately 1.98e308. + let scale = 3.2e76; + let m = black_box( + Matrix::<4>::try_from_rows([ + [4.0 * scale, scale, scale, scale], + [scale, 4.0 * scale, scale, scale], + [scale, scale, 4.0 * scale, scale], + [scale, scale, scale, 4.0 * scale], + ]) + .unwrap(), + ); + let expected = LaError::non_finite_computation_scalar(ArithmeticOperation::Determinant); + + assert_eq!(m.det_direct(), Err(expected)); + assert_eq!(m.det(), Err(expected)); + } + #[test] fn det_direct_d4_skips_zero_coefficient_cofactors_that_would_overflow() { let m = black_box( diff --git a/tests/exact_bench_config.rs b/tests/exact_bench_config.rs index a1e419c..501c32f 100644 --- a/tests/exact_bench_config.rs +++ b/tests/exact_bench_config.rs @@ -37,6 +37,16 @@ fn validate_baseline_and_random_corpus() { } } +/// Report whether `det_sign_exact` can certify this fixture through its direct filter. +/// +/// Both a missing bound and a direct-path error select the exact fallback. +fn fast_filter_is_conclusive(input: ExactInput) -> bool { + match input.matrix.det_direct_with_errbound() { + Ok(Some(estimate)) => estimate.determinant().abs() > estimate.absolute_error_bound(), + Ok(None) | Err(_) => false, + } +} + #[test] fn i16_range_rejects_unordered_bounds() { let Err(err) = I16Range::try_new(5, 4) else { @@ -138,6 +148,19 @@ fn exact_adversarial_benchmark_fixtures_are_correct() { let _ = validate_exact_fixture(hilbert_input::<5>()); } +#[test] +fn det_sign_exact_benchmark_filter_paths_are_stable() { + assert!(fast_filter_is_conclusive(baseline_input::<2>())); + assert!(fast_filter_is_conclusive(baseline_input::<3>())); + assert!(fast_filter_is_conclusive(baseline_input::<4>())); + assert!(!fast_filter_is_conclusive(baseline_input::<5>())); + + assert!(!fast_filter_is_conclusive(near_singular_3x3_input())); + assert!(!fast_filter_is_conclusive(large_entries_3x3_input())); + assert!(fast_filter_is_conclusive(hilbert_input::<4>())); + assert!(!fast_filter_is_conclusive(hilbert_input::<5>())); +} + #[test] fn validated_fixture_exposes_only_checked_inputs() { let raw = baseline_input::<3>();