From d45572a05488dbcb73b66a4b36a0e0cfc4656ca9 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Fri, 6 Mar 2026 20:36:04 -0800 Subject: [PATCH 1/5] feat: add det_sign_exact() behind "exact" feature flag (#33) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two-stage adaptive-precision determinant sign for Matrix: - Stage 1 (D ≤ 4): f64 fast filter via det_direct() + Shewchuk-style error bounds against the matrix permanent. Resolves sign without allocating for well-conditioned inputs. - Stage 2: exact Bareiss algorithm (fraction-free Gaussian elimination) in BigRational for guaranteed-correct sign. Closes #33 --- AGENTS.md | 4 + Cargo.lock | 4 +- Cargo.toml | 5 +- README.md | 32 +++- REFERENCES.md | 17 ++ justfile | 12 +- src/exact.rs | 404 ++++++++++++++++++++++++++++++++++++++++ src/lib.rs | 2 + tests/proptest_exact.rs | 72 +++++++ typos.toml | 1 + 10 files changed, 547 insertions(+), 6 deletions(-) create mode 100644 src/exact.rs create mode 100644 tests/proptest_exact.rs diff --git a/AGENTS.md b/AGENTS.md index 8f77f26..971b249 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -24,10 +24,12 @@ When making changes in this repo, prioritize (in order): - Format: `cargo fmt` (or `just fmt`) - Integration tests: `just test-integration` - Lint (Clippy): `cargo clippy --all-targets --all-features -- -D warnings` (or `just clippy`) +- Lint (Clippy, exact feature): `cargo clippy --features exact --all-targets -- -D warnings` (or `just clippy-exact`) - Lint/validate: `just check` - Pre-commit validation: `just ci` - Python tests: `just test-python` - Run a single test (by name filter): `cargo test solve_2x2_basic` (or the full path: `cargo test lu::tests::solve_2x2_basic`) +- Run exact-feature tests: `cargo test --features exact --verbose` (or `just test-exact`) - Run examples: `just examples` (or `cargo run --example det_5x5` / `cargo run --example solve_5x5` / `cargo run --example const_det_4x4`) - Spell check: `just spell-check` (uses `typos.toml` at repo root; add false positives to `[default.extend-words]`) @@ -40,6 +42,8 @@ When making changes in this repo, prioritize (in order): - `src/matrix.rs`: `Matrix` (`[[f64; D]; D]`) + helpers (`get`, `set`, `inf_norm`, `det`, `det_direct`) - `src/lu.rs`: `Lu` factorization with partial pivoting (`solve_vec`, `det`) - `src/ldlt.rs`: `Ldlt` factorization without pivoting for symmetric SPD/PSD matrices (`solve_vec`, `det`) + - `src/exact.rs`: `det_sign_exact()` — adaptive-precision determinant sign + (Shewchuk-style f64 filter + Bareiss in `BigRational`); `features = ["exact"]` - A minimal `justfile` exists for common workflows (see `just --list`). - The public API re-exports these items from `src/lib.rs`. - Dev-only benchmarks live in `benches/vs_linalg.rs` (Criterion + nalgebra/faer comparison). diff --git a/Cargo.lock b/Cargo.lock index 9950998..5f9ba00 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -675,12 +675,14 @@ dependencies = [ [[package]] name = "la-stack" -version = "0.1.3" +version = "0.2.0" dependencies = [ "approx", "criterion", "faer", "nalgebra", + "num-bigint", + "num-rational", "pastey", "proptest", ] diff --git a/Cargo.toml b/Cargo.toml index 4c2af4a..8f36a48 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "la-stack" -version = "0.1.3" +version = "0.2.0" edition = "2024" rust-version = "1.94" license = "BSD-3-Clause" @@ -17,6 +17,8 @@ keywords = [ "linear-algebra", "geometry", "const-generics" ] criterion = { version = "0.8.2", features = [ "html_reports" ], optional = true } faer = { version = "0.24.0", default-features = false, features = [ "std", "linalg" ], optional = true } nalgebra = { version = "0.34.1", optional = true } +num-bigint = { version = "0.4", optional = true } +num-rational = { version = "0.4", features = [ "num-bigint-std" ], optional = true } [dev-dependencies] approx = "0.5.1" @@ -26,6 +28,7 @@ proptest = "1.10.0" [features] default = [ ] bench = [ "criterion", "faer", "nalgebra" ] +exact = [ "num-bigint", "num-rational" ] [[bench]] name = "vs_linalg" diff --git a/README.md b/README.md index da85d72..6d4e4be 100644 --- a/README.md +++ b/README.md @@ -53,7 +53,7 @@ Add this to your `Cargo.toml`: ```toml [dependencies] -la-stack = "0.1" +la-stack = "0.2" ``` Solve a 5×5 system via LU: @@ -130,17 +130,45 @@ assert_eq!(DET, Some(30.0)); The public `det()` method automatically dispatches through the closed-form path for D ≤ 4 and falls back to LU for D ≥ 5 — no API change needed. +## 🔬 Exact determinant sign (`"exact"` feature) + +Enable the `exact` Cargo feature for `det_sign_exact()`, which returns the +provably correct sign (−1, 0, or +1) of the determinant using adaptive-precision +arithmetic: + +```toml +[dependencies] +la-stack = { version = "0.2", features = ["exact"] } +``` + +```rust,ignore +use la_stack::prelude::*; + +let m = Matrix::<3>::from_rows([ + [1.0, 2.0, 3.0], + [4.0, 5.0, 6.0], + [7.0, 8.0, 9.0], +]); +assert_eq!(m.det_sign_exact(), 0); // exactly singular +``` + +For D ≤ 4, a fast f64 filter (error-bounded `det_direct()`) resolves the sign +without allocating. Only near-degenerate or large (D ≥ 5) matrices fall through +to the exact Bareiss algorithm in `BigRational`. + ## 🧩 API at a glance | Type | Storage | Purpose | Key methods | |---|---|---|---| | `Vector` | `[f64; D]` | Fixed-length vector | `new`, `zero`, `dot`, `norm2_sq` | -| `Matrix` | `[[f64; D]; D]` | Fixed-size square matrix | `from_rows`, `zero`, `identity`, `lu`, `ldlt`, `det`, `det_direct` | +| `Matrix` | `[[f64; D]; D]` | Fixed-size square matrix | `from_rows`, `zero`, `identity`, `lu`, `ldlt`, `det`, `det_direct`, `det_sign_exact`¹ | | `Lu` | `Matrix` + pivot array | Factorization for solves/det | `solve_vec`, `det` | | `Ldlt` | `Matrix` | Factorization for symmetric SPD/PSD solves/det | `solve_vec`, `det` | Storage shown above reflects the current `f64` implementation. +¹ Requires `features = ["exact"]`. + ## 📋 Examples The `examples/` directory contains small, runnable programs: diff --git a/REFERENCES.md b/REFERENCES.md index 89d4a8c..05de553 100644 --- a/REFERENCES.md +++ b/REFERENCES.md @@ -26,6 +26,11 @@ 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) + +`det_sign_exact()` uses a Shewchuk-style f64 error-bound filter [8] backed by exact Bareiss +elimination [7] in `BigRational`. See `src/exact.rs` for the full architecture description. + ## References ### LU / Gaussian elimination @@ -53,3 +58,15 @@ matrices, see [6]. 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 sign (`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) · + [PDF](https://www.ams.org/journals/mcom/1968-22-103/S0025-5718-1968-0226829-0/S0025-5718-1968-0226829-0.pdf) +8. Shewchuk, Jonathan Richard. "Adaptive Precision Floating-Point Arithmetic and Fast + Robust Geometric Predicates." *Discrete & Computational Geometry* 18.3 (1997): 305–363. + [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. diff --git a/justfile b/justfile index b238f7c..5cae377 100644 --- a/justfile +++ b/justfile @@ -127,9 +127,13 @@ check-fast: # CI simulation: comprehensive validation (matches CI expectations) # Runs: checks + all tests (Rust + Python) + examples + bench compile -ci: check bench-compile test-all examples +ci: check bench-compile test-all examples clippy-exact @echo "🎯 CI checks complete!" +# Clippy for the "exact" feature (catches feature-gated lint issues) +clippy-exact: + cargo clippy --features exact --all-targets -- -D warnings -W clippy::pedantic + # Clean build artifacts clean: cargo clean @@ -500,9 +504,13 @@ test: cargo test --lib --verbose cargo test --doc --verbose -test-all: test test-integration test-python +test-all: test test-integration test-exact test-python @echo "✅ All tests passed" +# Tests for the "exact" feature (det_sign_exact + BigRational Bareiss) +test-exact: + cargo test --features exact --verbose + test-integration: cargo test --tests --verbose diff --git a/src/exact.rs b/src/exact.rs new file mode 100644 index 0000000..a10c955 --- /dev/null +++ b/src/exact.rs @@ -0,0 +1,404 @@ +//! Exact determinant sign via arbitrary-precision rational arithmetic. +//! +//! This module is only compiled when the `"exact"` Cargo feature is enabled. +//! +//! # Architecture +//! +//! `det_sign_exact` uses a two-stage adaptive-precision approach inspired by +//! Shewchuk's robust geometric predicates: +//! +//! 1. **Fast filter (D ≤ 4)**: compute `det_direct()` and a conservative error +//! bound. If `|det| > bound`, the f64 sign is provably correct — return +//! immediately without allocating. +//! 2. **Exact fallback**: convert entries to `BigRational` and run the Bareiss +//! algorithm (fraction-free Gaussian elimination) for a guaranteed-correct +//! sign. + +use num_bigint::{BigInt, Sign}; +use num_rational::BigRational; + +use crate::matrix::Matrix; + +// --------------------------------------------------------------------------- +// Error-bound constants for the f64 fast filter. +// +// Each constant bounds 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; +// over-estimating just means we fall through to Bareiss more often. +// --------------------------------------------------------------------------- + +const EPS: f64 = f64::EPSILON; // 2^-52 + +/// D=2: one f64 multiply + one FMA → 2 rounding events. +const ERR_COEFF_2: f64 = 3.0 * EPS + 16.0 * EPS * EPS; + +/// D=3: three 2×2 FMA minors + nested FMA combination. +const ERR_COEFF_3: f64 = 8.0 * EPS + 64.0 * EPS * EPS; + +/// D=4: six hoisted 2×2 minors → four 3×3 cofactors → FMA row combination. +const ERR_COEFF_4: f64 = 12.0 * EPS + 128.0 * EPS * EPS; + +/// Conservative absolute error bound for `det_direct()`, or `None` for D ≥ 5. +/// +/// Returns `Some(bound)` such that `|det_direct() - det_exact| ≤ bound`. +fn det_errbound(m: &Matrix) -> Option { + match D { + 0 | 1 => Some(0.0), // No arithmetic — result is exact. + 2 => { + let r = &m.rows; + let permanent = (r[0][0] * r[1][1]).abs() + (r[0][1] * r[1][0]).abs(); + Some(ERR_COEFF_2 * permanent) + } + 3 => { + let r = &m.rows; + let pm00 = (r[1][1] * r[2][2]).abs() + (r[1][2] * r[2][1]).abs(); + let pm01 = (r[1][0] * r[2][2]).abs() + (r[1][2] * r[2][0]).abs(); + let pm02 = (r[1][0] * r[2][1]).abs() + (r[1][1] * r[2][0]).abs(); + let permanent = r[0][2] + .abs() + .mul_add(pm02, r[0][1].abs().mul_add(pm01, r[0][0].abs() * pm00)); + Some(ERR_COEFF_3 * permanent) + } + 4 => { + let r = &m.rows; + // 2×2 minor permanents from rows 2–3. + 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(); + // 3×3 cofactor permanents from row 1. + 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)); + let permanent = 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)), + ); + Some(ERR_COEFF_4 * permanent) + } + _ => None, + } +} + +/// Convert an `f64` to an exact `BigRational`. +/// +/// Every finite `f64` is exactly representable as a rational number (`m × 2^e`), +/// so this conversion is lossless. +/// +/// # Panics +/// Panics if `x` is NaN or infinite. +fn f64_to_bigrational(x: f64) -> BigRational { + BigRational::from_float(x).expect("non-finite matrix entry in det_sign_exact") +} + +/// Compute the exact determinant of a `D×D` matrix using the Bareiss algorithm +/// (fraction-free Gaussian elimination) in `BigRational` arithmetic. +/// +/// Returns the determinant as an exact `BigRational`. +fn bareiss_det(m: &Matrix) -> BigRational { + if D == 0 { + return BigRational::from_integer(BigInt::from(1)); + } + if D == 1 { + return f64_to_bigrational(m.rows[0][0]); + } + + // Convert f64 entries to exact BigRational. + let mut a: Vec> = Vec::with_capacity(D); + for r in 0..D { + let mut row = Vec::with_capacity(D); + for c in 0..D { + row.push(f64_to_bigrational(m.rows[r][c])); + } + a.push(row); + } + + let zero = BigRational::from_integer(BigInt::from(0)); + let mut prev_pivot = BigRational::from_integer(BigInt::from(1)); + let mut sign: i8 = 1; + + for k in 0..D { + // Partial pivoting: find non-zero entry in column k at or below row k. + if a[k][k] == zero { + let mut found = false; + for i in (k + 1)..D { + if a[i][k] != zero { + a.swap(k, i); + sign = -sign; + found = true; + break; + } + } + if !found { + // Entire column below (and including) diagonal is zero → singular. + return zero; + } + } + + // Bareiss elimination for rows below k. + for i in (k + 1)..D { + for j in (k + 1)..D { + // a[i][j] = (a[k][k] * a[i][j] - a[i][k] * a[k][j]) / prev_pivot + a[i][j] = (&a[k][k] * &a[i][j] - &a[i][k] * &a[k][j]) / &prev_pivot; + } + // Zero out the eliminated column entry (not needed for det, but keeps + // the matrix consistent for potential debugging). + a[i][k] = zero.clone(); + } + + prev_pivot = a[k][k].clone(); + } + + let det = &a[D - 1][D - 1]; + if sign < 0 { -det } else { det.clone() } +} + +impl Matrix { + /// Exact determinant sign using adaptive-precision arithmetic. + /// + /// Returns `1` if `det > 0`, `-1` if `det < 0`, and `0` if `det == 0` (singular). + /// + /// For D ≤ 4, a fast f64 filter is tried first: `det_direct()` is compared + /// against a conservative error bound derived from the matrix permanent. + /// If the f64 result clearly exceeds the bound, the sign is returned + /// immediately without allocating. Otherwise (and always for D ≥ 5), the + /// Bareiss algorithm runs in exact [`BigRational`] arithmetic. + /// + /// # When to use + /// + /// Use this when the sign of the determinant must be correct regardless of + /// floating-point conditioning (e.g. geometric predicates on near-degenerate + /// configurations). For well-conditioned matrices the fast filter resolves + /// the sign without touching `BigRational`, so the overhead is minimal. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// let m = Matrix::<3>::from_rows([ + /// [1.0, 2.0, 3.0], + /// [4.0, 5.0, 6.0], + /// [7.0, 8.0, 9.0], + /// ]); + /// // This matrix is singular (row 3 = row 1 + row 2 in exact arithmetic). + /// assert_eq!(m.det_sign_exact(), 0); + /// + /// assert_eq!(Matrix::<3>::identity().det_sign_exact(), 1); + /// ``` + /// + /// # Panics + /// Panics if any matrix entry is NaN or infinite. + #[must_use] + pub fn det_sign_exact(&self) -> i8 { + // Stage 1: f64 fast filter for D ≤ 4. + if let (Some(det_f64), Some(err)) = (self.det_direct(), det_errbound(self)) { + if det_f64 > err { + return 1; + } + if det_f64 < -err { + return -1; + } + } + + // Stage 2: exact Bareiss fallback. + let det = bareiss_det(self); + match det.numer().sign() { + Sign::Plus => 1, + Sign::Minus => -1, + Sign::NoSign => 0, + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::DEFAULT_PIVOT_TOL; + + #[test] + fn det_sign_exact_d0_is_positive() { + assert_eq!(Matrix::<0>::zero().det_sign_exact(), 1); + } + + #[test] + fn det_sign_exact_d1_positive() { + let m = Matrix::<1>::from_rows([[42.0]]); + assert_eq!(m.det_sign_exact(), 1); + } + + #[test] + fn det_sign_exact_d1_negative() { + let m = Matrix::<1>::from_rows([[-3.5]]); + assert_eq!(m.det_sign_exact(), -1); + } + + #[test] + fn det_sign_exact_d1_zero() { + let m = Matrix::<1>::from_rows([[0.0]]); + assert_eq!(m.det_sign_exact(), 0); + } + + #[test] + fn det_sign_exact_identity_2d() { + assert_eq!(Matrix::<2>::identity().det_sign_exact(), 1); + } + + #[test] + fn det_sign_exact_identity_3d() { + assert_eq!(Matrix::<3>::identity().det_sign_exact(), 1); + } + + #[test] + fn det_sign_exact_identity_4d() { + assert_eq!(Matrix::<4>::identity().det_sign_exact(), 1); + } + + #[test] + fn det_sign_exact_identity_5d() { + assert_eq!(Matrix::<5>::identity().det_sign_exact(), 1); + } + + #[test] + fn det_sign_exact_singular_duplicate_rows() { + let m = Matrix::<3>::from_rows([ + [1.0, 2.0, 3.0], + [4.0, 5.0, 6.0], + [1.0, 2.0, 3.0], // duplicate of row 0 + ]); + assert_eq!(m.det_sign_exact(), 0); + } + + #[test] + fn det_sign_exact_singular_linear_combination() { + // Row 2 = row 0 + row 1 in exact arithmetic. + let m = Matrix::<3>::from_rows([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [5.0, 7.0, 9.0]]); + assert_eq!(m.det_sign_exact(), 0); + } + + #[test] + fn det_sign_exact_negative_det_row_swap() { + // Swapping two rows of the identity negates the determinant. + let m = Matrix::<3>::from_rows([[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]); + assert_eq!(m.det_sign_exact(), -1); + } + + #[test] + fn det_sign_exact_negative_det_known() { + // det([[1,2],[3,4]]) = 1*4 - 2*3 = -2 + let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); + assert_eq!(m.det_sign_exact(), -1); + } + + #[test] + fn det_sign_exact_agrees_with_det_for_spd() { + // SPD matrix → positive determinant. + let m = Matrix::<3>::from_rows([[4.0, 2.0, 0.0], [2.0, 5.0, 1.0], [0.0, 1.0, 3.0]]); + assert_eq!(m.det_sign_exact(), 1); + assert!(m.det(DEFAULT_PIVOT_TOL).unwrap() > 0.0); + } + + /// Near-singular matrix with an exact perturbation. + /// + /// The base matrix `[[1,2,3],[4,5,6],[7,8,9]]` is exactly singular (rows in + /// arithmetic progression). Adding `2^-50` to entry (0,0) makes + /// `det = 2^-50 × cofactor(0,0) = 2^-50 × (5×9 − 6×8) = −3 × 2^-50 < 0`. + /// Both f64 `det_direct()` and `det_sign_exact()` should agree here. + #[test] + fn det_sign_exact_near_singular_perturbation() { + let perturbation = f64::from_bits(0x3CD0_0000_0000_0000); // 2^-50 + let m = Matrix::<3>::from_rows([ + [1.0 + perturbation, 2.0, 3.0], + [4.0, 5.0, 6.0], + [7.0, 8.0, 9.0], + ]); + // Exact: det = perturbation × (5×9 − 6×8) = perturbation × (−3) < 0. + assert_eq!(m.det_sign_exact(), -1); + } + + /// For D ≤ 4, well-conditioned matrices should hit the fast filter + /// and never allocate `BigRational`. We can't directly observe this, + /// but we verify correctness for a range of known signs. + #[test] + fn det_sign_exact_fast_filter_positive_4x4() { + let m = Matrix::<4>::from_rows([ + [2.0, 1.0, 0.0, 0.0], + [1.0, 3.0, 1.0, 0.0], + [0.0, 1.0, 4.0, 1.0], + [0.0, 0.0, 1.0, 5.0], + ]); + // SPD tridiagonal → positive det. + assert_eq!(m.det_sign_exact(), 1); + } + + #[test] + fn det_sign_exact_fast_filter_negative_4x4() { + // Swap rows 0 and 1 of the above → negate det. + let m = Matrix::<4>::from_rows([ + [1.0, 3.0, 1.0, 0.0], + [2.0, 1.0, 0.0, 0.0], + [0.0, 1.0, 4.0, 1.0], + [0.0, 0.0, 1.0, 5.0], + ]); + assert_eq!(m.det_sign_exact(), -1); + } + + #[test] + fn det_sign_exact_subnormal_entries() { + // Subnormal f64 values should convert losslessly. + let tiny = 5e-324_f64; // smallest positive subnormal + assert!(tiny.is_subnormal()); + + let m = Matrix::<2>::from_rows([[tiny, 0.0], [0.0, tiny]]); + // det = tiny^2 > 0 + assert_eq!(m.det_sign_exact(), 1); + } + + #[test] + #[should_panic(expected = "non-finite matrix entry")] + fn det_sign_exact_panics_on_nan() { + let m = Matrix::<2>::from_rows([[f64::NAN, 0.0], [0.0, 1.0]]); + let _ = m.det_sign_exact(); + } + + #[test] + #[should_panic(expected = "non-finite matrix entry")] + fn det_sign_exact_panics_on_infinity() { + let m = Matrix::<2>::from_rows([[f64::INFINITY, 0.0], [0.0, 1.0]]); + let _ = m.det_sign_exact(); + } + + #[test] + fn det_sign_exact_pivot_needed_in_first_column() { + // First column has zero on diagonal but non-zero below → needs pivoting. + let m = Matrix::<3>::from_rows([[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]); + assert_eq!(m.det_sign_exact(), -1); + } + + #[test] + fn det_sign_exact_5x5_known() { + // det of a permutation matrix with two swaps = +1 (even permutation). + let m = Matrix::<5>::from_rows([ + [0.0, 1.0, 0.0, 0.0, 0.0], + [1.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 1.0, 0.0], + [0.0, 0.0, 1.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 1.0], + ]); + // Two transpositions → even permutation → det = +1 + assert_eq!(m.det_sign_exact(), 1); + } +} diff --git a/src/lib.rs b/src/lib.rs index 3210188..46e72e5 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -48,6 +48,8 @@ mod readme_doctests { fn det_5x5_ldlt_example() {} } +#[cfg(feature = "exact")] +mod exact; mod ldlt; mod lu; mod matrix; diff --git a/tests/proptest_exact.rs b/tests/proptest_exact.rs new file mode 100644 index 0000000..1ef388c --- /dev/null +++ b/tests/proptest_exact.rs @@ -0,0 +1,72 @@ +//! Property-based tests for the `det_sign_exact` API (requires `exact` feature). + +#![cfg(feature = "exact")] + +use pastey::paste; +use proptest::prelude::*; + +use la_stack::prelude::*; + +fn small_nonzero_f64() -> impl Strategy { + prop_oneof![(-1000i16..=-1i16), (1i16..=1000i16)].prop_map(|x| f64::from(x) / 10.0) +} + +macro_rules! gen_det_sign_exact_proptests { + ($d:literal) => { + paste! { + proptest! { + #![proptest_config(ProptestConfig::with_cases(64))] + + #[test] + fn []( + diag in proptest::array::[](small_nonzero_f64()), + ) { + // Diagonal matrix: determinant sign = product of diagonal signs. + let mut rows = [[0.0f64; $d]; $d]; + for i in 0..$d { + rows[i][i] = diag[i]; + } + let m = Matrix::<$d>::from_rows(rows); + + let exact_sign = m.det_sign_exact(); + + // Expected sign from the product of diagonal entries. + let neg_count = diag.iter().filter(|&&x| x < 0.0).count(); + let expected_sign: i8 = if neg_count % 2 == 0 { 1 } else { -1 }; + + prop_assert_eq!(exact_sign, expected_sign); + } + + #[test] + fn []( + diag in proptest::array::[](small_nonzero_f64()), + ) { + // For well-conditioned diagonal matrices, det().signum() + // should agree with det_sign_exact(). + let mut rows = [[0.0f64; $d]; $d]; + for i in 0..$d { + rows[i][i] = diag[i]; + } + let m = Matrix::<$d>::from_rows(rows); + + let exact_sign = m.det_sign_exact(); + let fp_det = m.det(DEFAULT_PIVOT_TOL).unwrap(); + let fp_sign: i8 = if fp_det > 0.0 { + 1 + } else if fp_det < 0.0 { + -1 + } else { + 0 + }; + + prop_assert_eq!(exact_sign, fp_sign); + } + } + } + }; +} + +gen_det_sign_exact_proptests!(2); +gen_det_sign_exact_proptests!(3); +gen_det_sign_exact_proptests!(4); +gen_det_sign_exact_proptests!(5); diff --git a/typos.toml b/typos.toml index 2a83eba..0af273f 100644 --- a/typos.toml +++ b/typos.toml @@ -29,3 +29,4 @@ Trefethen = "Trefethen" faer = "faer" frhs = "frhs" nrhs = "nrhs" +numer = "numer" From abfa297d1e8e771a60d1cc9a2a9227011876a965 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Sat, 7 Mar 2026 06:42:04 -0800 Subject: [PATCH 2/5] =?UTF-8?q?refactor:=20polish=20exact=20feature=20?= =?UTF-8?q?=E2=80=94=20docs,=20CI,=20example,=20review=20fixes?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Add exact_sign_3x3 example (near-singular 3×3 with 2^-50 perturbation) - Add docs.rs metadata to build published docs with the exact feature - Add non-finite assert guard in det_sign_exact fast filter - Add --features exact to Windows CI job and doc-check recipe - Add --features exact to tarpaulin coverage runs - Update README: zero-dep-by-default wording, scalar types section, example descriptions, design goals, references wording - Update Cargo.toml deps comment - Remove redundant clippy-exact from ci target --- .github/workflows/ci.yml | 1 + Cargo.toml | 9 +++++++- README.md | 24 +++++++++++++------- examples/exact_sign_3x3.rs | 45 ++++++++++++++++++++++++++++++++++++++ justfile | 8 ++++--- src/exact.rs | 4 ++++ 6 files changed, 79 insertions(+), 12 deletions(-) create mode 100644 examples/exact_sign_3x3.rs diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 7f01214..f799acc 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -232,3 +232,4 @@ jobs: cargo test --lib --verbose cargo test --doc --verbose cargo test --tests --verbose + cargo test --features exact --verbose diff --git a/Cargo.toml b/Cargo.toml index 8f36a48..0df0e97 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -13,7 +13,7 @@ categories = [ "mathematics", "science" ] keywords = [ "linear-algebra", "geometry", "const-generics" ] [dependencies] -# Intentionally empty (bench-only deps are optional features below) +# All runtime deps are optional; see [features] below. criterion = { version = "0.8.2", features = [ "html_reports" ], optional = true } faer = { version = "0.24.0", default-features = false, features = [ "std", "linalg" ], optional = true } nalgebra = { version = "0.34.1", optional = true } @@ -30,6 +30,10 @@ default = [ ] bench = [ "criterion", "faer", "nalgebra" ] exact = [ "num-bigint", "num-rational" ] +[[example]] +name = "exact_sign_3x3" +required-features = [ "exact" ] + [[bench]] name = "vs_linalg" harness = false @@ -39,6 +43,9 @@ required-features = [ "bench" ] lto = "fat" codegen-units = 1 +[package.metadata.docs.rs] +features = [ "exact" ] + [lints.rust] unsafe_code = "forbid" missing_docs = "warn" diff --git a/README.md b/README.md index 6d4e4be..376dfa2 100644 --- a/README.md +++ b/README.md @@ -31,7 +31,8 @@ while keeping the API intentionally small and explicit. - ✅ Const-generic dimensions (no dynamic sizes) - ✅ `const fn` where possible (compile-time evaluation of determinants, dot products, etc.) - ✅ Explicit algorithms (LU, solve, determinant) -- ✅ No runtime dependencies (dev-dependencies are for contributors only) +- ✅ Robust geometric predicates via optional exact arithmetic (`det_sign_exact`) +- ✅ No runtime dependencies by default (optional features may add deps) - ✅ Stack storage only (no heap allocation in core types) - ✅ `unsafe` forbidden @@ -44,8 +45,8 @@ while keeping the API intentionally small and explicit. ## 🔢 Scalar types Today, the core types are implemented for `f64`. The intent is to support `f32` and `f64` -(and `f128` if/when Rust gains a stable primitive for it). Longer term, we may add optional -arbitrary-precision support (e.g. via `rug`) depending on performance. +(and `f128` if/when Rust gains a stable primitive for it). Arbitrary-precision arithmetic +is available via the optional `"exact"` feature (see below). ## 🚀 Quickstart @@ -132,9 +133,10 @@ for D ≤ 4 and falls back to LU for D ≥ 5 — no API change needed. ## 🔬 Exact determinant sign (`"exact"` feature) -Enable the `exact` Cargo feature for `det_sign_exact()`, which returns the -provably correct sign (−1, 0, or +1) of the determinant using adaptive-precision -arithmetic: +The default build has **zero runtime dependencies**. Enable the optional +`exact` Cargo feature to add `det_sign_exact()`, which returns the provably +correct sign (−1, 0, or +1) of the determinant using adaptive-precision +arithmetic (this pulls in `num-bigint` and `num-rational` for `BigRational`): ```toml [dependencies] @@ -173,12 +175,18 @@ Storage shown above reflects the current `f64` implementation. The `examples/` directory contains small, runnable programs: +- **`solve_5x5`** — solve a 5×5 system via LU with partial pivoting +- **`det_5x5`** — determinant of a 5×5 matrix via LU +- **`const_det_4x4`** — compile-time 4×4 determinant via `det_direct()` +- **`exact_sign_3x3`** — exact determinant sign of a near-singular 3×3 matrix (requires `exact` feature) + ```bash just examples -# or: +# or individually: cargo run --example solve_5x5 cargo run --example det_5x5 cargo run --example const_det_4x4 +cargo run --features exact --example exact_sign_3x3 ``` ## 🤝 Contributing @@ -202,7 +210,7 @@ If you use this library in academic work, please cite it using [CITATION.cff](CI ## 📚 References -For canonical references to LU / `LDL^T` algorithms used by this crate, see [REFERENCES.md](REFERENCES.md). +For canonical references to the algorithms used by this crate, see [REFERENCES.md](REFERENCES.md). ## 📊 Benchmarks (vs nalgebra/faer) diff --git a/examples/exact_sign_3x3.rs b/examples/exact_sign_3x3.rs new file mode 100644 index 0000000..4a9a078 --- /dev/null +++ b/examples/exact_sign_3x3.rs @@ -0,0 +1,45 @@ +//! Exact determinant sign for a near-singular 3×3 matrix. +//! +//! This example demonstrates `det_sign_exact()`, which uses adaptive-precision +//! arithmetic to return the provably correct sign of the determinant — even when +//! the matrix is so close to singular that f64 rounding could give the wrong answer. +//! +//! Run with: `cargo run --features exact --example exact_sign_3x3` + +use la_stack::prelude::*; + +fn main() { + // Base matrix: rows in arithmetic progression → exactly singular (det = 0). + // [[1, 2, 3], + // [4, 5, 6], + // [7, 8, 9]] + // + // Perturb entry (0,0) by 2^-50 ≈ 8.9e-16. + // Exact det = 2^-50 × cofactor(0,0) = 2^-50 × (5×9 − 6×8) = −3 × 2^-50 < 0. + let perturbation = f64::from_bits(0x3CD0_0000_0000_0000); // 2^-50 + let m = Matrix::<3>::from_rows([ + [1.0 + perturbation, 2.0, 3.0], + [4.0, 5.0, 6.0], + [7.0, 8.0, 9.0], + ]); + + let sign = m.det_sign_exact(); + let det_f64 = m.det(DEFAULT_PIVOT_TOL).unwrap(); + + println!("Near-singular 3×3 matrix (perturbation = 2^-50 ≈ {perturbation:.2e}):"); + for r in 0..3 { + print!(" ["); + for c in 0..3 { + if c > 0 { + print!(", "); + } + print!("{:22.18}", m.get(r, c).unwrap()); + } + println!("]"); + } + println!(); + println!("f64 det() = {det_f64:+.6e}"); + println!("det_sign_exact() = {sign}"); + println!(); + println!("The exact sign is −1 (negative), matching the analytical result."); +} diff --git a/justfile b/justfile index 5cae377..4f6e578 100644 --- a/justfile +++ b/justfile @@ -127,7 +127,7 @@ check-fast: # CI simulation: comprehensive validation (matches CI expectations) # Runs: checks + all tests (Rust + Python) + examples + bench compile -ci: check bench-compile test-all examples clippy-exact +ci: check bench-compile test-all examples @echo "🎯 CI checks complete!" # Clippy for the "exact" feature (catches feature-gated lint issues) @@ -149,6 +149,7 @@ clippy: # Common tarpaulin arguments for all coverage runs # Note: -t 300 sets per-test timeout to 5 minutes (needed for slow CI environments) _coverage_base_args := '''--exclude-files 'benches/*' --exclude-files 'examples/*' \ + --features exact \ --workspace --lib --tests \ -t 300 --verbose --implicit-test-threads''' @@ -183,15 +184,16 @@ coverage-ci: default: @just --list -# Documentation build check +# Documentation build check (includes exact feature for full API coverage) doc-check: - RUSTDOCFLAGS='-D warnings' cargo doc --no-deps + RUSTDOCFLAGS='-D warnings' cargo doc --no-deps --features exact # Examples examples: cargo run --quiet --example det_5x5 cargo run --quiet --example solve_5x5 cargo run --quiet --example const_det_4x4 + cargo run --quiet --features exact --example exact_sign_3x3 # Fix (mutating): apply formatters/auto-fixes fix: toml-fmt fmt python-fix shell-fmt markdown-fix yaml-fix diff --git a/src/exact.rs b/src/exact.rs index a10c955..796ebaf 100644 --- a/src/exact.rs +++ b/src/exact.rs @@ -206,6 +206,10 @@ impl Matrix { pub fn det_sign_exact(&self) -> i8 { // Stage 1: f64 fast filter for D ≤ 4. if let (Some(det_f64), Some(err)) = (self.det_direct(), det_errbound(self)) { + assert!( + det_f64.is_finite(), + "non-finite matrix entry in det_sign_exact" + ); if det_f64 > err { return 1; } From ecd0a30d1e580ae47f215b811a60110b8f94223a Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Sat, 7 Mar 2026 06:58:00 -0800 Subject: [PATCH 3/5] Changed: make Codacy analysis non-blocking on pull requests (internal) Allow the Codacy analysis step to fail without blocking CI on pull requests, mitigating the impact of transient service outages. Ensures that subsequent SARIF processing steps always execute and surfaces a warning when the analysis fails. Refs: #33 --- .github/workflows/codacy.yml | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/.github/workflows/codacy.yml b/.github/workflows/codacy.yml index dfac01c..3b7135a 100644 --- a/.github/workflows/codacy.yml +++ b/.github/workflows/codacy.yml @@ -77,6 +77,7 @@ jobs: # Execute Codacy Analysis CLI and generate a SARIF output with # the security issues identified during the analysis - name: Run Codacy Analysis CLI + id: codacy_analysis uses: codacy/codacy-analysis-cli-action@562ee3e92b8e92df8b67e0a5ff8aa8e261919c08 with: # Check https://github.com/codacy/codacy-analysis-cli#project-token @@ -94,9 +95,18 @@ jobs: # Force 0 exit code to allow SARIF file generation # This will handover control about PR rejection to the GitHub side max-allowed-issues: 2147483647 + # Codacy can fail transiently on PRs (e.g. remote config/tools service outages). + # Keep PR checks non-blocking and continue to SARIF fallback/upload. + continue-on-error: ${{ github.event_name == 'pull_request' }} + + - name: Warn when Codacy analysis fails on PR + if: ${{ always() && github.event_name == 'pull_request' && steps.codacy_analysis.outcome == 'failure' }} + run: | + echo "::warning::Codacy Analysis CLI failed on this pull_request run; continuing with SARIF fallback." # Process SARIF file to split by tool - name: Split SARIF by tool + if: always() run: | # Fail fast and surface errors clearly set -euo pipefail @@ -124,6 +134,7 @@ jobs: # Select SARIF file for upload - name: Select SARIF file for upload + if: always() run: | set -euo pipefail # Honor preselected SARIF_FILE from earlier steps (e.g., empty SARIF case) From d911fcc51c02807e91372e1db94aaaa196d6ccb9 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Sat, 7 Mar 2026 07:27:30 -0800 Subject: [PATCH 4/5] Changed: harden Codacy CI and expand exact feature tests (internal) Update the Codacy workflow to handle missing project tokens gracefully on pull requests and clarify SARIF validation steps. Add comprehensive unit tests for the internal `det_errbound` and `bareiss_det` functions to ensure logic correctness across various matrix dimensions and edge cases like singular matrices. Refs: #33 --- .github/workflows/codacy.yml | 17 ++++++--- src/exact.rs | 71 ++++++++++++++++++++++++++++++++++++ 2 files changed, 83 insertions(+), 5 deletions(-) diff --git a/.github/workflows/codacy.yml b/.github/workflows/codacy.yml index 3b7135a..d8a157a 100644 --- a/.github/workflows/codacy.yml +++ b/.github/workflows/codacy.yml @@ -41,6 +41,8 @@ jobs: # only required for a private repository by # github/codeql-action/upload-sarif to get the Action run status actions: read + env: + CODACY_PROJECT_TOKEN: ${{ secrets.CODACY_PROJECT_TOKEN }} name: Codacy Security Scan runs-on: ubuntu-latest timeout-minutes: 30 @@ -77,14 +79,13 @@ jobs: # Execute Codacy Analysis CLI and generate a SARIF output with # the security issues identified during the analysis - name: Run Codacy Analysis CLI + if: ${{ env.CODACY_PROJECT_TOKEN != '' }} id: codacy_analysis uses: codacy/codacy-analysis-cli-action@562ee3e92b8e92df8b67e0a5ff8aa8e261919c08 with: # Check https://github.com/codacy/codacy-analysis-cli#project-token # to get your project token from your Codacy repository. - # You can also omit the token and run the tools that support - # default configurations - project-token: ${{ secrets.CODACY_PROJECT_TOKEN }} + project-token: ${{ env.CODACY_PROJECT_TOKEN }} verbose: true directory: ${{ env.CODACY_WORKDIR }} output: ${{ env.CODACY_SARIF }} @@ -99,13 +100,19 @@ jobs: # Keep PR checks non-blocking and continue to SARIF fallback/upload. continue-on-error: ${{ github.event_name == 'pull_request' }} + - name: Warn when Codacy token is unavailable on PR + if: ${{ github.event_name == 'pull_request' && env.CODACY_PROJECT_TOKEN == '' }} + run: | + echo "::warning::CODACY_PROJECT_TOKEN is unavailable for this pull_request." + echo "::warning::Skipping Codacy Analysis CLI and using SARIF fallback." + - name: Warn when Codacy analysis fails on PR if: ${{ always() && github.event_name == 'pull_request' && steps.codacy_analysis.outcome == 'failure' }} run: | echo "::warning::Codacy Analysis CLI failed on this pull_request run; continuing with SARIF fallback." - # Process SARIF file to split by tool - - name: Split SARIF by tool + # Validate SARIF output or create an empty fallback for upload + - name: Validate or create SARIF if: always() run: | # Fail fast and surface errors clearly diff --git a/src/exact.rs b/src/exact.rs index 796ebaf..cbc5fe4 100644 --- a/src/exact.rs +++ b/src/exact.rs @@ -405,4 +405,75 @@ mod tests { // Two transpositions → even permutation → det = +1 assert_eq!(m.det_sign_exact(), 1); } + + // ----------------------------------------------------------------------- + // Direct tests for internal helpers (coverage of private functions) + // ----------------------------------------------------------------------- + + #[test] + fn det_errbound_d0_is_zero() { + assert_eq!(det_errbound(&Matrix::<0>::zero()), Some(0.0)); + } + + #[test] + fn det_errbound_d1_is_zero() { + assert_eq!(det_errbound(&Matrix::<1>::from_rows([[42.0]])), Some(0.0)); + } + + #[test] + fn det_errbound_d2_positive() { + let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); + let bound = det_errbound(&m).unwrap(); + assert!(bound > 0.0); + // bound = ERR_COEFF_2 * (|1*4| + |2*3|) = ERR_COEFF_2 * 10 + assert!(ERR_COEFF_2.mul_add(-10.0, bound).abs() < 1e-30); + } + + #[test] + fn det_errbound_d3_positive() { + let m = Matrix::<3>::identity(); + let bound = det_errbound(&m).unwrap(); + assert!(bound > 0.0); + } + + #[test] + fn det_errbound_d4_positive() { + let m = Matrix::<4>::identity(); + let bound = det_errbound(&m).unwrap(); + assert!(bound > 0.0); + } + + #[test] + fn det_errbound_d5_is_none() { + assert_eq!(det_errbound(&Matrix::<5>::identity()), None); + } + + #[test] + fn bareiss_det_d0_is_one() { + let det = bareiss_det(&Matrix::<0>::zero()); + assert_eq!(det, BigRational::from_integer(BigInt::from(1))); + } + + #[test] + fn bareiss_det_d1_returns_entry() { + let det = bareiss_det(&Matrix::<1>::from_rows([[7.0]])); + assert_eq!(det, f64_to_bigrational(7.0)); + } + + #[test] + fn bareiss_det_d3_with_pivoting() { + // First column has zero on diagonal → exercises pivot swap + break. + let m = Matrix::<3>::from_rows([[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]); + let det = bareiss_det(&m); + // det of this permutation matrix = -1 + assert_eq!(det, BigRational::from_integer(BigInt::from(-1))); + } + + #[test] + fn bareiss_det_singular_all_zeros_in_column() { + // Column 1 is all zeros below diagonal after elimination → singular. + let m = Matrix::<3>::from_rows([[1.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]); + let det = bareiss_det(&m); + assert_eq!(det, BigRational::from_integer(BigInt::from(0))); + } } From 2205f12a94743c49333c6f6a7228edde21a0d74f Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Sat, 7 Mar 2026 07:43:15 -0800 Subject: [PATCH 5/5] Changed: expand det_sign_exact tests for large matrices (internal) Add unit tests for dimensions $D \ge 5$ to exercise the Bareiss algorithm path. Ensure non-finite entries trigger expected panics and verify pivoting logic when the fast filter is bypassed. Refs: #33 --- src/exact.rs | 30 +++++++++++++++++++++++++++--- 1 file changed, 27 insertions(+), 3 deletions(-) diff --git a/src/exact.rs b/src/exact.rs index cbc5fe4..97b53d0 100644 --- a/src/exact.rs +++ b/src/exact.rs @@ -386,9 +386,33 @@ mod tests { } #[test] - fn det_sign_exact_pivot_needed_in_first_column() { - // First column has zero on diagonal but non-zero below → needs pivoting. - let m = Matrix::<3>::from_rows([[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]); + #[should_panic(expected = "non-finite matrix entry")] + fn det_sign_exact_panics_on_nan_5x5() { + // D ≥ 5 bypasses the fast filter, exercising the bareiss_det path. + let mut m = Matrix::<5>::identity(); + m.set(2, 3, f64::NAN); + let _ = m.det_sign_exact(); + } + + #[test] + #[should_panic(expected = "non-finite matrix entry")] + fn det_sign_exact_panics_on_infinity_5x5() { + let mut m = Matrix::<5>::identity(); + m.set(0, 0, f64::INFINITY); + let _ = m.det_sign_exact(); + } + + #[test] + fn det_sign_exact_pivot_needed_5x5() { + // D ≥ 5 skips the fast filter → exercises Bareiss pivoting. + // Permutation matrix with a single swap (rows 0↔1) → det = −1. + let m = Matrix::<5>::from_rows([ + [0.0, 1.0, 0.0, 0.0, 0.0], + [1.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 1.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 1.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 1.0], + ]); assert_eq!(m.det_sign_exact(), -1); }