feat: exact linear system solve (solve_exact, solve_exact_f64)

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 <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 +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<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 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`.

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,19 +137,26 @@ 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::*;
 
@@ -204,12 +212,15 @@ 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.
 
+`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"]`.
 
 ## 📋 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

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