From 81ecb35bdaf159f1f44d1eb24274ecf82c6567d5 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Mon, 20 Apr 2026 14:11:22 -0700 Subject: [PATCH 1/3] feat: const-ify Lu/Ldlt det + solve_vec and Matrix inf_norm + det_errbound MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Make six public methods `const fn`, completing const-evaluation parity with `Matrix::det_direct` and `Vector::dot` now that MSRV 1.95 exposes `f64::mul_add` as const (1.94) and `core::hint::cold_path` as const (1.95): - `Lu::det`, `Lu::solve_vec` - `Ldlt::det`, `Ldlt::solve_vec` - `Matrix::inf_norm`, `Matrix::det_errbound` Iterator chains (`.iter().map().sum()`, `.enumerate().take(i)`, `.enumerate().skip(i + 1)`, `.iter_mut().enumerate().take(D)`) were rewritten as `while` loops since they are not const-stable. Fix error-variant correctness in both solve_vec paths: - A corrupt stored `U` / `D` diagonal at `(i, i)` now surfaces as `LaError::NonFinite { row: Some(i), col: i }`, matching the convention used by `Matrix::det`, `Lu::factor`, and `Ldlt::factor`. - A computed-intermediate overflow keeps `row: None, col: i`. - Previously both were conflated into `row: None, col: i`, defeating debuggability for callers who construct factorizations directly. Sharpen `LaError::NonFinite` variant docs and the `# Errors` sections of `Lu::solve_vec` / `Ldlt::solve_vec` to spell out the `(row, col)` contract for each failure mode. Cleanup: - Import `ERR_COEFF_{2,3,4}` at `src/matrix.rs` top; drop `crate::` prefixes inside `Matrix::det_errbound`. - Rename intermediate bindings `q` / `v` → `quotient` to stay under `clippy::many_single_char_names`. - Rename const-evaluability tests from `*_is_const_evaluable_*` to `*_const_eval_*` for a concise, consistent style. Tests: - Added `{lu,ldlt,inf_norm,det_errbound}_const_eval_*` tests that force compile-time evaluation inside `const { … }` initializers. - Added `solve_vec_defensive_non_finite_diagonal_{2,3,4,5}d` in `src/lu.rs` covering the new split error path (previously only the `Singular` defensive path was exercised). - Updated matching `Ldlt` defensive tests to assert the new `row: Some(D - 1), col: D - 1` payload. Validation: `just ci` passes — 164 unit tests, 324 integration/feature tests, 28 doctests; clippy pedantic + nursery as errors; all examples. Co-Authored-By: Oz --- src/ldlt.rs | 130 +++++++++++++++++++++++++++++++----- src/lib.rs | 17 ++++- src/lu.rs | 181 ++++++++++++++++++++++++++++++++++++++++++-------- src/matrix.rs | 70 ++++++++++++++++--- 4 files changed, 342 insertions(+), 56 deletions(-) diff --git a/src/ldlt.rs b/src/ldlt.rs index 965ab2d..d2bbf84 100644 --- a/src/ldlt.rs +++ b/src/ldlt.rs @@ -124,10 +124,12 @@ impl Ldlt { /// ``` #[inline] #[must_use] - pub fn det(&self) -> f64 { + pub const fn det(&self) -> f64 { let mut det = 1.0; - for i in 0..D { + let mut i = 0; + while i < D { det *= self.factors.rows[i][i]; + i += 1; } det } @@ -154,57 +156,82 @@ impl Ldlt { /// # Errors /// Returns [`LaError::Singular`] if a diagonal entry `d = D[i,i]` satisfies `d <= tol` /// (non-positive or too small), where `tol` is the tolerance that was used during factorization. - /// Returns [`LaError::NonFinite`] if NaN/∞ is detected. + /// + /// Returns [`LaError::NonFinite`] if NaN/∞ is detected. The `row`/`col` coordinates + /// follow the convention documented on [`LaError::NonFinite`]: + /// + /// - `row: Some(i), col: i` — the stored `D` diagonal at `(i, i)` is non-finite + /// (only reachable via direct `Ldlt` construction; [`Matrix::ldlt`](crate::Matrix::ldlt) + /// rejects such factorizations). + /// - `row: None, col: i` — a computed intermediate (forward/back-substitution + /// accumulator or the quotient `x[i] / diag`) overflowed to NaN/∞ at step `i`. #[inline] - pub fn solve_vec(&self, b: Vector) -> Result, LaError> { + pub const fn solve_vec(&self, b: Vector) -> Result, LaError> { let mut x = b.data; // Forward substitution: L y = b (L has unit diagonal). - for i in 0..D { + let mut i = 0; + while i < D { let mut sum = x[i]; let row = self.factors.rows[i]; - for (j, x_j) in x.iter().enumerate().take(i) { - sum = (-row[j]).mul_add(*x_j, sum); + let mut j = 0; + while j < i { + sum = (-row[j]).mul_add(x[j], sum); + j += 1; } if !sum.is_finite() { cold_path(); return Err(LaError::NonFinite { row: None, col: i }); } x[i] = sum; + i += 1; } // Diagonal solve: D z = y. - for (i, x_i) in x.iter_mut().enumerate().take(D) { + let mut i = 0; + while i < D { let diag = self.factors.rows[i][i]; + // A corrupt stored diagonal is a specific matrix cell (i, i), + // distinct from a computed overflow — report it with + // `row: Some(i)` per the `LaError::NonFinite` convention used by + // `Matrix::det`, `Lu::factor`, and `Ldlt::factor`. if !diag.is_finite() { cold_path(); - return Err(LaError::NonFinite { row: None, col: i }); + return Err(LaError::NonFinite { + row: Some(i), + col: i, + }); } if diag <= self.tol { cold_path(); return Err(LaError::Singular { pivot_col: i }); } - let v = *x_i / diag; - if !v.is_finite() { + let quotient = x[i] / diag; + if !quotient.is_finite() { cold_path(); return Err(LaError::NonFinite { row: None, col: i }); } - *x_i = v; + x[i] = quotient; + i += 1; } // Back substitution: Lᵀ x = z. - for ii in 0..D { + let mut ii = 0; + while ii < D { let i = D - 1 - ii; let mut sum = x[i]; - for (j, x_j) in x.iter().enumerate().skip(i + 1) { - sum = (-self.factors.rows[j][i]).mul_add(*x_j, sum); + let mut j = i + 1; + while j < D { + sum = (-self.factors.rows[j][i]).mul_add(x[j], sum); + j += 1; } if !sum.is_finite() { cold_path(); return Err(LaError::NonFinite { row: None, col: i }); } x[i] = sum; + ii += 1; } Ok(Vector::new(x)) @@ -487,7 +514,9 @@ mod tests { paste! { /// `solve_vec` must surface `NonFinite` when a stored /// diagonal is NaN, even though `factor` cannot produce - /// such a factorization. + /// such a factorization. The error must pinpoint the + /// corrupt cell at `(D-1, D-1)` per the + /// [`LaError::NonFinite`] convention. #[test] fn []() { let mut factors = Matrix::<$d>::identity(); @@ -498,7 +527,13 @@ mod tests { }; let b = Vector::<$d>::new([1.0; $d]); let err = ldlt.solve_vec(b).unwrap_err(); - assert_eq!(err, LaError::NonFinite { row: None, col: $d - 1 }); + assert_eq!( + err, + LaError::NonFinite { + row: Some($d - 1), + col: $d - 1, + } + ); } /// `solve_vec` must surface `Singular` when a stored @@ -524,4 +559,65 @@ mod tests { gen_solve_vec_defensive_tests!(3); gen_solve_vec_defensive_tests!(4); gen_solve_vec_defensive_tests!(5); + + // ----------------------------------------------------------------------- + // Const-evaluability tests. + // + // These prove that `Ldlt::det` and `Ldlt::solve_vec` are truly `const fn` + // by forcing the compiler to evaluate them inside a `const` initializer. + // `Ldlt::factor` is not (yet) `const fn` because the rank-1 update loop + // uses array indexing patterns that still require non-const helpers on + // some toolchains; we therefore construct `Ldlt` directly. + // ----------------------------------------------------------------------- + + #[test] + fn ldlt_det_const_eval_d2() { + const DET: f64 = { + // Diagonal D = [4.0, 0.25] ⇒ det = 1.0. + let mut factors = Matrix::<2>::identity(); + factors.rows[0][0] = 4.0; + factors.rows[1][1] = 0.25; + let ldlt = Ldlt::<2> { + factors, + tol: DEFAULT_SINGULAR_TOL, + }; + ldlt.det() + }; + assert!((DET - 1.0).abs() <= 1e-12); + } + + #[test] + fn ldlt_det_const_eval_d3() { + const DET: f64 = { + // Diagonal D = [2.0, 3.0, 5.0] ⇒ det = 30.0. + let mut factors = Matrix::<3>::identity(); + factors.rows[0][0] = 2.0; + factors.rows[1][1] = 3.0; + factors.rows[2][2] = 5.0; + let ldlt = Ldlt::<3> { + factors, + tol: DEFAULT_SINGULAR_TOL, + }; + ldlt.det() + }; + assert!((DET - 30.0).abs() <= 1e-12); + } + + #[test] + fn ldlt_solve_vec_const_eval_d2() { + // Identity factors ⇒ solve_vec returns the RHS untouched. + const X: [f64; 2] = { + let ldlt = Ldlt::<2> { + factors: Matrix::<2>::identity(), + tol: DEFAULT_SINGULAR_TOL, + }; + let b = Vector::<2>::new([1.0, 2.0]); + match ldlt.solve_vec(b) { + Ok(v) => v.into_array(), + Err(_) => [0.0, 0.0], + } + }; + assert!((X[0] - 1.0).abs() <= 1e-12); + assert!((X[1] - 2.0).abs() <= 1e-12); + } } diff --git a/src/lib.rs b/src/lib.rs index ae0ea84..070cf76 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -119,11 +119,22 @@ pub enum LaError { pivot_col: usize, }, /// A non-finite value (NaN/∞) was encountered. + /// + /// The `(row, col)` coordinate follows a consistent convention across the crate: + /// + /// - `row: Some(r), col: c` — a *stored* matrix cell at `(r, c)` is non-finite. + /// Used by `Matrix::det`, `Lu::factor`, `Ldlt::factor`, and the `solve_vec` + /// paths when they detect a corrupt stored factor (only reachable via + /// direct struct construction; `factor` itself rejects such inputs). + /// - `row: None, col: c` — the non-finite value is either a *vector input* + /// entry at index `c`, or a *computed intermediate* at step `c` + /// (e.g. an accumulator that overflowed during forward/back substitution). NonFinite { - /// Row of the non-finite entry (for matrix inputs), or `None` when - /// the error originates from a vector input or a computed intermediate. + /// Row of the non-finite entry for a stored matrix cell, or `None` for + /// a vector-input entry or a computed intermediate. See the variant + /// docs for the full convention. row: Option, - /// Column index (for matrix inputs), vector index, or factorization + /// Column index (stored cell), vector index, or factorization/solve /// step where the non-finite value was detected. col: usize, }, diff --git a/src/lu.rs b/src/lu.rs index 62d1817..92bce28 100644 --- a/src/lu.rs +++ b/src/lu.rs @@ -123,39 +123,66 @@ impl Lu { /// # Errors /// Returns [`LaError::Singular`] if a diagonal entry of `U` satisfies `|u_ii| <= tol`, where /// `tol` is the tolerance that was used during factorization. - /// Returns [`LaError::NonFinite`] if NaN/∞ is detected. + /// + /// Returns [`LaError::NonFinite`] if NaN/∞ is detected. The `row`/`col` coordinates + /// follow the convention documented on [`LaError::NonFinite`]: + /// + /// - `row: Some(i), col: i` — the stored `U` diagonal at `(i, i)` is non-finite + /// (only reachable via direct `Lu` construction; [`Matrix::lu`](crate::Matrix::lu) + /// rejects such factorizations). + /// - `row: None, col: i` — a computed intermediate (forward/back-substitution + /// accumulator or the quotient `sum / diag`) overflowed to NaN/∞ at step `i`. #[inline] - pub fn solve_vec(&self, b: Vector) -> Result, LaError> { + pub const fn solve_vec(&self, b: Vector) -> Result, LaError> { let mut x = [0.0; D]; - for (i, x_i) in x.iter_mut().enumerate() { - *x_i = b.data[self.piv[i]]; + let mut i = 0; + while i < D { + x[i] = b.data[self.piv[i]]; + i += 1; } // Forward substitution for L (unit diagonal). - for i in 0..D { + let mut i = 0; + while i < D { let mut sum = x[i]; let row = self.factors.rows[i]; - for (j, x_j) in x.iter().enumerate().take(i) { - sum = (-row[j]).mul_add(*x_j, sum); + let mut j = 0; + while j < i { + sum = (-row[j]).mul_add(x[j], sum); + j += 1; } if !sum.is_finite() { cold_path(); return Err(LaError::NonFinite { row: None, col: i }); } x[i] = sum; + i += 1; } // Back substitution for U. - for ii in 0..D { + let mut ii = 0; + while ii < D { let i = D - 1 - ii; let mut sum = x[i]; let row = self.factors.rows[i]; - for (j, x_j) in x.iter().enumerate().skip(i + 1) { - sum = (-row[j]).mul_add(*x_j, sum); + let mut j = i + 1; + while j < D { + sum = (-row[j]).mul_add(x[j], sum); + j += 1; } let diag = row[i]; - if !diag.is_finite() || !sum.is_finite() { + // Distinguish a corrupt stored pivot (row: Some(i), col: i) from + // a computed intermediate overflow (row: None, col: i) so callers + // can diagnose the failure source without inspecting internals. + if !diag.is_finite() { + cold_path(); + return Err(LaError::NonFinite { + row: Some(i), + col: i, + }); + } + if !sum.is_finite() { cold_path(); return Err(LaError::NonFinite { row: None, col: i }); } @@ -164,12 +191,13 @@ impl Lu { return Err(LaError::Singular { pivot_col: i }); } - let q = sum / diag; - if !q.is_finite() { + let quotient = sum / diag; + if !quotient.is_finite() { cold_path(); return Err(LaError::NonFinite { row: None, col: i }); } - x[i] = q; + x[i] = quotient; + ii += 1; } Ok(Vector::new(x)) @@ -192,10 +220,12 @@ impl Lu { /// ``` #[inline] #[must_use] - pub fn det(&self) -> f64 { + pub const fn det(&self) -> f64 { let mut det = self.piv_sign; - for i in 0..D { + let mut i = 0; + while i < D { det *= self.factors.rows[i][i]; + i += 1; } det } @@ -505,15 +535,17 @@ mod tests { // ----------------------------------------------------------------------- // Defensive-path coverage for `solve_vec`. // - // `Lu::factor` guarantees that every stored U diagonal satisfies - // `|U[i,i]| > tol`. `solve_vec` still re-checks during back-substitution - // as a safety net (see the `diag.abs() <= self.tol` guard). That branch - // is unreachable through the public API, so the only way to exercise it - // is to construct `Lu` directly with a corrupt U. The tests below - // document and verify that the safety net returns `Singular`. + // `Lu::factor` guarantees that every stored U diagonal is finite and + // satisfies `|U[i,i]| > tol`. `solve_vec` still re-checks during + // back-substitution as a safety net (see the `!diag.is_finite()` and + // `diag.abs() <= self.tol` guards). Those branches are unreachable + // through the public API, so the only way to exercise them is to + // construct `Lu` directly with a corrupt U. The tests below document + // and verify that the safety nets return the documented error variants + // with coordinates that locate the offending stored cell. // ----------------------------------------------------------------------- - macro_rules! gen_solve_vec_defensive_singular_tests { + macro_rules! gen_solve_vec_defensive_tests { ($d:literal) => { paste! { /// `solve_vec` must surface `Singular` when a stored U @@ -539,12 +571,107 @@ mod tests { let err = lu.solve_vec(b).unwrap_err(); assert_eq!(err, LaError::Singular { pivot_col: $d - 1 }); } + + /// `solve_vec` must surface `NonFinite` with the corrupt + /// cell's coordinates when a stored U diagonal is NaN, + /// even though `factor` cannot produce such a + /// factorization. The error must pinpoint `(D-1, D-1)` + /// per the [`LaError::NonFinite`] convention. + #[test] + fn []() { + let mut factors = Matrix::<$d>::identity(); + factors.rows[$d - 1][$d - 1] = f64::NAN; + + let mut piv = [0usize; $d]; + for (i, p) in piv.iter_mut().enumerate() { + *p = i; + } + + let lu = Lu::<$d> { + factors, + piv, + piv_sign: 1.0, + tol: DEFAULT_PIVOT_TOL, + }; + let b = Vector::<$d>::new([1.0; $d]); + let err = lu.solve_vec(b).unwrap_err(); + assert_eq!( + err, + LaError::NonFinite { + row: Some($d - 1), + col: $d - 1, + } + ); + } } }; } - gen_solve_vec_defensive_singular_tests!(2); - gen_solve_vec_defensive_singular_tests!(3); - gen_solve_vec_defensive_singular_tests!(4); - gen_solve_vec_defensive_singular_tests!(5); + gen_solve_vec_defensive_tests!(2); + gen_solve_vec_defensive_tests!(3); + gen_solve_vec_defensive_tests!(4); + gen_solve_vec_defensive_tests!(5); + + // ----------------------------------------------------------------------- + // Const-evaluability tests. + // + // These prove that `Lu::det` and `Lu::solve_vec` are truly `const fn` by + // forcing the compiler to evaluate them inside a `const` initializer. + // `Lu::factor` is not (yet) `const fn` because it relies on `<[T]>::swap`, + // which is not const-stable; we therefore construct `Lu` directly. + // ----------------------------------------------------------------------- + + #[test] + fn lu_det_const_eval_d2() { + const DET: f64 = { + // Triangular factors with diag [2.0, 3.0] and no row swaps. + let mut factors = Matrix::<2>::identity(); + factors.rows[0][0] = 2.0; + factors.rows[1][1] = 3.0; + let lu = Lu::<2> { + factors, + piv: [0, 1], + piv_sign: 1.0, + tol: DEFAULT_PIVOT_TOL, + }; + lu.det() + }; + assert!((DET - 6.0).abs() <= 1e-12); + } + + #[test] + fn lu_det_const_eval_d3_row_swap() { + const DET: f64 = { + // Identity factors but `piv_sign = -1.0` encoding a single row swap; + // the determinant magnitude is 1 but the sign flips. + let lu = Lu::<3> { + factors: Matrix::<3>::identity(), + piv: [1, 0, 2], + piv_sign: -1.0, + tol: DEFAULT_PIVOT_TOL, + }; + lu.det() + }; + assert!((DET - (-1.0)).abs() <= 1e-12); + } + + #[test] + fn lu_solve_vec_const_eval_d2() { + // Identity LU ⇒ solve_vec returns the permuted RHS untouched. + const X: [f64; 2] = { + let lu = Lu::<2> { + factors: Matrix::<2>::identity(), + piv: [0, 1], + piv_sign: 1.0, + tol: DEFAULT_PIVOT_TOL, + }; + let b = Vector::<2>::new([1.0, 2.0]); + match lu.solve_vec(b) { + Ok(v) => v.into_array(), + Err(_) => [0.0, 0.0], + } + }; + assert!((X[0] - 1.0).abs() <= 1e-12); + assert!((X[1] - 2.0).abs() <= 1e-12); + } } diff --git a/src/matrix.rs b/src/matrix.rs index ee54844..f0af519 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -5,6 +5,7 @@ use core::hint::cold_path; use crate::LaError; use crate::ldlt::Ldlt; use crate::lu::Lu; +use crate::{ERR_COEFF_2, ERR_COEFF_3, ERR_COEFF_4}; /// Fixed-size square matrix `D×D`, stored inline. #[must_use] @@ -132,11 +133,21 @@ impl Matrix { /// ``` #[inline] #[must_use] - pub fn inf_norm(&self) -> f64 { + pub const fn inf_norm(&self) -> f64 { let mut max_row_sum: f64 = 0.0; - for row in &self.rows { - let row_sum: f64 = row.iter().map(|&x| x.abs()).sum(); + let mut r = 0; + while r < D { + // Iterator chains like `row.iter().map(|x| x.abs()).sum()` are + // not yet const-stable, so accumulate the absolute row sum with + // a manual `while` loop. + let row = &self.rows[r]; + let mut row_sum: f64 = 0.0; + let mut c = 0; + while c < D { + row_sum += row[c].abs(); + c += 1; + } // Propagate NaN explicitly: `f64::max` drops NaN (IEEE 754 `maxNum`) // and `f64::maximum` (IEEE 754-2019 `maximum`) is still unstable, // so we short-circuit on NaN instead. @@ -147,6 +158,7 @@ impl Matrix { if row_sum > max_row_sum { max_row_sum = row_sum; } + r += 1; } max_row_sum @@ -510,13 +522,13 @@ impl Matrix { /// ``` #[must_use] #[inline] - pub fn det_errbound(&self) -> Option { + pub const fn det_errbound(&self) -> Option { match D { 0 | 1 => Some(0.0), // No arithmetic — result is exact. 2 => { let r = &self.rows; let permanent = (r[0][0] * r[1][1]).abs() + (r[0][1] * r[1][0]).abs(); - Some(crate::ERR_COEFF_2 * permanent) + Some(ERR_COEFF_2 * permanent) } 3 => { let r = &self.rows; @@ -526,7 +538,7 @@ impl Matrix { let permanent = r[0][2] .abs() .mul_add(pm02, r[0][1].abs().mul_add(pm01, r[0][0].abs() * pm00)); - Some(crate::ERR_COEFF_3 * permanent) + Some(ERR_COEFF_3 * permanent) } 4 => { let r = &self.rows; @@ -556,7 +568,7 @@ impl Matrix { .abs() .mul_add(pc2, r[0][1].abs().mul_add(pc1, r[0][0].abs() * pc0)), ); - Some(crate::ERR_COEFF_4 * permanent) + Some(ERR_COEFF_4 * permanent) } _ => None, } @@ -861,7 +873,7 @@ mod tests { } #[test] - fn det_direct_is_const_evaluable_d2() { + fn det_direct_const_eval_d2() { // Const evaluation proves the function is truly const fn. const DET: Option = { let m = Matrix::<2>::from_rows([[1.0, 0.0], [0.0, 1.0]]); @@ -871,7 +883,7 @@ mod tests { } #[test] - fn det_direct_is_const_evaluable_d3() { + fn det_direct_const_eval_d3() { const DET: Option = { let m = Matrix::<3>::from_rows([[2.0, 0.0, 0.0], [0.0, 3.0, 0.0], [0.0, 0.0, 5.0]]); m.det_direct() @@ -896,6 +908,46 @@ mod tests { assert_eq!(Matrix::<5>::identity().det_errbound(), None); } + #[test] + fn det_errbound_const_eval_d2() { + // Const evaluation proves the function is truly const fn. + const BOUND: Option = { + let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); + m.det_errbound() + }; + assert!(BOUND.is_some()); + assert!(BOUND.unwrap() > 0.0); + } + + #[test] + fn det_errbound_const_eval_d3() { + const BOUND: Option = { + let m = Matrix::<3>::from_rows([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]); + m.det_errbound() + }; + assert!(BOUND.is_some()); + assert!(BOUND.unwrap() > 0.0); + } + + #[test] + fn inf_norm_const_eval_d2() { + // Maximum absolute row sum: row 0 = 3, row 1 = 7 ⇒ 7. + const NORM: f64 = { + let m = Matrix::<2>::from_rows([[1.0, -2.0], [3.0, 4.0]]); + m.inf_norm() + }; + assert!((NORM - 7.0).abs() <= 1e-12); + } + + #[test] + fn inf_norm_const_eval_d3() { + const NORM: f64 = { + let m = Matrix::<3>::identity(); + m.inf_norm() + }; + assert!((NORM - 1.0).abs() <= 1e-12); + } + // === inf_norm NaN / Inf propagation (regression tests for #85) === macro_rules! gen_inf_norm_nonfinite_tests { From 0b98d3f6dbdd74699c318c4744a2b2f9a1b78481 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Mon, 20 Apr 2026 15:12:00 -0700 Subject: [PATCH 2/3] Changed: finalize documentation, benchmarks, and error handling Synchronize the CHANGELOG, README, and REFERENCES for the upcoming release. This includes documenting the new absolute error bound coefficients for determinants, adding a dedicated section for AI agent governance, and refining internal error helpers to use a more consistent coordinate convention for non-finite values. Refs: #86 --- CHANGELOG.md | 35 +++++++++++++ README.md | 53 +++++++++++++------ REFERENCES.md | 72 +++++++++++++++----------- src/exact.rs | 7 +-- src/ldlt.rs | 26 +++------- src/lib.rs | 139 ++++++++++++++++++++++++++++++++++++++++++++------ src/lu.rs | 31 +++-------- src/matrix.rs | 17 +++--- typos.toml | 3 ++ 9 files changed, 266 insertions(+), 117 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 19994aa..b4f5dd9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,39 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] + +### Added + +- Regression tests for solver and determinant overflow handling [`f763b11`](https://github.com/acgetchell/la-stack/commit/f763b119bcc57276b83f370b0bf7abce654c7eb8) +- Defensive-path test coverage for LU and LDLT solve_vec [`87d426f`](https://github.com/acgetchell/la-stack/commit/87d426fca1627445b804fd26b62fc7d9d4f0ae48) +- Const-ify Lu/Ldlt det + solve_vec and Matrix inf_norm + det_errbound [`81ecb35`](https://github.com/acgetchell/la-stack/commit/81ecb35bdaf159f1f44d1eb24274ecf82c6567d5) + +### Changed + +- Report infinite vs finite off-diagonal pairs as asymmetric [`1805779`](https://github.com/acgetchell/la-stack/commit/1805779dbca49183fbfa95c68ec00984966aa551) + +### Documentation + +- Strengthen LDLT symmetry precondition and add is_symmetric API [`1693307`](https://github.com/acgetchell/la-stack/commit/1693307c58e30a804e9b79abc783af99be8dc947) + +### Fixed + +- Propagate NaN in Matrix::inf_norm [#85](https://github.com/acgetchell/la-stack/pull/85) [`16ffa45`](https://github.com/acgetchell/la-stack/commit/16ffa45ade11b21a179cad8fcecc51d802086a1d) + +### Maintenance + +- Bump taiki-e/install-action from 2.73.0 to 2.75.7 [`a1d1c1e`](https://github.com/acgetchell/la-stack/commit/a1d1c1edba7d5bf6928de4e8cab86ba853685430) +- Bump actions/download-artifact from 4.3.0 to 8.0.1 [`8772ea2`](https://github.com/acgetchell/la-stack/commit/8772ea2f18cf1b94a21e269b4bc0c8e0a3b650eb) +- Bump rand from 0.9.2 to 0.9.4 [`b91670a`](https://github.com/acgetchell/la-stack/commit/b91670a46f82306dde5433ca7c406a8757e3fc59) +- Bump actions/upload-artifact from 7.0.0 to 7.0.1 [`7579274`](https://github.com/acgetchell/la-stack/commit/7579274d56ab8c46eae6906c07a3c0d0fd41c84f) +- Bump actions/github-script from 8.0.0 to 9.0.0 [`a522abd`](https://github.com/acgetchell/la-stack/commit/a522abdabf211d47f0f52cbfc2a349fa72ce87e2) +- Bump actions/cache from 5.0.4 to 5.0.5 [`b12d95c`](https://github.com/acgetchell/la-stack/commit/b12d95c8ce242ea6014c6a6aaa7845ff501d2e4d) +- Bump actions-rust-lang/setup-rust-toolchain [`c728774`](https://github.com/acgetchell/la-stack/commit/c72877460efba3d598f64ff253d2df8f3695acb1) +- Bump astral-sh/setup-uv from 8.0.0 to 8.1.0 [`7a18733`](https://github.com/acgetchell/la-stack/commit/7a18733ba0f39102678bbe30772fcde10f20e9d8) +- Bump taiki-e/install-action from 2.75.7 to 2.75.18 [`433cfc1`](https://github.com/acgetchell/la-stack/commit/433cfc1b01c9128e93c82cb553aa63d4091bace3) +- Bump MSRV to Rust 1.95 and adopt new stable features [`0ab3c33`](https://github.com/acgetchell/la-stack/commit/0ab3c336074f2b866256fbe5db8a8ec5306d580a) + ## [0.4.0] - 2026-04-11 ### Added @@ -35,6 +68,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Bump taiki-e/install-action from 2.70.1 to 2.73.0 [`7de720d`](https://github.com/acgetchell/la-stack/commit/7de720dfb8328d01843cfabd15d23086ee98832b) - Bump astral-sh/setup-uv from 7.6.0 to 8.0.0 [`af40753`](https://github.com/acgetchell/la-stack/commit/af40753130fac56d48da7fce2f18b11dc391ebe6) - Add performance regression detection for exact-arithmetic benchmarks [`44bce99`](https://github.com/acgetchell/la-stack/commit/44bce99bcae8f852dbf7500e71fa182730e08bca) +- Release v0.4.0 [`843bb3c`](https://github.com/acgetchell/la-stack/commit/843bb3cfcc88e114bdd5e00d242b605af9b2b434) ### Performance @@ -228,6 +262,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Add tarpaulin coverage upload [`7486dfd`](https://github.com/acgetchell/la-stack/commit/7486dfd54e16a6dbde41575c3f35a1acb65f57d2) +[unreleased]: https://github.com/acgetchell/la-stack/compare/v0.4.0..HEAD [0.4.0]: https://github.com/acgetchell/la-stack/compare/v0.3.0..v0.4.0 [0.3.0]: https://github.com/acgetchell/la-stack/compare/v0.2.2..v0.3.0 [0.2.2]: https://github.com/acgetchell/la-stack/compare/v0.2.1..v0.2.2 diff --git a/README.md b/README.md index d7be54f..b4a1507 100644 --- a/README.md +++ b/README.md @@ -242,6 +242,27 @@ Storage shown above reflects the `f64` implementation. ¹ Requires `features = ["exact"]`. +## 📊 Benchmarks (vs nalgebra/faer) + +![LU solve (factor + solve): median time vs dimension](docs/assets/bench/vs_linalg_lu_solve_median.svg) + +Raw data: [docs/assets/bench/vs_linalg_lu_solve_median.csv](docs/assets/bench/vs_linalg_lu_solve_median.csv) + +Summary (median time; lower is better). The “la-stack vs nalgebra/faer” columns show the % time reduction relative to each baseline (positive = la-stack faster): + + +| D | la-stack median (ns) | nalgebra median (ns) | faer median (ns) | la-stack vs nalgebra | la-stack vs faer | +|---:|--------------------:|--------------------:|----------------:|---------------------:|----------------:| +| 2 | 2.309 | 4.365 | 140.156 | +47.1% | +98.4% | +| 3 | 18.331 | 22.706 | 181.074 | +19.3% | +89.9% | +| 4 | 27.430 | 51.372 | 210.451 | +46.6% | +87.0% | +| 5 | 53.819 | 70.722 | 276.064 | +23.9% | +80.5% | +| 8 | 143.611 | 160.309 | 356.960 | +10.4% | +59.8% | +| 16 | 611.393 | 580.793 | 871.704 | -5.3% | +29.9% | +| 32 | 2,631.241 | 2,733.946 | 2,832.816 | +3.8% | +7.1% | +| 64 | 17,233.345 | 14,112.678 | 12,164.571 | -22.1% | -41.7% | + + ## 📋 Examples The `examples/` directory contains small, runnable programs: @@ -289,26 +310,26 @@ If you use this library in academic work, please cite it using [CITATION.cff](CI For canonical references to the algorithms used by this crate, see [REFERENCES.md](REFERENCES.md). -## 📊 Benchmarks (vs nalgebra/faer) +## 🤖 AI Agents -![LU solve (factor + solve): median time vs dimension](docs/assets/bench/vs_linalg_lu_solve_median.svg) +This repository contains an `AGENTS.md` file, which defines the canonical rules and invariants +for all AI coding assistants and autonomous agents working on this codebase. -Raw data: [docs/assets/bench/vs_linalg_lu_solve_median.csv](docs/assets/bench/vs_linalg_lu_solve_median.csv) +AI tools (including `ChatGPT`, `Claude`, `CodeRabbit`, `KiloCode`, and `WARP`) are expected to read +and follow `AGENTS.md` when proposing or applying changes. -Summary (median time; lower is better). The “la-stack vs nalgebra/faer” columns show the % time reduction relative to each baseline (positive = la-stack faster): +Portions of this library were developed with the assistance of these AI tools: - -| D | la-stack median (ns) | nalgebra median (ns) | faer median (ns) | la-stack vs nalgebra | la-stack vs faer | -|---:|--------------------:|--------------------:|----------------:|---------------------:|----------------:| -| 2 | 2.309 | 4.365 | 140.156 | +47.1% | +98.4% | -| 3 | 18.331 | 22.706 | 181.074 | +19.3% | +89.9% | -| 4 | 27.430 | 51.372 | 210.451 | +46.6% | +87.0% | -| 5 | 53.819 | 70.722 | 276.064 | +23.9% | +80.5% | -| 8 | 143.611 | 160.309 | 356.960 | +10.4% | +59.8% | -| 16 | 611.393 | 580.793 | 871.704 | -5.3% | +29.9% | -| 32 | 2,631.241 | 2,733.946 | 2,832.816 | +3.8% | +7.1% | -| 64 | 17,233.345 | 14,112.678 | 12,164.571 | -22.1% | -41.7% | - +- [ChatGPT](https://openai.com/chatgpt) +- [Claude](https://www.anthropic.com/claude) +- [CodeRabbit](https://coderabbit.ai/) +- [KiloCode](https://kilocode.ai/) +- [WARP](https://www.warp.dev) + +> All code was written and/or reviewed and validated by the author. + +For full tool citation metadata, see the [AI-Assisted Development Tools](REFERENCES.md#ai-assisted-development-tools) +section of `REFERENCES.md`. ## 📄 License diff --git a/REFERENCES.md b/REFERENCES.md index 3d394ba..bd850c7 100644 --- a/REFERENCES.md +++ b/REFERENCES.md @@ -8,14 +8,46 @@ processed by GitHub and other platforms. A Zenodo DOI will be added for tagged releases. +## AI-Assisted Development Tools + +- Anthropic. "Claude." . +- CodeRabbit AI, Inc. "CodeRabbit." . +- KiloCode. "KiloCode AI Engineering Assistant." . +- OpenAI. "ChatGPT." . +- Warp Dev, Inc. "WARP." . + +All AI-generated output was reviewed and/or edited by the maintainer. +No generated content was used without human oversight. + ## Linear algebra algorithms -### LU decomposition (Gaussian elimination with partial pivoting) +### Absolute error bound for closed-form determinants -The LU implementation in `la-stack` follows the standard Gaussian elimination / LU factorization -approach with partial pivoting for numerical stability. +`Matrix::det_errbound()` returns a conservative Shewchuk-style absolute error bound [8] +for `Matrix::det_direct()` in dimensions 2–4. The bound has the form +`ERR_COEFF_D · p(|A|)`, where `p(|A|)` is the absolute Leibniz sum (the cofactor-expansion +tree with `|·|` at every leaf) and `ERR_COEFF_D ∈ {ERR_COEFF_2, ERR_COEFF_3, ERR_COEFF_4}` +is a dimension-specific constant derived from the rounding-event count of `det_direct`. +The same bound is used internally by `det_sign_exact()`'s fast filter, but +`det_errbound()` itself is available without the `exact` feature, so downstream crates +can build custom adaptive-precision logic with pure f64 arithmetic. -See references [1-3] below. +### Exact determinant sign (adaptive-precision Bareiss) + +`det_sign_exact()` uses a Shewchuk-style f64 error-bound filter [8] (the same bound exposed +by `det_errbound()` above) backed by integer-only Bareiss elimination [7] in `BigInt`. Each +f64 entry is decomposed into `mantissa × 2^exponent`, scaled to a common integer base, and +eliminated without any `BigRational` or GCD overhead. +See `src/exact.rs` for the full architecture description. + +### f64 → BigRational conversion (`f64_to_bigrational`) + +`f64_to_bigrational` converts an f64 to an exact `BigRational` by decomposing the IEEE 754 +binary64 bit representation [9] into its sign, exponent, and significand fields. Because +every finite f64 is exactly `±m × 2^e` (where `m` is an integer), the rational can be +constructed directly via `BigRational::new_raw` without GCD normalization — trailing zeros +in the significand are stripped first so the fraction is already in lowest terms. See +Goldberg [10] for background on IEEE 754 representation and exact rational reconstruction. ### LDL^T factorization (symmetric SPD/PSD) @@ -26,16 +58,14 @@ pivoting. For background on the SPD/PSD setting, see [4-5]. For pivoted variants used for symmetric *indefinite* matrices, see [6]. -### Exact determinant sign (adaptive-precision Bareiss) +### LU decomposition (Gaussian elimination with partial pivoting) -`det_sign_exact()` uses a Shewchuk-style f64 error-bound filter [8] backed by integer-only -Bareiss elimination [7] in `BigInt`. Each f64 entry is decomposed into `mantissa × 2^exponent`, -scaled to a common integer base, and eliminated without any `BigRational` or GCD overhead. -See `src/exact.rs` for the full architecture description. +The LU implementation in `la-stack` follows the standard Gaussian elimination / LU factorization +approach with partial pivoting for numerical stability. -## References +See references [1-3] below. -### LU / Gaussian elimination +## References 1. Trefethen, Lloyd N., and Robert S. Schreiber. "Average-case stability of Gaussian elimination." *SIAM Journal on Matrix Analysis and Applications* 11.3 (1990): 335–360. @@ -46,23 +76,14 @@ See `src/exact.rs` for the full architecture description. 3. Huang, Han, and K. Tikhomirov. "Average-case analysis of the Gaussian elimination with partial pivoting." *Probability Theory and Related Fields* 189 (2024): 501–567. [Open-access PDF](https://link.springer.com/article/10.1007/s00440-024-01276-2) (also: [arXiv:2206.01726](https://arxiv.org/abs/2206.01726)) - -### LDL^T / Cholesky (symmetric SPD/PSD) - 4. Cholesky, Andre-Louis. "On the numerical solution of systems of linear equations" (manuscript dated 2 Dec 1910; published 2005). Scan + English analysis: [BibNum](https://www.bibnum.education.fr/mathematiques/algebre/sur-la-resolution-numerique-des-systemes-d-equations-lineaires) 5. Brezinski, Claude. "La methode de Cholesky." (2005). [PDF](https://eudml.org/doc/252115) - -### Pivoted LDL^T (symmetric indefinite) - 6. Bunch, J. R., L. Kaufman, and B. N. Parlett. "Decomposition of a Symmetric Matrix." *Numerische Mathematik* 27 (1976/1977): 95–110. [Full text](https://eudml.org/doc/132435) - -### Exact determinant (`det_exact`, `det_exact_f64`, `det_sign_exact`) - 7. Bareiss, Erwin H. "Sylvester's Identity and Multistep Integer-Preserving Gaussian Elimination." *Mathematics of Computation* 22.103 (1968): 565–578. [DOI](https://doi.org/10.1090/S0025-5718-1968-0226829-0) · @@ -72,17 +93,6 @@ See `src/exact.rs` for the full architecture description. [DOI](https://doi.org/10.1007/PL00009321) · [PDF](https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf) Also: Technical Report CMU-CS-96-140, Carnegie Mellon University, May 1996. - -### f64 → BigRational conversion (`f64_to_bigrational`) - -`f64_to_bigrational` converts an f64 to an exact `BigRational` by decomposing the IEEE 754 -binary64 bit representation into its sign, exponent, and significand fields. Because every -finite f64 is exactly `±m × 2^e` (where `m` is an integer), the rational can be constructed -directly via `BigRational::new_raw` without GCD normalization — trailing zeros in the -significand are stripped first so the fraction is already in lowest terms. - -See references [9-10] below. - 9. IEEE Computer Society. "IEEE Standard for Floating-Point Arithmetic." *IEEE Std 754-2019* (Revision of IEEE 754-2008), 2019. [DOI](https://doi.org/10.1109/IEEESTD.2019.8766229) diff --git a/src/exact.rs b/src/exact.rs index 6bd81bf..83d9562 100644 --- a/src/exact.rs +++ b/src/exact.rs @@ -63,10 +63,7 @@ fn validate_finite(m: &Matrix) -> Result<(), LaError> { for c in 0..D { if !m.rows[r][c].is_finite() { cold_path(); - return Err(LaError::NonFinite { - row: Some(r), - col: c, - }); + return Err(LaError::non_finite_cell(r, c)); } } } @@ -81,7 +78,7 @@ fn validate_finite_vec(v: &Vector) -> Result<(), LaError> { for (i, &x) in v.data.iter().enumerate() { if !x.is_finite() { cold_path(); - return Err(LaError::NonFinite { row: None, col: i }); + return Err(LaError::non_finite_at(i)); } } Ok(()) diff --git a/src/ldlt.rs b/src/ldlt.rs index d2bbf84..a2542e2 100644 --- a/src/ldlt.rs +++ b/src/ldlt.rs @@ -62,10 +62,7 @@ impl Ldlt { let d = f.rows[j][j]; if !d.is_finite() { cold_path(); - return Err(LaError::NonFinite { - row: Some(j), - col: j, - }); + return Err(LaError::non_finite_cell(j, j)); } if d <= tol { cold_path(); @@ -77,10 +74,7 @@ impl Ldlt { let l = f.rows[i][j] / d; if !l.is_finite() { cold_path(); - return Err(LaError::NonFinite { - row: Some(i), - col: j, - }); + return Err(LaError::non_finite_cell(i, j)); } f.rows[i][j] = l; } @@ -95,10 +89,7 @@ impl Ldlt { let new_val = (-l_i_d).mul_add(l_k, f.rows[i][k]); if !new_val.is_finite() { cold_path(); - return Err(LaError::NonFinite { - row: Some(i), - col: k, - }); + return Err(LaError::non_finite_cell(i, k)); } f.rows[i][k] = new_val; } @@ -181,7 +172,7 @@ impl Ldlt { } if !sum.is_finite() { cold_path(); - return Err(LaError::NonFinite { row: None, col: i }); + return Err(LaError::non_finite_at(i)); } x[i] = sum; i += 1; @@ -197,10 +188,7 @@ impl Ldlt { // `Matrix::det`, `Lu::factor`, and `Ldlt::factor`. if !diag.is_finite() { cold_path(); - return Err(LaError::NonFinite { - row: Some(i), - col: i, - }); + return Err(LaError::non_finite_cell(i, i)); } if diag <= self.tol { cold_path(); @@ -210,7 +198,7 @@ impl Ldlt { let quotient = x[i] / diag; if !quotient.is_finite() { cold_path(); - return Err(LaError::NonFinite { row: None, col: i }); + return Err(LaError::non_finite_at(i)); } x[i] = quotient; i += 1; @@ -228,7 +216,7 @@ impl Ldlt { } if !sum.is_finite() { cold_path(); - return Err(LaError::NonFinite { row: None, col: i }); + return Err(LaError::non_finite_at(i)); } x[i] = sum; ii += 1; diff --git a/src/lib.rs b/src/lib.rs index 070cf76..8b9eb38 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -61,34 +61,111 @@ mod vector; use core::fmt; // --------------------------------------------------------------------------- -// Error-bound constants for determinant error analysis. +// Error-bound constants for `Matrix::det_errbound()`. // -// These constants bound the absolute error of `det_direct()` relative to the -// *permanent* (sum of absolute products in the Leibniz expansion). The -// constants are conservative over-estimates following Shewchuk's methodology. +// For `D ∈ {2, 3, 4}`, `Matrix::det_direct()` evaluates the Leibniz expansion +// of the determinant as a tree of f64 multiplies and fused multiply-adds +// (FMAs). Following Shewchuk's error-analysis methodology (REFERENCES.md +// [8]), the absolute error of that computation is bounded by // -// These are NOT feature-gated because they use pure f64 arithmetic and are -// useful for adaptive-precision logic even without the `exact` feature. +// |det_direct(A) - det_exact(A)| ≤ ERR_COEFF_D · p(|A|) +// +// where `p(|A|)` is the **absolute Leibniz sum** +// +// p(|A|) = Σ_σ ∏ᵢ |A[i, σ(i)]|, +// +// i.e. the same cofactor-expansion tree as `det_direct` but with each +// entry replaced by its magnitude. Note that `p(|A|)` is *not* the +// combinatorial matrix permanent — the name "permanent" appears in the +// source for brevity and to match the cited literature. +// +// Each constant has the shape `a · EPS + b · EPS²`: the linear term bounds +// the first-order rounding and the quadratic term absorbs the interaction +// of errors in nested FMAs. The coefficients `a` and `b` are conservative +// over-estimates derived from the longest dependency chain of `det_direct` +// at that dimension. +// +// These constants are NOT feature-gated — they rely only on f64 arithmetic +// and are useful for adaptive-precision logic even without the `exact` +// feature. Most callers should prefer `Matrix::det_errbound()`, which +// applies these constants to the actual matrix; the raw constants are +// exposed for advanced use cases (composing the bound with a pre-reduced +// permanent, rolling a custom adaptive filter, etc.). See +// `Matrix::det_sign_exact()` (behind the `exact` feature) for the +// reference adaptive-filter that consumes these internally. // --------------------------------------------------------------------------- const EPS: f64 = f64::EPSILON; // 2^-52 -/// Error coefficient for D=2 determinant error bound. +/// Absolute error coefficient for [`Matrix::<2>::det_direct`](crate::Matrix::det_direct). +/// +/// For any 2×2 matrix `A = [[a, b], [c, d]]` with finite f64 entries, /// -/// Accounts for one f64 multiply + one FMA → 2 rounding events. -/// Used in computing the absolute error bound for 2×2 determinants. +/// ```text +/// |A.det_direct() - det_exact(A)| ≤ ERR_COEFF_2 · (|a·d| + |b·c|) +/// ``` +/// +/// `det_direct` evaluates `a·d - b·c` as one multiply followed by one FMA +/// (2 rounding events); the linear `3·EPS` term bounds those roundings +/// and the quadratic `16·EPS²` term is a conservative cushion for their +/// interaction. Derivation follows Shewchuk's framework; see +/// `REFERENCES.md` \[8\]. +/// +/// Prefer [`Matrix::det_errbound`](crate::Matrix::det_errbound) unless +/// you already have the absolute-Leibniz sum available; see +/// [`Matrix::det_sign_exact`](crate::Matrix::det_sign_exact) (under the +/// `exact` feature) for the reference adaptive-precision filter. +/// +/// # Example +/// ``` +/// use la_stack::prelude::*; +/// +/// let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); +/// let det = m.det_direct().unwrap(); +/// // Compute the bound from the raw constant for illustration; most +/// // callers would just use `m.det_errbound().unwrap()` instead. +/// let p = (1.0_f64 * 4.0).abs() + (2.0_f64 * 3.0).abs(); +/// let bound = ERR_COEFF_2 * p; +/// if det.abs() > bound { +/// // The f64 sign is provably correct without exact arithmetic. +/// } +/// ``` pub const ERR_COEFF_2: f64 = 3.0 * EPS + 16.0 * EPS * EPS; -/// Error coefficient for D=3 determinant error bound. +/// Absolute error coefficient for [`Matrix::<3>::det_direct`](crate::Matrix::det_direct). +/// +/// For any 3×3 matrix `A` with finite f64 entries, +/// +/// ```text +/// |A.det_direct() - det_exact(A)| ≤ ERR_COEFF_3 · p(|A|) +/// ``` /// -/// Accounts for three 2×2 FMA minors + nested FMA combination. -/// Used in computing the absolute error bound for 3×3 determinants. +/// where `p(|A|)` is the absolute Leibniz sum (the same cofactor +/// expansion as `det_direct` but with `|·|` at every leaf). +/// `det_direct` for D=3 uses three 2×2 FMA minors combined by a nested +/// FMA, yielding the `8·EPS + 64·EPS²` bound. See `REFERENCES.md` +/// \[8\] for the Shewchuk framework these bounds follow. +/// +/// Prefer [`Matrix::det_errbound`](crate::Matrix::det_errbound) over this +/// constant for typical use; see [`ERR_COEFF_2`] for a worked example. pub const ERR_COEFF_3: f64 = 8.0 * EPS + 64.0 * EPS * EPS; -/// Error coefficient for D=4 determinant error bound. +/// Absolute error coefficient for [`Matrix::<4>::det_direct`](crate::Matrix::det_direct). +/// +/// For any 4×4 matrix `A` with finite f64 entries, +/// +/// ```text +/// |A.det_direct() - det_exact(A)| ≤ ERR_COEFF_4 · p(|A|) +/// ``` /// -/// Accounts for six hoisted 2×2 minors → four 3×3 cofactors → FMA row combination. -/// Used in computing the absolute error bound for 4×4 determinants. +/// where `p(|A|)` is the absolute Leibniz sum. `det_direct` for D=4 +/// hoists six 2×2 minors, combines them into four 3×3 cofactors, then +/// reduces those with an FMA row combination, yielding the +/// `12·EPS + 128·EPS²` bound. See `REFERENCES.md` \[8\] for the +/// Shewchuk framework these bounds follow. +/// +/// Prefer [`Matrix::det_errbound`](crate::Matrix::det_errbound) over this +/// constant for typical use; see [`ERR_COEFF_2`] for a worked example. pub const ERR_COEFF_4: f64 = 12.0 * EPS + 128.0 * EPS * EPS; /// Default absolute threshold used for singularity/degeneracy detection. @@ -150,6 +227,38 @@ pub enum LaError { }, } +impl LaError { + /// Construct a [`LaError::NonFinite`] pinpointing a stored matrix cell at `(row, col)`. + /// + /// Use this for non-finite values read from a stored `Matrix` entry or + /// factorization cell. The resulting error has `row: Some(row), col`, + /// matching the stored-cell convention documented on + /// [`NonFinite`](Self::NonFinite). For vector-input entries or computed + /// intermediates, use [`non_finite_at`](Self::non_finite_at). + #[inline] + #[must_use] + pub const fn non_finite_cell(row: usize, col: usize) -> Self { + Self::NonFinite { + row: Some(row), + col, + } + } + + /// Construct a [`LaError::NonFinite`] pinpointing a vector-input entry or + /// computed-intermediate step at index `col`. + /// + /// Use this for non-finite values in a `Vector` input or an accumulator + /// that overflowed during forward/back substitution. The resulting error + /// has `row: None, col`, matching the vector/intermediate convention + /// documented on [`NonFinite`](Self::NonFinite). For stored matrix cells, + /// use [`non_finite_cell`](Self::non_finite_cell). + #[inline] + #[must_use] + pub const fn non_finite_at(col: usize) -> Self { + Self::NonFinite { row: None, col } + } +} + impl fmt::Display for LaError { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { match *self { diff --git a/src/lu.rs b/src/lu.rs index 92bce28..6fbf23f 100644 --- a/src/lu.rs +++ b/src/lu.rs @@ -34,20 +34,14 @@ impl Lu { let mut pivot_abs = lu.rows[k][k].abs(); if !pivot_abs.is_finite() { cold_path(); - return Err(LaError::NonFinite { - row: Some(k), - col: k, - }); + return Err(LaError::non_finite_cell(k, k)); } for r in (k + 1)..D { let v = lu.rows[r][k].abs(); if !v.is_finite() { cold_path(); - return Err(LaError::NonFinite { - row: Some(r), - col: k, - }); + return Err(LaError::non_finite_cell(r, k)); } if v > pivot_abs { pivot_abs = v; @@ -69,10 +63,7 @@ impl Lu { let pivot = lu.rows[k][k]; if !pivot.is_finite() { cold_path(); - return Err(LaError::NonFinite { - row: Some(k), - col: k, - }); + return Err(LaError::non_finite_cell(k, k)); } // Eliminate below pivot. @@ -80,10 +71,7 @@ impl Lu { let mult = lu.rows[r][k] / pivot; if !mult.is_finite() { cold_path(); - return Err(LaError::NonFinite { - row: Some(r), - col: k, - }); + return Err(LaError::non_finite_cell(r, k)); } lu.rows[r][k] = mult; @@ -153,7 +141,7 @@ impl Lu { } if !sum.is_finite() { cold_path(); - return Err(LaError::NonFinite { row: None, col: i }); + return Err(LaError::non_finite_at(i)); } x[i] = sum; i += 1; @@ -177,14 +165,11 @@ impl Lu { // can diagnose the failure source without inspecting internals. if !diag.is_finite() { cold_path(); - return Err(LaError::NonFinite { - row: Some(i), - col: i, - }); + return Err(LaError::non_finite_cell(i, i)); } if !sum.is_finite() { cold_path(); - return Err(LaError::NonFinite { row: None, col: i }); + return Err(LaError::non_finite_at(i)); } if diag.abs() <= self.tol { cold_path(); @@ -194,7 +179,7 @@ impl Lu { let quotient = sum / diag; if !quotient.is_finite() { cold_path(); - return Err(LaError::NonFinite { row: None, col: i }); + return Err(LaError::non_finite_at(i)); } x[i] = quotient; ii += 1; diff --git a/src/matrix.rs b/src/matrix.rs index f0af519..a28b977 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -456,15 +456,12 @@ impl Matrix { for r in 0..D { for c in 0..D { if !self.rows[r][c].is_finite() { - return Err(LaError::NonFinite { - row: Some(r), - col: c, - }); + return Err(LaError::non_finite_cell(r, c)); } } } // All entries are finite but the determinant overflowed. - Err(LaError::NonFinite { row: None, col: 0 }) + Err(LaError::non_finite_at(0)) }; } self.lu(tol).map(|lu| lu.det()) @@ -475,9 +472,13 @@ impl Matrix { /// Returns `Some(bound)` such that `|det_direct() - det_exact| ≤ bound`, /// or `None` for D ≥ 5 where no fast bound is available. /// - /// For D ≤ 4, the bound is derived from the matrix permanent using - /// Shewchuk-style error analysis. For D = 0 or 1, returns `Some(0.0)` - /// since the determinant computation is exact (no arithmetic). + /// For D ≤ 4, the bound is derived from the absolute Leibniz sum using + /// Shewchuk-style error analysis (see `REFERENCES.md` \[8\] and the + /// per-constant docs on [`ERR_COEFF_2`](crate::ERR_COEFF_2), + /// [`ERR_COEFF_3`](crate::ERR_COEFF_3), and + /// [`ERR_COEFF_4`](crate::ERR_COEFF_4)). For D = 0 or 1, returns + /// `Some(0.0)` since the determinant computation is exact (no + /// arithmetic). /// /// This method does NOT require the `exact` feature — the bounds use /// pure f64 arithmetic and are useful for custom adaptive-precision logic. diff --git a/typos.toml b/typos.toml index 545032f..6c98b70 100644 --- a/typos.toml +++ b/typos.toml @@ -14,6 +14,9 @@ extend-exclude = [ ".pytest_cache/**", # Contains regex fragments (e.g. [Ss]upport, [Dd]elete) that trigger false positives. "cliff.toml", + # Auto-generated by `just changelog` (git-cliff); contains git commit + # short-SHAs like `a522abd` that can accidentally match dictionary words. + "CHANGELOG.md", ] [default.extend-words] From f8d80a01e9913f87e1b19b2ad5ffbc0994e2bfdb Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Mon, 20 Apr 2026 15:22:56 -0700 Subject: [PATCH 3/3] Changed: consolidate and expand const-evaluability tests via macros Refactor manual const-evaluation tests in `Ldlt` and `Matrix` into macros to standardize coverage across matrix dimensions 2 through 5. This ensures all determinant and norm variants are fully exercised at compile time. Refs: #86 --- src/ldlt.rs | 101 ++++++++++++++++++++++++--------------------- src/matrix.rs | 111 +++++++++++++++++++++++++++++++------------------- 2 files changed, 123 insertions(+), 89 deletions(-) diff --git a/src/ldlt.rs b/src/ldlt.rs index a2542e2..0dce894 100644 --- a/src/ldlt.rs +++ b/src/ldlt.rs @@ -558,54 +558,63 @@ mod tests { // some toolchains; we therefore construct `Ldlt` directly. // ----------------------------------------------------------------------- - #[test] - fn ldlt_det_const_eval_d2() { - const DET: f64 = { - // Diagonal D = [4.0, 0.25] ⇒ det = 1.0. - let mut factors = Matrix::<2>::identity(); - factors.rows[0][0] = 4.0; - factors.rows[1][1] = 0.25; - let ldlt = Ldlt::<2> { - factors, - tol: DEFAULT_SINGULAR_TOL, - }; - ldlt.det() - }; - assert!((DET - 1.0).abs() <= 1e-12); - } - - #[test] - fn ldlt_det_const_eval_d3() { - const DET: f64 = { - // Diagonal D = [2.0, 3.0, 5.0] ⇒ det = 30.0. - let mut factors = Matrix::<3>::identity(); - factors.rows[0][0] = 2.0; - factors.rows[1][1] = 3.0; - factors.rows[2][2] = 5.0; - let ldlt = Ldlt::<3> { - factors, - tol: DEFAULT_SINGULAR_TOL, - }; - ldlt.det() - }; - assert!((DET - 30.0).abs() <= 1e-12); - } + macro_rules! gen_ldlt_const_eval_tests { + ($d:literal) => { + paste! { + /// `Ldlt::det` must be fully const-evaluable. Setting + /// `factors[0][0] = 2.0` and leaving the remaining identity + /// diagonals at `1.0` gives `det = 2.0` for every `D ≥ 1`, + /// exercising the multiply-accumulate loop at each dimension. + #[test] + fn []() { + const DET: f64 = { + let mut factors = Matrix::<$d>::identity(); + factors.rows[0][0] = 2.0; + let ldlt = Ldlt::<$d> { + factors, + tol: DEFAULT_SINGULAR_TOL, + }; + ldlt.det() + }; + assert!((DET - 2.0).abs() <= 1e-12); + } - #[test] - fn ldlt_solve_vec_const_eval_d2() { - // Identity factors ⇒ solve_vec returns the RHS untouched. - const X: [f64; 2] = { - let ldlt = Ldlt::<2> { - factors: Matrix::<2>::identity(), - tol: DEFAULT_SINGULAR_TOL, - }; - let b = Vector::<2>::new([1.0, 2.0]); - match ldlt.solve_vec(b) { - Ok(v) => v.into_array(), - Err(_) => [0.0, 0.0], + /// `Ldlt::solve_vec` must be fully const-evaluable. Identity + /// factors with RHS `b = [1.0, 2.0, …, D]` round-trips `b` + /// unchanged, exercising the full forward sub / diagonal solve + /// / back sub pipeline inside a `const { … }` initializer. + #[test] + fn []() { + #[allow(clippy::cast_precision_loss)] + const X: [f64; $d] = { + let ldlt = Ldlt::<$d> { + factors: Matrix::<$d>::identity(), + tol: DEFAULT_SINGULAR_TOL, + }; + let mut b_arr = [0.0f64; $d]; + let mut i = 0; + while i < $d { + b_arr[i] = i as f64 + 1.0; + i += 1; + } + let b = Vector::<$d>::new(b_arr); + match ldlt.solve_vec(b) { + Ok(v) => v.into_array(), + Err(_) => [0.0f64; $d], + } + }; + #[allow(clippy::cast_precision_loss)] + for i in 0..$d { + let expected = i as f64 + 1.0; + assert!((X[i] - expected).abs() <= 1e-12); + } + } } }; - assert!((X[0] - 1.0).abs() <= 1e-12); - assert!((X[1] - 2.0).abs() <= 1e-12); } + + gen_ldlt_const_eval_tests!(2); + gen_ldlt_const_eval_tests!(3); + gen_ldlt_const_eval_tests!(4); + gen_ldlt_const_eval_tests!(5); } diff --git a/src/matrix.rs b/src/matrix.rs index a28b977..b38bf76 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -873,23 +873,36 @@ mod tests { ); } - #[test] - fn det_direct_const_eval_d2() { - // Const evaluation proves the function is truly const fn. - const DET: Option = { - let m = Matrix::<2>::from_rows([[1.0, 0.0], [0.0, 1.0]]); - m.det_direct() + // === det_direct const-evaluability tests (D = 2..=5) === + // + // Every dimension hits a distinct arm of the `match D { … }` body inside + // `det_direct`, so exercising each at compile time is the tightest + // const-fn proof available. + + macro_rules! gen_det_direct_const_eval_tests { + ($d:literal) => { + paste! { + /// `Matrix::::det_direct()` on the identity must const-evaluate + /// to `Some(1.0)` for every closed-form dimension `D ∈ {1, 2, 3, 4}`. + #[test] + fn []() { + const DET: Option = Matrix::<$d>::identity().det_direct(); + assert_eq!(DET, Some(1.0)); + } + } }; - assert_eq!(DET, Some(1.0)); } + gen_det_direct_const_eval_tests!(2); + gen_det_direct_const_eval_tests!(3); + gen_det_direct_const_eval_tests!(4); + #[test] - fn det_direct_const_eval_d3() { - const DET: Option = { - let m = Matrix::<3>::from_rows([[2.0, 0.0, 0.0], [0.0, 3.0, 0.0], [0.0, 0.0, 5.0]]); - m.det_direct() - }; - assert_eq!(DET, Some(30.0)); + fn det_direct_const_eval_d5_is_none() { + // D ≥ 5 has no closed-form arm; `det_direct` returns `None`. Verify + // that the wildcard arm is reachable in a `const { … }` context. + const DET: Option = Matrix::<5>::identity().det_direct(); + assert_eq!(DET, None); } // === det_errbound tests (no `exact` feature required) === @@ -909,46 +922,58 @@ mod tests { assert_eq!(Matrix::<5>::identity().det_errbound(), None); } - #[test] - fn det_errbound_const_eval_d2() { - // Const evaluation proves the function is truly const fn. - const BOUND: Option = { - let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); - m.det_errbound() - }; - assert!(BOUND.is_some()); - assert!(BOUND.unwrap() > 0.0); - } + // === det_errbound const-evaluability tests (D = 2..=5) === - #[test] - fn det_errbound_const_eval_d3() { - const BOUND: Option = { - let m = Matrix::<3>::from_rows([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]); - m.det_errbound() + macro_rules! gen_det_errbound_const_eval_tests { + ($d:literal) => { + paste! { + /// `Matrix::::det_errbound()` on the identity must const-evaluate + /// to `Some(bound)` with `bound > 0` for every closed-form dimension + /// `D ∈ {2, 3, 4}`. Each dimension hits a distinct arm of + /// `det_errbound` with a dimension-specific permanent computation. + #[test] + fn []() { + const BOUND: Option = Matrix::<$d>::identity().det_errbound(); + assert!(BOUND.is_some()); + assert!(BOUND.unwrap() > 0.0); + } + } }; - assert!(BOUND.is_some()); - assert!(BOUND.unwrap() > 0.0); } + gen_det_errbound_const_eval_tests!(2); + gen_det_errbound_const_eval_tests!(3); + gen_det_errbound_const_eval_tests!(4); + #[test] - fn inf_norm_const_eval_d2() { - // Maximum absolute row sum: row 0 = 3, row 1 = 7 ⇒ 7. - const NORM: f64 = { - let m = Matrix::<2>::from_rows([[1.0, -2.0], [3.0, 4.0]]); - m.inf_norm() - }; - assert!((NORM - 7.0).abs() <= 1e-12); + fn det_errbound_const_eval_d5_is_none() { + // D ≥ 5 has no fast-filter bound; `det_errbound` returns `None`. + const BOUND: Option = Matrix::<5>::identity().det_errbound(); + assert_eq!(BOUND, None); } - #[test] - fn inf_norm_const_eval_d3() { - const NORM: f64 = { - let m = Matrix::<3>::identity(); - m.inf_norm() + // === inf_norm const-evaluability tests (D = 2..=5) === + + macro_rules! gen_inf_norm_const_eval_tests { + ($d:literal) => { + paste! { + /// `Matrix::::inf_norm()` on the identity must const-evaluate + /// to `1.0` for every `D ≥ 1` — each row has a single `1.0` + /// entry, so the max absolute row sum is exactly `1.0`. + #[test] + fn []() { + const NORM: f64 = Matrix::<$d>::identity().inf_norm(); + assert!((NORM - 1.0).abs() <= 1e-12); + } + } }; - assert!((NORM - 1.0).abs() <= 1e-12); } + gen_inf_norm_const_eval_tests!(2); + gen_inf_norm_const_eval_tests!(3); + gen_inf_norm_const_eval_tests!(4); + gen_inf_norm_const_eval_tests!(5); + // === inf_norm NaN / Inf propagation (regression tests for #85) === macro_rules! gen_inf_norm_nonfinite_tests {