From f6d68855a314fc4802c559080d9faca85db81bb9 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Sat, 11 Jul 2026 13:11:57 -0700 Subject: [PATCH 1/2] perf(matrix): restore det_direct throughput for D=2..4 - Reintroduce branch-free dense expansions for D=2 and D=3. - Share D=4 minors while retaining guarded sparse evaluation. - Preserve non-finite handling for mathematically inactive terms. Closes #158 --- src/lib.rs | 5 ++- src/matrix.rs | 110 ++++++++++++++++++++++++++++++++++++++++++-------- 2 files changed, 97 insertions(+), 18 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 3a6e1ae..b5dbe11 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -341,8 +341,9 @@ pub const ERR_COEFF_3: f64 = 8.0 * EPS + 64.0 * EPS * EPS; /// ``` /// /// where `p(|A|)` is the absolute Leibniz sum. `det_direct` for D=4 -/// evaluates four nested 3×3 cofactors and reduces them with an FMA row -/// combination, yielding the +/// evaluates four nested 3×3 cofactors, sharing their six 2×2 minors when +/// every coefficient in the first two rows is non-zero, and reduces them with +/// an FMA row combination, yielding the /// `12·EPS + 128·EPS²` bound. See `REFERENCES.md` \[8\] for the /// Shewchuk framework these bounds follow. /// diff --git a/src/matrix.rs b/src/matrix.rs index 48bb4cf..b6276c8 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -861,15 +861,10 @@ impl Matrix { let b = self.rows[0][1]; let c = self.rows[1][0]; let d = self.rows[1][1]; - if b == 0.0 { - Some(FilterArithmetic::::multiply(a, d)) - } else { - let subtrahend = FilterArithmetic::::multiply(b, c); - let mut det = - FilterArithmetic::::mul_add(a, d, -subtrahend.value); - det.underflow_safe &= subtrahend.underflow_safe; - Some(det) - } + let subtrahend = FilterArithmetic::::multiply(b, c); + let mut det = FilterArithmetic::::mul_add(a, d, -subtrahend.value); + det.underflow_safe &= subtrahend.underflow_safe; + Some(det) } 3 => Some(Self::det3_elements::( [self.rows[0][0], self.rows[0][1], self.rows[0][2]], @@ -879,6 +874,26 @@ impl Matrix { 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 { + return Some(FilterArithmetic { + value: Self::det4_dense_elements(r), + underflow_safe: true, + }); + } + let mut det = if r[0][3] == 0.0 { FilterArithmetic { value: 0.0, @@ -937,6 +952,36 @@ 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. + #[expect( + clippy::inline_always, + 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 { + 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])); + let s03 = r[2][0].mul_add(r[3][3], -(r[2][3] * r[3][0])); + let s02 = r[2][0].mul_add(r[3][2], -(r[2][2] * r[3][0])); + let s01 = r[2][0].mul_add(r[3][1], -(r[2][1] * r[3][0])); + + let c00 = r[1][1].mul_add(s23, (-r[1][2]).mul_add(s13, r[1][3] * s12)); + let c01 = r[1][0].mul_add(s23, (-r[1][2]).mul_add(s03, r[1][3] * s02)); + let c02 = r[1][0].mul_add(s13, (-r[1][1]).mul_add(s03, r[1][3] * s01)); + let c03 = r[1][0].mul_add(s12, (-r[1][1]).mul_add(s02, r[1][2] * s01)); + + r[0][0].mul_add( + c00, + (-r[0][1]).mul_add(c01, r[0][2].mul_add(c02, -(r[0][3] * c03))), + ) + } + /// Floating-point determinant, using closed-form formulas for D ≤ 4 and /// LU decomposition for D ≥ 5. /// @@ -1283,12 +1328,14 @@ impl Matrix { } } - /// Evaluate a 3×3 determinant expansion while skipping zero coefficients. + /// Evaluate a 3×3 determinant expansion with a guarded sparse fallback. /// - /// This helper protects the public [`det_direct`](Self::det_direct) - /// contract: a mathematically absent term must not compute an overflowing - /// minor and poison the determinant with `0.0 * inf == NaN`. Nonzero terms - /// keep the same fused multiply-add ordering as the closed-form expansion. + /// When all three first-row coefficients are non-zero, one branch-free + /// closed form is used. The sparse fallback protects the public + /// [`det_direct`](Self::det_direct) contract: a mathematically absent term + /// must not compute an overflowing minor and poison the determinant with + /// `0.0 * inf == NaN`. Nonzero terms keep the same fused multiply-add + /// ordering as the closed-form expansion. #[expect( clippy::inline_always, reason = "det_direct callers must eliminate unused filter-safety bookkeeping" @@ -1299,6 +1346,17 @@ 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 m00 = r1[1].mul_add(r2[2], -(r1[2] * r2[1])); + let m01 = r1[0].mul_add(r2[2], -(r1[2] * r2[0])); + let m02 = r1[0].mul_add(r2[1], -(r1[1] * r2[0])); + return FilterArithmetic { + value: r0[0].mul_add(m00, (-r0[1]).mul_add(m01, r0[2] * m02)), + underflow_safe: true, + }; + } + let mut det = if r0[2] == 0.0 { FilterArithmetic { value: 0.0, @@ -1731,13 +1789,17 @@ mod tests { let m = black_box( Matrix::<4>::try_from_rows([ [4.0, 1.0, 3.0, 2.0], - [0.0, 5.0, 2.0, 1.0], + [1.0, 5.0, 2.0, 1.0], [7.0, 2.0, 6.0, 3.0], [1.0, 8.0, 4.0, 9.0], ]) .unwrap(), ); - assert_abs_diff_eq!(m.det_direct().unwrap().unwrap(), 92.0, epsilon = 1e-12); + let direct = m.det_direct().unwrap().unwrap(); + let paired = m.det_direct_with_errbound().unwrap().unwrap(); + + assert_abs_diff_eq!(direct, 112.0, epsilon = 1e-12); + assert_eq!(paired.determinant().to_bits(), direct.to_bits()); } #[test] @@ -1754,6 +1816,22 @@ mod tests { assert_eq!(m.det_direct(), Ok(Some(0.0))); } + #[test] + fn det_direct_d4_sparse_second_row_skips_inactive_overflowing_minors() { + let m = black_box( + Matrix::<4>::try_from_rows([ + [1.0e-300, 1.0, 1.0, 1.0], + [0.0, 1.0, 0.0, 0.0], + [0.0, 1.0e300, 1.0, 1.0e300], + [0.0, 1.0e300, 0.0, 1.0e300], + ]) + .unwrap(), + ); + + assert_eq!(m.det_direct(), Ok(Some(1.0))); + assert_eq!(m.det(), Ok(1.0)); + } + #[test] fn det_direct_d5_returns_none() { assert_eq!(Matrix::<5>::identity().det_direct(), Ok(None)); From 5e15d9250455fd1467e47a735412f65449f30b56 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Sat, 11 Jul 2026 13:44:06 -0700 Subject: [PATCH 2/2] ci(coderabbit): avoid duplicate review checks - Use the legacy required status while preserving automatic approvals. - Add explicit dense D3 determinant coverage and document the D4 error bound. --- .coderabbit.yaml | 12 ++++++++++++ src/matrix.rs | 20 ++++++++++++++++++++ 2 files changed, 32 insertions(+) diff --git a/.coderabbit.yaml b/.coderabbit.yaml index cf4fc40..a6fac5c 100644 --- a/.coderabbit.yaml +++ b/.coderabbit.yaml @@ -1,4 +1,16 @@ # yaml-language-server: $schema=https://coderabbit.ai/integrations/schema.v2.json --- +language: en-US +early_access: false reviews: + profile: chill request_changes_workflow: true + # Publish only the legacy commit status consumed by the main-branch ruleset. + # The review-progress path creates a second GitHub check run for the same review. + review_progress: false + commit_status: true + fail_commit_status: false + tools: + github-checks: + enabled: true + timeout_ms: 90000 diff --git a/src/matrix.rs b/src/matrix.rs index b6276c8..b43e9a7 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -958,6 +958,12 @@ impl Matrix { /// 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. + /// + /// When no intermediate undergoes gradual underflow, the rounding error is + /// bounded by `ERR_COEFF_4 · p(|A|)`, where `p(|A|)` is the absolute Leibniz + /// sum. This helper returns only the determinant; use + /// [`Self::det_errbound`] or [`Self::det_direct_with_errbound`] to obtain the + /// certified bound. #[expect( clippy::inline_always, reason = "the D=4 determinant hot path must inline its shared-minor expansion" @@ -1749,6 +1755,20 @@ mod tests { assert_abs_diff_eq!(m.det_direct().unwrap().unwrap(), 0.0, epsilon = 1e-12); } + #[test] + fn det_direct_d3_dense_known_value() { + // det = 1*(5*8 - 7*6) - 2*(4*8 - 7*2) + 3*(4*6 - 5*2) = 4 + let m = black_box( + Matrix::<3>::try_from_rows([[1.0, 2.0, 3.0], [4.0, 5.0, 7.0], [2.0, 6.0, 8.0]]) + .unwrap(), + ); + let direct = m.det_direct().unwrap().unwrap(); + let paired = m.det_direct_with_errbound().unwrap().unwrap(); + + assert_abs_diff_eq!(direct, 4.0, epsilon = 1e-12); + assert_eq!(paired.determinant().to_bits(), direct.to_bits()); + } + #[test] fn det_direct_d3_nonsingular() { // [[2,1,0],[0,3,1],[1,0,2]] → det = 2*(6-0) - 1*(0-1) + 0 = 13