Skip to content

feat!(api): enforce fallible numeric invariants #127

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
acgetchell merged 1 commit into from
Jun 3, 2026
Merged

feat!(api): enforce fallible numeric invariants #127

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions AGENTS.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,6 +111,9 @@ invariant over the convenient edit.
`git --no-pager log`, `git --no-pager show`, `git --no-pager blame`) to inspect changes/history
- **ALWAYS** use `git --no-pager` when reading git output
- Suggest git commands that modify version control state for the user to run manually
- When suggesting branch names, prefer `{type}/{issue}-descriptor-or-two`, e.g. `fix/307-topology-validation`,
`perf/315-bench-profile`, or `doc/329-branch-guidance`. If an environment requires an owner/tool prefix,
keep this structure after the prefix, e.g. `codex/fix/307-topology-validation`.

### Commit Messages

Expand Down Expand Up @@ -240,6 +243,9 @@ just examples # Run all examples
- Python setup: `uv sync --group dev` (or `just python-sync`)
- 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`)
Cargo accepts only one positional test filter. To run multiple focused
filters, run separate `cargo test <filter>` commands rather than passing
multiple filter arguments.
- 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 ldlt_solve_3x3` / `cargo run --example const_det_4x4` /
Expand Down
9 changes: 9 additions & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0

## [Unreleased]

### Added

- Feat!(matrix): enforce fallible matrix invariants [`e26c283`](https://github.com/acgetchell/la-stack/commit/e26c28358b2358100353b2895441b68892e92cd7)

### Changed

- Remove redundant cache restore keys for cargo-llvm-cov [`f75a01c`](https://github.com/acgetchell/la-stack/commit/f75a01c99c8dbcc8b6ffc36ae9f94ba968a2f111)
Expand Down Expand Up @@ -37,6 +41,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
### Documentation

- Document feature requirement for exact APIs [`19b10d5`](https://github.com/acgetchell/la-stack/commit/19b10d552e83b6a7f9e91695b4850b8fab3f4550)
- Document scalar scope and release roadmap [`bfb0393`](https://github.com/acgetchell/la-stack/commit/bfb039386588f94b95561c610181ca6d486acd6e)

- Clarify that la-stack intentionally supports f64 floating-point APIs plus optional exact rationals, not alternate scalar families.
- Add a roadmap covering the v0.4.x stable-Rust issue sequence and the v0.5.0 generic_const_exprs anchor.
- Refresh generated changelog entries and archived changelog grouping.

### Maintenance

Expand Down
4 changes: 2 additions & 2 deletions Cargo.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,8 @@ proptest = "1.11.0"

[features]
default = [ ]
bench = [ "criterion", "faer", "nalgebra" ]
exact = [ "num-bigint", "num-rational", "num-traits" ]
bench = [ "dep:criterion", "dep:faer", "dep:nalgebra" ]
exact = [ "dep:num-bigint", "dep:num-rational", "dep:num-traits" ]

[[example]]
name = "exact_det_3x3"
Expand Down
139 changes: 78 additions & 61 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,25 +76,29 @@ Solve a 5×5 system via LU:
```rust
use la_stack::prelude::*;

// This system requires pivoting (a[0][0] = 0), so it's a good LU demo.
// A = J - I: zeros on diagonal, ones elsewhere.
let a = Matrix::<5>::from_rows([
[0.0, 1.0, 1.0, 1.0, 1.0],
[1.0, 0.0, 1.0, 1.0, 1.0],
[1.0, 1.0, 0.0, 1.0, 1.0],
[1.0, 1.0, 1.0, 0.0, 1.0],
[1.0, 1.0, 1.0, 1.0, 0.0],
]);

let b = Vector::<5>::new([14.0, 13.0, 12.0, 11.0, 10.0]);

let lu = a.lu(DEFAULT_PIVOT_TOL).unwrap();
let x = lu.solve_vec(b).unwrap().into_array();

// Floating-point rounding is expected; compare with a tolerance.
let expected = [1.0, 2.0, 3.0, 4.0, 5.0];
for (x_i, e_i) in x.iter().zip(expected.iter()) {
assert!((*x_i - *e_i).abs() <= 1e-12);
fn main() -> Result<(), LaError> {
// This system requires pivoting (a[0][0] = 0), so it's a good LU demo.
// A = J - I: zeros on diagonal, ones elsewhere.
let a = Matrix::<5>::from_rows([
[0.0, 1.0, 1.0, 1.0, 1.0],
[1.0, 0.0, 1.0, 1.0, 1.0],
[1.0, 1.0, 0.0, 1.0, 1.0],
[1.0, 1.0, 1.0, 0.0, 1.0],
[1.0, 1.0, 1.0, 1.0, 0.0],
]);

let b = Vector::<5>::new([14.0, 13.0, 12.0, 11.0, 10.0]);

let lu = a.lu(DEFAULT_PIVOT_TOL)?;
let x = lu.solve_vec(b)?.into_array();

// Floating-point rounding is expected; compare with a tolerance.
let expected = [1.0, 2.0, 3.0, 4.0, 5.0];
for (x_i, e_i) in x.iter().zip(expected.iter()) {
assert!((*x_i - *e_i).abs() <= 1e-12);
}

Ok(())
}
```

Expand All @@ -106,17 +110,21 @@ For symmetric positive-definite matrices, `LDL^T` is essentially a square-root-f
```rust
use la_stack::prelude::*;

// This matrix is symmetric positive-definite (A = L*L^T) so LDLT works without pivoting.
let a = Matrix::<5>::from_rows([
[1.0, 1.0, 0.0, 0.0, 0.0],
[1.0, 2.0, 1.0, 0.0, 0.0],
[0.0, 1.0, 2.0, 1.0, 0.0],
[0.0, 0.0, 1.0, 2.0, 1.0],
[0.0, 0.0, 0.0, 1.0, 2.0],
]);

let det = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap().det();
assert!((det - 1.0).abs() <= 1e-12);
fn main() -> Result<(), LaError> {
// This matrix is symmetric positive-definite (A = L*L^T) so LDLT works without pivoting.
let a = Matrix::<5>::from_rows([
[1.0, 1.0, 0.0, 0.0, 0.0],
[1.0, 2.0, 1.0, 0.0, 0.0],
[0.0, 1.0, 2.0, 1.0, 0.0],
[0.0, 0.0, 1.0, 2.0, 1.0],
[0.0, 0.0, 0.0, 1.0, 2.0],
]);

let det = a.ldlt(DEFAULT_SINGULAR_TOL)?.det()?;
assert!((det - 1.0).abs() <= 1e-12);

Ok(())
}
```

> ⚠️ **LDLT invariant:** The input matrix must be **symmetric**. Asymmetric
Expand All @@ -133,23 +141,23 @@ assert!((det - 1.0).abs() <= 1e-12);

`det_direct()` is a `const fn` providing closed-form determinants for D=0–4,
using fused multiply-add where applicable. `Matrix::<0>::zero().det_direct()`
returns `Some(1.0)` (the empty-product convention). For D=1–4, cofactor
returns `Ok(Some(1.0))` (the empty-product convention). For D=1–4, cofactor
expansion bypasses LU factorization entirely. This enables compile-time
evaluation when inputs are known:

```rust
use la_stack::prelude::*;

// Evaluated entirely at compile time — no runtime cost.
const DET: Option<f64> = {
const DET: Result<Option<f64>, LaError> = {
let m = Matrix::<3>::from_rows([
[2.0, 0.0, 0.0],
[0.0, 3.0, 0.0],
[0.0, 0.0, 5.0],
]);
m.det_direct()
};
assert_eq!(DET, Some(30.0));
assert_eq!(DET, Ok(Some(30.0)));
```

The public `det()` method automatically dispatches through the closed-form path
Expand Down Expand Up @@ -181,23 +189,27 @@ la-stack = { version = "0.4.1", 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],
[7.0, 8.0, 9.0],
]);
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);
fn main() -> Result<(), LaError> {
// Exact determinant
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

let det = m.det_exact()?;
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)?.into_array();
assert!((x[0] - 1.0).abs() <= f64::EPSILON);
assert!((x[1] - 2.0).abs() <= f64::EPSILON);

Ok(())
}
```

With the `exact` feature enabled, `BigInt` and `BigRational` are re-exported
Expand All @@ -222,19 +234,24 @@ adaptive-precision logic for geometric predicates:
```rust,ignore
use la_stack::prelude::*;

let m = Matrix::<3>::identity();
if let Some(bound) = m.det_errbound() {
let det = m.det_direct().unwrap();
if det.abs() > bound {
// f64 sign is guaranteed correct
let sign = det.signum() as i8;
fn main() -> Result<(), LaError> {
let m = Matrix::<3>::identity();
if let Some(bound) = m.det_errbound()? {
if let Some(det) = m.det_direct()? {
if det.abs() > bound {
// f64 sign is guaranteed correct
let sign = det.signum() as i8;
} else {
// Fall back to exact arithmetic (requires `exact` feature)
let sign = m.det_sign_exact()?;
}
}
} else {
// Fall back to exact arithmetic (requires `exact` feature)
let sign = m.det_sign_exact().unwrap();
// D ≥ 5: no fast filter, use exact directly (requires `exact` feature)
let sign = m.det_sign_exact()?;
}
} else {
// D ≥ 5: no fast filter, use exact directly (requires `exact` feature)
let sign = m.det_sign_exact().unwrap();

Ok(())
}
```

Expand Down
4 changes: 2 additions & 2 deletions benches/exact.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@
//! empirical evidence for `docs/PERFORMANCE.md`.

use criterion::{BenchmarkGroup, Criterion, measurement::WallTime};
use la_stack::{Matrix, Vector};
use la_stack::{DEFAULT_PIVOT_TOL, Matrix, Vector};
use pastey::paste;
use std::hint::black_box;

Expand Down Expand Up @@ -179,7 +179,7 @@ macro_rules! gen_exact_benches_for_dim {
[<group_d $d>].bench_function("det", |bencher| {
bencher.iter(|| {
let det = black_box(a)
.det(la_stack::DEFAULT_PIVOT_TOL)
.det(DEFAULT_PIVOT_TOL)
.expect("diagonally dominant matrix is non-singular");
black_box(det);
});
Expand Down
39 changes: 23 additions & 16 deletions benches/vs_linalg.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,9 @@
//! - Matrix infinity norm is the maximum absolute row sum on all sides.

use criterion::Criterion;
use faer::linalg::solvers::Solve;
use faer::linalg::solvers::{PartialPivLu, Solve};
use faer::perm::PermRef;
use la_stack::{DEFAULT_PIVOT_TOL, Matrix, Vector};
use pastey::paste;
use std::hint::black_box;

Expand Down Expand Up @@ -43,7 +44,7 @@ fn faer_perm_sign(p: PermRef<'_, usize>) -> f64 {
}
}

fn faer_det_from_partial_piv_lu(lu: &faer::linalg::solvers::PartialPivLu<f64>) -> f64 {
fn faer_det_from_partial_piv_lu(lu: &PartialPivLu<f64>) -> f64 {
// For PA = LU with unit-lower L, det(A) = det(P) * det(U).
let u = lu.U();
let mut det = 1.0;
Expand Down Expand Up @@ -128,10 +129,10 @@ macro_rules! gen_vs_linalg_benches_for_dim {
paste! {{
// Isolate each dimension's inputs to keep types and captures clean.
{
let a = la_stack::Matrix::<$d>::from_rows(make_matrix_rows::<$d>());
let rhs = la_stack::Vector::<$d>::new(make_vector_array::<$d>(0.0));
let v1 = la_stack::Vector::<$d>::new(make_vector_array::<$d>(0.0));
let v2 = la_stack::Vector::<$d>::new(make_vector_array::<$d>(1.0));
let a = Matrix::<$d>::from_rows(make_matrix_rows::<$d>());
let rhs = Vector::<$d>::new(make_vector_array::<$d>(0.0));
let v1 = Vector::<$d>::new(make_vector_array::<$d>(0.0));
let v2 = Vector::<$d>::new(make_vector_array::<$d>(1.0));

let na = nalgebra::SMatrix::<f64, $d, $d>::from_fn(|r, c| matrix_entry::<$d>(r, c));
let nrhs = nalgebra::SVector::<f64, $d>::from_fn(|i, _| vector_entry(i, 0.0));
Expand All @@ -145,7 +146,7 @@ macro_rules! gen_vs_linalg_benches_for_dim {

// Precompute LU once for solve-only / det-only benchmarks.
let a_lu = a
.lu(la_stack::DEFAULT_PIVOT_TOL)
.lu(DEFAULT_PIVOT_TOL)
.expect("matrix should be non-singular");
let na_lu = na.clone().lu();
let fa_lu = fa.partial_piv_lu();
Expand All @@ -156,9 +157,12 @@ macro_rules! gen_vs_linalg_benches_for_dim {
[<group_d $d>].bench_function("la_stack_det_via_lu", |bencher| {
bencher.iter(|| {
let lu = black_box(a)
.lu(la_stack::DEFAULT_PIVOT_TOL)
.lu(DEFAULT_PIVOT_TOL)
.expect("matrix should be non-singular");
let det = lu.det();
let det = match lu.det() {
Ok(det) => det,
Err(err) => panic!("finite benchmark matrix determinant failed: {err}"),
};
black_box(det);
});
});
Expand All @@ -183,7 +187,7 @@ macro_rules! gen_vs_linalg_benches_for_dim {
[<group_d $d>].bench_function("la_stack_det", |bencher| {
bencher.iter(|| {
let det = black_box(a)
.det(la_stack::DEFAULT_PIVOT_TOL)
.det(DEFAULT_PIVOT_TOL)
.expect("matrix should be non-singular");
black_box(det);
});
Expand All @@ -193,7 +197,7 @@ macro_rules! gen_vs_linalg_benches_for_dim {
[<group_d $d>].bench_function("la_stack_lu", |bencher| {
bencher.iter(|| {
let lu = black_box(a)
.lu(la_stack::DEFAULT_PIVOT_TOL)
.lu(DEFAULT_PIVOT_TOL)
.expect("matrix should be non-singular");
let _ = black_box(lu);
});
Expand All @@ -217,7 +221,7 @@ macro_rules! gen_vs_linalg_benches_for_dim {
[<group_d $d>].bench_function("la_stack_lu_solve", |bencher| {
bencher.iter(|| {
let lu = black_box(a)
.lu(la_stack::DEFAULT_PIVOT_TOL)
.lu(DEFAULT_PIVOT_TOL)
.expect("matrix should be non-singular");
let x = lu
.solve_vec(black_box(rhs))
Expand Down Expand Up @@ -273,7 +277,10 @@ macro_rules! gen_vs_linalg_benches_for_dim {
// === Determinant from a precomputed LU ===
[<group_d $d>].bench_function("la_stack_det_from_lu", |bencher| {
bencher.iter(|| {
let det = a_lu.det();
let det = match a_lu.det() {
Ok(det) => det,
Err(err) => panic!("finite benchmark matrix determinant failed: {err}"),
};
black_box(det);
});
});
Expand All @@ -295,7 +302,7 @@ macro_rules! gen_vs_linalg_benches_for_dim {
// === Vector dot product ===
[<group_d $d>].bench_function("la_stack_dot", |bencher| {
bencher.iter(|| {
let result = black_box(v1).dot(black_box(v2));
let result = black_box(v1).dot(black_box(v2)).unwrap();
black_box(result);
});
});
Expand All @@ -322,7 +329,7 @@ macro_rules! gen_vs_linalg_benches_for_dim {
// === Vector norm squared ===
[<group_d $d>].bench_function("la_stack_norm2_sq", |bencher| {
bencher.iter(|| {
let result = black_box(v1).norm2_sq();
let result = black_box(v1).norm2_sq().unwrap();
black_box(result);
});
});
Expand All @@ -349,7 +356,7 @@ macro_rules! gen_vs_linalg_benches_for_dim {
// === Matrix infinity norm (max absolute row sum) ===
[<group_d $d>].bench_function("la_stack_inf_norm", |bencher| {
bencher.iter(|| {
let result = black_box(a).inf_norm();
let result = black_box(a).inf_norm().unwrap();
black_box(result);
});
});
Expand Down
Loading
Loading