feat!(api): enforce fallible numeric invariants

- Introduce validated tolerance values for factorization and symmetry APIs.
- Return typed errors for non-finite matrix, vector, determinant, and norm results.
- Make determinant and vector operations propagate overflow/non-finite failures.
- Update examples, docs, benchmarks, and property tests for the fallible API.
- Document the v0.4.2 roadmap order for finite Matrix/Vector proof-type work.

BREAKING CHANGE: tolerance arguments now use Tolerance instead of raw f64, Matrix::set returns Option<()>, determinant helpers return Result, Lu::det and Ldlt::det return Result<f64, LaError>, and Vector::dot / Vector::norm2_sq return Result<f64, LaError>.

Closes #83

Author: acgetchell

AGENTS.md

@@ -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
 
@@ -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` /

CHANGELOG.md

@@ -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)
@@ -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
 

Cargo.toml

@@ -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"

README.md

@@ -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(())
 }
 ```
 
@@ -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
@@ -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
@@ -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
@@ -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(())
 }
 ```
 

benches/exact.rs

@@ -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;
 
@@ -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);
                 });

benches/vs_linalg.rs

@@ -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;
 
@@ -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;
@@ -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));
@@ -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();
@@ -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);
                     });
                 });
@@ -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);
                     });
@@ -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);
                     });
@@ -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))
@@ -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);
                     });
                 });
@@ -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);
                     });
                 });
@@ -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);
                     });
                 });
@@ -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);
                     });
                 });