From d04fcd3b24843b7377a34849dbb2e6469bf9fc56 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Thu, 12 Mar 2026 12:40:14 -0700 Subject: [PATCH 1/4] feat: exact linear system solve (solve_exact, solve_exact_f64) - Add `solve_exact()` returning `[BigRational; D]` via Gaussian elimination with partial pivoting in BigRational - Add `solve_exact_f64()` converting exact result to `Vector` - Add `validate_finite_vec()` helper for vector input validation - Add `exact_solve_3x3` example (near-singular system: f64 vs exact) - Macro-generate solve tests for D=2..5 - Update Cargo.toml: keywords (exact-arithmetic, robust-predicates), description ("Fast" not "Small") - Simplify scalar types section (remove f32/f128 speculation) - Add `gh` CLI usage guidance to AGENTS.md - Update README, AGENTS.md with solve_exact documentation Closes #42 --- AGENTS.md | 28 ++- Cargo.toml | 8 +- README.md | 25 +- examples/exact_solve_3x3.rs | 65 +++++ src/exact.rs | 462 +++++++++++++++++++++++++++++++++++- 5 files changed, 570 insertions(+), 18 deletions(-) create mode 100644 examples/exact_solve_3x3.rs diff --git a/AGENTS.md b/AGENTS.md index 4d0bc15..9bb3f62 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -93,7 +93,7 @@ testable. #### Reference examples - `src/matrix.rs` — `gen_public_api_matrix_tests!` -- `src/exact.rs` — `gen_det_exact_tests!`, `gen_det_exact_f64_tests!` +- `src/exact.rs` — `gen_det_exact_tests!`, `gen_det_exact_f64_tests!`, `gen_solve_exact_tests!`, `gen_solve_exact_f64_tests!` #### When single-dimension tests are acceptable @@ -144,7 +144,8 @@ just examples # Run all examples - 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` / `cargo run --features exact --example exact_det_3x3` / - `cargo run --features exact --example exact_sign_3x3`) + `cargo run --features exact --example exact_sign_3x3` / + `cargo run --features exact --example exact_solve_3x3`) - Spell check: `just spell-check` (uses `typos.toml` at repo root; add false positives to `[default.extend-words]`) ### Changelog @@ -155,6 +156,15 @@ just examples # Run all examples ### GitHub Issues +Use the `gh` CLI to read, create, and edit issues: + +- **Read**: `gh issue view ` (or `--json title,body,labels,milestone` for structured data) +- **List**: `gh issue list` (add `--label enhancement`, `--milestone v0.3.0`, etc. to filter) +- **Create**: `gh issue create --title "..." --body "..." --label enhancement --label rust` +- **Edit**: `gh issue edit --add-label "..."`, `--milestone "..."`, `--title "..."` +- **Comment**: `gh issue comment --body "..."` +- **Close**: `gh issue close ` (with optional `--reason completed` or `--reason "not planned"`) + When creating or updating issues: - **Labels**: Use appropriate labels: `enhancement`, `bug`, `performance`, `documentation`, `rust`, `python`, etc. @@ -173,10 +183,10 @@ When creating or updating issues: ## Feature flags -- `exact` — enables exact determinant methods via `BigRational` arithmetic: - `det_exact()`, `det_exact_f64()`, and `det_sign_exact()`. +- `exact` — enables exact arithmetic methods via `BigRational`: + `det_exact()`, `det_exact_f64()`, `det_sign_exact()`, `solve_exact()`, and `solve_exact_f64()`. Also re-exports `BigRational` from the crate root and prelude. - Gates `src/exact.rs`, additional tests, and the `exact_det_3x3`/`exact_sign_3x3` examples. + Gates `src/exact.rs`, additional tests, and the `exact_det_3x3`/`exact_sign_3x3`/`exact_solve_3x3` examples. Clippy, doc builds, and test commands have dedicated `--features exact` variants. ## Code structure (big picture) @@ -188,9 +198,11 @@ When creating or updating issues: - `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_exact()`, `det_exact_f64()`, `det_sign_exact()` — exact determinant - value and sign via Bareiss in `BigRational`; `det_sign_exact()` adds a Shewchuk-style f64 - filter for fast sign resolution; `features = ["exact"]` + - `src/exact.rs`: exact arithmetic behind `features = ["exact"]`: + - Determinants: `det_exact()`, `det_exact_f64()`, `det_sign_exact()` via Bareiss in + `BigRational`; `det_sign_exact()` adds a Shewchuk-style f64 filter for fast sign resolution + - Linear system solve: `solve_exact()`, `solve_exact_f64()` via Gaussian elimination + with partial pivoting in `BigRational` - Rust tests are inline `#[cfg(test)]` modules in each `src/*.rs` file. - Python tests live in `scripts/tests/` and run via `just test-python` (`uv run pytest`). - The public API re-exports these items from `src/lib.rs`. diff --git a/Cargo.toml b/Cargo.toml index d397461..dd9b22e 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -4,13 +4,13 @@ version = "0.2.2" edition = "2024" rust-version = "1.94" license = "BSD-3-Clause" -description = "Small, stack-allocated linear algebra for fixed dimensions" +description = "Fast, stack-allocated linear algebra for fixed dimensions" readme = "README.md" documentation = "https://docs.rs/la-stack" repository = "https://github.com/acgetchell/la-stack" homepage = "https://github.com/acgetchell/la-stack" categories = [ "mathematics", "science" ] -keywords = [ "linear-algebra", "geometry", "const-generics" ] +keywords = [ "linear-algebra", "geometry", "const-generics", "exact-arithmetic", "robust-predicates" ] [dependencies] # All runtime deps are optional; see [features] below. @@ -39,6 +39,10 @@ required-features = [ "exact" ] name = "exact_sign_3x3" required-features = [ "exact" ] +[[example]] +name = "exact_solve_3x3" +required-features = [ "exact" ] + [[bench]] name = "vs_linalg" harness = false diff --git a/README.md b/README.md index 942f91f..aafadf9 100644 --- a/README.md +++ b/README.md @@ -32,6 +32,7 @@ while keeping the API intentionally small and explicit. - ✅ `const fn` where possible (compile-time evaluation of determinants, dot products, etc.) - ✅ Explicit algorithms (LU, solve, determinant) - ✅ Robust geometric predicates via optional exact arithmetic (`det_sign_exact`, `det_errbound`) +- ✅ Exact linear system solve via optional arbitrary-precision arithmetic (`solve_exact`, `solve_exact_f64`) - ✅ No runtime dependencies by default (optional features may add deps) - ✅ Stack storage only (no heap allocation in core types) - ✅ `unsafe` forbidden @@ -46,9 +47,9 @@ See [CHANGELOG.md](CHANGELOG.md) for details. ## 🔢 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). Arbitrary-precision arithmetic -is available via the optional `"exact"` feature (see below). +The core types use `f64`. When f64 precision is insufficient (e.g. near-degenerate +geometric configurations), the optional `"exact"` feature provides arbitrary-precision +arithmetic via `BigRational` (see below). ## 🚀 Quickstart @@ -136,8 +137,8 @@ for D ≤ 4 and falls back to LU for D ≥ 5 — no API change needed. ## 🔬 Exact arithmetic (`"exact"` feature) The default build has **zero runtime dependencies**. Enable the optional -`exact` Cargo feature to add exact determinant methods using adaptive-precision -arithmetic (this pulls in `num-bigint`, `num-rational`, and `num-traits` for +`exact` Cargo feature to add exact arithmetic methods using arbitrary-precision +rationals (this pulls in `num-bigint`, `num-rational`, and `num-traits` for `BigRational`): ```toml @@ -145,10 +146,17 @@ arithmetic (this pulls in `num-bigint`, `num-rational`, and `num-traits` for la-stack = { version = "0.2.2", features = ["exact"] } ``` +**Determinants:** + - **`det_exact()`** — returns the exact determinant as a `BigRational` - **`det_exact_f64()`** — returns the exact determinant converted to the nearest `f64` - **`det_sign_exact()`** — returns the provably correct sign (−1, 0, or +1) +**Linear system solve:** + +- **`solve_exact(b)`** — solves `Ax = b` exactly, returning `[BigRational; D]` +- **`solve_exact_f64(b)`** — solves `Ax = b` exactly, converting the result to `Vector` (f64) + ```rust,ignore use la_stack::prelude::*; @@ -204,12 +212,15 @@ exposed for advanced use cases. | Type | Storage | Purpose | Key methods | |---|---|---|---| | `Vector` | `[f64; D]` | Fixed-length vector | `new`, `zero`, `dot`, `norm2_sq` | -| `Matrix` | `[[f64; D]; D]` | Square matrix | `lu`, `ldlt`, `det`, `det_direct`, `det_errbound`, `det_exact`¹, `det_exact_f64`¹, `det_sign_exact`¹ | +| `Matrix` | `[[f64; D]; D]` | Square matrix | See below | | `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. +`Matrix` key methods: `lu`, `ldlt`, `det`, `det_direct`, `det_errbound`, +`det_exact`¹, `det_exact_f64`¹, `det_sign_exact`¹, `solve_exact`¹, `solve_exact_f64`¹. + ¹ Requires `features = ["exact"]`. ## 📋 Examples @@ -221,6 +232,7 @@ The `examples/` directory contains small, runnable programs: - **`const_det_4x4`** — compile-time 4×4 determinant via `det_direct()` - **`exact_det_3x3`** — exact determinant value of a near-singular 3×3 matrix (requires `exact` feature) - **`exact_sign_3x3`** — exact determinant sign of a near-singular 3×3 matrix (requires `exact` feature) +- **`exact_solve_3x3`** — exact solve of a near-singular 3×3 system vs f64 LU (requires `exact` feature) ```bash just examples @@ -230,6 +242,7 @@ cargo run --example det_5x5 cargo run --example const_det_4x4 cargo run --features exact --example exact_det_3x3 cargo run --features exact --example exact_sign_3x3 +cargo run --features exact --example exact_solve_3x3 ``` ## 🤝 Contributing diff --git a/examples/exact_solve_3x3.rs b/examples/exact_solve_3x3.rs new file mode 100644 index 0000000..045b059 --- /dev/null +++ b/examples/exact_solve_3x3.rs @@ -0,0 +1,65 @@ +//! Exact linear system solve for a near-singular 3×3 system. +//! +//! This example demonstrates `solve_exact()` and `solve_exact_f64()`, which use +//! arbitrary-precision rational arithmetic to compute a provably correct +//! solution — even when the matrix is so close to singular that the f64 LU +//! solve produces a wildly inaccurate result. +//! +//! Run with: `cargo run --features exact --example exact_solve_3x3` + +use la_stack::prelude::*; + +fn main() { + // Near-singular 3×3 system. + // + // The base matrix [[1,2,3],[4,5,6],[7,8,9]] is exactly singular (rows in + // arithmetic progression). Perturbing entry (0,0) by 2^-50 ≈ 8.9e-16 + // makes it invertible but extremely ill-conditioned. + let perturbation = f64::from_bits(0x3CD0_0000_0000_0000); // 2^-50 + let a = Matrix::<3>::from_rows([ + [1.0 + perturbation, 2.0, 3.0], + [4.0, 5.0, 6.0], + [7.0, 8.0, 9.0], + ]); + + let b = Vector::<3>::new([1.0, 2.0, 3.0]); + + // f64 LU solve (using zero pivot tolerance since the matrix is nearly singular + // and would be rejected by DEFAULT_PIVOT_TOL). + let lu_x = a.lu(0.0).unwrap().solve_vec(b).unwrap().into_array(); + + // Exact solve. + let exact_x = a.solve_exact(b).unwrap(); + let exact_x_f64 = a.solve_exact_f64(b).unwrap().into_array(); + + println!("Near-singular 3×3 system (perturbation = 2^-50 ≈ {perturbation:.2e}):"); + for r in 0..3 { + print!(" ["); + for c in 0..3 { + if c > 0 { + print!(", "); + } + print!("{:22.18}", a.get(r, c).unwrap()); + } + println!("]"); + } + println!( + "b = [{}, {}, {}]", + b.as_array()[0], + b.as_array()[1], + b.as_array()[2] + ); + println!(); + println!( + "f64 LU solve: x = [{:+.6e}, {:+.6e}, {:+.6e}]", + lu_x[0], lu_x[1], lu_x[2] + ); + println!( + "solve_exact(): x = [{}, {}, {}]", + exact_x[0], exact_x[1], exact_x[2] + ); + println!( + "solve_exact_f64(): x = [{:+.6e}, {:+.6e}, {:+.6e}]", + exact_x_f64[0], exact_x_f64[1], exact_x_f64[2] + ); +} diff --git a/src/exact.rs b/src/exact.rs index 8d1351e..631b03f 100644 --- a/src/exact.rs +++ b/src/exact.rs @@ -1,10 +1,12 @@ -//! Exact determinant value and sign via arbitrary-precision rational arithmetic. +//! Exact arithmetic operations via arbitrary-precision rational numbers. //! //! This module is only compiled when the `"exact"` Cargo feature is enabled. //! //! # Architecture //! -//! All three public methods (`det_exact`, `det_exact_f64`, `det_sign_exact`) +//! ## Determinants +//! +//! All three determinant methods (`det_exact`, `det_exact_f64`, `det_sign_exact`) //! share the same core: entries are converted losslessly to `BigRational` and //! the Bareiss algorithm (fraction-free Gaussian elimination) computes the //! exact determinant. @@ -17,6 +19,13 @@ //! immediately without allocating. //! 2. **Exact fallback**: run the Bareiss algorithm for a guaranteed-correct //! sign. +//! +//! ## Linear system solve +//! +//! `solve_exact` and `solve_exact_f64` solve `A x = b` using standard Gaussian +//! elimination with partial pivoting in `BigRational` arithmetic. Every finite +//! `f64` is exactly representable as a rational, so the result is provably +//! correct. use num_bigint::{BigInt, Sign}; use num_rational::BigRational; @@ -24,6 +33,7 @@ use num_traits::ToPrimitive; use crate::LaError; use crate::matrix::Matrix; +use crate::vector::Vector; /// Validate that all entries in a `D×D` matrix are finite (not NaN or infinite). /// @@ -40,6 +50,19 @@ fn validate_finite(m: &Matrix) -> Result<(), LaError> { Ok(()) } +/// Validate that all entries in a length-`D` vector are finite. +/// +/// Returns `Ok(())` if all entries are finite, or `Err(LaError::NonFinite)` with +/// the index of the first non-finite entry found. +fn validate_finite_vec(v: &Vector) -> Result<(), LaError> { + for (i, &x) in v.data.iter().enumerate() { + if !x.is_finite() { + return Err(LaError::NonFinite { col: i }); + } + } + Ok(()) +} + /// Convert an `f64` to an exact `BigRational`. /// /// Every finite `f64` is exactly representable as a rational number (`m × 2^e`), @@ -113,6 +136,78 @@ fn bareiss_det(m: &Matrix) -> BigRational { if sign < 0 { -det } else { det.clone() } } +/// Solve `A x = b` using Gaussian elimination with partial pivoting in +/// `BigRational` arithmetic. +/// +/// Returns the exact solution as a `Vec` of length `D`. +/// Returns `Err(LaError::Singular)` if the matrix is exactly singular. +fn gauss_solve(m: &Matrix, b: &Vector) -> Result, LaError> { + if D == 0 { + return Ok(Vec::new()); + } + + let zero = BigRational::from_integer(BigInt::from(0)); + + // Build augmented matrix [A | b] as D × (D+1). + let mut aug: Vec> = Vec::with_capacity(D); + for r in 0..D { + let mut row = Vec::with_capacity(D + 1); + for c in 0..D { + row.push(f64_to_bigrational(m.rows[r][c])); + } + row.push(f64_to_bigrational(b.data[r])); + aug.push(row); + } + + // Forward elimination with partial pivoting. + for k in 0..D { + // Find first non-zero pivot in column k at or below row k. + if aug[k][k] == zero { + let mut found = false; + for i in (k + 1)..D { + if aug[i][k] != zero { + aug.swap(k, i); + found = true; + break; + } + } + if !found { + return Err(LaError::Singular { pivot_col: k }); + } + } + + // Eliminate below pivot. + let pivot = aug[k][k].clone(); + for i in (k + 1)..D { + if aug[i][k] == zero { + continue; + } + let factor = &aug[i][k] / &pivot; + // We need index `j` to read aug[k][j] and write aug[i][j] + // (two distinct rows) — iterators can't borrow both. + #[allow(clippy::needless_range_loop)] + for j in (k + 1)..=D { + let term = &factor * &aug[k][j]; + aug[i][j] -= term; + } + aug[i][k] = zero.clone(); + } + } + + // Back-substitution. + let mut x: Vec = vec![zero; D]; + for ii in 0..D { + let i = D - 1 - ii; + let mut sum = aug[i][D].clone(); + for j in (i + 1)..D { + sum -= &aug[i][j] * &x[j]; + } + x[i] = sum / &aug[i][i]; + } + + Ok(x) +} + impl Matrix { /// Exact determinant using arbitrary-precision rational arithmetic. /// @@ -176,6 +271,81 @@ impl Matrix { } } + /// Exact linear system solve using arbitrary-precision rational arithmetic. + /// + /// Solves `A x = b` where `A` is `self` and `b` is the given vector. + /// Returns the exact solution as `[BigRational; D]`. Every finite `f64` + /// is exactly representable as a rational, so the conversion is lossless + /// and the result is provably correct. + /// + /// # When to use + /// + /// Use this when you need a provably correct solution — for example, + /// exact circumcenter computation for near-degenerate simplices where + /// f64 arithmetic may produce wildly wrong results. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// // A x = b where A = [[1,2],[3,4]], b = [5, 11] → x = [1, 2] + /// let a = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); + /// let b = Vector::<2>::new([5.0, 11.0]); + /// let x = a.solve_exact(b).unwrap(); + /// assert_eq!(x[0], BigRational::from_integer(1.into())); + /// assert_eq!(x[1], BigRational::from_integer(2.into())); + /// ``` + /// + /// # Errors + /// Returns [`LaError::NonFinite`] if any matrix or vector entry is NaN or + /// infinite. + /// Returns [`LaError::Singular`] if the matrix is exactly singular. + #[inline] + pub fn solve_exact(&self, b: Vector) -> Result<[BigRational; D], LaError> { + validate_finite(self)?; + validate_finite_vec(&b)?; + let solution = gauss_solve(self, &b)?; + Ok(std::array::from_fn(|i| solution[i].clone())) + } + + /// Exact linear system solve converted to `f64`. + /// + /// Computes the exact [`BigRational`] solution via + /// [`solve_exact`](Self::solve_exact) and converts each component to the + /// nearest `f64`. This is useful when you want the most accurate f64 + /// solution possible without committing to `BigRational` downstream. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// let a = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); + /// let b = Vector::<2>::new([5.0, 11.0]); + /// let x = a.solve_exact_f64(b).unwrap().into_array(); + /// assert!((x[0] - 1.0).abs() <= f64::EPSILON); + /// assert!((x[1] - 2.0).abs() <= f64::EPSILON); + /// ``` + /// + /// # Errors + /// Returns [`LaError::NonFinite`] if any matrix or vector entry is NaN or + /// infinite. + /// Returns [`LaError::Singular`] if the matrix is exactly singular. + /// Returns [`LaError::Overflow`] if any component of the exact solution is + /// too large to represent as a finite `f64`. + #[inline] + pub fn solve_exact_f64(&self, b: Vector) -> Result, LaError> { + let exact = self.solve_exact(b)?; + let mut result = [0.0f64; D]; + for (i, val) in exact.iter().enumerate() { + let f = val.to_f64().unwrap_or(f64::INFINITY); + if !f.is_finite() { + return Err(LaError::Overflow); + } + result[i] = f; + } + Ok(Vector::new(result)) + } + /// Exact determinant sign using adaptive-precision arithmetic. /// /// Returns `1` if `det > 0`, `-1` if `det < 0`, and `0` if `det == 0` (singular). @@ -711,6 +881,294 @@ mod tests { assert_eq!(m.det_exact_f64(), Err(LaError::Overflow)); } + // ----------------------------------------------------------------------- + // solve_exact: macro-generated per-dimension tests (D=2..5) + // ----------------------------------------------------------------------- + + /// Helper: build an arbitrary RHS vector for dimension `$d`. + fn arbitrary_rhs() -> Vector { + let values = [1.0, -2.5, 3.0, 0.25, -4.0]; + let mut arr = [0.0f64; D]; + for (dst, src) in arr.iter_mut().zip(values.iter()) { + *dst = *src; + } + Vector::::new(arr) + } + + macro_rules! gen_solve_exact_tests { + ($d:literal) => { + paste! { + #[test] + fn []() { + let a = Matrix::<$d>::identity(); + let b = arbitrary_rhs::<$d>(); + let x = a.solve_exact(b).unwrap(); + for (i, xi) in x.iter().enumerate() { + assert_eq!(*xi, f64_to_bigrational(b.data[i])); + } + } + + #[test] + fn []() { + let mut a = Matrix::<$d>::identity(); + a.set(0, 0, f64::NAN); + let b = arbitrary_rhs::<$d>(); + assert_eq!(a.solve_exact(b), Err(LaError::NonFinite { col: 0 })); + } + + #[test] + fn []() { + let mut a = Matrix::<$d>::identity(); + a.set(0, 0, f64::INFINITY); + let b = arbitrary_rhs::<$d>(); + assert_eq!(a.solve_exact(b), Err(LaError::NonFinite { col: 0 })); + } + + #[test] + fn []() { + let a = Matrix::<$d>::identity(); + let mut b_arr = [1.0f64; $d]; + b_arr[0] = f64::NAN; + let b = Vector::<$d>::new(b_arr); + assert_eq!(a.solve_exact(b), Err(LaError::NonFinite { col: 0 })); + } + + #[test] + fn []() { + let a = Matrix::<$d>::identity(); + let mut b_arr = [1.0f64; $d]; + b_arr[$d - 1] = f64::INFINITY; + let b = Vector::<$d>::new(b_arr); + assert_eq!(a.solve_exact(b), Err(LaError::NonFinite { col: $d - 1 })); + } + + #[test] + fn []() { + // Zero matrix is singular. + let a = Matrix::<$d>::zero(); + let b = arbitrary_rhs::<$d>(); + assert_eq!(a.solve_exact(b), Err(LaError::Singular { pivot_col: 0 })); + } + } + }; + } + + gen_solve_exact_tests!(2); + gen_solve_exact_tests!(3); + gen_solve_exact_tests!(4); + gen_solve_exact_tests!(5); + + macro_rules! gen_solve_exact_f64_tests { + ($d:literal) => { + paste! { + #[test] + fn []() { + let a = Matrix::<$d>::identity(); + let b = arbitrary_rhs::<$d>(); + let x = a.solve_exact_f64(b).unwrap().into_array(); + for i in 0..$d { + assert!((x[i] - b.data[i]).abs() <= f64::EPSILON); + } + } + + #[test] + fn []() { + let mut a = Matrix::<$d>::identity(); + a.set(0, 0, f64::NAN); + let b = arbitrary_rhs::<$d>(); + assert_eq!(a.solve_exact_f64(b), Err(LaError::NonFinite { col: 0 })); + } + } + }; + } + + gen_solve_exact_f64_tests!(2); + gen_solve_exact_f64_tests!(3); + gen_solve_exact_f64_tests!(4); + gen_solve_exact_f64_tests!(5); + + /// For D ≤ 4, `solve_exact_f64` should agree with `Lu::solve_vec` on + /// well-conditioned matrices. + macro_rules! gen_solve_exact_f64_agrees_with_lu { + ($d:literal) => { + paste! { + #[test] + #[allow(clippy::cast_precision_loss)] + fn []() { + // Diagonally dominant → well-conditioned. + let mut rows = [[0.0f64; $d]; $d]; + for r in 0..$d { + for c in 0..$d { + rows[r][c] = if r == c { + (r as f64) + f64::from($d) + 1.0 + } else { + 0.1 / ((r + c + 1) as f64) + }; + } + } + let a = Matrix::<$d>::from_rows(rows); + let b = arbitrary_rhs::<$d>(); + let exact = a.solve_exact_f64(b).unwrap().into_array(); + let lu_sol = a.lu(DEFAULT_PIVOT_TOL).unwrap() + .solve_vec(b).unwrap().into_array(); + for i in 0..$d { + let eps = lu_sol[i].abs().mul_add(1e-12, 1e-12); + assert!((exact[i] - lu_sol[i]).abs() <= eps); + } + } + } + }; + } + + gen_solve_exact_f64_agrees_with_lu!(2); + gen_solve_exact_f64_agrees_with_lu!(3); + gen_solve_exact_f64_agrees_with_lu!(4); + gen_solve_exact_f64_agrees_with_lu!(5); + + // ----------------------------------------------------------------------- + // solve_exact: dimension-specific tests + // ----------------------------------------------------------------------- + + #[test] + fn solve_exact_d0_returns_empty() { + let a = Matrix::<0>::zero(); + let b = Vector::<0>::zero(); + let x = a.solve_exact(b).unwrap(); + assert!(x.is_empty()); + } + + #[test] + fn solve_exact_known_2x2() { + // [[1,2],[3,4]] x = [5, 11] → x = [1, 2] + let a = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); + let b = Vector::<2>::new([5.0, 11.0]); + let x = a.solve_exact(b).unwrap(); + assert_eq!(x[0], BigRational::from_integer(BigInt::from(1))); + assert_eq!(x[1], BigRational::from_integer(BigInt::from(2))); + } + + #[test] + fn solve_exact_pivoting_needed() { + // First column has zero on diagonal → pivot swap required. + let a = Matrix::<3>::from_rows([[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]); + let b = Vector::<3>::new([2.0, 3.0, 4.0]); + let x = a.solve_exact(b).unwrap(); + // x = [3, 2, 4] + assert_eq!(x[0], f64_to_bigrational(3.0)); + assert_eq!(x[1], f64_to_bigrational(2.0)); + assert_eq!(x[2], f64_to_bigrational(4.0)); + } + + #[test] + fn solve_exact_fractional_result() { + // [[2, 1], [1, 3]] x = [1, 1] → x = [2/5, 1/5] + let a = Matrix::<2>::from_rows([[2.0, 1.0], [1.0, 3.0]]); + let b = Vector::<2>::new([1.0, 1.0]); + let x = a.solve_exact(b).unwrap(); + assert_eq!(x[0], BigRational::new(BigInt::from(2), BigInt::from(5))); + assert_eq!(x[1], BigRational::new(BigInt::from(1), BigInt::from(5))); + } + + #[test] + fn solve_exact_singular_duplicate_rows() { + let a = Matrix::<3>::from_rows([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [1.0, 2.0, 3.0]]); + let b = Vector::<3>::new([1.0, 2.0, 3.0]); + assert!(matches!(a.solve_exact(b), Err(LaError::Singular { .. }))); + } + + #[test] + fn solve_exact_5x5_permutation() { + // Permutation matrix (swap rows 0↔1): P x = b → x = P^T b. + let a = 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], + ]); + let b = Vector::<5>::new([10.0, 20.0, 30.0, 40.0, 50.0]); + let x = a.solve_exact(b).unwrap(); + assert_eq!(x[0], f64_to_bigrational(20.0)); + assert_eq!(x[1], f64_to_bigrational(10.0)); + assert_eq!(x[2], f64_to_bigrational(30.0)); + assert_eq!(x[3], f64_to_bigrational(40.0)); + assert_eq!(x[4], f64_to_bigrational(50.0)); + } + + // ----------------------------------------------------------------------- + // solve_exact_f64: dimension-specific tests + // ----------------------------------------------------------------------- + + #[test] + fn solve_exact_f64_known_2x2() { + let a = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); + let b = Vector::<2>::new([5.0, 11.0]); + let x = a.solve_exact_f64(b).unwrap().into_array(); + assert!((x[0] - 1.0).abs() <= f64::EPSILON); + assert!((x[1] - 2.0).abs() <= f64::EPSILON); + } + + #[test] + fn solve_exact_f64_overflow_returns_err() { + // [[1/big, 0], [0, 1/big]] x = [big, big] → x = [big², big²], + // which overflows f64. + let big = f64::MAX / 2.0; + let a = Matrix::<2>::from_rows([[1.0 / big, 0.0], [0.0, 1.0 / big]]); + let b = Vector::<2>::new([big, big]); + assert_eq!(a.solve_exact_f64(b), Err(LaError::Overflow)); + } + + // ----------------------------------------------------------------------- + // gauss_solve: internal helper tests + // ----------------------------------------------------------------------- + + #[test] + fn gauss_solve_d0_returns_empty() { + let a = Matrix::<0>::zero(); + let b = Vector::<0>::zero(); + assert_eq!(gauss_solve(&a, &b).unwrap().len(), 0); + } + + #[test] + fn gauss_solve_d1() { + let a = Matrix::<1>::from_rows([[2.0]]); + let b = Vector::<1>::new([6.0]); + let x = gauss_solve(&a, &b).unwrap(); + assert_eq!(x[0], f64_to_bigrational(3.0)); + } + + #[test] + fn gauss_solve_singular_column_all_zero() { + let a = Matrix::<3>::from_rows([[1.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]); + let b = Vector::<3>::new([1.0, 2.0, 3.0]); + assert_eq!(gauss_solve(&a, &b), Err(LaError::Singular { pivot_col: 1 })); + } + + // ----------------------------------------------------------------------- + // validate_finite_vec tests + // ----------------------------------------------------------------------- + + #[test] + fn validate_finite_vec_ok() { + assert!(validate_finite_vec(&Vector::<3>::new([1.0, 2.0, 3.0])).is_ok()); + } + + #[test] + fn validate_finite_vec_err_on_nan() { + assert_eq!( + validate_finite_vec(&Vector::<2>::new([f64::NAN, 1.0])), + Err(LaError::NonFinite { col: 0 }) + ); + } + + #[test] + fn validate_finite_vec_err_on_inf() { + assert_eq!( + validate_finite_vec(&Vector::<2>::new([1.0, f64::NEG_INFINITY])), + Err(LaError::NonFinite { col: 1 }) + ); + } + // ----------------------------------------------------------------------- // validate_finite tests // ----------------------------------------------------------------------- From 125b259608a2c1965aba486d2594ff9fa328609c Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Thu, 12 Mar 2026 12:57:54 -0700 Subject: [PATCH 2/4] refactor: restructure gauss_solve to avoid break/continue Replace `found` boolean + `break` with idiomatic `.find()` + `if let`, and `if zero { continue }` with `if not_zero { ... }` block. Eliminates coverage instrumentation gaps on these keywords. --- src/exact.rs | 33 +++++++++++++-------------------- 1 file changed, 13 insertions(+), 20 deletions(-) diff --git a/src/exact.rs b/src/exact.rs index 631b03f..4410fad 100644 --- a/src/exact.rs +++ b/src/exact.rs @@ -163,15 +163,9 @@ fn gauss_solve(m: &Matrix, b: &Vector) -> Result(m: &Matrix, b: &Vector) -> Result Date: Thu, 12 Mar 2026 13:32:14 -0700 Subject: [PATCH 3/4] docs: clarify pivoting strategy and use idiomatic rev() in exact.rs - Replace "partial pivoting" with "first-non-zero pivoting" in module doc, function doc, and inline comments for both bareiss_det and gauss_solve - Add rationale: exact BigRational arithmetic means any non-zero pivot gives a correct result (no numerical stability concern) - Use idiomatic (0..D).rev() for back-substitution loop --- src/exact.rs | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/src/exact.rs b/src/exact.rs index 4410fad..9f0c28d 100644 --- a/src/exact.rs +++ b/src/exact.rs @@ -22,10 +22,11 @@ //! //! ## Linear system solve //! -//! `solve_exact` and `solve_exact_f64` solve `A x = b` using standard Gaussian -//! elimination with partial pivoting in `BigRational` arithmetic. Every finite -//! `f64` is exactly representable as a rational, so the result is provably -//! correct. +//! `solve_exact` and `solve_exact_f64` solve `A x = b` using Gaussian +//! elimination with first-non-zero pivoting in `BigRational` arithmetic. +//! Since all arithmetic is exact, any non-zero pivot gives the correct result +//! (there is no numerical stability concern). Every finite `f64` is exactly +//! representable as a rational, so the result is provably correct. use num_bigint::{BigInt, Sign}; use num_rational::BigRational; @@ -101,7 +102,7 @@ fn bareiss_det(m: &Matrix) -> BigRational { let mut sign: i8 = 1; for k in 0..D { - // Partial pivoting: find non-zero entry in column k at or below row k. + // Find first 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 { @@ -136,8 +137,12 @@ fn bareiss_det(m: &Matrix) -> BigRational { if sign < 0 { -det } else { det.clone() } } -/// Solve `A x = b` using Gaussian elimination with partial pivoting in -/// `BigRational` arithmetic. +/// Solve `A x = b` using Gaussian elimination with first-non-zero pivoting +/// in `BigRational` arithmetic. +/// +/// Since all arithmetic is exact, any non-zero pivot gives the correct result +/// (no numerical stability concern). This matches the pivoting strategy used +/// by `bareiss_det`. /// /// Returns the exact solution as a `Vec` of length `D`. /// Returns `Err(LaError::Singular)` if the matrix is exactly singular. @@ -159,7 +164,7 @@ fn gauss_solve(m: &Matrix, b: &Vector) -> Result(m: &Matrix, b: &Vector) -> Result = vec![zero; D]; - for ii in 0..D { - let i = D - 1 - ii; + for i in (0..D).rev() { let mut sum = aug[i][D].clone(); for j in (i + 1)..D { sum -= &aug[i][j] * &x[j]; From 2bddccf3f20156fadb4fee6a4093e76b1b631de1 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Thu, 12 Mar 2026 14:18:45 -0700 Subject: [PATCH 4/4] docs: add LDLT example, solve_exact README snippet, and sync examples MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Add examples/ldlt_solve_3x3.rs: 3×3 SPD tridiagonal LDLT solve + det - Add solve_exact_f64 code example to README exact arithmetic section - Add missing examples to justfile (ldlt_solve_3x3, exact_det_3x3, exact_solve_3x3) - Add lu.rs and ldlt.rs test macro references to AGENTS.md - Fix stale "partial pivoting" → "first-non-zero pivoting" in AGENTS.md code structure - Drop "current" from README storage note --- AGENTS.md | 7 ++++-- README.md | 12 ++++++++++- examples/ldlt_solve_3x3.rs | 44 ++++++++++++++++++++++++++++++++++++++ justfile | 3 +++ 4 files changed, 63 insertions(+), 3 deletions(-) create mode 100644 examples/ldlt_solve_3x3.rs diff --git a/AGENTS.md b/AGENTS.md index 9bb3f62..67d1947 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -93,6 +93,8 @@ testable. #### Reference examples - `src/matrix.rs` — `gen_public_api_matrix_tests!` +- `src/lu.rs` — `gen_public_api_pivoting_solve_vec_and_det_tests!`, `gen_public_api_tridiagonal_smoke_solve_vec_and_det_tests!` +- `src/ldlt.rs` — `gen_public_api_ldlt_identity_tests!`, `gen_public_api_ldlt_diagonal_tests!` - `src/exact.rs` — `gen_det_exact_tests!`, `gen_det_exact_f64_tests!`, `gen_solve_exact_tests!`, `gen_solve_exact_f64_tests!` #### When single-dimension tests are acceptable @@ -143,7 +145,8 @@ just examples # Run all examples - 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` / `cargo run --features exact --example exact_det_3x3` / + `cargo run --example ldlt_solve_3x3` / `cargo run --example const_det_4x4` / + `cargo run --features exact --example exact_det_3x3` / `cargo run --features exact --example exact_sign_3x3` / `cargo run --features exact --example exact_solve_3x3`) - Spell check: `just spell-check` (uses `typos.toml` at repo root; add false positives to `[default.extend-words]`) @@ -202,7 +205,7 @@ When creating or updating issues: - Determinants: `det_exact()`, `det_exact_f64()`, `det_sign_exact()` via Bareiss in `BigRational`; `det_sign_exact()` adds a Shewchuk-style f64 filter for fast sign resolution - Linear system solve: `solve_exact()`, `solve_exact_f64()` via Gaussian elimination - with partial pivoting in `BigRational` + with first-non-zero pivoting in `BigRational` - Rust tests are inline `#[cfg(test)]` modules in each `src/*.rs` file. - Python tests live in `scripts/tests/` and run via `just test-python` (`uv run pytest`). - The public API re-exports these items from `src/lib.rs`. diff --git a/README.md b/README.md index aafadf9..767292d 100644 --- a/README.md +++ b/README.md @@ -160,6 +160,7 @@ la-stack = { version = "0.2.2", features = ["exact"] } ```rust,ignore use la_stack::prelude::*; +// Exact determinant let m = Matrix::<3>::from_rows([ [1.0, 2.0, 3.0], [4.0, 5.0, 6.0], @@ -169,6 +170,13 @@ assert_eq!(m.det_sign_exact().unwrap(), 0); // exactly singular let det = m.det_exact().unwrap(); assert_eq!(det, BigRational::from_integer(0.into())); // exact zero + +// Exact linear system solve +let a = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); +let b = Vector::<2>::new([5.0, 11.0]); +let x = a.solve_exact_f64(b).unwrap().into_array(); +assert!((x[0] - 1.0).abs() <= f64::EPSILON); +assert!((x[1] - 2.0).abs() <= f64::EPSILON); ``` `BigRational` is re-exported from the crate root and prelude when the `exact` @@ -216,7 +224,7 @@ exposed for advanced use cases. | `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. +Storage shown above reflects the `f64` implementation. `Matrix` key methods: `lu`, `ldlt`, `det`, `det_direct`, `det_errbound`, `det_exact`¹, `det_exact_f64`¹, `det_sign_exact`¹, `solve_exact`¹, `solve_exact_f64`¹. @@ -229,6 +237,7 @@ 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 +- **`ldlt_solve_3x3`** — solve a 3×3 symmetric positive definite system via LDLT - **`const_det_4x4`** — compile-time 4×4 determinant via `det_direct()` - **`exact_det_3x3`** — exact determinant value of a near-singular 3×3 matrix (requires `exact` feature) - **`exact_sign_3x3`** — exact determinant sign of a near-singular 3×3 matrix (requires `exact` feature) @@ -239,6 +248,7 @@ just examples # or individually: cargo run --example solve_5x5 cargo run --example det_5x5 +cargo run --example ldlt_solve_3x3 cargo run --example const_det_4x4 cargo run --features exact --example exact_det_3x3 cargo run --features exact --example exact_sign_3x3 diff --git a/examples/ldlt_solve_3x3.rs b/examples/ldlt_solve_3x3.rs new file mode 100644 index 0000000..6fc5073 --- /dev/null +++ b/examples/ldlt_solve_3x3.rs @@ -0,0 +1,44 @@ +//! Solve a 3×3 symmetric positive definite system via LDLT factorization. +//! +//! LDLT is the natural choice for SPD matrices (e.g. Gram matrices, covariance +//! matrices, stiffness matrices). It avoids the row swaps of LU and exploits +//! symmetry, using roughly half the work. +//! +//! Run with: `cargo run --example ldlt_solve_3x3` + +use la_stack::prelude::*; + +fn main() -> Result<(), LaError> { + // Symmetric positive definite 3×3 matrix (classic SPD tridiagonal). + let a = Matrix::<3>::from_rows([[4.0, -1.0, 0.0], [-1.0, 4.0, -1.0], [0.0, -1.0, 4.0]]); + + // Choose x = [1, 2, 3]. Then b = A x = [2, 4, 10]. + let b = Vector::<3>::new([2.0, 4.0, 10.0]); + + let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL)?; + let x = ldlt.solve_vec(b)?.into_array(); + let det = ldlt.det(); + + println!("A (3×3 SPD tridiagonal):"); + for r in 0..3 { + print!(" ["); + for c in 0..3 { + if c > 0 { + print!(", "); + } + print!("{:5.1}", a.get(r, c).unwrap()); + } + println!("]"); + } + println!( + "b = [{}, {}, {}]", + b.as_array()[0], + b.as_array()[1], + b.as_array()[2] + ); + println!(); + println!("x = [{:.6}, {:.6}, {:.6}]", x[0], x[1], x[2]); + println!("det = {det}"); + + Ok(()) +} diff --git a/justfile b/justfile index 01428e6..7fb8c78 100644 --- a/justfile +++ b/justfile @@ -211,8 +211,11 @@ doc-check: examples: cargo run --quiet --example det_5x5 cargo run --quiet --example solve_5x5 + cargo run --quiet --example ldlt_solve_3x3 cargo run --quiet --example const_det_4x4 + cargo run --quiet --features exact --example exact_det_3x3 cargo run --quiet --features exact --example exact_sign_3x3 + cargo run --quiet --features exact --example exact_solve_3x3 # Fix (mutating): apply formatters/auto-fixes fix: toml-fmt fmt python-fix shell-fmt markdown-fix yaml-fix