Merge pull request #48 from acgetchell/feat/solve-exact

feat: exact linear system solve (solve_exact, solve_exact_f64)

Author: acgetchell

AGENTS.md

@@ -93,7 +93,9 @@ testable.
 #### Reference examples
 
 - `src/matrix.rs` — `gen_public_api_matrix_tests!`
-- `src/exact.rs` — `gen_det_exact_tests!`, `gen_det_exact_f64_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,8 +145,10 @@ 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 --features exact --example exact_sign_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]`)
 
 ### Changelog
@@ -155,6 +159,15 @@ just examples         # Run all examples
 
 ### GitHub Issues
 
+Use the `gh` CLI to read, create, and edit issues:
+
+- **Read**: `gh issue view <number>` (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 <number> --add-label "..."`, `--milestone "..."`, `--title "..."`
+- **Comment**: `gh issue comment <number> --body "..."`
+- **Close**: `gh issue close <number>` (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 +186,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 +201,11 @@ When creating or updating issues:
   - `src/matrix.rs`: `Matrix<const D: usize>` (`[[f64; D]; D]`) + helpers (`get`, `set`, `inf_norm`, `det`, `det_direct`)
   - `src/lu.rs`: `Lu<const D: usize>` factorization with partial pivoting (`solve_vec`, `det`)
   - `src/ldlt.rs`: `Ldlt<const D: usize>` 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 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`.

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

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,22 +137,30 @@ 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
 [dependencies]
 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<D>` (f64)
+
 ```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],
@@ -161,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`
@@ -204,11 +220,14 @@ exposed for advanced use cases.
 | Type | Storage | Purpose | Key methods |
 |---|---|---|---|
 | `Vector<D>` | `[f64; D]` | Fixed-length vector | `new`, `zero`, `dot`, `norm2_sq` |
-| `Matrix<D>` | `[[f64; D]; D]` | Square matrix | `lu`, `ldlt`, `det`, `det_direct`, `det_errbound`, `det_exact`¹, `det_exact_f64`¹, `det_sign_exact`¹ |
+| `Matrix<D>` | `[[f64; D]; D]` | Square matrix | See below |
 | `Lu<D>` | `Matrix<D>` + pivot array | Factorization for solves/det | `solve_vec`, `det` |
 | `Ldlt<D>` | `Matrix<D>` | 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<D>` 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"]`.
 
@@ -218,18 +237,22 @@ 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)
+- **`exact_solve_3x3`** — exact solve of a near-singular 3×3 system vs f64 LU (requires `exact` feature)
 
 ```bash
 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
+cargo run --features exact --example exact_solve_3x3
 ```
 
 ## 🤝 Contributing

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]
+    );
+}

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(())
+}

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