diff --git a/.github/workflows/benchmarks.yml b/.github/workflows/benchmarks.yml index 357e380..1858803 100644 --- a/.github/workflows/benchmarks.yml +++ b/.github/workflows/benchmarks.yml @@ -58,7 +58,9 @@ env: jobs: bench: runs-on: ubuntu-latest - timeout-minutes: 30 + # The full exact suite has 260+ cases: Criterion's default 3 s warm-up + # plus 5 s measurement already needs 35+ minutes, before build/analysis. + timeout-minutes: 90 steps: - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 @@ -109,8 +111,7 @@ jobs: repo: context.repo.repo, workflow_id: 'benchmarks.yml', branch: 'main', - status: 'completed', - conclusion: 'success', + status: 'success', per_page: 5, }); @@ -158,15 +159,16 @@ jobs: echo "::notice::Baseline found — comparing against main" echo "comparison_available=true" >> "$GITHUB_OUTPUT" # --baseline-lenient rather than --baseline: benches added on the - # PR branch that don't yet exist in the main baseline get a - # "no baseline data" notice instead of aborting the whole run. + # PR branch that don't yet exist in the main baseline run without + # a comparison instead of aborting the whole run. + # Keep normal sampling; plots/HTML are unused by the summary. cargo bench --locked --features bench,exact --bench exact \ - -- --baseline-lenient main 2>&1 | tee bench-output.txt + -- --noplot --baseline-lenient main 2>&1 | tee bench-output.txt else echo "::notice::No baseline found — running without comparison" echo "comparison_available=false" >> "$GITHUB_OUTPUT" cargo bench --locked --features bench,exact --bench exact \ - 2>&1 | tee bench-output.txt + -- --noplot 2>&1 | tee bench-output.txt fi if grep -q "Performance has regressed" bench-output.txt; then @@ -182,11 +184,11 @@ jobs: github.ref == 'refs/heads/main' run: > cargo bench --locked --features bench,exact --bench exact - -- --save-baseline main + -- --noplot --save-baseline main - name: Run benchmarks (manual ref) if: github.event_name == 'workflow_dispatch' && github.ref != 'refs/heads/main' - run: cargo bench --locked --features bench,exact --bench exact + run: cargo bench --locked --features bench,exact --bench exact -- --noplot - name: Upload baseline artifact if: > @@ -208,12 +210,24 @@ jobs: comparison_available="${BENCH_COMPARISON_AVAILABLE:-}" regression="${BENCH_REGRESSION:-}" + # Lenient comparisons silently skip missing baselines. Require a + # change report for every analyzed case before claiming coverage. + if [ -f bench-output.txt ]; then + analyzed="$(grep -c '^Benchmarking .*: Analyzing$' bench-output.txt || true)" + compared="$(grep -c '^[[:space:]]*change:' bench-output.txt || true)" + if [ "$analyzed" = "0" ] || [ "$analyzed" != "$compared" ]; then + comparison_available=false + fi + else + comparison_available=false + fi + if [ "$comparison_available" != "true" ] || [ -z "$regression" ]; then { echo "### ❓ Benchmark Comparison Unavailable" echo "" echo "No usable comparison against the main baseline was produced." - echo "The benchmark still ran, but no regression claim can be made." + echo "Benchmarks may be incomplete or lack a baseline; no regression claim can be made." } >> "$GITHUB_STEP_SUMMARY" echo "::warning::Benchmark comparison unavailable" elif [ "$regression" = "true" ]; then diff --git a/AGENTS.md b/AGENTS.md index e1a6b66..36d6414 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -230,12 +230,32 @@ When user requests commit message generation: - `src/lib.rs` includes `README.md` with `#![doc = include_str!("../README.md")]`, so README examples are the docs.rs landing page examples. +- Keep the README quickstart and brief capability descriptions discoverable; + put fuller worked API examples and caller contracts in the documentation-only + `guide` module in `src/lib.rs`. Preserve useful detail in the linked guide when + shortening README sections, including numerical limitations and error semantics. +- Link API references and worked API guides to docs.rs. Keep repository-owned + mathematical background, benchmark reports, roadmap, contributing, and release + instructions on GitHub, using absolute URLs to the intended repository revision. +- README links and Contents anchors must work both on GitHub and in generated + rustdoc. Repository-relative file links can resolve incorrectly from rustdoc; + verify destinations in both contexts and use intra-doc links within Rust docs + where the referenced item is available under the selected features. +- Choose docs.rs `latest` links for intentionally current guidance and explicit + versions for release-specific contracts. docs.rs builds published crates, so + merging changes does not publish new guide pages. Local rendering checks do + not establish availability on the published site. - When changing Rust examples in `README.md`, mirror executable versions in the private `readme_doctests` module in `src/lib.rs`. Keep mirrors hidden/private so they do not duplicate the docs.rs landing page, but make them runnable by `cargo test --doc`. - README examples that require optional features may remain `rust,ignore` in README for default-feature doctest compatibility, but must have a `#[cfg(feature = "...")]` hidden doctest mirror in `src/lib.rs` and be verified with the matching feature set (for example, `cargo test --features exact --doc`). +- Guide examples run directly as doctests; gate feature-dependent guides with + the matching feature and remove obsolete private mirrors when moving examples + out of README. Run `just doc-check` for changed guide docs and inspect generated + pages and anchors for explicit links, which rustdoc does not validate. Validate + changed executable examples with the default and matching feature doctest recipes. - When intentionally updating package versions or dependency snippets, keep README `la-stack` dependency examples in sync with the package `version` in `Cargo.toml`. Do not perform version bumps unless explicitly requested by the maintainer; see **Public-API stability** above. diff --git a/README.md b/README.md index 652e64d..1eb1399 100644 --- a/README.md +++ b/README.md @@ -17,6 +17,34 @@ Fast, stack-allocated linear algebra for fixed dimensions in Rust. This crate grew from the need to support [`delaunay`](https://crates.io/crates/delaunay) with fast, stack-allocated linear algebra primitives and algorithms while keeping the API intentionally small and explicit. +## Contents + +- [Introduction](#-introduction) +- [Quickstart](#-quickstart) +- [Mathematical basis](#-mathematical-basis) +- [Design goals](#-design-goals) +- [Anti-goals](#-anti-goals) +- [Use this crate when](#-use-this-crate-when) +- [Scalar and bounded-value types](#-scalar-and-bounded-value-types) +- [Features](#-features) + - [Adaptive determinant filtering (D ≤ 4)](#adaptive-determinant-filtering-d--4) + - [Certified dot products and affine differences](#certified-dot-products-and-affine-differences) + - [Compile-time determinants (D ≤ 4)](#compile-time-determinants-d--4) + - [Exact arithmetic](#exact-arithmetic-exact-feature) + - [LDLT determinant](#ldlt-determinant) + - [LU solve](#lu-solve) + - [Outward-rounded interval determinants](#outward-rounded-interval-determinants) + - [Overflow-safe Euclidean norms](#overflow-safe-euclidean-norms) +- [API at a glance](#-api-at-a-glance) +- [Documentation Map](#documentation-map) +- [Benchmarks](#-benchmarks-vs-nalgebrafaer) +- [Examples](#-examples) +- [Contributing](#-contributing) +- [Citation](#-citation) +- [References](#-references) +- [AI Agents](#-ai-agents) +- [License](#-license) + ## 📐 Introduction `la-stack` provides a handful of const-generic, stack-backed building blocks: @@ -42,6 +70,56 @@ while keeping the API intentionally small and explicit. - `Ldlt` for no-pivot factorization intended for exactly symmetric positive-definite matrices (solve + det; typed pivot diagnostics) +## 🚀 Quickstart + +The minimum supported Rust version (MSRV) is 1.98.1. + +Add this to your `Cargo.toml`: + +```toml +[dependencies] +la-stack = "0.4.5" +``` + +### Solve a 5×5 system + +This system has solution `[1, 2, 3, 4, 5]` and requires partial pivoting: + +```rust +use la_stack::prelude::*; + +fn main() -> Result<(), LaError> { + // The zero leading entry requires LU pivoting. + let a = Matrix::<5>::try_from_rows([ + [0.0, 2.0, -1.0, 1.0, 3.0], + [4.0, -1.0, 2.0, 0.0, 1.0], + [1.0, 3.0, 5.0, -2.0, 0.0], + [2.0, 0.0, -1.0, 4.0, 1.0], + [-1.0, 2.0, 0.0, 1.0, 6.0], + ])?; + let b = Vector::try_new([20.0, 13.0, 14.0, 20.0, 37.0])?; + let lu = a.lu(DEFAULT_SINGULAR_TOL)?; + let x = lu.solve(b)?; + + for (&actual, expected) in x.as_array().iter().zip([1.0, 2.0, 3.0, 4.0, 5.0]) { + assert!((actual - expected).abs() <= 1e-12); + } + Ok(()) +} +``` + +The assertion tolerance is suitable for this known example; LU does not +provide a certified solution error bound. + +### Feature flags + +- `default`: no runtime dependencies; includes outward-rounded `Interval` and + `IntervalMatrix` APIs +- `exact`: exact determinant signs, determinant values, and solves over stored + `f64` values or caller-supplied `BigRational` inputs +- `bench`: repository-development gate used only by benchmark targets and + benchmark-input tests; application crates should not enable it + ## 🧮 Mathematical basis `la-stack` operates on finite IEEE 754 binary64 values in small, fixed @@ -70,7 +148,7 @@ See the [mathematical basis](https://github.com/acgetchell/la-stack/blob/v0.4.5/docs/mathematical_basis.md) for the algorithms, validity boundaries, and supporting references. -## ✨ Design goals +## 🎯 Design goals - ✅ `const fn` where possible (compile-time evaluation of determinants, dot products, etc.) - ✅ Const-generic storage (no dynamically sized matrix or vector representation) @@ -150,447 +228,76 @@ Lower-precision `f32` / `f16` throughput-oriented workloads are outside the crate's scope; they usually indicate large-matrix or accelerator-oriented use cases better served by broader linear-algebra libraries. -## 🚀 Quickstart - -The minimum supported Rust version (MSRV) is 1.98.1. - -Add this to your `Cargo.toml`: +## ✨ Features -```toml -[dependencies] -la-stack = "0.4.5" -``` +### Adaptive determinant filtering (D ≤ 4) -### Feature flags +`det_direct_with_errbound()` pairs a determinant with its certified absolute +bound, without optional dependencies. Resolve the sign when `|det| > bound`; +otherwise an exact fallback is needed. With `exact`, `det_sign_exact()` handles +filtering and fallback automatically. +[Worked examples: the floating-point filter and exact fallback][guide-adaptive]. -- `default`: no runtime dependencies; includes outward-rounded `Interval` and - `IntervalMatrix` APIs -- `exact`: exact determinant signs, determinant values, and solves over stored - `f64` values or caller-supplied `BigRational` inputs -- `bench`: repository-development gate used only by benchmark targets and - benchmark-input tests; application crates should not enable it +### Certified dot products and affine differences -### LU solve +`dot_with_errbound()` and `dot_difference_with_errbound()` return certified +bounds for dot products and `axis · (left - right)` over the original stored +coordinates. Their endpoints support sign and threshold proofs; a bound that +straddles the threshold or an unavailable certificate is inconclusive. +[Worked examples: dot-product signs and affine threshold tests][guide-certified]. -Solve a 5×5 system via LU: +### Compile-time determinants (D ≤ 4) -```rust -use la_stack::prelude::*; +`det_direct()` evaluates closed-form determinants in `const` contexts through +D=4. `det()` selects those formulas automatically and uses zero-tolerance LU +for larger dimensions; a failed numerical pivot remains `LaError::Singular`. +[Compile-time example and dimension contracts][guide-compile-time]. -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>::try_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], - ])?; +### Exact arithmetic (`"exact"` feature) - let b = Vector::<5>::try_new([14.0, 13.0, 12.0, 11.0, 10.0])?; - - let lu = a.lu(DEFAULT_SINGULAR_TOL)?; - let x = lu.solve(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(()) -} -``` - -### LDLT determinant - -Compute a determinant for a symmetric positive-definite matrix via LDLT (no -pivoting). - -For these matrices, `LDLᵀ` is a square-root-free Cholesky form. Multiplying each -column of `L` by the square root of the corresponding diagonal entry yields a -Cholesky factor: - -```rust -use la_stack::prelude::*; - -fn main() -> Result<(), LaError> { - // This matrix is symmetric positive-definite (A = L*L^T) so LDLT works without pivoting. - let a = Matrix::<5>::try_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 ldlt = match a.ldlt(DEFAULT_SINGULAR_TOL) { - Ok(ldlt) => ldlt, - Err(err @ LaError::Asymmetric { - row, - col, - upper, - lower, - allowed_abs_diff, - .. - }) => { - eprintln!( - "LDLT mismatch at ({row}, {col}): {upper} vs {lower} (allowed {allowed_abs_diff})" - ); - return Err(err); - } - Err(err) => return Err(err), - }; - - let det = ldlt.det()?; - assert!((det - 1.0).abs() <= 1e-12); - - Ok(()) -} -``` - -> ⚠️ **LDLT invariant:** The input matrix must be **exactly symmetric**: every -> mirrored pair must compare equal (`+0.0 == -0.0` is accepted). Asymmetric -> inputs passed to -> [`Matrix::ldlt`](https://docs.rs/la-stack/latest/la_stack/struct.Matrix.html#method.ldlt) -> return a typed `LaError::Asymmetric` containing both observed values and the -> required allowed difference of zero. The tolerance-based -> [`Matrix::first_asymmetry`](https://docs.rs/la-stack/latest/la_stack/struct.Matrix.html#method.first_asymmetry) -> and `Matrix::is_symmetric` methods remain useful diagnostics, but do not prove -> the exact precondition required by LDLT. Use `lu()` when exact symmetry or -> positive definiteness is not guaranteed. A negative LDLT diagonal or a zero -> diagonal with nonzero remaining coupling returns -> `LaError::NotPositiveSemidefinite` with a typed -> `PositiveSemidefiniteViolation`. An uncoupled zero or positive pivot -> at or below the caller's tolerance returns `LaError::Singular` with a -> numerical `SingularityReason`. Because these pivots are computed in binary64, -> success is not an exact positive-definiteness certificate for the stored -> matrix. - -## ⚡ Compile-time determinants (D ≤ 4) - -`det_direct()` is a `const fn` providing closed-form determinants for D=0–4, -using fused multiply-add where applicable. It returns `Ok(Some(det))` for those -dimensions and `Ok(None)` for D ≥ 5. `Matrix::<0>::zero().det_direct()` returns -`Ok(Some(1.0))` (the empty-product convention). For D=1–4, direct formulas -bypass 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: Result, LaError> = match Matrix::<4>::try_from_rows([ - [2.0, 0.0, 0.0, 0.0], - [0.0, 3.0, 0.0, 0.0], - [0.0, 0.0, 5.0, 0.0], - [0.0, 0.0, 0.0, 7.0], -]) { - Ok(matrix) => matrix.det_direct(), - Err(err) => Err(err), -}; - -fn main() -> Result<(), LaError> { - assert_eq!(DET?, Some(210.0)); - Ok(()) -} -``` - -The public `det()` method automatically dispatches through the closed-form path -for D ≤ 4 and falls back to zero-tolerance LU for D ≥ 5. Tiny nonzero -determinants are not flattened by a configured pivot tolerance. The LU fallback -returns `LaError::Singular` when floating-point elimination cannot produce a -non-zero pivot; it does not misreport that numerical failure as an exact zero. -Use `lu()` directly when you need a different tolerance policy, and use the -exact determinant APIs when exact singularity classification matters. - -## 📦 Outward-rounded interval determinants - -`Interval` encloses expression construction that has not yet been reduced to a -single stored `f64`. Point intervals preserve finite binary64 values exactly; -`try_from_subtraction`, `try_add`, `try_mul`, `negate`, and `try_square` enclose -the corresponding exact-real operations. `IntervalMatrix::det_sign()` then -uses a division-free subset expansion through D=7, returning positive, -negative, zero, or inconclusive evidence. - -```rust -use la_stack::prelude::*; - -fn main() -> Result<(), LaError> { - // Relative coordinates and the lifted norm retain their construction error. - let x = Interval::try_from_subtraction(0.1, 0.0)?; - let y = Interval::try_from_subtraction(0.1, 0.0)?; - let z = Interval::try_from_subtraction(0.1, 0.0)?; - let lifted = x - .try_square()? - .try_add(&y.try_square()?)? - .try_add(&z.try_square()?)?; - - let matrix = IntervalMatrix::<4>::from_rows([ - [Interval::ONE, Interval::ZERO, Interval::ZERO, Interval::ONE], - [Interval::ZERO, Interval::ONE, Interval::ZERO, Interval::ONE], - [Interval::ZERO, Interval::ZERO, Interval::ONE, Interval::ONE], - [x, y, z, lifted], - ]); - assert_eq!( - matrix.det_sign()?, - IntervalDeterminantSign::Negative, - ); - Ok(()) -} -``` - -Every successful interval keeps finite ordered endpoints. Subnormal bounds are -preserved, both signed zeros are treated as real zero and canonicalized to -`+0.0`, and underflowed nonzero products widen toward the least subnormal value. -If an exact result range cannot fit between finite binary64 endpoints, the -operation returns `LaError::IntervalRangeExhausted` with its interval operation -recorded in `ArithmeticOperation`. - -`Positive`, `Negative`, and `Zero` are proofs. `Inconclusive` only means that -the determinant enclosure overlaps zero; it must not be converted to equality -or singularity. A filtered-exact caller should rebuild the same derived -expression with `RationalMatrix` and call `det_sign()` when the interval result -is inconclusive or reports range failure. Lifting a finished `Matrix` with -`IntervalMatrix::from_matrix` encloses its stored entries, but cannot recover -rounding that occurred while those entries were assembled. - -## 🔬 Exact arithmetic (`"exact"` feature) - -The default build has **zero runtime dependencies**. Enable the optional -`exact` Cargo feature to add exact arithmetic methods using arbitrary-precision -rationals (this pulls in `num-bigint`, `num-rational`, and `num-traits` for -`BigRational`): +Enable exact determinant signs, determinant values, and solves: ```toml [dependencies] la-stack = { version = "0.4.5", features = ["exact"] } ``` -The feature exposes two deliberate input domains: - -- `Matrix` / `Vector` store finite binary64 inputs. Their exact methods - treat each stored bit pattern as its exact rational value, so the determinant - or solve stage introduces no further roundoff. They cannot recover information - already lost before construction. -- `RationalMatrix` / `RationalVector` accept coefficients already - assembled as `BigRational`. They preserve derived differences, squared norms, - affine coefficients, and other rational expressions without an intermediate - `f64` conversion. - -**Determinants:** - -- **`det_exact()`** — returns the exact determinant as a `BigRational` -- **`det_exact_f64()`** — returns the exact determinant as `f64` only when - it is exactly representable (or `LaError::Unrepresentable` otherwise) -- **`det_exact_rounded_f64()`** — returns the exact determinant rounded to a - finite `f64` using IEEE 754 round-to-nearest, ties-to-even -- **`det_sign_exact()`** — infallibly returns the provably correct - `DeterminantSign` variant (`Negative`, `Zero`, or `Positive`) - -**Linear system solve:** - -- **`solve_exact(b)`** — solves `Ax = b` exactly, returning a - `RationalVector` -- **`solve_exact_f64(b)`** — solves `Ax = b` exactly, returning `Vector` only when - every component is exactly representable as `f64` -- **`solve_exact_rounded_f64(b)`** — solves `Ax = b` exactly, returning each - component rounded to finite `f64` using IEEE 754 round-to-nearest, - ties-to-even -- **`ExactF64Conversion`** — converts an existing exact determinant or solution - under the strict or rounded contract without repeating exact elimination - -**Already-exact rational input:** - -- **`RationalMatrix::det_sign()`** — returns the exact sign without constructing - a rational determinant -- **`RationalMatrix::det()`** — returns the exact `BigRational` determinant -- **`RationalMatrix::solve(&rhs)`** — returns a `RationalVector` exact - solution -- **`try_with_rational_matrix!`** — dispatches a runtime-selected dimension - through D=8 to a const-generic rational matrix on stable Rust - -The `Matrix::det_exact*` value and conversion methods return -`LaError::DeterminantScaleOverflow` if their aggregate power-of-two scaling -exceeds the internal exponent representation. `RationalMatrix::det()` is -infallible because it clears rational row denominators without an exponent-scale -conversion. The exact solve methods for both input domains return -`LaError::Singular` with `SingularityReason::Exact` when the stored matrix is -exactly singular. - -For exact-to-f64 output, strict conversions use -`UnrepresentableReason::RequiresRounding` when explicit rounding can produce a -finite value and `UnrepresentableReason::NotFinite` otherwise. Rounded -conversions opt into nearest-even rounding but still report `NotFinite` when no -finite `f64` exists. - -The following 5×5 system has exact determinant 2^-60. Its exact rational inputs -therefore produce a unique solution through the general Bareiss path. Supplying -the same coefficients as `f64` inputs loses the `2^-60` perturbation at `1.0`, -making the leading rows identical and the binary64 system singular. - -```rust,ignore -use core::assert_matches; +`Matrix` / `Vector` exact methods preserve stored `f64` values; +`RationalMatrix` / `RationalVector` also preserve rational expressions before +any `f64` rounding. Keep exact results or explicitly choose strict versus +rounded conversion with `ExactF64Conversion`. +[Worked examples: rational inputs, exact solves, and output conversion][api-exact]. -use la_stack::prelude::*; +### LDLT determinant -fn main() -> Result<(), LaError> { - // This is far below one binary64 ULP at 1.0, so 1.0 + 2^-60 rounds to 1.0. - let epsilon = BigRational::new(1.into(), (1_u64 << 60).into()); - let one = BigRational::from_integer(1.into()); - let zero = BigRational::from_integer(0.into()); - - // The leading block is [[1, 1], [1, 1 + 2^-60]]. The remaining diagonal - // extends the example to D=5, where the general Bareiss path is used. - let matrix = RationalMatrix::<5>::try_from_fn(|row, col| match (row, col) { - (0, 0 | 1) | (1, 0) => one.clone(), - (1, 1) => &one + &epsilon, - _ if row == col => one.clone(), - _ => zero.clone(), - })?; - assert_eq!(matrix.det_sign(), DeterminantSign::Positive); - assert_eq!(matrix.det(), epsilon); - - let rhs = RationalVector::try_new([ - zero, - -&epsilon, - BigRational::from_integer(2.into()), - BigRational::from_integer(3.into()), - BigRational::from_integer(4.into()), - ])?; - let exact_solution = matrix.solve(&rhs)?; - assert_eq!( - exact_solution.as_array(), - &[ - BigRational::from_integer(1.into()), - BigRational::from_integer((-1).into()), - BigRational::from_integer(2.into()), - BigRational::from_integer(3.into()), - BigRational::from_integer(4.into()), - ] - ); - - // Supplying the same coefficients as f64 inputs destroys the perturbation - // and makes the matrix singular, even though the exact solution is integral. - let epsilon_f64 = epsilon.try_to_f64()?; - assert_eq!((1.0 + epsilon_f64).to_bits(), 1.0_f64.to_bits()); - let f64_matrix = Matrix::<5>::try_from_rows([ - [1.0, 1.0, 0.0, 0.0, 0.0], - [1.0, 1.0 + epsilon_f64, 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 f64_rhs = Vector::<5>::try_new([0.0, -epsilon_f64, 2.0, 3.0, 4.0])?; - let f64_solve = f64_matrix - .lu(DEFAULT_SINGULAR_TOL) - .and_then(|lu| lu.solve(f64_rhs)); - assert_matches!( - f64_solve, - Err(LaError::Singular { .. }) - ); - Ok(()) -} -``` +`Matrix::ldlt()` provides a square-root-free factorization for exactly symmetric +positive-definite matrices, supporting determinants and solves without pivoting. +Approximate symmetry is not sufficient, and floating-point success is not an +exact positive-definiteness certificate. +[Worked example and typed pivot diagnostics][guide-ldlt]. -```rust,ignore -use la_stack::prelude::*; +### LU solve -fn main() -> Result<(), LaError> { - // Exact determinant - let m = Matrix::<3>::try_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(), DeterminantSign::Zero); // exactly singular - - let det = m.det_exact()?; - assert_eq!(det, BigRational::from_integer(0.into())); // exact zero - let det_f64 = det.try_to_f64()?; - assert_eq!(det_f64, 0.0); - - // If strict exact-to-f64 conversion would require rounding, opt in - // explicitly with the rounded API. - let inexact = Matrix::<2>::try_from_rows([ - [1.0 + f64::EPSILON, 0.0], - [0.0, 1.0 - f64::EPSILON], - ])?; - let exact_det = inexact.det_exact()?; - let rounded_det = match exact_det.try_to_f64() { - Ok(det) => det, - Err(err) if err.requires_rounding() => exact_det.to_rounded_f64()?, - Err(err) => return Err(err), - }; - assert_eq!(rounded_det.to_bits(), 1.0f64.to_bits()); - - // If the exact determinant cannot fit in f64, keep the BigRational value. - let big = f64::MAX / 2.0; - let huge = Matrix::<3>::try_from_rows([ - [0.0, 0.0, 1.0], - [big, 0.0, 1.0], - [0.0, big, 1.0], - ])?; - let huge_det = huge.det_exact()?; - assert_eq!( - huge_det - .try_to_f64() - .err() - .and_then(|err| err.unrepresentable_reason()), - Some(UnrepresentableReason::NotFinite) - ); - println!("exact determinant = {huge_det}"); - - // Exact linear system solve - let a = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?; - let b = Vector::<2>::try_new([5.0, 11.0])?; - let exact_x = a.solve_exact(b)?; - let x = exact_x.try_to_f64()?.into_array(); - assert!((x[0] - 1.0).abs() <= f64::EPSILON); - assert!((x[1] - 2.0).abs() <= f64::EPSILON); +`Matrix::lu()` uses partial pivoting for general square systems. Reuse one +`Lu` factorization for multiple right-hand sides or a determinant; pivot +tolerances control rejection, not solution accuracy. +[Worked example: solving and reusing factors][guide-lu]. - Ok(()) -} -``` +### Outward-rounded interval determinants + +`Interval` preserves bounds while assembling differences, squares, and other +expressions. `IntervalMatrix::det_sign()` certifies determinant signs through +D=7: an enclosure separated from zero proves its sign, and `[0, 0]` proves +exact zero. Other overlaps with zero are inconclusive and may need exact fallback. +[Worked example: lifted coordinates, range errors, and fallback][guide-intervals]. + +### Overflow-safe Euclidean norms -With the `exact` feature enabled, `RationalMatrix`, `RationalVector`, -`DeterminantSign`, `ExactF64Conversion`, `BigInt`, and `BigRational` are -re-exported from the crate root and prelude, -alongside the most commonly needed `num-traits` items (`FromPrimitive`, -`ToPrimitive`, `Signed`). This lets consumers construct exact values -(`BigRational::from_f64`, `from_i64`), query sign (`is_positive` / -`is_negative`), and convert back to `f64` (`try_to_f64`, `to_rounded_f64`, or -the raw `to_f64`) with a single -`use la_stack::prelude::*;` — no need to add `num-bigint`, `num-rational`, -or `num-traits` to their own `Cargo.toml`. Use -`DeterminantSign::as_i8()` only when numeric −1/0/+1 interoperability is -required. - -For `det_sign_exact()`, D ≤ 4 matrices first use a fast f64 filter -(error-bounded `det_direct()`) when its rounded intermediates stay in the normal -range or are exact structural zeros. An inconclusive filter falls back to the -same direct determinant expansion in `BigInt`. D ≥ 5 skips the closed-form -filter and uses fraction-free Bareiss elimination in `BigInt`. -Because `Matrix` stores only finite entries, arithmetic range failures in the -filter are inconclusive rather than errors and the exact fallback is total. - -## 📏 Overflow-safe Euclidean norms - -`Vector::norm()` computes the Euclidean norm with a deterministic scaled -sum-of-squares recurrence, so large or subnormal finite coordinates do not fail -merely because their raw squares overflow or underflow. It returns positive zero -for empty and all-zero vectors and reports `LaError::NonFinite` with -`ArithmeticOperation::VectorNorm` only when the exact norm rounds to infinity. -Near the upper range, a fixed-size stack accumulator sums squares exactly and -compares squared rounding midpoints to prevent false or hidden overflow. This -fallback needs no optional dependencies. The general binary64 result remains -approximate and has no certified error bound. - -`Vector::norm_squared()` remains the direct left-to-right FMA sum of squares for -callers that need the squared norm. Its distinct contract deliberately reports -overflow when that square is not finite, even when `norm()` can return a finite -norm. +`Vector::norm()` avoids unnecessary overflow and underflow from squaring +coordinates. `norm_squared()` computes the squared norm and can overflow even +when the norm is finite. Both results remain approximate, without a certified +error bound. +[Worked example and range contracts][guide-norms]. **v0.4.6 migration:** `Vector::norm2_sq()` is renamed to `Vector::norm_squared()`, the unreleased `Vector::norm2()` API is named `Vector::norm()`, and @@ -598,217 +305,68 @@ the unreleased `Vector::norm2()` API is named `Vector::norm()`, and are removed; their numerical behavior and error contracts are unchanged by the renames. `Matrix::norm_inf()` remains the maximum absolute row sum. -## 🎯 Certified dot products and affine differences - -`Vector::dot_with_errbound()` evaluates the same left-to-right FMA tree as -`Vector::dot()` and returns its estimate together with a certified absolute -roundoff bound. `Vector::dot_difference_with_errbound()` directly evaluates - -```text -Σᵢ axis[i] × (left[i] - right[i]) -``` - -as two FMAs per coordinate. It does not first round `left - right` into a new -`Vector`, so the certificate covers the intended expression over the original -stored binary64 coordinates. - -The opaque `ScalarWithErrorBound` exposes the estimate, absolute error bound, -and finite outward-rounded lower and upper bounds. Those endpoints support -positive, negative, and caller-selected threshold proofs: - -```rust -use la_stack::prelude::*; - -fn is_separated( - axis: &Vector, - left: &Vector, - right: &Vector, - threshold: f64, -) -> Result, LaError> { - let Some(value) = axis.dot_difference_with_errbound(left, right)? else { - return Ok(None); - }; - if value.lower_bound() > threshold { - Ok(Some(true)) - } else if value.upper_bound() <= threshold { - Ok(Some(false)) - } else { - Ok(None) - } -} - -# fn main() -> Result<(), LaError> { -let axis = Vector::<2>::try_new([2.0, -1.0])?; -let left = Vector::<2>::try_new([4.0, 1.0])?; -let right = Vector::<2>::try_new([1.0, 3.0])?; -assert_eq!(is_separated(&axis, &left, &right, 1.0)?, Some(true)); -# Ok(()) -# } -``` - -An interval that overlaps the threshold is inconclusive, not equal. Likewise, -`Ok(None)` means gradual underflow or proof-only range exhaustion prevented a -certificate. A filtered-exact caller should rebuild the same dot or affine -expression in `BigRational` (available through the `exact` feature) or another -exact backend. A `LaError::NonFinite` instead reports that the specified FMA -estimate itself overflowed. These certified bounds describe roundoff in a fixed -arithmetic tree; they are distinct from user-selected numerical tolerances. - -## 🛡️ Adaptive determinant filtering (D ≤ 4) - -`det_direct_with_errbound()` returns a closed-form determinant together with -the conservative absolute error bound used by the fast filter, computed from -one call that evaluates the determinant once and computes its matching bound. -It returns `None` when a D ≤ 4 computation may be affected by gradual -underflow, as well as for unsupported D ≥ 5 dimensions. -It returns `LaError::NonFinite` if the determinant or bound computation -overflows to NaN or infinity. -This method does NOT require the `exact` feature — it uses pure f64 arithmetic -and is available by default. Use `det_errbound()` when only the bound is needed. -The paired API enables custom adaptive-precision logic for geometric predicates: - -```rust,ignore -use la_stack::prelude::*; - -fn adaptive_det_sign( - matrix: &Matrix, -) -> DeterminantSign { - if let Ok(Some(estimate)) = matrix.det_direct_with_errbound() { - if estimate.determinant().abs() > estimate.absolute_error_bound() { - return if estimate.determinant() > 0.0 { - DeterminantSign::Positive - } else { - DeterminantSign::Negative - }; - } - } - - matrix.det_sign_exact() -} - -fn main() -> Result<(), LaError> { - let identity = Matrix::<3>::identity(); - assert_eq!( - adaptive_det_sign(&identity), - DeterminantSign::Positive - ); - - // A zero determinant cannot pass the f64 sign filter, so this exercises - // the exact fallback. - let singular = Matrix::<3>::try_from_rows([ - [1.0, 2.0, 3.0], - [4.0, 5.0, 6.0], - [7.0, 8.0, 9.0], - ])?; - assert_eq!(adaptive_det_sign(&singular), DeterminantSign::Zero); - - // The f64 filter overflows for this finite matrix, but the exact fallback - // still resolves its positive determinant sign. - let big = f64::MAX / 2.0; - let overflowing = Matrix::<3>::try_from_rows([ - [0.0, 0.0, 1.0], - [big, 0.0, 1.0], - [0.0, big, 1.0], - ])?; - assert_eq!( - adaptive_det_sign(&overflowing), - DeterminantSign::Positive - ); - - Ok(()) -} -``` - -The error coefficients (`ERR_COEFF_2`, `ERR_COEFF_3`, `ERR_COEFF_4`) are -conservative, dimension-specific constants, not caller-tunable tolerances. The -[mathematical basis](https://github.com/acgetchell/la-stack/blob/v0.4.5/docs/mathematical_basis.md#determinants-and-certified-sign-filtering) -documents the bound and states its range preconditions. The constants are explicit -crate-root exports for advanced users who want to compose the same bound: -`use la_stack::{ERR_COEFF_2, ERR_COEFF_3, ERR_COEFF_4};`. They intentionally stay -out of the common prelude. +[guide-lu]: https://docs.rs/la-stack/latest/la_stack/guide/index.html#solving-and-reusing-factors +[guide-ldlt]: https://docs.rs/la-stack/latest/la_stack/guide/ldlt/index.html +[guide-compile-time]: https://docs.rs/la-stack/latest/la_stack/guide/compile_time/index.html +[guide-intervals]: https://docs.rs/la-stack/latest/la_stack/guide/intervals/index.html +[guide-norms]: https://docs.rs/la-stack/latest/la_stack/guide/norms/index.html +[guide-certified]: https://docs.rs/la-stack/latest/la_stack/guide/certified/index.html +[guide-adaptive]: https://docs.rs/la-stack/latest/la_stack/guide/adaptive/index.html ## 🧩 API at a glance -| Type | Storage | Purpose | Key methods | -|---|---|---|---| -| `Vector` | `[f64; D]` | Finite fixed-length vector for input and computation | `try_new`, `as_array`, `into_array`, `dot`, `dot_with_errbound`, `dot_difference_with_errbound`, `norm`, `norm_squared` | -| `Matrix` | `[[f64; D]; D]` | Finite square matrix for input and computation | See below | -| `Interval` | Two finite ordered `f64` bounds | Outward-rounded exact-real enclosure | `try_new`, `point`, `try_from_subtraction`, `try_add`, `try_mul`, `negate`, `try_square` | -| `IntervalMatrix` | `[[Interval; D]; D]` | Division-free determinant enclosure and sign proof through D=7 | `from_rows`, `try_from_point_rows`, `from_matrix`, `det`, `det_sign` | -| `IntervalDeterminantSign` | enum | Positive, negative, zero, or inconclusive determinant evidence | — | -| `RationalVector`¹ | `[BigRational; D]` | Exact rational right-hand side and solution | `try_new`, `try_from_fn`, `as_array`, `into_array`, `get` | -| `RationalMatrix`¹ | `[[BigRational; D]; D]` | Exact rational matrix for determinant and solve operations | `try_from_rows`, `try_from_fn`, `as_rows`, `det_sign`, `det`, `solve` | -| `DeterminantWithErrorBound` | Opaque validated pair | Paired direct determinant and certified absolute bound | `determinant`, `absolute_error_bound` | -| `ScalarWithErrorBound` | Opaque validated certificate | Paired scalar estimate, absolute bound, and outward endpoints | `estimate`, `absolute_error_bound`, `lower_bound`, `upper_bound` | -| `Lu` | Inline factors + permutation | Factorization for solves/det | `solve`, `det` | -| `Ldlt` | Inline factors | No-pivot SPD factorization for solves/det | `solve`, `det` | -| `Tolerance` | finite non-negative `f64` | Validated numerical threshold | `try_new`, `get` | -| `LaError` | typed variants and reasons | Structured, actionable failure reporting | See error semantics below | -| `DeterminantSign`¹ | enum | Exact determinant sign | `as_i8` | -| `ExactF64Conversion`¹ | trait | Strict or explicitly rounded conversion of exact results to `f64` | `try_to_f64`, `to_rounded_f64` | - -`Matrix` and `Vector` use the intentional inline `f64` scalar model. -`IntervalMatrix` retains inline fixed-size storage and uses a fixed 128-entry -stack workspace for its supported determinant dimensions. The exact-feature -rational types retain fixed-size outer arrays while their `BigRational` scalars -use arbitrary-precision integer storage. - -For a runtime-selected interval dimension from 0 through -`MAX_INTERVAL_MATRIX_DIM` (7), `try_with_interval_matrix!` dispatches to a -concrete `IntervalMatrix`. This is useful when stable Rust cannot express a -derived const dimension such as `D + 1`. - -For a runtime dimension from 0 through `MAX_STACK_MATRIX_DISPATCH_DIM` (7), -`try_with_stack_matrix!` dispatches to a concrete `Matrix` while preserving -inline stack storage. Larger dimensions produce -`LaError::UnsupportedDimension`, converted through `From` into the -closure's declared `Result` error type; the macro does not introduce a -dynamically sized matrix representation. - -With the exact feature, a runtime dimension from 0 through -`MAX_RATIONAL_MATRIX_DISPATCH_DIM` (8) can similarly be dispatched with -`try_with_rational_matrix!` to a concrete `RationalMatrix`. Larger -dimensions produce `LaError::UnsupportedDimension`; the rational macro also -preserves the const-generic representation rather than introducing a -dynamically sized matrix type. - -`Matrix` key methods: `as_rows`, `into_rows`, `lu`, `ldlt`, `det`, -`det_direct`, `det_direct_with_errbound`, `det_errbound`, -`det_exact`¹, `det_exact_f64`¹, `det_exact_rounded_f64`¹, `det_sign_exact`¹, -`solve_exact`¹, `solve_exact_f64`¹, `solve_exact_rounded_f64`¹. -Matrix and vector constructors validate non-finite inputs at public API -boundaries. After construction, `Matrix` and `Vector` carry that -finite-storage invariant directly, so factorization kernels do not repeat an -O(D²) input scan. Computed factor matrices are still checked before they become -observable results. - -`Matrix::as_rows` and `Vector::as_array` borrow their validated backing arrays; -`Matrix::into_rows` and `Vector::into_array` consume the value and return the -owned fixed-size arrays. - -`Matrix::get(row, col)` returns `None` for out-of-bounds coordinates; -`Matrix::try_get` instead returns a structured `LaError` preserving those -coordinates. The single fallible `Matrix::set` validates both coordinates and -finiteness before mutating the matrix. - -`LaError` and its reason/location enums are non-exhaustive. Numerical -singularity records the `FactorizationKind`, -observed pivot magnitude, and tolerance, while exact-arithmetic singularity is -identified separately. `LaError::NonFinite` retains the crate-wide non-finite -contract but uses `NonFiniteOrigin`, `NonFiniteLocation`, and -`ArithmeticOperation` to distinguish invalid inputs from computed overflow. -`LaError::InvertedInterval` preserves rejected finite bounds when the lower -endpoint exceeds the upper endpoint, and `LaError::IntervalRangeExhausted` -distinguishes finite-input interval range loss from a non-finite value. -`InvalidToleranceReason` distinguishes negative from non-finite tolerances, and -`PositiveSemidefiniteViolation` distinguishes negative LDLT pivots from a zero -pivot with nonzero coupling. Match these public enums with a wildcard and use -`..` for struct-style variants so future error context can be added without -breaking callers. +Start with the capability you need; the [API reference][api-reference] lists +the complete public surface, and the [worked examples][api-guide] show how +to combine operations. + +| Capability | Main entry points | +|---|---| +| Vector operations and norms | [`Vector`][api-vector] | +| Floating-point determinants and solves | [`Matrix`][api-matrix], [`Lu`][api-lu], [`Ldlt`][api-ldlt] | +| Gram matrix construction | [`gram_matrix`][api-gram] | +| Certified dot, affine-difference, and determinant estimates | [`ScalarWithErrorBound`][api-scalar-bound], [`DeterminantWithErrorBound`][api-det-bound] | +| Interval expressions and determinant signs | [`Interval`][api-interval], [`IntervalMatrix`][api-interval-matrix] | +| Exact signs, determinants, solves, and output conversion¹ | [Exact arithmetic examples][api-exact] | +| Runtime selection of a const-generic matrix dimension | [Dimension dispatch examples][api-dispatch] | + +[`Tolerance`][api-tolerance] validates numerical rejection thresholds. +[`LaError`][api-error] and its reason/location enums preserve structured +failure details; match non-exhaustive enums with a wildcard and struct-style +variants with `..`. See the [storage, access, and error guide][api-contracts] +for the full contracts. ¹ Requires `features = ["exact"]`. -## 📊 Benchmarks (vs nalgebra/faer) +[api-reference]: https://docs.rs/la-stack/latest/la_stack/ +[api-guide]: https://docs.rs/la-stack/latest/la_stack/guide/index.html +[api-vector]: https://docs.rs/la-stack/latest/la_stack/struct.Vector.html +[api-matrix]: https://docs.rs/la-stack/latest/la_stack/struct.Matrix.html +[api-lu]: https://docs.rs/la-stack/latest/la_stack/struct.Lu.html +[api-ldlt]: https://docs.rs/la-stack/latest/la_stack/struct.Ldlt.html +[api-gram]: https://docs.rs/la-stack/latest/la_stack/fn.gram_matrix.html +[api-scalar-bound]: https://docs.rs/la-stack/latest/la_stack/struct.ScalarWithErrorBound.html +[api-det-bound]: https://docs.rs/la-stack/latest/la_stack/struct.DeterminantWithErrorBound.html +[api-interval]: https://docs.rs/la-stack/latest/la_stack/struct.Interval.html +[api-interval-matrix]: https://docs.rs/la-stack/latest/la_stack/struct.IntervalMatrix.html +[api-exact]: https://docs.rs/la-stack/latest/la_stack/guide/exact/index.html +[api-dispatch]: https://docs.rs/la-stack/latest/la_stack/guide/index.html#dimension-dispatch +[api-tolerance]: https://docs.rs/la-stack/latest/la_stack/struct.Tolerance.html +[api-error]: https://docs.rs/la-stack/latest/la_stack/enum.LaError.html +[api-contracts]: https://docs.rs/la-stack/latest/la_stack/guide/index.html#storage-access-and-errors + + + +## 🗺️ Documentation Map + +- [API guide][api-guide] — worked examples, API selection, storage, and error contracts. +- [Mathematical basis](https://github.com/acgetchell/la-stack/blob/v0.4.5/docs/mathematical_basis.md) — algorithms, numerical guarantees, and limitations. +- [Benchmarking](https://github.com/acgetchell/la-stack/blob/v0.4.5/docs/BENCHMARKING.md) — benchmark suites, comparison workflows, and measurement methodology. +- [Performance reports](https://github.com/acgetchell/la-stack/blob/v0.4.5/docs/PERFORMANCE.md) — release-to-release measurement results and provenance. +- [Coverage](https://github.com/acgetchell/la-stack/blob/v0.4.5/docs/COVERAGE.md) — local and CI coverage commands and report locations. +- [Roadmap](https://github.com/acgetchell/la-stack/blob/v0.4.5/docs/roadmap.md) — release planning, future directions, and non-goals. +- [Releasing](https://github.com/acgetchell/la-stack/blob/v0.4.5/docs/RELEASING.md) — release preparation, validation, and publication. + +## 📈 Benchmarks (vs nalgebra/faer) ![LU solve (factor + solve): median time vs dimension][lu-solve-benchmark] @@ -928,7 +486,7 @@ For coverage commands and report locations, see For the full contributor workflow, see [CONTRIBUTING.md](https://github.com/acgetchell/la-stack/blob/v0.4.5/CONTRIBUTING.md). -## 📝 Citation +## 📚 Citation If you use this library in academic work, please cite it using [CITATION.cff](https://github.com/acgetchell/la-stack/blob/v0.4.5/CITATION.cff) @@ -936,7 +494,7 @@ If you use this library in academic work, please cite it using Zenodo under the [all-versions concept DOI](https://doi.org/10.5281/zenodo.18158926). -## 📚 References +## 🔎 References For canonical references to the algorithms used by this crate, see [REFERENCES.md](https://github.com/acgetchell/la-stack/blob/v0.4.5/REFERENCES.md). @@ -949,7 +507,7 @@ before proposing or applying changes. See [CONTRIBUTING.md](https://github.com/acgetchell/la-stack/blob/v0.4.5/CONTRIBUTING.md) for the repository's AI-assisted development note. -## 📄 License +## 📜 License BSD 3-Clause License. See [LICENSE](https://github.com/acgetchell/la-stack/blob/v0.4.5/LICENSE). diff --git a/SECURITY.md b/SECURITY.md index a365ccc..b56485a 100644 --- a/SECURITY.md +++ b/SECURITY.md @@ -101,3 +101,36 @@ advisories or release notes unless anonymity is requested. This project uses GitHub CodeQL, Dependabot security updates, secret scanning with push protection, `cargo audit`, zizmor, Clippy SARIF analysis, and repository-owned Semgrep rules. + +### Numerical logging false positives + +CodeQL's `rust/cleartext-logging` query uses name-based heuristics to identify +potentially sensitive data. Here, `Matrix::certified_error_bound` and the +benchmark helper's `certified_bound` refer to numerical rounding-error bounds. +They contain no credentials, cryptographic certificates, or personal information. + +Alerts #171–#178 were reviewed against commit +`bd80cc05df3ebf409d8db9a471b671a5737bf0c4`: their sinks are assertion or panic +diagnostics for deterministic test and benchmark fixtures. These individual +alerts were dismissed as false positives with that rationale. Keep the query +enabled and preserve diagnostic values needed to investigate numerical failures. +Review new alerts on their own data flow. + +### Benchmark dependency maintenance + +As of September 7, 2026, `paste` 1.0.15 enters through the development dependency +`faer` 0.24.4, via `gemm` 0.19.0 and `pulp` 0.22.3. Those are the latest +published upstream versions checked on that date. Repository-owned Rust code +already uses `pastey`; changing that direct dependency cannot replace upstream +uses of `paste`. + +[RUSTSEC-2024-0436](https://rustsec.org/advisories/RUSTSEC-2024-0436.html) is an +unmaintained-package advisory with no patched version. It is a genuine +maintenance concern, separate from the logging false positives. `paste` is absent +from the library's normal and build dependency graph, including with `exact` +enabled, but its procedural macro executes when building development targets. + +Keep this advisory visible in `cargo audit`. Recheck the dependency path with +`cargo tree --locked --all-features -i paste` when updating `faer`, `gemm`, or +`pulp`, and remove the dependency through a maintained upstream release when +available. diff --git a/justfile b/justfile index 209c204..affc2dc 100644 --- a/justfile +++ b/justfile @@ -18,7 +18,7 @@ _coverage_base_args := '''--features exact \ --workspace --lib --tests \ --verbose''' cargo_edit_version := "0.13.13" -cargo_llvm_cov_version := "0.9.0" +cargo_llvm_cov_version := "0.9.1" cargo_machete_version := "0.9.2" cargo_nextest_version := "0.9.143" cargo_update_version := "22.1.1" @@ -26,7 +26,7 @@ clippy_sarif_version := "0.8.0" dprint_version := "0.57.4" git_cliff_version := "2.14.1" just_version := "1.58.0" -rumdl_version := "0.2.66" +rumdl_version := "0.2.67" sarif_fmt_version := "0.8.0" taplo_version := "0.10.0" typos_version := "1.50.1" @@ -571,7 +571,7 @@ lint-config: json-check toml-ci yaml-ci github-actions-check justfile-fmt-check lint-docs: markdown-ci docs-version-check # Markdown -markdown-check: _ensure-rumdl +markdown-check: _ensure-rumdl _ensure-uv #!/usr/bin/env bash set -euo pipefail files=() @@ -585,24 +585,7 @@ markdown-check: _ensure-rumdl done < <(git ls-files -co --exclude-standard -z -- '*.md') if [ "${#files[@]}" -gt 0 ]; then printf '%s\0' "${files[@]}" | xargs -0 -n100 rumdl check - violations=0 - for file in "${files[@]}"; do - line_number=0 - while IFS= read -r line || [[ -n "$line" ]]; do - line_number=$((line_number + 1)) - case "$line" in - '|'*) continue ;; - esac - if [ "${#line}" -gt 160 ]; then - printf '%s:%d: line length %d exceeds 160\n' "$file" "$line_number" "${#line}" >&2 - violations=$((violations + 1)) - fi - done < "$file" - done - if [ "$violations" -gt 0 ]; then - echo "Markdown raw line-length check failed." >&2 - exit 1 - fi + uv run --locked scripts/check_markdown_lines.py "${files[@]}" else echo "No markdown files found to check." fi diff --git a/scripts/check_markdown_lines.py b/scripts/check_markdown_lines.py new file mode 100644 index 0000000..94a5744 --- /dev/null +++ b/scripts/check_markdown_lines.py @@ -0,0 +1,34 @@ +"""Check raw Markdown line lengths using Unicode characters in every locale.""" + +import argparse +import sys +from pathlib import Path + +MAX_LINE_LENGTH = 160 + + +def main() -> int: + """Check UTF-8 files, preserving the Markdown recipe's table exemption.""" + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("files", nargs="+", type=Path) + args = parser.parse_args() + failed = False + for path in args.files: + try: + with path.open(encoding="utf-8") as source: + for line_number, line in enumerate(source, start=1): + # Text mode normalizes CRLF; whitespace still counts toward the limit. + line = line.removesuffix("\n") + if not line.startswith("|") and len(line) > MAX_LINE_LENGTH: + print(f"{path}:{line_number}: line length {len(line)} exceeds {MAX_LINE_LENGTH}", file=sys.stderr) + failed = True + except (OSError, UnicodeError) as error: + print(f"{path}: {error}", file=sys.stderr) + failed = True + if failed: + print("Markdown raw line-length check failed.", file=sys.stderr) + return int(failed) + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/scripts/tests/test_check_markdown_lines.py b/scripts/tests/test_check_markdown_lines.py new file mode 100644 index 0000000..4787218 --- /dev/null +++ b/scripts/tests/test_check_markdown_lines.py @@ -0,0 +1,51 @@ +"""Regression coverage for locale-independent Markdown line-length checks.""" + +import os +import subprocess +import sys +from pathlib import Path + +import pytest + +CHECKER = Path(__file__).resolve().parents[1] / "check_markdown_lines.py" + + +@pytest.mark.parametrize( + ("content", "diagnostic"), + [ + (("a" * 159 + "—\n").encode(), ""), + (("a" * 159 + "🗺\r\n").encode(), ""), + (("a" * 159 + "—").encode(), ""), + (("a" * 160 + "—\n").encode(), ":1: line length 161 exceeds 160"), + (b"a" * 161, ":1: line length 161 exceeds 160"), + (b"a" * 159 + b" \n", ":1: line length 161 exceeds 160"), + (b"|" + b"a" * 200 + b"\n", ""), + (b"short\r\n" + b"a" * 161 + b"\r\n", ":2: line length 161 exceeds 160"), + (b"invalid UTF-8: \xff\n", "decode"), + ], + ids=["unicode-limit", "emoji-crlf", "no-final-newline", "unicode-too-long", "ascii-too-long", "trailing-spaces", "table", "line-number", "invalid-utf8"], +) +def test_raw_line_limit_under_c_locale(tmp_path: Path, content: bytes, diagnostic: str) -> None: + """Count UTF-8 characters under the locale that exposed the Windows failure.""" + markdown = tmp_path / "document with spaces.md" + markdown.write_bytes(content) + environment = os.environ.copy() + environment["LC_ALL"] = "C" + result = subprocess.run( # noqa: S603 - fixed local checker and test-owned input. + [sys.executable, str(CHECKER), str(markdown)], + check=False, + capture_output=True, + encoding="utf-8", + env=environment, + timeout=30, + ) + + assert result.stdout == "" + if diagnostic: + assert result.returncode == 1 + assert str(markdown) in result.stderr + assert diagnostic in result.stderr + assert "Markdown raw line-length check failed." in result.stderr + else: + assert result.returncode == 0 + assert result.stderr == "" diff --git a/src/lib.rs b/src/lib.rs index 87b8cb9..24dd85b 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -5,315 +5,904 @@ #[cfg(doc)] mod readme_doctests { - //! Executable versions of README examples. + //! Executable version of the README quickstart. /// ```rust /// use la_stack::prelude::*; /// - /// # fn main() -> Result<(), LaError> { - /// // This system requires pivoting (a[0][0] = 0), so it's a good LU demo. - /// let a = Matrix::<5>::try_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], - /// ])?; + /// fn main() -> Result<(), LaError> { + /// // The zero leading entry requires LU pivoting. + /// let a = Matrix::<5>::try_from_rows([ + /// [0.0, 2.0, -1.0, 1.0, 3.0], + /// [4.0, -1.0, 2.0, 0.0, 1.0], + /// [1.0, 3.0, 5.0, -2.0, 0.0], + /// [2.0, 0.0, -1.0, 4.0, 1.0], + /// [-1.0, 2.0, 0.0, 1.0, 6.0], + /// ])?; + /// let b = Vector::try_new([20.0, 13.0, 14.0, 20.0, 37.0])?; + /// let lu = a.lu(DEFAULT_SINGULAR_TOL)?; + /// let x = lu.solve(b)?; /// - /// let b = Vector::<5>::try_new([14.0, 13.0, 12.0, 11.0, 10.0])?; - /// - /// let lu = a.lu(DEFAULT_SINGULAR_TOL)?; - /// let x = lu.solve(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); + /// for (&actual, expected) in x.as_array().iter().zip([1.0, 2.0, 3.0, 4.0, 5.0]) { + /// assert!((actual - expected).abs() <= 1e-12); + /// } + /// Ok(()) /// } - /// # Ok(()) - /// # } /// ``` fn solve_5x5_example() {} +} - /// ```rust - /// use la_stack::prelude::*; - /// - /// # fn main() -> Result<(), LaError> { - /// // This matrix is symmetric positive-definite (A = L*L^T) so LDLT works without pivoting. - /// let a = Matrix::<5>::try_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 ldlt = match a.ldlt(DEFAULT_SINGULAR_TOL) { - /// Ok(ldlt) => ldlt, - /// Err(err @ LaError::Asymmetric { row, col, .. }) => { - /// eprintln!("LDLT requires symmetry; first mismatch at ({row}, {col})"); - /// return Err(err); - /// } - /// Err(err) => return Err(err), - /// }; - /// - /// let det = ldlt.det()?; - /// assert!((det - 1.0).abs() <= 1e-12); - /// # Ok(()) - /// # } - /// ``` - fn det_5x5_ldlt_example() {} +// Documentation-only workflows keep the README overview compact without adding +// a runtime API. Their examples are executed by the default/exact doctest gates. +#[cfg(doc)] +pub mod guide { + //! Worked examples and contracts for choosing and combining APIs. + //! + //! Start with [`prelude`](crate::prelude) for common imports. The crate's + //! generated reference lists the complete public surface; these examples + //! demonstrate how the pieces fit together. + //! + //! Explore the topic guides for more workflows: + //! + //! - [LDLT determinants and symmetry](ldlt) + //! - [Compile-time determinants](compile_time) + //! - [Outward-rounded interval determinants](intervals) + //! - [Overflow-safe Euclidean norms](norms) + //! - [Certified dot products and affine differences](certified) + //! - [Adaptive determinant filtering](adaptive) + //! + //! Enabling `exact` also adds the exact-arithmetic guide to the module list. + //! + //! # Solving and reusing factors + //! + //! [`Matrix::lu`](crate::Matrix::lu) computes a partially pivoted + //! factorization. Keep the resulting [`Lu`](crate::Lu) to solve multiple + //! right-hand sides without repeating factorization. This 5×5 system has a + //! zero leading entry, so the first elimination step requires pivoting. + //! + //! ```rust + //! use la_stack::prelude::*; + //! + //! # fn main() -> Result<(), LaError> { + //! let a = Matrix::<5>::try_from_rows([ + //! [0.0, 2.0, -1.0, 1.0, 3.0], + //! [4.0, -1.0, 2.0, 0.0, 1.0], + //! [1.0, 3.0, 5.0, -2.0, 0.0], + //! [2.0, 0.0, -1.0, 4.0, 1.0], + //! [-1.0, 2.0, 0.0, 1.0, 6.0], + //! ])?; + //! let lu = a.lu(DEFAULT_SINGULAR_TOL)?; + //! let systems = [ + //! ([20.0, 13.0, 14.0, 20.0, 37.0], [1.0, 2.0, 3.0, 4.0, 5.0]), + //! ([5.0, 6.0, 7.0, 6.0, 8.0], [1.0; 5]), + //! ]; + //! for (rhs, expected) in systems { + //! let solution = lu.solve(Vector::try_new(rhs)?)?; + //! for (&actual, expected) in solution.as_array().iter().zip(expected) { + //! assert!((actual - expected).abs() <= 1e-12); + //! } + //! } + //! # Ok(()) + //! # } + //! ``` + //! + //! The assertions use a tolerance suitable for this known fixture; they do + //! not establish a general error bound for LU. Factorization tolerances + //! reject small pivots and are not accuracy guarantees. + //! [`Ldlt`](crate::Ldlt) offers the same solve/determinant workflow for + //! exactly symmetric positive-definite input, without pivoting. Approximate + //! symmetry from [`Matrix::is_symmetric`](crate::Matrix::is_symmetric) or + //! [`Matrix::first_asymmetry`](crate::Matrix::first_asymmetry) does not prove + //! the exact symmetry required by [`Matrix::ldlt`](crate::Matrix::ldlt). + //! + //! # Gram matrices + //! + //! [`gram_matrix`](crate::gram_matrix) accepts `M` vectors of dimension `N` + //! and returns a `Matrix` of pairwise inner products. Here five vectors + //! in six dimensions produce a 5×5 matrix. Each dot product is computed once + //! and mirrored, so the result is bit-for-bit symmetric. + //! + //! ```rust + //! use la_stack::prelude::*; + //! + //! # fn main() -> Result<(), LaError> { + //! let vectors = [ + //! Vector::try_new([1.0, 1.0, 0.0, 0.0, 0.0, 0.0])?, + //! Vector::try_new([0.0, 1.0, 1.0, 0.0, 0.0, 0.0])?, + //! Vector::try_new([0.0, 0.0, 1.0, 1.0, 0.0, 0.0])?, + //! Vector::try_new([0.0, 0.0, 0.0, 1.0, 1.0, 0.0])?, + //! Vector::try_new([0.0, 0.0, 0.0, 0.0, 1.0, 1.0])?, + //! ]; + //! let gram = gram_matrix(&vectors)?; + //! assert_eq!(gram.norm_inf()?, 4.0); + //! assert!(gram.is_symmetric(Tolerance::try_new(0.0)?)?); + //! + //! // This fixture's exact Gram matrix is tridiagonal: 2 on the diagonal, + //! // 1 immediately above/below it. Its 5×5 determinant is 6. + //! let determinant = gram.ldlt(DEFAULT_SINGULAR_TOL)?.det()?; + //! assert!((determinant - 6.0).abs() <= 1e-12); + //! # Ok(()) + //! # } + //! ``` + //! + //! Gram construction provides no certified rounding-error bound and does + //! not prove rank or positive definiteness. The generated function + //! documentation explains its conditioning and geometric interpretation. + //! + //! # Dimension dispatch + //! + //! [`try_with_stack_matrix!`](crate::try_with_stack_matrix) selects a + //! concrete `Matrix` for runtime dimensions 0 through + //! [`MAX_STACK_MATRIX_DISPATCH_DIM`](crate::MAX_STACK_MATRIX_DISPATCH_DIM) + //! (7). The closure receives a zero matrix and returns its declared result. + //! + //! ```rust + //! use core::assert_matches; + //! + //! use la_stack::prelude::*; + //! + //! # fn main() -> Result<(), LaError> { + //! let requested = 5usize; + //! let determinant = try_with_stack_matrix!(requested, |mut matrix| -> Result { + //! for row in 0..requested { + //! matrix.set(row, row, 2.0)?; + //! if row + 1 < requested { + //! matrix.set(row, row + 1, 1.0)?; + //! matrix.set(row + 1, row, 1.0)?; + //! } + //! } + //! matrix.det() + //! })?; + //! assert!((determinant - 6.0).abs() <= 1e-12); + //! + //! let unsupported = try_with_stack_matrix!(8, |matrix| -> Result { + //! matrix.det() + //! }); + //! assert_matches!( + //! unsupported, + //! Err(LaError::UnsupportedDimension { requested: 8, max: 7, .. }) + //! ); + //! # Ok(()) + //! # } + //! ``` + //! + //! [`try_with_interval_matrix!`](crate::try_with_interval_matrix) similarly + //! dispatches dimensions 0 through + //! [`MAX_INTERVAL_MATRIX_DIM`](crate::MAX_INTERVAL_MATRIX_DIM) (7) to an + //! [`IntervalMatrix`](crate::IntervalMatrix). These macros are useful when + //! stable Rust cannot express a derived const dimension such as `D + 1`. + //! Larger dimensions produce [`LaError::UnsupportedDimension`](crate::LaError::UnsupportedDimension), + //! converted through `From` into the closure's declared error type. + //! Dispatch preserves const-generic storage; it does not create a dynamically + //! sized matrix representation or limit dimensions chosen directly at compile time. + //! + //! # Storage, access, and errors + //! + //! [`Matrix`](crate::Matrix) and [`Vector`](crate::Vector) store + //! `[[f64; D]; D]` and `[f64; D]` inline. Constructors validate non-finite + //! inputs, and the types preserve that finite-storage invariant. Factorization + //! kernels therefore avoid a repeated O(D²) input scan; computed factor + //! matrices are still checked before becoming observable results. + //! + //! [`Matrix::as_rows`](crate::Matrix::as_rows) and + //! [`Vector::as_array`](crate::Vector::as_array) borrow validated backing + //! arrays. [`Matrix::into_rows`](crate::Matrix::into_rows) and + //! [`Vector::into_array`](crate::Vector::into_array) consume the value and + //! return owned fixed-size arrays. + //! [`Matrix::get`](crate::Matrix::get) returns `None` for invalid coordinates; + //! [`Matrix::try_get`](crate::Matrix::try_get) preserves them in a typed error. + //! [`Matrix::set`](crate::Matrix::set) checks coordinates and finiteness + //! before mutation. [`Matrix::norm_inf`](crate::Matrix::norm_inf) computes + //! the maximum absolute row sum. + //! + //! [`Vector::dot`](crate::Vector::dot), [`Vector::norm`](crate::Vector::norm), + //! and [`Vector::norm_squared`](crate::Vector::norm_squared) provide ordinary + //! vector reductions. [`ScalarWithErrorBound`](crate::ScalarWithErrorBound) + //! is the opaque result of the certified dot and affine-difference methods; + //! it exposes the estimate, absolute bound, and outward-rounded endpoints. + //! [`DeterminantWithErrorBound`](crate::DeterminantWithErrorBound) pairs a + //! direct determinant with its certified absolute bound. Use + //! [`Matrix::det_errbound`](crate::Matrix::det_errbound) for the bound alone. + //! + //! [`Interval`](crate::Interval) stores two finite ordered bounds and supports + //! point construction, outward-rounded subtraction, addition, multiplication, + //! negation, and square. [`IntervalMatrix`](crate::IntervalMatrix) stores + //! `[[Interval; D]; D]` inline and uses a fixed 128-entry stack workspace for + //! supported determinant dimensions. [`IntervalDeterminantSign`](crate::IntervalDeterminantSign) + //! distinguishes positive, negative, exact zero, and inconclusive evidence. + //! + //! Parse numerical thresholds through [`Tolerance::try_new`](crate::Tolerance::try_new). + //! [`LaError`](crate::LaError) and its reason/location enums are non-exhaustive; + //! use wildcard match arms and `..` for struct-style variants. In particular: + //! + //! - [`SingularityReason`](crate::SingularityReason) separates exact singularity + //! from numerical rejection, retaining the [`FactorizationKind`](crate::FactorizationKind), + //! observed pivot magnitude, and tolerance for the latter. + //! - [`NonFiniteOrigin`](crate::NonFiniteOrigin), [`NonFiniteLocation`](crate::NonFiniteLocation), + //! and [`ArithmeticOperation`](crate::ArithmeticOperation) distinguish invalid + //! inputs from computed non-finite values. + //! - [`LaError::InvertedInterval`](crate::LaError::InvertedInterval) preserves + //! rejected finite endpoints; [`LaError::IntervalRangeExhausted`](crate::LaError::IntervalRangeExhausted) + //! distinguishes finite-input interval range loss from a non-finite value. + //! - [`InvalidToleranceReason`](crate::InvalidToleranceReason) distinguishes + //! negative and non-finite tolerances. + //! - [`PositiveSemidefiniteViolation`](crate::PositiveSemidefiniteViolation) + //! distinguishes a negative LDLT pivot from a zero pivot with nonzero coupling. - /// ```rust - /// use la_stack::prelude::*; - /// - /// // Evaluated entirely at compile time — no runtime cost. - /// const DET: Result, LaError> = match Matrix::<4>::try_from_rows([ - /// [2.0, 0.0, 0.0, 0.0], - /// [0.0, 3.0, 0.0, 0.0], - /// [0.0, 0.0, 5.0, 0.0], - /// [0.0, 0.0, 0.0, 7.0], - /// ]) { - /// Ok(matrix) => matrix.det_direct(), - /// Err(err) => Err(err), - /// }; - /// - /// # fn main() -> Result<(), LaError> { - /// assert_eq!(DET?, Some(210.0)); - /// # Ok(()) - /// # } - /// ``` - fn det_direct_4x4_const_example() {} + pub mod ldlt { + //! LDLT determinants and exact symmetry. + //! + //! Compute a determinant for a symmetric positive-definite matrix via LDLT (no + //! pivoting). + //! + //! For these matrices, `LDLᵀ` is a square-root-free Cholesky form. Multiplying each + //! column of `L` by the square root of the corresponding diagonal entry yields a + //! Cholesky factor: + //! + //! ```rust + //! use la_stack::prelude::*; + //! + //! fn main() -> Result<(), LaError> { + //! // This matrix is symmetric positive-definite (A = L*L^T) so LDLT works without pivoting. + //! let a = Matrix::<5>::try_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 ldlt = match a.ldlt(DEFAULT_SINGULAR_TOL) { + //! Ok(ldlt) => ldlt, + //! Err(err @ LaError::Asymmetric { + //! row, + //! col, + //! upper, + //! lower, + //! allowed_abs_diff, + //! .. + //! }) => { + //! eprintln!( + //! "LDLT mismatch at ({row}, {col}): {upper} vs {lower} (allowed {allowed_abs_diff})" + //! ); + //! return Err(err); + //! } + //! Err(err) => return Err(err), + //! }; + //! + //! let det = ldlt.det()?; + //! assert!((det - 1.0).abs() <= 1e-12); + //! + //! Ok(()) + //! } + //! ``` + //! + //! > ⚠️ **LDLT invariant:** The input matrix must be **exactly symmetric**: every + //! > mirrored pair must compare equal (`+0.0 == -0.0` is accepted). Asymmetric + //! > inputs passed to + //! > [`Matrix::ldlt`](crate::Matrix::ldlt) + //! > return a typed `LaError::Asymmetric` containing both observed values and the + //! > required allowed difference of zero. The tolerance-based + //! > [`Matrix::first_asymmetry`](crate::Matrix::first_asymmetry) + //! > and [`Matrix::is_symmetric`](crate::Matrix::is_symmetric) methods remain useful diagnostics, but do not prove + //! > the exact precondition required by LDLT. Use `lu()` when exact symmetry or + //! > positive definiteness is not guaranteed. A negative LDLT diagonal or a zero + //! > diagonal with nonzero remaining coupling returns + //! > `LaError::NotPositiveSemidefinite` with a typed + //! > `PositiveSemidefiniteViolation`. An uncoupled zero or positive pivot + //! > at or below the caller's tolerance returns `LaError::Singular` with a + //! > numerical `SingularityReason`. Because these pivots are computed in binary64, + //! > success is not an exact positive-definiteness certificate for the stored + //! > matrix. + } - /// ```rust - /// use la_stack::prelude::*; - /// - /// # fn main() -> Result<(), LaError> { - /// let x = Interval::try_from_subtraction(0.1, 0.0)?; - /// let y = Interval::try_from_subtraction(0.1, 0.0)?; - /// let z = Interval::try_from_subtraction(0.1, 0.0)?; - /// let lifted = x - /// .try_square()? - /// .try_add(&y.try_square()?)? - /// .try_add(&z.try_square()?)?; - /// - /// let matrix = IntervalMatrix::<4>::from_rows([ - /// [Interval::ONE, Interval::ZERO, Interval::ZERO, Interval::ONE], - /// [Interval::ZERO, Interval::ONE, Interval::ZERO, Interval::ONE], - /// [Interval::ZERO, Interval::ZERO, Interval::ONE, Interval::ONE], - /// [x, y, z, lifted], - /// ]); - /// assert_eq!( - /// matrix.det_sign()?, - /// IntervalDeterminantSign::Negative, - /// ); - /// # Ok(()) - /// # } - /// ``` - fn interval_determinant_example() {} + pub mod compile_time { + //! Compile-time determinants and dimension dispatch. + //! + //! [`det_direct()`](crate::Matrix::det_direct) is a `const fn` providing closed-form determinants for D=0–4, + //! using fused multiply-add where applicable. It returns `Ok(Some(det))` for those + //! dimensions and `Ok(None)` for D ≥ 5. `Matrix::<0>::zero().det_direct()` returns + //! `Ok(Some(1.0))` (the empty-product convention). For D=1–4, direct formulas + //! bypass 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: Result, LaError> = match Matrix::<4>::try_from_rows([ + //! [2.0, 0.0, 0.0, 0.0], + //! [0.0, 3.0, 0.0, 0.0], + //! [0.0, 0.0, 5.0, 0.0], + //! [0.0, 0.0, 0.0, 7.0], + //! ]) { + //! Ok(matrix) => matrix.det_direct(), + //! Err(err) => Err(err), + //! }; + //! + //! fn main() -> Result<(), LaError> { + //! assert_eq!(DET?, Some(210.0)); + //! Ok(()) + //! } + //! ``` + //! + //! The public `det()` method automatically dispatches through the closed-form path + //! for D ≤ 4 and falls back to zero-tolerance LU for D ≥ 5. Tiny nonzero + //! determinants are not flattened by a configured pivot tolerance. The LU fallback + //! returns `LaError::Singular` when floating-point elimination cannot produce a + //! non-zero pivot; it does not misreport that numerical failure as an exact zero. + //! Use `lu()` directly when you need a different tolerance policy, and use the + //! exact determinant APIs when exact singularity classification matters. + } - /// ```rust - /// use la_stack::prelude::*; - /// - /// fn is_separated( - /// axis: &Vector, - /// left: &Vector, - /// right: &Vector, - /// threshold: f64, - /// ) -> Result, LaError> { - /// let Some(value) = axis.dot_difference_with_errbound(left, right)? else { - /// return Ok(None); - /// }; - /// if value.lower_bound() > threshold { - /// Ok(Some(true)) - /// } else if value.upper_bound() <= threshold { - /// Ok(Some(false)) - /// } else { - /// Ok(None) - /// } - /// } - /// - /// # fn main() -> Result<(), LaError> { - /// let axis = Vector::<2>::try_new([2.0, -1.0])?; - /// let left = Vector::<2>::try_new([4.0, 1.0])?; - /// let right = Vector::<2>::try_new([1.0, 3.0])?; - /// assert_eq!(is_separated(&axis, &left, &right, 1.0)?, Some(true)); - /// # Ok(()) - /// # } - /// ``` - fn certified_linear_form_example() {} + pub mod intervals { + //! Outward-rounded interval expressions and determinant signs. + //! + //! `Interval` encloses expression construction that has not yet been reduced to a + //! single stored `f64`. Point intervals preserve finite binary64 values exactly; + //! `try_from_subtraction`, `try_add`, `try_mul`, `negate`, and `try_square` enclose + //! the corresponding exact-real operations. [`IntervalMatrix::det_sign()`](crate::IntervalMatrix::det_sign) then + //! uses a division-free subset expansion through D=7, returning positive, + //! negative, zero, or inconclusive evidence. + //! + //! ```rust + //! use la_stack::prelude::*; + //! + //! fn main() -> Result<(), LaError> { + //! // Relative coordinates and the lifted norm retain their construction error. + //! let x = Interval::try_from_subtraction(0.1, 0.0)?; + //! let y = Interval::try_from_subtraction(0.1, 0.0)?; + //! let z = Interval::try_from_subtraction(0.1, 0.0)?; + //! let lifted = x + //! .try_square()? + //! .try_add(&y.try_square()?)? + //! .try_add(&z.try_square()?)?; + //! + //! let matrix = IntervalMatrix::<4>::from_rows([ + //! [Interval::ONE, Interval::ZERO, Interval::ZERO, Interval::ONE], + //! [Interval::ZERO, Interval::ONE, Interval::ZERO, Interval::ONE], + //! [Interval::ZERO, Interval::ZERO, Interval::ONE, Interval::ONE], + //! [x, y, z, lifted], + //! ]); + //! assert_eq!( + //! matrix.det_sign()?, + //! IntervalDeterminantSign::Negative, + //! ); + //! Ok(()) + //! } + //! ``` + //! + //! Every successful interval keeps finite ordered endpoints. Subnormal bounds are + //! preserved, both signed zeros are treated as real zero and canonicalized to + //! `+0.0`, and underflowed nonzero products widen toward the least subnormal value. + //! If an exact result range cannot fit between finite binary64 endpoints, the + //! operation returns `LaError::IntervalRangeExhausted` with its interval operation + //! recorded in `ArithmeticOperation`. + //! + //! `Positive`, `Negative`, and `Zero` are proofs. `Inconclusive` only means that + //! the determinant enclosure overlaps zero; it must not be converted to equality + //! or singularity. A filtered-exact caller should rebuild the same derived + //! expression with `RationalMatrix` and call `det_sign()` when the interval result + //! is inconclusive or reports range failure. Lifting a finished `Matrix` with + //! [`IntervalMatrix::from_matrix`](crate::IntervalMatrix::from_matrix) encloses its stored entries, but cannot recover + //! rounding that occurred while those entries were assembled. + } - #[cfg(feature = "exact")] - /// ```rust - /// use la_stack::prelude::*; - /// - /// # fn main() -> Result<(), LaError> { - /// // Exact determinant - /// let m = Matrix::<3>::try_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(), DeterminantSign::Zero); // exactly singular - /// - /// let det = m.det_exact()?; - /// assert_eq!(det, BigRational::from_integer(0.into())); // exact zero - /// let det_f64 = det.try_to_f64()?; - /// assert_eq!(det_f64, 0.0); - /// - /// // If strict exact-to-f64 conversion would require rounding, opt in - /// // explicitly with the rounded API. - /// let inexact = Matrix::<2>::try_from_rows([ - /// [1.0 + f64::EPSILON, 0.0], - /// [0.0, 1.0 - f64::EPSILON], - /// ])?; - /// let exact_det = inexact.det_exact()?; - /// let rounded_det = match exact_det.try_to_f64() { - /// Ok(det) => det, - /// Err(err) if err.requires_rounding() => exact_det.to_rounded_f64()?, - /// Err(err) => return Err(err), - /// }; - /// assert_eq!(rounded_det.to_bits(), 1.0f64.to_bits()); - /// - /// // If the exact determinant cannot fit in f64, keep the BigRational value. - /// let big = f64::MAX / 2.0; - /// let huge = Matrix::<3>::try_from_rows([ - /// [0.0, 0.0, 1.0], - /// [big, 0.0, 1.0], - /// [0.0, big, 1.0], - /// ])?; - /// let huge_det = huge.det_exact()?; - /// assert_eq!( - /// huge_det - /// .try_to_f64() - /// .err() - /// .and_then(|err| err.unrepresentable_reason()), - /// Some(UnrepresentableReason::NotFinite) - /// ); - /// println!("exact determinant = {huge_det}"); - /// - /// // Exact linear system solve - /// let a = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?; - /// let b = Vector::<2>::try_new([5.0, 11.0])?; - /// let exact_x = a.solve_exact(b)?; - /// let x = exact_x.try_to_f64()?.into_array(); - /// assert!((x[0] - 1.0).abs() <= f64::EPSILON); - /// assert!((x[1] - 2.0).abs() <= f64::EPSILON); - /// # Ok(()) - /// # } - /// ``` - fn exact_arithmetic_example() {} + pub mod norms { + //! Overflow-safe Euclidean norms and squared norms. + //! + //! `Vector::norm()` computes the Euclidean norm with a deterministic scaled + //! sum-of-squares recurrence, so large or subnormal finite coordinates do not fail + //! merely because their raw squares overflow or underflow. It returns positive zero + //! for empty and all-zero vectors and reports `LaError::NonFinite` with + //! `ArithmeticOperation::VectorNorm` only when the exact norm rounds to infinity. + //! Near the upper range, a fixed-size stack accumulator sums squares exactly and + //! compares squared rounding midpoints to prevent false or hidden overflow. This + //! fallback needs no optional dependencies. The general binary64 result remains + //! approximate and has no certified error bound. + //! + //! `Vector::norm_squared()` remains the direct left-to-right FMA sum of squares for + //! callers that need the squared norm. Its distinct contract deliberately reports + //! overflow when that square is not finite, even when `norm()` can return a finite + //! norm. + //! + //! # Large and subnormal coordinates + //! + //! ```rust + //! use core::assert_matches; + //! + //! use la_stack::prelude::*; + //! + //! # fn main() -> Result<(), LaError> { + //! let large = Vector::<5>::try_new([3e200, 4e200, 0.0, 0.0, 0.0])?; + //! assert!((large.norm()? / 1e200 - 5.0).abs() <= 1e-12); + //! assert_matches!( + //! large.norm_squared(), + //! Err(LaError::NonFinite { + //! origin: NonFiniteOrigin::Computation { + //! operation: ArithmeticOperation::VectorSquaredNorm, .. + //! }, + //! location: NonFiniteLocation::Step { index: 0, .. }, + //! .. + //! }) + //! ); + //! + //! let tiny = f64::from_bits(16); + //! let small = Vector::<5>::try_new([3.0 * tiny, 4.0 * tiny, 0.0, 0.0, 0.0])?; + //! assert_eq!(small.norm()?, 5.0 * tiny); + //! assert_eq!(small.norm_squared()?, 0.0); // The raw squares underflow. + //! # Ok(()) + //! # } + //! ``` + //! + //! The large-vector assertion uses a tolerance for this fixture, not a certified + //! error bound. The small vector shows why taking the square root of + //! `norm_squared()` can lose a representable nonzero norm. + } - #[cfg(feature = "exact")] - /// ```rust - /// use core::assert_matches; - /// - /// use la_stack::prelude::*; - /// - /// # fn main() -> Result<(), LaError> { - /// let epsilon = BigRational::new(1.into(), (1_u64 << 60).into()); - /// let one = BigRational::from_integer(1.into()); - /// let zero = BigRational::from_integer(0.into()); - /// - /// let matrix = RationalMatrix::<5>::try_from_fn(|row, col| match (row, col) { - /// (0, 0 | 1) | (1, 0) => one.clone(), - /// (1, 1) => &one + &epsilon, - /// _ if row == col => one.clone(), - /// _ => zero.clone(), - /// })?; - /// assert_eq!(matrix.det_sign(), DeterminantSign::Positive); - /// assert_eq!(matrix.det(), epsilon); - /// - /// let rhs = RationalVector::try_new([ - /// zero, - /// -&epsilon, - /// BigRational::from_integer(2.into()), - /// BigRational::from_integer(3.into()), - /// BigRational::from_integer(4.into()), - /// ])?; - /// let exact_solution = matrix.solve(&rhs)?; - /// assert_eq!( - /// exact_solution.as_array(), - /// &[ - /// BigRational::from_integer(1.into()), - /// BigRational::from_integer((-1).into()), - /// BigRational::from_integer(2.into()), - /// BigRational::from_integer(3.into()), - /// BigRational::from_integer(4.into()), - /// ] - /// ); - /// - /// let epsilon_f64 = epsilon.try_to_f64()?; - /// assert_eq!((1.0 + epsilon_f64).to_bits(), 1.0_f64.to_bits()); - /// let f64_matrix = Matrix::<5>::try_from_rows([ - /// [1.0, 1.0, 0.0, 0.0, 0.0], - /// [1.0, 1.0 + epsilon_f64, 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 f64_rhs = Vector::<5>::try_new([0.0, -epsilon_f64, 2.0, 3.0, 4.0])?; - /// let f64_solve = f64_matrix - /// .lu(DEFAULT_SINGULAR_TOL) - /// .and_then(|lu| lu.solve(f64_rhs)); - /// assert_matches!( - /// f64_solve, - /// Err(LaError::Singular { .. }) - /// ); - /// # Ok(()) - /// # } - /// ``` - fn rational_input_example() {} + pub mod certified { + //! Certified dot products and affine differences. + //! + //! `Vector::dot_with_errbound()` evaluates the same left-to-right FMA tree as + //! `Vector::dot()` and returns its estimate together with a certified absolute + //! roundoff bound. `Vector::dot_difference_with_errbound()` directly evaluates + //! + //! ```text + //! Σᵢ axis[i] × (left[i] - right[i]) + //! ``` + //! + //! as two FMAs per coordinate. It does not first round `left - right` into a new + //! `Vector`, so the certificate covers the intended expression over the original + //! stored binary64 coordinates. + //! + //! The opaque [`ScalarWithErrorBound`](crate::ScalarWithErrorBound) exposes the estimate, absolute error bound, + //! and finite outward-rounded lower and upper bounds. Those endpoints support + //! positive, negative, and caller-selected threshold proofs: + //! + //! # A certified dot-product sign + //! + //! ```rust + //! use la_stack::prelude::*; + //! + //! # fn main() -> Result<(), LaError> { + //! let axis = Vector::<5>::try_new([2.0, -1.0, 3.0, 1.0, -2.0])?; + //! let point = Vector::<5>::try_new([4.0, 1.0, 2.0, 3.0, 1.0])?; + //! let positive = axis.dot_with_errbound(&point)?.and_then(|value| { + //! if value.lower_bound() > 0.0 { + //! Some(true) + //! } else if value.upper_bound() <= 0.0 { + //! Some(false) + //! } else { + //! None // An enclosure overlapping zero is inconclusive. + //! } + //! }); + //! assert_eq!(positive, Some(true)); + //! # Ok(()) + //! # } + //! ``` + //! + //! # An affine threshold test + //! + //! ```rust + //! use la_stack::prelude::*; + //! + //! fn is_separated( + //! axis: &Vector, + //! left: &Vector, + //! right: &Vector, + //! threshold: f64, + //! ) -> Result, LaError> { + //! let Some(value) = axis.dot_difference_with_errbound(left, right)? else { + //! return Ok(None); + //! }; + //! if value.lower_bound() > threshold { + //! Ok(Some(true)) + //! } else if value.upper_bound() <= threshold { + //! Ok(Some(false)) + //! } else { + //! Ok(None) + //! } + //! } + //! + //! # fn main() -> Result<(), LaError> { + //! let axis = Vector::<2>::try_new([2.0, -1.0])?; + //! let left = Vector::<2>::try_new([4.0, 1.0])?; + //! let right = Vector::<2>::try_new([1.0, 3.0])?; + //! assert_eq!(is_separated(&axis, &left, &right, 1.0)?, Some(true)); + //! # Ok(()) + //! # } + //! ``` + //! + //! An interval that overlaps the threshold is inconclusive, not equal. Likewise, + //! `Ok(None)` means gradual underflow or proof-only range exhaustion prevented a + //! certificate. A filtered-exact caller should rebuild the same dot or affine + //! expression in `BigRational` (available through the `exact` feature) or another + //! exact backend. A `LaError::NonFinite` instead reports that the specified FMA + //! estimate itself overflowed. These certified bounds describe roundoff in a fixed + //! arithmetic tree; they are distinct from user-selected numerical tolerances. + } + + pub mod adaptive { + //! Adaptive determinant filtering with certified bounds. + //! + //! [`det_direct_with_errbound()`](crate::Matrix::det_direct_with_errbound) returns a closed-form determinant together with + //! the conservative absolute error bound used by the fast filter, computed from + //! one call that evaluates the determinant once and computes its matching bound. + //! It returns `None` when a D ≤ 4 computation may be affected by gradual + //! underflow, as well as for unsupported D ≥ 5 dimensions. + //! It returns `LaError::NonFinite` if the determinant or bound computation + //! overflows to NaN or infinity. + //! This method does NOT require the `exact` feature — it uses pure f64 arithmetic + //! and is available by default. Use [`det_errbound()`](crate::Matrix::det_errbound) when only the bound is needed. + //! The paired API enables custom adaptive-precision logic for geometric predicates: + //! + //! ```rust + //! use la_stack::prelude::*; + //! + //! # fn main() -> Result<(), LaError> { + //! let matrix = Matrix::<3>::identity(); + //! let sign = matrix.det_direct_with_errbound()?.and_then(|value| { + //! if value.determinant() > value.absolute_error_bound() { + //! Some(1) + //! } else if -value.determinant() > value.absolute_error_bound() { + //! Some(-1) + //! } else { + //! None // The bound cannot establish a sign. + //! } + //! }); + //! assert_eq!(sign, Some(1)); + //! # Ok(()) + //! # } + //! ``` + //! + //! With the `exact` feature, `Matrix::det_sign_exact()` + //! already handles filtering and exact fallback. The + //! [custom adaptive example](https://docs.rs/la-stack/latest/la_stack/guide/exact/index.html#adaptive-determinant-filtering) + //! shows positive, singular, and overflowing filter cases. It requires `exact`. + //! + //! The error coefficients (`ERR_COEFF_2`, `ERR_COEFF_3`, `ERR_COEFF_4`) are + //! conservative, dimension-specific constants, not caller-tunable tolerances. The + //! [mathematical basis](https://github.com/acgetchell/la-stack/blob/v0.4.5/docs/mathematical_basis.md#determinants-and-certified-sign-filtering) + //! documents the bound and states its range preconditions. The constants are explicit + //! crate-root exports for advanced users who want to compose the same bound: + //! `use la_stack::{ERR_COEFF_2, ERR_COEFF_3, ERR_COEFF_4};`. They intentionally stay + //! out of the common prelude. + } #[cfg(feature = "exact")] - /// ```rust - /// use la_stack::prelude::*; - /// - /// fn adaptive_det_sign( - /// matrix: &Matrix, - /// ) -> DeterminantSign { - /// if let Ok(Some(estimate)) = matrix.det_direct_with_errbound() { - /// if estimate.determinant().abs() > estimate.absolute_error_bound() { - /// return if estimate.determinant() > 0.0 { - /// DeterminantSign::Positive - /// } else { - /// DeterminantSign::Negative - /// }; - /// } - /// } - /// - /// matrix.det_sign_exact() - /// } - /// - /// # fn main() -> Result<(), LaError> { - /// let identity = Matrix::<3>::identity(); - /// assert_eq!( - /// adaptive_det_sign(&identity), - /// DeterminantSign::Positive - /// ); - /// - /// let singular = Matrix::<3>::try_from_rows([ - /// [1.0, 2.0, 3.0], - /// [4.0, 5.0, 6.0], - /// [7.0, 8.0, 9.0], - /// ])?; - /// assert_eq!(adaptive_det_sign(&singular), DeterminantSign::Zero); - /// - /// let big = f64::MAX / 2.0; - /// let overflowing = Matrix::<3>::try_from_rows([ - /// [0.0, 0.0, 1.0], - /// [big, 0.0, 1.0], - /// [0.0, big, 1.0], - /// ])?; - /// assert_eq!( - /// adaptive_det_sign(&overflowing), - /// DeterminantSign::Positive - /// ); - /// # Ok(()) - /// # } - /// ``` - fn adaptive_precision_example() {} + pub mod exact { + //! Exact arithmetic over stored binary64 and rational inputs. + //! + //! The default build has **zero runtime dependencies**. Enable the optional + //! `exact` Cargo feature to add exact arithmetic methods using arbitrary-precision + //! rationals (this pulls in `num-bigint`, `num-rational`, and `num-traits` for + //! `BigRational`): + //! + //! See the [crate-level installation instructions](crate) for Cargo configuration. + //! + //! The feature exposes two deliberate input domains: + //! + //! - `Matrix` / `Vector` store finite binary64 inputs. Their exact methods + //! treat each stored bit pattern as its exact rational value, so the determinant + //! or solve stage introduces no further roundoff. They cannot recover information + //! already lost before construction. + //! - `RationalMatrix` / `RationalVector` accept coefficients already + //! assembled as `BigRational`. They preserve derived differences, squared norms, + //! affine coefficients, and other rational expressions without an intermediate + //! `f64` conversion. + //! + //! **Determinants:** + //! + //! - **`det_exact()`** — returns the exact determinant as a `BigRational` + //! - **`det_exact_f64()`** — returns the exact determinant as `f64` only when + //! it is exactly representable (or `LaError::Unrepresentable` otherwise) + //! - **`det_exact_rounded_f64()`** — returns the exact determinant rounded to a + //! finite `f64` using IEEE 754 round-to-nearest, ties-to-even + //! - **`det_sign_exact()`** — infallibly returns the provably correct + //! `DeterminantSign` variant (`Negative`, `Zero`, or `Positive`) + //! + //! **Linear system solve:** + //! + //! - **`solve_exact(b)`** — solves `Ax = b` exactly, returning a + //! `RationalVector` + //! - **`solve_exact_f64(b)`** — solves `Ax = b` exactly, returning `Vector` only when + //! every component is exactly representable as `f64` + //! - **`solve_exact_rounded_f64(b)`** — solves `Ax = b` exactly, returning each + //! component rounded to finite `f64` using IEEE 754 round-to-nearest, + //! ties-to-even + //! - **`ExactF64Conversion`** — converts an existing exact determinant or solution + //! under the strict or rounded contract without repeating exact elimination + //! + //! **Already-exact rational input:** + //! + //! - **`RationalMatrix::det_sign()`** — returns the exact sign without constructing + //! a rational determinant + //! - **`RationalMatrix::det()`** — returns the exact `BigRational` determinant + //! - **`RationalMatrix::solve(&rhs)`** — returns a `RationalVector` exact + //! solution + //! - **`try_with_rational_matrix!`** — dispatches a runtime-selected dimension + //! through D=8 to a const-generic rational matrix on stable Rust + //! + //! The `Matrix::det_exact*` value and conversion methods return + //! `LaError::DeterminantScaleOverflow` if their aggregate power-of-two scaling + //! exceeds the internal exponent representation. `RationalMatrix::det()` is + //! infallible because it clears rational row denominators without an exponent-scale + //! conversion. The exact solve methods for both input domains return + //! `LaError::Singular` with `SingularityReason::Exact` when the stored matrix is + //! exactly singular. + //! + //! For exact-to-f64 output, strict conversions use + //! `UnrepresentableReason::RequiresRounding` when explicit rounding can produce a + //! finite value and `UnrepresentableReason::NotFinite` otherwise. Rounded + //! conversions opt into nearest-even rounding but still report `NotFinite` when no + //! finite `f64` exists. + //! + //! # Preserving rational inputs + //! + //! The following 5×5 system has exact determinant 2^-60. Its exact rational inputs + //! therefore produce a unique solution through the general Bareiss path. Supplying + //! the same coefficients as `f64` inputs loses the `2^-60` perturbation at `1.0`, + //! making the leading rows identical and the binary64 system singular. + //! + //! ```rust + //! use core::assert_matches; + //! + //! use la_stack::prelude::*; + //! + //! fn main() -> Result<(), LaError> { + //! // This is far below one binary64 ULP at 1.0, so 1.0 + 2^-60 rounds to 1.0. + //! let epsilon = BigRational::new(1.into(), (1_u64 << 60).into()); + //! let one = BigRational::from_integer(1.into()); + //! let zero = BigRational::from_integer(0.into()); + //! + //! // The leading block is [[1, 1], [1, 1 + 2^-60]]. The remaining diagonal + //! // extends the example to D=5, where the general Bareiss path is used. + //! let matrix = RationalMatrix::<5>::try_from_fn(|row, col| match (row, col) { + //! (0, 0 | 1) | (1, 0) => one.clone(), + //! (1, 1) => &one + &epsilon, + //! _ if row == col => one.clone(), + //! _ => zero.clone(), + //! })?; + //! assert_eq!(matrix.det_sign(), DeterminantSign::Positive); + //! assert_eq!(matrix.det(), epsilon); + //! + //! let rhs = RationalVector::try_new([ + //! zero, + //! -&epsilon, + //! BigRational::from_integer(2.into()), + //! BigRational::from_integer(3.into()), + //! BigRational::from_integer(4.into()), + //! ])?; + //! let exact_solution = matrix.solve(&rhs)?; + //! assert_eq!( + //! exact_solution.as_array(), + //! &[ + //! BigRational::from_integer(1.into()), + //! BigRational::from_integer((-1).into()), + //! BigRational::from_integer(2.into()), + //! BigRational::from_integer(3.into()), + //! BigRational::from_integer(4.into()), + //! ] + //! ); + //! + //! // Supplying the same coefficients as f64 inputs destroys the perturbation + //! // and makes the matrix singular, even though the exact solution is integral. + //! let epsilon_f64 = epsilon.try_to_f64()?; + //! assert_eq!((1.0 + epsilon_f64).to_bits(), 1.0_f64.to_bits()); + //! let f64_matrix = Matrix::<5>::try_from_rows([ + //! [1.0, 1.0, 0.0, 0.0, 0.0], + //! [1.0, 1.0 + epsilon_f64, 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 f64_rhs = Vector::<5>::try_new([0.0, -epsilon_f64, 2.0, 3.0, 4.0])?; + //! let f64_solve = f64_matrix + //! .lu(DEFAULT_SINGULAR_TOL) + //! .and_then(|lu| lu.solve(f64_rhs)); + //! assert_matches!( + //! f64_solve, + //! Err(LaError::Singular { .. }) + //! ); + //! Ok(()) + //! } + //! ``` + //! + //! # Stored binary64 inputs and output conversion + //! + //! ```rust + //! use la_stack::prelude::*; + //! + //! fn main() -> Result<(), LaError> { + //! // Exact determinant + //! let m = Matrix::<3>::try_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(), DeterminantSign::Zero); // exactly singular + //! + //! let det = m.det_exact()?; + //! assert_eq!(det, BigRational::from_integer(0.into())); // exact zero + //! let det_f64 = det.try_to_f64()?; + //! assert_eq!(det_f64, 0.0); + //! + //! // If strict exact-to-f64 conversion would require rounding, opt in + //! // explicitly with the rounded API. + //! let inexact = Matrix::<2>::try_from_rows([ + //! [1.0 + f64::EPSILON, 0.0], + //! [0.0, 1.0 - f64::EPSILON], + //! ])?; + //! let exact_det = inexact.det_exact()?; + //! let rounded_det = match exact_det.try_to_f64() { + //! Ok(det) => det, + //! Err(err) if err.requires_rounding() => exact_det.to_rounded_f64()?, + //! Err(err) => return Err(err), + //! }; + //! assert_eq!(rounded_det.to_bits(), 1.0f64.to_bits()); + //! + //! // If the exact determinant cannot fit in f64, keep the BigRational value. + //! let big = f64::MAX / 2.0; + //! let huge = Matrix::<3>::try_from_rows([ + //! [0.0, 0.0, 1.0], + //! [big, 0.0, 1.0], + //! [0.0, big, 1.0], + //! ])?; + //! let huge_det = huge.det_exact()?; + //! assert_eq!( + //! huge_det + //! .try_to_f64() + //! .err() + //! .and_then(|err| err.unrepresentable_reason()), + //! Some(UnrepresentableReason::NotFinite) + //! ); + //! println!("exact determinant = {huge_det}"); + //! + //! // Exact linear system solve + //! let a = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?; + //! let b = Vector::<2>::try_new([5.0, 11.0])?; + //! let exact_x = a.solve_exact(b)?; + //! let x = exact_x.try_to_f64()?.into_array(); + //! assert!((x[0] - 1.0).abs() <= f64::EPSILON); + //! assert!((x[1] - 2.0).abs() <= f64::EPSILON); + //! + //! Ok(()) + //! } + //! ``` + //! + //! With the `exact` feature enabled, `RationalMatrix`, `RationalVector`, + //! `DeterminantSign`, `ExactF64Conversion`, `BigInt`, and `BigRational` are + //! re-exported from the crate root and prelude, + //! alongside the most commonly needed `num-traits` items (`FromPrimitive`, + //! `ToPrimitive`, `Signed`). This lets consumers construct exact values + //! (`BigRational::from_f64`, `from_i64`), query sign (`is_positive` / + //! `is_negative`), and convert back to `f64` (`try_to_f64`, `to_rounded_f64`, or + //! the raw `to_f64`) with a single + //! `use la_stack::prelude::*;` — no need to add `num-bigint`, `num-rational`, + //! or `num-traits` to their own `Cargo.toml`. Use + //! `DeterminantSign::as_i8()` only when numeric −1/0/+1 interoperability is + //! required. + //! + //! For `det_sign_exact()`, D ≤ 4 matrices first use a fast f64 filter + //! (error-bounded [`det_direct_with_errbound()`](crate::Matrix::det_direct_with_errbound)) when its rounded intermediates stay in the normal + //! range or are exact structural zeros. An inconclusive filter falls back to the + //! same direct determinant expansion in `BigInt`. D ≥ 5 skips the closed-form + //! filter and uses fraction-free Bareiss elimination in `BigInt`. + //! Because `Matrix` stores only finite entries, arithmetic range failures in the + //! filter are inconclusive rather than errors and the exact fallback is total. + //! + //! # A five-dimensional rational solve + //! + //! ```rust + //! use core::assert_matches; + //! + //! use la_stack::prelude::*; + //! + //! # fn main() -> Result<(), LaError> { + //! // A tridiagonal exact matrix with determinant 6. + //! let matrix = RationalMatrix::<5>::try_from_fn(|row, col| { + //! BigRational::from_integer(if row == col { + //! 2.into() + //! } else if row.abs_diff(col) == 1 { + //! 1.into() + //! } else { + //! 0.into() + //! }) + //! })?; + //! let numerators = [4, 8, 12, 16, 14]; + //! let rhs = RationalVector::try_from_fn(|row| { + //! BigRational::new(numerators[row].into(), 3.into()) + //! })?; + //! let solution = matrix.solve(&rhs)?; + //! let expected = [1, 2, 3, 4, 5].map(|n| BigRational::new(n.into(), 3.into())); + //! assert_eq!(solution.as_array(), &expected); + //! + //! // Keep the exact solution until the caller explicitly opts into rounding. + //! assert_matches!( + //! solution.try_to_f64(), + //! Err(LaError::Unrepresentable { + //! index: Some(0), + //! reason: UnrepresentableReason::RequiresRounding, + //! .. + //! }) + //! ); + //! let rounded = solution.to_rounded_f64()?; + //! assert_eq!(rounded.as_array(), &[1.0 / 3.0, 2.0 / 3.0, 1.0, 4.0 / 3.0, 5.0 / 3.0]); + //! # Ok(()) + //! # } + //! ``` + //! + //! # Rational dimension dispatch + //! + //! [`try_with_rational_matrix!`](crate::try_with_rational_matrix) selects a + //! concrete `RationalMatrix` for dimensions 0 through + //! [`MAX_RATIONAL_MATRIX_DISPATCH_DIM`](crate::MAX_RATIONAL_MATRIX_DISPATCH_DIM) + //! (8). Larger dimensions return [`LaError::UnsupportedDimension`](crate::LaError::UnsupportedDimension), + //! converted through `From` into the closure's declared error type. + //! The macro preserves the const-generic representation. + //! + //! Rational constructors include `try_from_rows` / `try_new` and + //! `try_from_fn`; `as_rows` / `as_array` and `get` borrow their stored + //! values, while `into_rows` / `into_array` return the owned arrays. + //! + //! # Adaptive determinant filtering + //! + //! This example requires `exact` and illustrates a custom filter with exact + //! fallback. Use [`Matrix::det_sign_exact`](crate::Matrix::det_sign_exact) directly + //! when no custom filtering policy is needed. + //! + //! ```rust + //! use la_stack::prelude::*; + //! + //! fn adaptive_det_sign( + //! matrix: &Matrix, + //! ) -> DeterminantSign { + //! if let Ok(Some(estimate)) = matrix.det_direct_with_errbound() { + //! if estimate.determinant().abs() > estimate.absolute_error_bound() { + //! return if estimate.determinant() > 0.0 { + //! DeterminantSign::Positive + //! } else { + //! DeterminantSign::Negative + //! }; + //! } + //! } + //! + //! matrix.det_sign_exact() + //! } + //! + //! fn main() -> Result<(), LaError> { + //! let identity = Matrix::<3>::identity(); + //! assert_eq!( + //! adaptive_det_sign(&identity), + //! DeterminantSign::Positive + //! ); + //! + //! // A zero determinant cannot pass the f64 sign filter, so this exercises + //! // the exact fallback. + //! let singular = Matrix::<3>::try_from_rows([ + //! [1.0, 2.0, 3.0], + //! [4.0, 5.0, 6.0], + //! [7.0, 8.0, 9.0], + //! ])?; + //! assert_eq!(adaptive_det_sign(&singular), DeterminantSign::Zero); + //! + //! // The f64 filter overflows for this finite matrix, but the exact fallback + //! // still resolves its positive determinant sign. + //! let big = f64::MAX / 2.0; + //! let overflowing = Matrix::<3>::try_from_rows([ + //! [0.0, 0.0, 1.0], + //! [big, 0.0, 1.0], + //! [0.0, big, 1.0], + //! ])?; + //! assert_eq!( + //! adaptive_det_sign(&overflowing), + //! DeterminantSign::Positive + //! ); + //! + //! Ok(()) + //! } + //! ``` + } } - mod error; #[cfg(feature = "exact")] mod exact;