diff --git a/CHANGELOG.md b/CHANGELOG.md index 34d851e..e826f64 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### ⚠️ Breaking Changes + +- Make exact f64 conversions strict + ### Added - Guard README dependency snippets [`7137fee`](https://github.com/acgetchell/la-stack/commit/7137fee16ab33e08f4dc6a60e02417e3e7c4e020) @@ -17,6 +21,40 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Refresh README determinant examples with explicit fallible handling and hidden doctest mirrors - Update CI uv pins to 0.11.19 +- Report determinant scale overflow precisely [`928f62b`](https://github.com/acgetchell/la-stack/commit/928f62bba0d837afe04cb8ccb3fbfed6b095d8f7) + - Add a typed LaError::DeterminantScaleOverflow path for exact determinant scale exponent failures + - Convert det_exact_f64 directly from the shared Bareiss integer/exponent pair while preserving Overflow for finite-f64 conversion failures + - Reuse vector finiteness scanning across raw and proof-bearing constructors + - Harden docs version sync checks for reordered inline-table dependency snippets and pruned Markdown traversal +- Add release performance comparison workflow [`53b5fde`](https://github.com/acgetchell/la-stack/commit/53b5fde13f3e4afbd3db80d324185a091203cb75) + - Extend vs_linalg with LDLT/Cholesky benchmark rows and shared deterministic inputs. + - Add smoke coverage that checks la-stack, nalgebra, and faer agree on benchmark inputs. + - Expand bench-compare to support latest-vs-last reports, suite/scope selection, peer baseline context, and clearer malformed Criterion diagnostics. + - Document the benchmark methodology, release baseline workflow, roadmap direction, and contributor guidance. +- Publish release benchmark baselines [`9497ca5`](https://github.com/acgetchell/la-stack/commit/9497ca5f88dd7800bdf3823123e2a295fcb5ced1) + - Add a release-only benchmark workflow that saves full Criterion baselines for published releases and attaches the archived baseline to the GitHub Release. + - Keep the regular benchmark workflow focused on PR and main-branch comparison runs. + - Document how to restore archived release baselines for future performance comparisons. +- Feat!(api): make Matrix and Vector finite by construction [`1fa2f55`](https://github.com/acgetchell/la-stack/commit/1fa2f55cfac6f249a7e2bf30922901539e580dd8) +- [**breaking**] Make exact f64 conversions strict [`8e33f1a`](https://github.com/acgetchell/la-stack/commit/8e33f1a8ec291bfcb6312375969efce076421e96) + - Add explicit rounded exact-to-f64 APIs for determinant and solve results + - Report exact conversion failures with typed Unrepresentable reasons + - Remove finite proof wrapper APIs now that Matrix and Vector carry finiteness directly + - Move error and tolerance contracts into first-class modules with prelude exports + - Update exact benchmarks to distinguish strict Result paths from rounded f64 paths + - Document and exercise the rounded fallback pattern for RequiresRounding errors + +### Changed + +- Cover determinant scale overflow boundaries [`532093a`](https://github.com/acgetchell/la-stack/commit/532093a1ed9f65ace159f80aac290907d570ea8a) + + - Extract determinant scale exponent calculation into a private helper + - Assert typed DeterminantScaleOverflow errors for dimension conversion and exponent product overflow +- Harden support script parsing [`87e1d00`](https://github.com/acgetchell/la-stack/commit/87e1d0042ad2de7888c1a065ab78524a63f4c045) + - Require Python 3.13 for support-script tooling and align Ruff/Ty with that baseline. + - Replace mypy with strict Ty checking in the Python workflow. + - Parse TOML, JSON, argparse, and Semgrep inputs into typed boundary objects before downstream use. + - Reject malformed Criterion estimates, non-finite timings, invalid confidence intervals, and malformed Semgrep result shapes. ### Documentation @@ -27,6 +65,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Add citation metadata validation to the release checklist and config lint flow. - Include CITATION.cff in YAML/CFF formatting checks. +### Fixed + +- Escape path regex in benchmark parser test [`1222c93`](https://github.com/acgetchell/la-stack/commit/1222c9325c7b7ea0e1a174e4ee6f0f08e1f7e94b) + + Use a literal regex pattern for the malformed Criterion JSON diagnostic so + Windows paths with backslashes do not break pytest's match expression. +- Align ty with Python 3.13 [`b9e0ba0`](https://github.com/acgetchell/la-stack/commit/b9e0ba08e54a15d8eddd5c5c53edc37bbc03939a) + ## [0.4.2] - 2026-06-04 ### Added @@ -378,7 +424,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Bump rust-version to 1.95 in Cargo.toml - Bump channel to 1.95.0 in rust-toolchain.toml - Add core::hint::cold_path() hints at cold/error branches: - - src/exact.rs: validate_finite, validate_finite_vec, gauss_solve singular return, det_exact_f64 / solve_exact_f64 overflow returns, det_sign_exact Stage 2 Bareiss fallback diff --git a/README.md b/README.md index ed3931b..7c92b06 100644 --- a/README.md +++ b/README.md @@ -34,7 +34,7 @@ while keeping the API intentionally small and explicit. - ✅ `const fn` where possible (compile-time evaluation of determinants, dot products, etc.) - ✅ Explicit algorithms (LU, solve, determinant) - ✅ Robust geometric predicates via optional exact arithmetic (`det_sign_exact`, `det_errbound`) -- ✅ Exact linear system solve via optional arbitrary-precision arithmetic (`solve_exact`, `solve_exact_f64`) +- ✅ Exact linear system solve via optional arbitrary-precision arithmetic (`solve_exact`, strict/rounded f64 conversions) - ✅ No runtime dependencies by default (optional features may add deps) - ✅ Stack storage only (no heap allocation in core types) - ✅ `unsafe` forbidden @@ -213,14 +213,19 @@ la-stack = { version = "0.4.2", features = ["exact"] } **Determinants:** - **`det_exact()`** — returns the exact determinant as a `BigRational` -- **`det_exact_f64()`** — returns the exact determinant converted to the nearest `f64` - (or `LaError::Overflow` when the exact value is unrepresentable) +- **`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` - **`det_sign_exact()`** — returns the provably correct sign (−1, 0, or +1) **Linear system solve:** - **`solve_exact(b)`** — solves `Ax = b` exactly, returning `[BigRational; D]` -- **`solve_exact_f64(b)`** — solves `Ax = b` exactly, converting the result to `Vector` (f64) +- **`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` ```rust,ignore use la_stack::prelude::*; @@ -239,6 +244,19 @@ fn main() -> Result<(), LaError> { let det_f64 = m.det_exact_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 rounded_det = match inexact.det_exact_f64() { + Ok(det) => det, + Err(err) if err.requires_rounding() => inexact.det_exact_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([ @@ -247,7 +265,12 @@ fn main() -> Result<(), LaError> { [0.0, big, 1.0], ])?; let huge_det = huge.det_exact()?; - assert_eq!(huge.det_exact_f64(), Err(LaError::Overflow { index: None })); + assert_eq!( + huge.det_exact_f64() + .err() + .and_then(|err| err.unrepresentable_reason()), + Some(UnrepresentableReason::NotFinite) + ); println!("exact determinant = {huge_det}"); // Exact linear system solve @@ -338,7 +361,8 @@ compose the same bound themselves. Storage shown above reflects the intentional `f64` scalar model. `Matrix` key methods: `lu`, `ldlt`, `det`, `det_direct`, `det_errbound`, -`det_exact`¹, `det_exact_f64`¹, `det_sign_exact`¹, `solve_exact`¹, `solve_exact_f64`¹. +`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 kernels do not revalidate stored entries. diff --git a/benches/exact.rs b/benches/exact.rs index 5348f51..6683333 100644 --- a/benches/exact.rs +++ b/benches/exact.rs @@ -19,6 +19,10 @@ //! a fixed-seed corpus of diagonally-dominant random matrices per //! dimension. Each operation is pre-timed across the corpus to select //! p50/p95/p99 cumulative input subsets, then measured with Criterion. +//! +//! Fallible exact-to-f64 conversions use a `_result` suffix. Those rows measure +//! the full `Result` path, including valid `Err(Unrepresentable)` outcomes for +//! inputs whose exact answer cannot be represented as finite binary64. use std::array; use std::cell::Cell; @@ -202,7 +206,8 @@ enum ExactRandomOperation { DetSignExact, DetExact, SolveExact, - SolveExactF64, + SolveExactF64Result, + SolveExactRoundedF64, } impl ExactRandomOperation { @@ -212,7 +217,8 @@ impl ExactRandomOperation { Self::DetSignExact => "det_sign_exact", Self::DetExact => "det_exact", Self::SolveExact => "solve_exact", - Self::SolveExactF64 => "solve_exact_f64", + Self::SolveExactF64Result => "solve_exact_f64_result", + Self::SolveExactRoundedF64 => "solve_exact_rounded_f64", } } } @@ -323,10 +329,14 @@ fn run_random_operation( ); let _ = black_box(x); } - ExactRandomOperation::SolveExactF64 => { + ExactRandomOperation::SolveExactF64Result => { + let x = black_box(input.matrix).solve_exact_f64(black_box(input.rhs)); + let _ = black_box(x); + } + ExactRandomOperation::SolveExactRoundedF64 => { let x = require_ok( - black_box(input.matrix).solve_exact_f64(black_box(input.rhs)), - "exact linear solve converted to f64", + black_box(input.matrix).solve_exact_rounded_f64(black_box(input.rhs)), + "exact linear solve rounded to f64", ); let _ = black_box(x); } @@ -488,9 +498,10 @@ fn hilbert() -> Matrix { ) } -/// Populate a Criterion group with the four headline exact-arithmetic +/// Populate a Criterion group with the five headline exact-arithmetic /// benches on a single `(matrix, rhs)` pair: `det_sign_exact`, -/// `det_exact`, `solve_exact`, `solve_exact_f64`. +/// `det_exact`, `solve_exact`, `solve_exact_f64_result`, and +/// `solve_exact_rounded_f64`. /// /// Used by every adversarial-input group so each one measures the same /// operations, making the resulting tables directly comparable. @@ -523,11 +534,18 @@ fn bench_extreme_group( }); }); - group.bench_function("solve_exact_f64", |bencher| { + group.bench_function("solve_exact_f64_result", |bencher| { + bencher.iter(|| { + let x = black_box(m).solve_exact_f64(black_box(rhs)); + let _ = black_box(x); + }); + }); + + group.bench_function("solve_exact_rounded_f64", |bencher| { bencher.iter(|| { let x = require_ok( - black_box(m).solve_exact_f64(black_box(rhs)), - "exact linear solve converted to f64", + black_box(m).solve_exact_rounded_f64(black_box(rhs)), + "exact linear solve rounded to f64", ); let _ = black_box(x); }); @@ -571,12 +589,20 @@ macro_rules! gen_exact_benches_for_dim { }); }); - // === det_exact_f64 (exact → f64) === - [].bench_function("det_exact_f64", |bencher| { + // === det_exact_f64 (fallible exact → f64 Result) === + [].bench_function("det_exact_f64_result", |bencher| { + bencher.iter(|| { + let det = black_box(a).det_exact_f64(); + black_box(det); + }); + }); + + // === det_exact_rounded_f64 (lossy exact → f64) === + [].bench_function("det_exact_rounded_f64", |bencher| { bencher.iter(|| { let det = require_ok( - black_box(a).det_exact_f64(), - "exact determinant converted to f64", + black_box(a).det_exact_rounded_f64(), + "exact determinant rounded to f64", ); black_box(det); }); @@ -601,12 +627,20 @@ macro_rules! gen_exact_benches_for_dim { }); }); - // === solve_exact_f64 (exact → f64) === - [].bench_function("solve_exact_f64", |bencher| { + // === solve_exact_f64 (fallible exact → f64 Result) === + [].bench_function("solve_exact_f64_result", |bencher| { + bencher.iter(|| { + let x = black_box(a).solve_exact_f64(black_box(rhs)); + black_box(x); + }); + }); + + // === solve_exact_rounded_f64 (lossy exact → f64) === + [].bench_function("solve_exact_rounded_f64", |bencher| { bencher.iter(|| { let x = require_ok( - black_box(a).solve_exact_f64(black_box(rhs)), - "exact linear solve converted to f64", + black_box(a).solve_exact_rounded_f64(black_box(rhs)), + "exact linear solve rounded to f64", ); black_box(x); }); @@ -642,7 +676,12 @@ macro_rules! gen_random_percentile_benches_for_dim { bench_random_percentile_operation( &mut [], &corpus, - ExactRandomOperation::SolveExactF64, + ExactRandomOperation::SolveExactF64Result, + ); + bench_random_percentile_operation( + &mut [], + &corpus, + ExactRandomOperation::SolveExactRoundedF64, ); [].finish(); @@ -677,8 +716,9 @@ fn main() { // === Adversarial / extreme-input groups === // - // Each group runs the same four exact-arithmetic benches - // (`det_sign_exact`, `det_exact`, `solve_exact`, `solve_exact_f64`) + // Each group runs the same five exact-arithmetic benches + // (`det_sign_exact`, `det_exact`, `solve_exact`, `solve_exact_f64_result`, + // `solve_exact_rounded_f64`) // via `bench_extreme_group`, so the resulting tables are directly // comparable across input classes. diff --git a/docs/BENCHMARKING.md b/docs/BENCHMARKING.md index 857b528..dfd27d7 100644 --- a/docs/BENCHMARKING.md +++ b/docs/BENCHMARKING.md @@ -15,7 +15,8 @@ la-stack has two Criterion benchmark suites: dependency version used here. - **`exact`** (`benches/exact.rs`) — measures exact-arithmetic methods - (`det_exact`, `solve_exact`, `det_sign_exact`, etc.) alongside f64 + (`det_exact`, `solve_exact`, `det_sign_exact`, strict `*_result` + conversions, and lossy `*_rounded_f64` conversions) alongside f64 baselines (`det`, `det_direct`) across D=2–5. Use this to understand the cost of exact arithmetic and track optimization progress. In addition to the fixed per-dimension groups (`exact_d{2..5}`), the @@ -35,10 +36,14 @@ la-stack has two Criterion benchmark suites: ill-conditioned matrices whose non-terminating-in-binary entries stress the `f64_decompose → BigInt` scaling path. - Each random percentile and adversarial group runs the same four + Each random percentile and adversarial group runs the same five exact-arithmetic benches (`det_sign_exact`, `det_exact`, `solve_exact`, - `solve_exact_f64`) so the resulting tables are directly comparable - across input classes. + `solve_exact_f64_result`, `solve_exact_rounded_f64`) so the resulting tables + are directly comparable across input classes. Rows with a `_result` suffix + measure the strict fallible conversion path, including valid + `Err(Unrepresentable)` outcomes when the exact answer is not + finite-binary64 representable. Rows with a `_rounded_f64` suffix measure the + intentionally lossy finite-binary64 conversion path. ## `vs_linalg` methodology @@ -188,6 +193,23 @@ local. The report includes per-dimension tables showing median times, percent change, speedup, and last-release nalgebra/faer context where a matching `vs_linalg` peer exists. +For exact-arithmetic comparisons against v0.4.2 or older baselines, rows such +as `det_exact_rounded_f64 (vs det_exact_f64)` mean the current rounded API is +being compared to the historical lossy `*_exact_f64` benchmark. Rows such as +`det_exact_f64_result (vs det_exact_f64)` intentionally show the overhead of +the new strict conversion contract against that same historical baseline. + +The default `release-signal` scope reports exact-arithmetic rows whose inputs +are fixed across versions: deterministic D=2..=5 cases plus adversarial fixed +matrices. Random percentile groups are exploratory tail probes; each benchmark +run selects p50/p95/p99 input sets by timing the implementation under test, so +those rows can measure different corpus subsets across versions. Include them +when investigating tails with: + +```bash +uv run bench-compare v0.4.2 --suite exact --scope all-benches +``` + To generate a current snapshot without a saved baseline: ```bash diff --git a/docs/archive/changelog/0.1.md b/docs/archive/changelog/0.1.md index 6545e28..473ce4a 100644 --- a/docs/archive/changelog/0.1.md +++ b/docs/archive/changelog/0.1.md @@ -19,7 +19,6 @@ Renames WARP.md to AGENTS.md for broader applicability. Refreshes dependencies in Cargo.lock and pyproject.toml, including faer. Adds a type checker (ty) for Python scripts. (internal) - - Enables benchmark feature for criterion [`d701b58`](https://github.com/acgetchell/la-stack/commit/d701b586e89dd3583a77cd7c96512a7e6d6343a9) Enables the "bench" feature for criterion and other benchmark @@ -105,7 +104,6 @@ nalgebra, providing a broader performance context. Adds faer implementations and benchmarks. Updates data and plots. - - Refactors criterion plotting script for clarity [`be24d03`](https://github.com/acgetchell/la-stack/commit/be24d033ec052fcd25b9598a067c774036675bc2) Refactors the criterion plotting script to improve readability and @@ -113,7 +111,6 @@ for structured data handling and introduces `ReadmeMarkerError` for better error context when handling README markers. The script now uses named attributes, enhancing code clarity. (internal) - - Updates benchmark data and plot [`f157436`](https://github.com/acgetchell/la-stack/commit/f157436439b684e28e87e7438eec9a0736d94b5b) Updates benchmark data for the linear algebra solver comparison. @@ -129,7 +126,6 @@ division by zero and displaying "n/a" instead. Also pins uv version for reproducible CI and relaxes whitespace check in codacy config. - - Correctly escapes backslashes in gnuplot strings [`732480a`](https://github.com/acgetchell/la-stack/commit/732480afe68a14b36ccb546d34ad969ef5b35f41) Fixes an issue where backslashes were not properly escaped when @@ -137,7 +133,6 @@ The gnuplot quoting function now escapes both backslashes and single quotes to ensure correct string interpretation by gnuplot. - - Improves error handling for README table updates [`cd6873f`](https://github.com/acgetchell/la-stack/commit/cd6873f74aaaff4cd272223ad7f0b41f889f5870) Refactors README table update error handling for better clarity @@ -152,7 +147,6 @@ - Implements CI/CD workflows and code quality tools [`1989c3d`](https://github.com/acgetchell/la-stack/commit/1989c3dc3603b3b7674c978826af2d4318da504a) Sets up CI/CD workflows for testing, linting, audit, and security analysis. - - Initial implementation of `la-stack` crate [`c3ab9c5`](https://github.com/acgetchell/la-stack/commit/c3ab9c59d48463a29992a6f3918d22b037059dfa) Adds the initial implementation of the `la-stack` crate, @@ -160,7 +154,6 @@ Includes `Vector`, `Matrix`, and `Lu` types, along with basic functionality and examples. Also sets up `justfile` for development workflow and a license file. - - Configures CI workflow for cross-platform builds [`ff96949`](https://github.com/acgetchell/la-stack/commit/ff9694913a52d6c54a643e9186a8d7481d4d8eda) Adds a comprehensive CI workflow using GitHub Actions to enable @@ -179,7 +172,6 @@ ### Changed - Initial commit [`e9a3553`](https://github.com/acgetchell/la-stack/commit/e9a3553e823fab736b603c97712f484cfa990f7c) - - Improves linting and validation workflows [`f874ff7`](https://github.com/acgetchell/la-stack/commit/f874ff78552e1f61333e30ca07b20a914ddcc98f) Adds markdown linting with configuration and JSON validation to the @@ -187,7 +179,6 @@ incorporate the new linting and validation steps. Updates cspell configuration to include markdownlint. Improves matrix tests. Updates WARP.md file. - - Overhauls public API tests and examples [`5699ac5`](https://github.com/acgetchell/la-stack/commit/5699ac59cfe923cd4db7f5fca0420c450b5dc27a) Modernizes tests by adding property-based tests and tests for edge @@ -195,7 +186,6 @@ This change focuses on the tests and examples, improving coverage and demonstrating common usage scenarios. - - Improves code quality with Codacy and stricter lints [`6e5b004`](https://github.com/acgetchell/la-stack/commit/6e5b004e79c07ed3ca02e89f0af77f911e791b57) Adds Codacy configuration for static analysis and security scanning, @@ -222,7 +212,6 @@ numerical stability and preventing false negatives. Also, refactors test setup in `lu.rs` for clarity. - - Improves spell check configuration [`0ad77f6`](https://github.com/acgetchell/la-stack/commit/0ad77f6d1383be1e6ddf885d66c1f8d2c5e8c252) Updates the spell check configuration to include gitignore files, diff --git a/docs/archive/changelog/0.2.md b/docs/archive/changelog/0.2.md index df56217..466298d 100644 --- a/docs/archive/changelog/0.2.md +++ b/docs/archive/changelog/0.2.md @@ -9,7 +9,6 @@ Expose the internal `det_errbound()` function and error coefficients as public API to enable consumer-side adaptive-precision logic for geometric predicates. - - Expose det_errbound without requiring exact feature [`d96f676`](https://github.com/acgetchell/la-stack/commit/d96f676bee8436bde622addcf3ac48c24111da75) Make `det_errbound()` and error coefficients available without the @@ -34,7 +33,6 @@ - Add tests for L-multiplier overflow and trailing submatrix overflow during LDLT factorization - - Add test for forward substitution overflow in LDLT solve_vec - Update README det_sign_exact example for Result API (.unwrap()) - Update changelog generation and release metadata (internal) @@ -66,14 +64,12 @@ - Stage 1 (D ≤ 4): f64 fast filter via det_direct() + Shewchuk-style error bounds against the matrix permanent. Resolves sign without allocating for well-conditioned inputs. - - Stage 2: exact Bareiss algorithm (fraction-free Gaussian elimination) in BigRational for guaranteed-correct sign. ### Changed - Bump the dependencies group with 2 updates [`ad4b505`](https://github.com/acgetchell/la-stack/commit/ad4b505ab5a5ec4c2c960a07f4d43553877828cb) - - Rename det benchmark and add black_box to tests [`40cc2b7`](https://github.com/acgetchell/la-stack/commit/40cc2b7bde49cdd0a2c7a247fb045447447d19df) Rename "la_stack_det_direct" benchmark to "la_stack_det" since it @@ -82,14 +78,12 @@ Add black_box to det_direct known-value tests to prevent the compiler from constant-folding the const fn at compile time, ensuring the function body is exercised at runtime for more robust coverage. - - Clarify det_direct D=0 support and fix example description [`1308612`](https://github.com/acgetchell/la-stack/commit/130861273dcd9b459a5c1891656e686a6760756c) Update the README to note that det_direct() supports D=0 using the empty-product convention (returning 1.0). Correct the doc comment in the 4x4 example matrix, which was inaccurately described as being Hilbert-like. These are non-functional documentation refinements. - - Polish exact feature — docs, CI, example, review fixes [`abfa297`](https://github.com/acgetchell/la-stack/commit/abfa297d1e8e771a60d1cc9a2a9227011876a965) - Add exact_sign_3x3 example (near-singular 3×3 with 2^-50 perturbation) @@ -99,10 +93,8 @@ - Add --features exact to tarpaulin coverage runs - Update README: zero-dep-by-default wording, scalar types section, example descriptions, design goals, references wording - - Update Cargo.toml deps comment - Remove redundant clippy-exact from ci target - - Make Codacy analysis non-blocking on pull requests (internal) [`ecd0a30`](https://github.com/acgetchell/la-stack/commit/ecd0a30d1e580ae47f215b811a60110b8f94223a) @@ -110,7 +102,6 @@ requests, mitigating the impact of transient service outages. Ensures that subsequent SARIF processing steps always execute and surfaces a warning when the analysis fails. - - Harden Codacy CI and expand exact feature tests (internal) [`d911fcc`](https://github.com/acgetchell/la-stack/commit/d911fcc51c02807e91372e1db94aaaa196d6ccb9) Update the Codacy workflow to handle missing project tokens gracefully @@ -118,25 +109,21 @@ unit tests for the internal `det_errbound` and `bareiss_det` functions to ensure logic correctness across various matrix dimensions and edge cases like singular matrices. - - Expand det_sign_exact tests for large matrices (internal) [`2205f12`](https://github.com/acgetchell/la-stack/commit/2205f12a94743c49333c6f6a7228edde21a0d74f) Add unit tests for dimensions $D \ge 5$ to exercise the Bareiss algorithm path. Ensure non-finite entries trigger expected panics and verify pivoting logic when the fast filter is bypassed. - - Improve release safety and consolidate error handling [`0ef2b94`](https://github.com/acgetchell/la-stack/commit/0ef2b94798f5fa530a09d92d0d117d42a3bc2b8d) Ensure cargo publish uses the lockfile for reproducible releases. Refactor the tag_release script to use a single exception handler for various process and validation failures. - - Update developer docs and optimize release tagging logic [`8aa129a`](https://github.com/acgetchell/la-stack/commit/8aa129ace22919510022a58def30b85c2639551d) Update AGENTS.md and scripts/README.md with details on feature flags, code structure, and release workflows. Improve the efficiency of changelog section extraction in tag_release.py by using slice-based trimming. - - Improve internal release scripts and expand test coverage [`5eb5064`](https://github.com/acgetchell/la-stack/commit/5eb5064b46b8cbe408ddc958b86bf01dfee1d62c) Refine SemVer regex to support digit-prefixed alphanumeric prerelease @@ -174,10 +161,8 @@ - _github_anchor now matches github-slugger (strips brackets, removes non-alphanumeric punctuation, replaces whitespace with hyphens) - - Tag deletion deferred until after changelog extraction succeeds, preventing orphaned tag loss on validation failure - - Force push instruction emitted when --force flag is used Add docs/RELEASING.md with step-by-step release workflow (version @@ -198,7 +183,6 @@ - Add _ensure-typos helper and typos-cli to setup-tools install/verify - Update CI to install typos-cli via taiki-e/install-action instead of npm cspell - - Add changelog tooling (git-cliff + tag script) [`5475651`](https://github.com/acgetchell/la-stack/commit/547565142a374d6a7fed432dd6dbaa8b24cf98a1) Add git-cliff-based changelog generation and a Python tag-release @@ -206,16 +190,12 @@ used in delaunay. New files: - - cliff.toml: Keep a Changelog template with commit parsers for mixed commit history, PR/commit link rewriting via preprocessors - - scripts/tag_release.py: CLI for annotated tag creation with semver validation, changelog section extraction, 125KB size limit handling - - scripts/subprocess_utils.py: safe subprocess wrappers (ported from delaunay) - - scripts/postprocess_changelog.py: strips trailing blank lines from generated output diff --git a/docs/archive/changelog/0.3.md b/docs/archive/changelog/0.3.md index 118719a..5c0ed16 100644 --- a/docs/archive/changelog/0.3.md +++ b/docs/archive/changelog/0.3.md @@ -17,12 +17,10 @@ - Add `num-traits` as explicit `exact` feature dependency - Add `examples/exact_det_3x3.rs` demonstrating exact value on near-singular matrix - - Add 22 tests for det_exact, det_exact_f64, and validate_finite - Update docs: README, AGENTS.md, REFERENCES.md, module doc, typos.toml - [**breaking**] Expose exact determinant value (det_exact, det_exact_f64) [`92ce476`](https://github.com/acgetchell/la-stack/commit/92ce476201d5759766d73221612cd27492bccbe5) - - Add `det_exact()` returning `BigRational` via Bareiss elimination - Add `det_exact_f64()` converting exact result to f64 - Add `LaError::Overflow` for f64 conversion overflow @@ -32,17 +30,14 @@ - Macro-generate exact tests for D=2..5 - Update README, AGENTS.md, REFERENCES.md, typos.toml - Exact linear system solve (solve_exact, solve_exact_f64) [`d04fcd3`](https://github.com/acgetchell/la-stack/commit/d04fcd3b24843b7377a34849dbb2e6469bf9fc56) - - Add `solve_exact()` returning `[BigRational; D]` via Gaussian elimination with partial pivoting in BigRational - - Add `solve_exact_f64()` converting exact result to `Vector<D>` - Add `validate_finite_vec()` helper for vector input validation - Add `exact_solve_3x3` example (near-singular system: f64 vs exact) - Macro-generate solve tests for D=2..5 - Update Cargo.toml: keywords (exact-arithmetic, robust-predicates), description ("Fast" not "Small") - - Simplify scalar types section (remove f32/f128 speculation) - Add `gh` CLI usage guidance to AGENTS.md - Update README, AGENTS.md with solve_exact documentation @@ -54,26 +49,21 @@ - Add LaError::Overflow variant for det_exact_f64 when the exact result exceeds f64 range (replaces misleading NonFinite { col: 0 }) - - Refactor det_exact/det_exact_f64 tests into macro-generated per-dimension suites (D=2..5) using pastey::paste - - Fix broken intra-doc links when exact feature is disabled - Add dimension coverage testing guidance to AGENTS.md - Clarify in AGENTS.md that only det_sign_exact uses the Shewchuk-style f64 filter - - Add #[non_exhaustive] to LaError [`d068c0e`](https://github.com/acgetchell/la-stack/commit/d068c0e62073a9746a36dab06761605635281481) Future error variants (e.g. from #42) won't require a breaking version bump. Simplify Overflow docstring accordingly. - - Use det_direct() in exact_det_3x3 example [`5b69fa8`](https://github.com/acgetchell/la-stack/commit/5b69fa8011cd5471af10be5db8fcde64ec9a9dfa) Replace det(DEFAULT_PIVOT_TOL) with det_direct() so the demo compares plain f64 cofactor expansion to the exact path without pivot-tolerance heuristics - - Restructure gauss_solve to avoid break/continue [`125b259`](https://github.com/acgetchell/la-stack/commit/125b259608a2c1965aba486d2594ff9fa328609c) Replace `found` boolean + `break` with idiomatic `.find()` + `if let`, @@ -87,23 +77,18 @@ - Replace "partial pivoting" with "first-non-zero pivoting" in module doc, function doc, and inline comments for both bareiss_det and gauss_solve - - Add rationale: exact BigRational arithmetic means any non-zero pivot gives a correct result (no numerical stability concern) - - Use idiomatic (0..D).rev() for back-substitution loop - Add LDLT example, solve_exact README snippet, and sync examples [`2bddccf`](https://github.com/acgetchell/la-stack/commit/2bddccf3f20156fadb4fee6a4093e76b1b631de1) - - Add examples/ldlt_solve_3x3.rs: 3×3 SPD tridiagonal LDLT solve + det - Add solve_exact_f64 code example to README exact arithmetic section - Add missing examples to justfile (ldlt_solve_3x3, exact_det_3x3, exact_solve_3x3) - - Add lu.rs and ldlt.rs test macro references to AGENTS.md - Fix stale "partial pivoting" → "first-non-zero pivoting" in AGENTS.md code structure - - Drop "current" from README storage note [0.3.0]: https://github.com/acgetchell/la-stack/compare/v0.2.2...v0.3.0 diff --git a/examples/exact_solve_3x3.rs b/examples/exact_solve_3x3.rs index 83c3c6b..2afe790 100644 --- a/examples/exact_solve_3x3.rs +++ b/examples/exact_solve_3x3.rs @@ -1,9 +1,10 @@ //! Exact linear system solve for a near-singular 3×3 system. //! -//! This example demonstrates `solve_exact()` and `solve_exact_f64()`, which use -//! arbitrary-precision rational arithmetic to compute a provably correct -//! solution — even when the matrix is so close to singular that the f64 LU -//! solve produces a wildly inaccurate result. +//! This example demonstrates `solve_exact()` and `solve_exact_f64()`. The exact +//! solve uses arbitrary-precision rational arithmetic to compute a provably +//! correct solution — even when the matrix is so close to singular that the f64 +//! LU solve produces a wildly inaccurate result. The `solve_exact_f64()` helper +//! only succeeds when every exact component is exactly representable as `f64`. //! //! Run with: `cargo run --features exact --example exact_solve_3x3` @@ -30,8 +31,6 @@ fn main() -> Result<(), LaError> { // Exact solve. let exact_x = a.solve_exact(b)?; - let exact_x_f64 = a.solve_exact_f64(b)?.into_array(); - println!("Near-singular 3×3 system (perturbation = 2^-50 ≈ {perturbation:.2e}):"); for r in 0..3 { print!(" ["); @@ -58,9 +57,23 @@ fn main() -> Result<(), LaError> { "solve_exact(): x = [{}, {}, {}]", exact_x[0], exact_x[1], exact_x[2] ); - println!( - "solve_exact_f64(): x = [{:+.6e}, {:+.6e}, {:+.6e}]", - exact_x_f64[0], exact_x_f64[1], exact_x_f64[2] - ); + match a.solve_exact_f64(b) { + Ok(x) => { + let x = x.into_array(); + println!( + "solve_exact_f64(): x = [{:+.6e}, {:+.6e}, {:+.6e}]", + x[0], x[1], x[2] + ); + } + Err(err) if err.requires_rounding() => { + println!("solve_exact_f64(): {err}"); + let x = a.solve_exact_rounded_f64(b)?.into_array(); + println!( + "rounded fallback: x = [{:+.6e}, {:+.6e}, {:+.6e}]", + x[0], x[1], x[2] + ); + } + Err(err) => return Err(err), + } Ok(()) } diff --git a/justfile b/justfile index 2d392c1..811f19b 100644 --- a/justfile +++ b/justfile @@ -169,6 +169,7 @@ bench: # This catches bench/release-profile-only warnings that won't show up in normal debug-profile runs. bench-compile: RUSTFLAGS='-D warnings' cargo bench --no-run --features bench + RUSTFLAGS='-D warnings' cargo bench --no-run --features bench,exact --bench exact # Run the cheaper latest measurements used for latest-vs-last reports. bench-latest: bench-vs-linalg-la-stack bench-exact diff --git a/scripts/bench_compare.py b/scripts/bench_compare.py index 5165de7..24dcea3 100644 --- a/scripts/bench_compare.py +++ b/scripts/bench_compare.py @@ -43,16 +43,35 @@ # # Mirrors the structure of `benches/exact.rs`: general-case per-dimension # groups (`exact_d{2..5}`), fixed-seed random percentile groups, plus -# adversarial/extreme-input groups that share a fixed four-bench layout -# (`det_sign_exact`, `det_exact`, `solve_exact`, `solve_exact_f64`). -_EXTREME_BENCHES: list[str] = ["det_sign_exact", "det_exact", "solve_exact", "solve_exact_f64"] +# adversarial/extreme-input groups that share a fixed five-bench layout +# (`det_sign_exact`, `det_exact`, `solve_exact`, `solve_exact_f64_result`, +# `solve_exact_rounded_f64`). +_EXTREME_BENCHES: list[str] = [ + "det_sign_exact", + "det_exact", + "solve_exact", + "solve_exact_f64_result", + "solve_exact_rounded_f64", +] _RANDOM_PERCENTILE_BENCHES: list[str] = [f"{operation}_{percentile}" for operation in _EXTREME_BENCHES for percentile in ("p50", "p95", "p99")] +_EXACT_DIMENSION_BENCHES: list[str] = [ + "det", + "det_direct", + "det_exact", + "det_exact_f64_result", + "det_exact_rounded_f64", + "det_sign_exact", + "solve_exact", + "solve_exact_f64_result", + "solve_exact_rounded_f64", +] + EXACT_GROUPS: dict[str, list[str]] = { - "exact_d2": ["det", "det_direct", "det_exact", "det_exact_f64", "det_sign_exact", "solve_exact", "solve_exact_f64"], - "exact_d3": ["det", "det_direct", "det_exact", "det_exact_f64", "det_sign_exact", "solve_exact", "solve_exact_f64"], - "exact_d4": ["det", "det_direct", "det_exact", "det_exact_f64", "det_sign_exact", "solve_exact", "solve_exact_f64"], - "exact_d5": ["det", "det_direct", "det_exact", "det_exact_f64", "det_sign_exact", "solve_exact", "solve_exact_f64"], + "exact_d2": _EXACT_DIMENSION_BENCHES, + "exact_d3": _EXACT_DIMENSION_BENCHES, + "exact_d4": _EXACT_DIMENSION_BENCHES, + "exact_d5": _EXACT_DIMENSION_BENCHES, "exact_random_percentile_d2": _RANDOM_PERCENTILE_BENCHES, "exact_random_percentile_d3": _RANDOM_PERCENTILE_BENCHES, "exact_random_percentile_d4": _RANDOM_PERCENTILE_BENCHES, @@ -63,6 +82,25 @@ "exact_hilbert_5x5": _EXTREME_BENCHES, } +EXACT_RELEASE_SIGNAL_GROUPS: frozenset[str] = frozenset(group for group in EXACT_GROUPS if not group.startswith("exact_random_percentile_d")) + +# v0.4.2 and earlier named the lossy exact-to-f64 benches after the public +# `*_exact_f64` API. Current benches split that behavior into strict `*_result` +# and lossy `*_rounded_f64` variants. Use the old baseline when present so +# release reports show both compatibility-successor performance and strict +# conversion overhead instead of silently dropping the new rows. +EXACT_LEGACY_BASELINE_BENCHES: dict[str, str] = { + "det_exact_f64_result": "det_exact_f64", + "det_exact_rounded_f64": "det_exact_f64", + "solve_exact_f64_result": "solve_exact_f64", + "solve_exact_rounded_f64": "solve_exact_f64", +} + +_EXACT_LEGACY_PREFIX_BASELINE_BENCHES: tuple[tuple[str, str], ...] = ( + ("solve_exact_f64_result_", "solve_exact_f64_"), + ("solve_exact_rounded_f64_", "solve_exact_f64_"), +) + VS_LINALG_LA_STACK_ONLY_BENCHES_BY_METRIC: dict[str, list[str]] = { "det_via_lu": ["la_stack_det"], } @@ -123,6 +161,7 @@ class Comparison: current_ns: float speedup: float # baseline / current (>1 = faster) pct_change: float # signed percent change (negative = faster) + baseline_bench: str | None = None baseline_nalgebra_ns: float | None = None baseline_faer_ns: float | None = None @@ -231,6 +270,36 @@ def _collect_exact_results(criterion_dir: Path, sample: str, stat: str) -> list[ return results +def _legacy_exact_baseline_bench(bench: str) -> str | None: + """Return the legacy exact benchmark name for renamed rows.""" + legacy_bench = EXACT_LEGACY_BASELINE_BENCHES.get(bench) + if legacy_bench is not None: + return legacy_bench + + for current_prefix, legacy_prefix in _EXACT_LEGACY_PREFIX_BASELINE_BENCHES: + if bench.startswith(current_prefix): + return f"{legacy_prefix}{bench.removeprefix(current_prefix)}" + + return None + + +def _exact_baseline_path(group_dir: Path, bench: str, baseline_name: str) -> tuple[str, Path]: + """Return the exact benchmark baseline path, falling back to legacy names.""" + base_path = group_dir / bench / baseline_name / "estimates.json" + if base_path.exists(): + return (bench, base_path) + + legacy_bench = _legacy_exact_baseline_bench(bench) + if legacy_bench is None: + return (bench, base_path) + + legacy_path = group_dir / legacy_bench / baseline_name / "estimates.json" + if legacy_path.exists(): + return (legacy_bench, legacy_path) + + return (bench, base_path) + + def _ordered_vs_linalg_benches(group_dir: Path, sample: str) -> list[str]: """Return present vs_linalg benches in a stable, metric-aware order.""" present = {child.name for child in group_dir.iterdir() if child.is_dir() and (child / sample / "estimates.json").exists()} @@ -275,18 +344,22 @@ def _collect_exact_comparisons( criterion_dir: Path, baseline_name: str, stat: str, + scope: str, ) -> list[Comparison]: """Compare current exact-arithmetic results against a named baseline.""" comparisons: list[Comparison] = [] for group, benches in EXACT_GROUPS.items(): + if scope == "release-signal" and group not in EXACT_RELEASE_SIGNAL_GROUPS: + continue + group_dir = criterion_dir / group if not group_dir.is_dir(): continue for bench in benches: new_path = group_dir / bench / "new" / "estimates.json" - base_path = group_dir / bench / baseline_name / "estimates.json" + baseline_bench, base_path = _exact_baseline_path(group_dir, bench, baseline_name) if not new_path.exists() or not base_path.exists(): continue @@ -306,6 +379,7 @@ def _collect_exact_comparisons( current_ns=new_point, speedup=speedup, pct_change=pct_change, + baseline_bench=baseline_bench if baseline_bench != bench else None, ) ) @@ -346,6 +420,13 @@ def _baseline_peer_times(group_dir: Path, bench: str, baseline_name: str, stat: return (nalgebra_ns, faer_ns) +def _comparison_bench_label(comparison: Comparison) -> str: + """Return the display label for a comparison table row.""" + if comparison.baseline_bench is None: + return comparison.bench + return f"{comparison.bench} (vs {comparison.baseline_bench})" + + def _collect_vs_linalg_comparisons( criterion_dir: Path, baseline_name: str, @@ -403,7 +484,7 @@ def _collect_comparisons( """Compare current (new) results against a named baseline.""" comparisons: list[Comparison] = [] if suite in ("all", "exact"): - comparisons.extend(_collect_exact_comparisons(criterion_dir, baseline_name, stat)) + comparisons.extend(_collect_exact_comparisons(criterion_dir, baseline_name, stat, scope)) if suite in ("all", "vs_linalg"): comparisons.extend(_collect_vs_linalg_comparisons(criterion_dir, baseline_name, stat, scope)) return comparisons @@ -542,7 +623,7 @@ def _comparison_tables(comparisons: list[Comparison], baseline_name: str) -> str ) for c in items: cells = [ - c.bench, + _comparison_bench_label(c), _format_time(c.baseline_ns), _format_time(c.current_ns), _format_pct(c.pct_change), diff --git a/scripts/tests/test_bench_compare.py b/scripts/tests/test_bench_compare.py index adb804e..e0744a2 100644 --- a/scripts/tests/test_bench_compare.py +++ b/scripts/tests/test_bench_compare.py @@ -34,12 +34,20 @@ def _build_criterion_tree(criterion_dir: Path, stat: str = "median") -> None: group = criterion_dir / f"exact_d{d}" _write_estimates(group / "det" / "new" / "estimates.json", stat, det) _write_estimates(group / "det_exact" / "new" / "estimates.json", stat, det_exact) + _write_estimates(group / "det_exact_f64_result" / "new" / "estimates.json", stat, det_exact * 1.1) + _write_estimates(group / "det_exact_rounded_f64" / "new" / "estimates.json", stat, det_exact * 1.2) + _write_estimates(group / "solve_exact_f64_result" / "new" / "estimates.json", stat, det_exact * 1.3) + _write_estimates(group / "solve_exact_rounded_f64" / "new" / "estimates.json", stat, det_exact * 1.4) ns_group = criterion_dir / "exact_near_singular_3x3" _write_estimates(ns_group / "det_sign_exact" / "new" / "estimates.json", stat, 12000.0) + _write_estimates(ns_group / "solve_exact_f64_result" / "new" / "estimates.json", stat, 51000.0) + _write_estimates(ns_group / "solve_exact_rounded_f64" / "new" / "estimates.json", stat, 52000.0) random_group = criterion_dir / "exact_random_percentile_d3" _write_estimates(random_group / "det_exact_p95" / "new" / "estimates.json", stat, 33000.0) + _write_estimates(random_group / "solve_exact_f64_result_p95" / "new" / "estimates.json", stat, 54000.0) + _write_estimates(random_group / "solve_exact_rounded_f64_p95" / "new" / "estimates.json", stat, 55000.0) def _build_vs_linalg_tree(criterion_dir: Path, stat: str = "median") -> None: @@ -196,7 +204,7 @@ def test_read_estimate_non_numeric_ci_bound_names_field(tmp_path: Path) -> None: def test_collect_results(tmp_path: Path) -> None: _build_criterion_tree(tmp_path) results = bench_compare._collect_results(tmp_path, "new", "median") - assert len(results) == 6 # 2 benches x 2 dims + 1 near-singular + 1 random percentile + assert len(results) == 18 # 6 benches x 2 dims + 3 near-singular + 3 random percentile groups = {r.group for r in results} assert "exact_d2" in groups assert "exact_d3" in groups @@ -216,11 +224,37 @@ def test_collect_comparisons(tmp_path: Path) -> None: group = tmp_path / f"exact_d{d}" _write_estimates(group / "det" / "v0.3.0" / "estimates.json", "median", det) _write_estimates(group / "det_exact" / "v0.3.0" / "estimates.json", "median", det_exact) + _write_estimates(group / "det_exact_f64" / "v0.3.0" / "estimates.json", "median", det_exact * 1.1) + _write_estimates(group / "solve_exact_f64" / "v0.3.0" / "estimates.json", "median", det_exact * 1.3) + random_group = tmp_path / "exact_random_percentile_d3" + _write_estimates(random_group / "det_exact_p95" / "v0.3.0" / "estimates.json", "median", 66000.0) + _write_estimates(random_group / "solve_exact_f64_p95" / "v0.3.0" / "estimates.json", "median", 108000.0) comparisons = bench_compare._collect_comparisons(tmp_path, "v0.3.0", "median") - assert len(comparisons) == 4 # 2 benches x 2 dims (near-singular has no baseline) + assert len(comparisons) == 12 # 6 benches x 2 dims (near-singular has no baseline) + assert {c.group for c in comparisons} == {"exact_d2", "exact_d3"} for c in comparisons: assert c.speedup == pytest.approx(c.baseline_ns / c.current_ns) + assert {(c.bench, c.baseline_bench) for c in comparisons if c.baseline_bench is not None} == { + ("det_exact_f64_result", "det_exact_f64"), + ("det_exact_rounded_f64", "det_exact_f64"), + ("solve_exact_f64_result", "solve_exact_f64"), + ("solve_exact_rounded_f64", "solve_exact_f64"), + } + + all_comparisons = bench_compare._collect_comparisons( + tmp_path, + "v0.3.0", + "median", + scope="all-benches", + ) + assert len(all_comparisons) == 15 + random_comparisons = [c for c in all_comparisons if c.group == "exact_random_percentile_d3"] + assert {(c.bench, c.baseline_bench) for c in random_comparisons} == { + ("det_exact_p95", None), + ("solve_exact_f64_result_p95", "solve_exact_f64_p95"), + ("solve_exact_rounded_f64_p95", "solve_exact_f64_p95"), + } def test_collect_comparisons_zero_current(tmp_path: Path) -> None: @@ -301,12 +335,16 @@ def test_comparison_tables_per_dimension(tmp_path: Path) -> None: group = tmp_path / f"exact_d{d}" _write_estimates(group / "det" / "v0.3.0" / "estimates.json", "median", det) _write_estimates(group / "det_exact" / "v0.3.0" / "estimates.json", "median", det_exact) + _write_estimates(group / "det_exact_f64" / "v0.3.0" / "estimates.json", "median", det_exact * 1.1) + _write_estimates(group / "solve_exact_f64" / "v0.3.0" / "estimates.json", "median", det_exact * 1.3) comparisons = bench_compare._collect_comparisons(tmp_path, "v0.3.0", "median") tables = bench_compare._comparison_tables(comparisons, "v0.3.0") assert "### D=2" in tables assert "### D=3" in tables assert "| Benchmark | v0.3.0 | Latest | Change | Speedup |" in tables + assert "det_exact_rounded_f64 (vs det_exact_f64)" in tables + assert "solve_exact_f64_result (vs solve_exact_f64)" in tables def test_comparison_tables_include_vs_linalg_peer_context(tmp_path: Path) -> None: diff --git a/semgrep.yaml b/semgrep.yaml index 793229f..f053965 100644 --- a/semgrep.yaml +++ b/semgrep.yaml @@ -110,22 +110,6 @@ rules: - "/tests/semgrep/src/project_rules/finite_api_contract.rs" pattern-regex: '(?ms)^\s*pub\s+struct\s+(?:Matrix|Vector)\s*(?:<[^>{]*>)?\s*\{(?:(?!^\s*\}).|\n)*^\s*pub(?:\([^)]*\))?\s+(?:rows|data)\s*:' - - id: la-stack.rust.no-public-finite-proof-reexports - languages: - - regex - severity: WARNING - message: "Do not publicly re-export FiniteMatrix/FiniteVector; Matrix/Vector are the public finite proof types." - metadata: - category: api - rationale: >- - Internal finite wrappers may exist to make algorithm boundaries explicit, - but downstream callers should not name or construct them directly. - paths: - include: - - "/src/**/*.rs" - - "/tests/semgrep/src/project_rules/finite_api_contract.rs" - pattern-regex: '(?m)^\s*pub\s+use\s+(?:crate::)?(?:matrix|vector)::(?:\{[^;\n]*(?:FiniteMatrix|FiniteVector)[^;\n]*\}|(?:FiniteMatrix|FiniteVector)(?:\s+as\s+[A-Za-z_][A-Za-z0-9_]*)?)\s*;' - - id: la-stack.rust.no-public-raw-linear-algebra-modules languages: - regex @@ -134,9 +118,9 @@ rules: metadata: category: api rationale: >- - The matrix and vector modules contain crate-internal proof wrappers and - unchecked helpers. Making the modules public would expose implementation - details that bypass the clean API surface. + The matrix and vector modules contain crate-internal unchecked helpers. + Making the modules public would expose implementation details that + bypass the clean API surface. paths: include: - "/src/**/*.rs" @@ -230,7 +214,7 @@ rules: include: - "/src/lib.rs" - "/tests/semgrep/src/project_rules/readme_doctest_mirrors.rs" - pattern-regex: '(?ms)^\s*mod\s+readme_doctests\s*\{(?:(?!^\}).|\n)*(?:\.unwrap\s*\(|\.expect\s*\()' + pattern-regex: '(?ms)^\s*mod\s+readme_doctests(?:_[A-Za-z0-9_]+)?\s*\{(?:(?!^\}).|\n)*(?:\.unwrap\s*\(|\.expect\s*\()' - id: la-stack.rust.no-unwrap-expect-in-benches-examples languages: diff --git a/src/error.rs b/src/error.rs new file mode 100644 index 0000000..41769dd --- /dev/null +++ b/src/error.rs @@ -0,0 +1,716 @@ +//! Error types and helpers for linear algebra operations. + +use core::fmt; + +use crate::Tolerance; + +/// Reason an exact result cannot satisfy an exact-to-`f64` conversion contract. +/// +/// `RequiresRounding` is recoverable when the caller is willing to opt into a +/// rounded exact-to-`f64` API. `NotFinite` means even the rounded result would +/// not be a finite `f64`. +/// +/// # Examples +/// ``` +/// use la_stack::prelude::*; +/// +/// let err = LaError::unrepresentable(None, UnrepresentableReason::RequiresRounding); +/// assert!(err.requires_rounding()); +/// +/// let err = LaError::unrepresentable(None, UnrepresentableReason::NotFinite); +/// assert_eq!( +/// err.unrepresentable_reason(), +/// Some(UnrepresentableReason::NotFinite) +/// ); +/// ``` +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +#[non_exhaustive] +pub enum UnrepresentableReason { + /// A finite `f64` exists only after rounding, but the requested conversion + /// requires an exact binary64 representation. + RequiresRounding, + /// The exact value would convert to NaN or infinity rather than a finite + /// `f64`. + NotFinite, +} + +/// Linear algebra errors. +/// +/// This enum is `#[non_exhaustive]` — downstream `match` arms must include a +/// wildcard (`_`) pattern to compile, allowing new variants to be added in +/// future minor releases without breaking existing code. +/// +/// # Examples +/// ``` +/// use la_stack::prelude::*; +/// +/// let err = LaError::unrepresentable(None, UnrepresentableReason::RequiresRounding); +/// assert!(err.requires_rounding()); +/// +/// match LaError::unsupported_dimension(8, MAX_STACK_MATRIX_DISPATCH_DIM) { +/// LaError::UnsupportedDimension { requested, max } => { +/// assert_eq!((requested, max), (8, MAX_STACK_MATRIX_DISPATCH_DIM)); +/// } +/// _ => unreachable!("constructor returns the requested variant"), +/// } +/// ``` +#[derive(Clone, Copy, Debug, PartialEq)] +#[non_exhaustive] +pub enum LaError { + /// The matrix is (numerically) singular. + Singular { + /// The factorization column/step where a suitable pivot/diagonal could not be found. + pivot_col: usize, + }, + /// A non-finite value (NaN/∞) was encountered. + /// + /// The `(row, col)` coordinate follows a consistent convention across the crate: + /// + /// - `row: Some(r), col: c` — the non-finite value is tied to a matrix/factor + /// cell at `(r, c)`, either because a stored input/factor cell is already + /// non-finite or because factorization computed a non-finite value for + /// that cell before storing it. + /// - `row: None, col: c` — the non-finite value is tied to a vector entry, + /// determinant product, solve accumulator, or other scalar/intermediate + /// that has no matrix row coordinate. + NonFinite { + /// Row of the non-finite entry for a stored matrix cell, or `None` for + /// a vector-input entry or a computed intermediate. See the variant + /// docs for the full convention. + row: Option, + /// Column index (stored cell), vector index, or factorization/solve + /// step where the non-finite value was detected. + col: usize, + }, + /// An exact result cannot satisfy the requested finite `f64` conversion. + /// + /// Returned by [`Matrix::det_exact_f64`](crate::Matrix::det_exact_f64) and + /// [`Matrix::solve_exact_f64`](crate::Matrix::solve_exact_f64) (requires the + /// `exact` feature) when the exact rational value is too large, too small, + /// or would require rounding in binary64. Also returned by the rounded + /// exact-to-`f64` APIs when the rounded result would be NaN or infinite. + Unrepresentable { + /// For vector results (e.g. `solve_exact_f64`), the index of the + /// component that failed conversion. `None` for scalar results. + index: Option, + /// Why the requested conversion cannot return a finite `f64`. + reason: UnrepresentableReason, + }, + /// Exact determinant scaling overflowed the internal exponent representation. + DeterminantScaleOverflow { + /// Matrix dimension `D`. + dim: usize, + /// Minimum decomposed f64 exponent among non-zero matrix entries. + min_exponent: i32, + }, + /// A requested runtime matrix dimension has no stack-dispatch arm. + UnsupportedDimension { + /// Runtime dimension requested by the caller. + requested: usize, + /// Largest runtime dimension supported by the dispatch helper. + max: usize, + }, + /// A matrix index is outside the `D×D` storage domain. + IndexOutOfBounds { + /// Requested row index. + row: usize, + /// Requested column index. + col: usize, + /// Matrix dimension `D`; valid row and column indices are `< dim`. + dim: usize, + }, + /// A tolerance value is not finite and non-negative. + InvalidTolerance { + /// Raw tolerance supplied by the caller. + value: f64, + }, + /// A matrix required to be symmetric has an asymmetric off-diagonal pair. + Asymmetric { + /// Row index of the first asymmetric pair. + row: usize, + /// Column index of the first asymmetric pair. + col: usize, + /// Matrix dimension `D`. + dim: usize, + }, + /// A symmetric matrix failed the positive-semidefinite LDLT domain check. + NotPositiveSemidefinite { + /// Factorization column/step where a negative LDLT diagonal was found. + pivot_col: usize, + /// Negative diagonal value observed at that step. + value: f64, + }, +} + +impl LaError { + /// Construct a [`LaError::NonFinite`] pinpointing a stored matrix cell at `(row, col)`. + /// + /// Use this for non-finite values read from a stored [`Matrix`](crate::Matrix) + /// entry or factorization cell, and for non-finite factorization updates + /// that would be stored at `(row, col)` if accepted. The resulting error has + /// `row: Some(row), col`, matching the matrix/factor-cell convention + /// documented on [`NonFinite`](Self::NonFinite). For vector-input entries + /// or scalar intermediates without a matrix row coordinate, use + /// [`non_finite_at`](Self::non_finite_at). + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// assert_eq!( + /// LaError::non_finite_cell(1, 2), + /// LaError::NonFinite { + /// row: Some(1), + /// col: 2, + /// } + /// ); + /// ``` + #[inline] + #[must_use] + pub const fn non_finite_cell(row: usize, col: usize) -> Self { + Self::NonFinite { + row: Some(row), + col, + } + } + + /// Construct a [`LaError::NonFinite`] pinpointing a vector-input entry or + /// computed scalar/intermediate at index `col`. + /// + /// Use this for non-finite values in a [`Vector`](crate::Vector) input, + /// determinant scalar, tolerance-scale accumulator, or solve accumulator + /// that overflowed during forward/back substitution. The resulting error + /// has `row: None, col`, matching the vector/scalar-intermediate convention + /// documented on [`NonFinite`](Self::NonFinite). For stored matrix cells or + /// computed factorization updates tied to a matrix cell, use + /// [`non_finite_cell`](Self::non_finite_cell). + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// assert_eq!( + /// LaError::non_finite_at(2), + /// LaError::NonFinite { row: None, col: 2 } + /// ); + /// ``` + #[inline] + #[must_use] + pub const fn non_finite_at(col: usize) -> Self { + Self::NonFinite { row: None, col } + } + + /// Construct a [`LaError::Unrepresentable`] for exact-to-`f64` conversion. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// assert_eq!( + /// LaError::unrepresentable(Some(2), UnrepresentableReason::RequiresRounding), + /// LaError::Unrepresentable { + /// index: Some(2), + /// reason: UnrepresentableReason::RequiresRounding, + /// } + /// ); + /// ``` + #[inline] + #[must_use] + pub const fn unrepresentable(index: Option, reason: UnrepresentableReason) -> Self { + Self::Unrepresentable { index, reason } + } + + /// Return the reason for an exact-to-`f64` conversion failure. + /// + /// This is a concise alternative to matching the full + /// [`LaError::Unrepresentable`] variant when callers only need the + /// conversion reason. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// let err = LaError::unrepresentable(None, UnrepresentableReason::RequiresRounding); + /// assert_eq!( + /// err.unrepresentable_reason(), + /// Some(UnrepresentableReason::RequiresRounding) + /// ); + /// assert_eq!(LaError::Singular { pivot_col: 0 }.unrepresentable_reason(), None); + /// ``` + #[inline] + #[must_use] + pub const fn unrepresentable_reason(&self) -> Option { + match self { + Self::Unrepresentable { reason, .. } => Some(*reason), + _ => None, + } + } + + /// Return `true` when strict exact-to-`f64` conversion only failed because + /// rounding would be required. + /// + /// This is useful at the call site that wants to retry with an explicit + /// rounded exact-to-`f64` API while still propagating non-finite conversion + /// failures. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// let err = LaError::unrepresentable(None, UnrepresentableReason::RequiresRounding); + /// assert!(err.requires_rounding()); + /// + /// let err = LaError::unrepresentable(None, UnrepresentableReason::NotFinite); + /// assert!(!err.requires_rounding()); + /// ``` + #[inline] + #[must_use] + pub const fn requires_rounding(&self) -> bool { + matches!( + self, + Self::Unrepresentable { + reason: UnrepresentableReason::RequiresRounding, + .. + } + ) + } + + /// Construct a [`LaError::DeterminantScaleOverflow`] for exact determinant scaling. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// assert_eq!( + /// LaError::determinant_scale_overflow(3, -1074), + /// LaError::DeterminantScaleOverflow { + /// dim: 3, + /// min_exponent: -1074, + /// } + /// ); + /// ``` + #[inline] + #[must_use] + pub const fn determinant_scale_overflow(dim: usize, min_exponent: i32) -> Self { + Self::DeterminantScaleOverflow { dim, min_exponent } + } + + /// Construct a [`LaError::UnsupportedDimension`] for runtime stack dispatch. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// assert_eq!( + /// LaError::unsupported_dimension(8, MAX_STACK_MATRIX_DISPATCH_DIM), + /// LaError::UnsupportedDimension { + /// requested: 8, + /// max: MAX_STACK_MATRIX_DISPATCH_DIM, + /// } + /// ); + /// ``` + #[inline] + #[must_use] + pub const fn unsupported_dimension(requested: usize, max: usize) -> Self { + Self::UnsupportedDimension { requested, max } + } + + /// Construct a [`LaError::IndexOutOfBounds`] for a `D×D` matrix index. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// assert_eq!( + /// LaError::index_out_of_bounds(2, 0, 2), + /// LaError::IndexOutOfBounds { + /// row: 2, + /// col: 0, + /// dim: 2, + /// } + /// ); + /// ``` + #[inline] + #[must_use] + pub const fn index_out_of_bounds(row: usize, col: usize, dim: usize) -> Self { + Self::IndexOutOfBounds { row, col, dim } + } + + /// Construct a [`LaError::InvalidTolerance`] for a raw tolerance value. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// assert_eq!( + /// LaError::invalid_tolerance(-1.0), + /// LaError::InvalidTolerance { value: -1.0 } + /// ); + /// ``` + #[inline] + #[must_use] + pub const fn invalid_tolerance(value: f64) -> Self { + Self::InvalidTolerance { value } + } + + /// Construct a [`LaError::Asymmetric`] for a `D×D` matrix. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// assert_eq!( + /// LaError::asymmetric(0, 1, 3), + /// LaError::Asymmetric { + /// row: 0, + /// col: 1, + /// dim: 3, + /// } + /// ); + /// ``` + #[inline] + #[must_use] + pub const fn asymmetric(row: usize, col: usize, dim: usize) -> Self { + Self::Asymmetric { row, col, dim } + } + + /// Construct a [`LaError::NotPositiveSemidefinite`] for LDLT factorization. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// assert_eq!( + /// LaError::not_positive_semidefinite(1, -3.0), + /// LaError::NotPositiveSemidefinite { + /// pivot_col: 1, + /// value: -3.0, + /// } + /// ); + /// ``` + #[inline] + #[must_use] + pub const fn not_positive_semidefinite(pivot_col: usize, value: f64) -> Self { + Self::NotPositiveSemidefinite { pivot_col, value } + } + + /// Parse a raw tolerance into a finite, non-negative [`Tolerance`]. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// assert_eq!(LaError::validate_tolerance(1e-12)?.get(), 1e-12); + /// + /// let raw = 0.0; + /// let tol = LaError::validate_tolerance(raw)?; + /// let _lu = Matrix::<2>::identity().lu(tol)?; + /// + /// assert_eq!( + /// LaError::validate_tolerance(-1.0), + /// Err(LaError::InvalidTolerance { value: -1.0 }) + /// ); + /// # Ok::<(), LaError>(()) + /// ``` + /// + /// # Errors + /// Returns [`LaError::InvalidTolerance`] when `value` is NaN, infinite, or + /// negative. + #[inline] + pub const fn validate_tolerance(value: f64) -> Result { + Tolerance::new(value) + } +} + +impl fmt::Display for LaError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match *self { + Self::Singular { pivot_col } => { + write!(f, "singular matrix at pivot column {pivot_col}") + } + Self::NonFinite { row: Some(r), col } => { + write!(f, "non-finite value at ({r}, {col})") + } + Self::NonFinite { row: None, col } => { + write!(f, "non-finite value at index {col}") + } + Self::Unrepresentable { + index: Some(i), + reason: UnrepresentableReason::RequiresRounding, + } => write!( + f, + "exact result requires rounding to fit finite f64 at index {i}" + ), + Self::Unrepresentable { + index: None, + reason: UnrepresentableReason::RequiresRounding, + } => write!(f, "exact result requires rounding to fit finite f64"), + Self::Unrepresentable { + index: Some(i), + reason: UnrepresentableReason::NotFinite, + } => write!(f, "exact result does not round to finite f64 at index {i}"), + Self::Unrepresentable { + index: None, + reason: UnrepresentableReason::NotFinite, + } => write!(f, "exact result does not round to finite f64"), + Self::DeterminantScaleOverflow { dim, min_exponent } => { + write!( + f, + "exact determinant scale exponent overflows for dimension {dim} with minimum entry exponent {min_exponent}" + ) + } + Self::UnsupportedDimension { requested, max } => { + write!( + f, + "unsupported matrix dimension {requested}; maximum stack-dispatch dimension is {max}" + ) + } + Self::IndexOutOfBounds { row, col, dim } => { + write!( + f, + "matrix index ({row}, {col}) is out of bounds for dimension {dim}" + ) + } + Self::InvalidTolerance { value } => { + write!(f, "invalid tolerance {value}; expected finite value >= 0") + } + Self::Asymmetric { row, col, dim } => { + write!( + f, + "matrix is not symmetric for dimension {dim}: asymmetric pair ({row}, {col})" + ) + } + Self::NotPositiveSemidefinite { pivot_col, value } => { + write!( + f, + "matrix is not positive semidefinite at LDLT pivot column {pivot_col}: diagonal value {value} < 0" + ) + } + } + } +} + +impl std::error::Error for LaError {} + +#[cfg(test)] +mod tests { + use super::*; + + use core::assert_matches; + + #[test] + fn laerror_display_formats_singular() { + let err = LaError::Singular { pivot_col: 3 }; + assert_eq!(err.to_string(), "singular matrix at pivot column 3"); + } + + #[test] + fn laerror_display_formats_nonfinite_with_row() { + let err = LaError::NonFinite { + row: Some(1), + col: 2, + }; + assert_eq!(err.to_string(), "non-finite value at (1, 2)"); + } + + #[test] + fn laerror_display_formats_nonfinite_without_row() { + let err = LaError::NonFinite { row: None, col: 3 }; + assert_eq!(err.to_string(), "non-finite value at index 3"); + } + + #[test] + fn laerror_display_formats_unrepresentable_requires_rounding() { + let err = LaError::Unrepresentable { + index: None, + reason: UnrepresentableReason::RequiresRounding, + }; + assert_eq!( + err.to_string(), + "exact result requires rounding to fit finite f64" + ); + } + + #[test] + fn laerror_display_formats_unrepresentable_requires_rounding_with_index() { + let err = LaError::Unrepresentable { + index: Some(2), + reason: UnrepresentableReason::RequiresRounding, + }; + assert_eq!( + err.to_string(), + "exact result requires rounding to fit finite f64 at index 2" + ); + } + + #[test] + fn laerror_display_formats_unrepresentable_not_finite() { + let err = LaError::Unrepresentable { + index: None, + reason: UnrepresentableReason::NotFinite, + }; + assert_eq!(err.to_string(), "exact result does not round to finite f64"); + } + + #[test] + fn laerror_display_formats_unrepresentable_not_finite_with_index() { + let err = LaError::Unrepresentable { + index: Some(2), + reason: UnrepresentableReason::NotFinite, + }; + assert_eq!( + err.to_string(), + "exact result does not round to finite f64 at index 2" + ); + } + + #[test] + fn laerror_unrepresentable_reason_reports_typed_reason() { + let rounding = LaError::Unrepresentable { + index: Some(2), + reason: UnrepresentableReason::RequiresRounding, + }; + let not_finite = LaError::Unrepresentable { + index: None, + reason: UnrepresentableReason::NotFinite, + }; + + assert_eq!( + rounding.unrepresentable_reason(), + Some(UnrepresentableReason::RequiresRounding) + ); + assert_eq!( + not_finite.unrepresentable_reason(), + Some(UnrepresentableReason::NotFinite) + ); + assert_eq!( + LaError::Singular { pivot_col: 0 }.unrepresentable_reason(), + None + ); + } + + #[test] + fn laerror_requires_rounding_only_matches_rounding_reason() { + assert!( + LaError::Unrepresentable { + index: Some(2), + reason: UnrepresentableReason::RequiresRounding, + } + .requires_rounding() + ); + assert!( + !LaError::Unrepresentable { + index: None, + reason: UnrepresentableReason::NotFinite, + } + .requires_rounding() + ); + assert!(!LaError::Singular { pivot_col: 0 }.requires_rounding()); + } + + #[test] + fn laerror_display_formats_determinant_scale_overflow() { + let err = LaError::DeterminantScaleOverflow { + dim: 3, + min_exponent: -1074, + }; + assert_eq!( + err.to_string(), + "exact determinant scale exponent overflows for dimension 3 with minimum entry exponent -1074" + ); + } + + #[test] + fn laerror_display_formats_unsupported_dimension() { + let err = LaError::UnsupportedDimension { + requested: 8, + max: crate::MAX_STACK_MATRIX_DISPATCH_DIM, + }; + assert_eq!( + err.to_string(), + "unsupported matrix dimension 8; maximum stack-dispatch dimension is 7" + ); + } + + #[test] + fn laerror_display_formats_index_out_of_bounds() { + let err = LaError::IndexOutOfBounds { + row: 3, + col: 0, + dim: 3, + }; + assert_eq!( + err.to_string(), + "matrix index (3, 0) is out of bounds for dimension 3" + ); + } + + #[test] + fn laerror_display_formats_invalid_tolerance() { + let err = LaError::InvalidTolerance { value: -1.0 }; + assert_eq!( + err.to_string(), + "invalid tolerance -1; expected finite value >= 0" + ); + } + + #[test] + fn validate_tolerance_matches_tolerance_new() { + for value in [0.0, 1e-12, f64::MAX] { + assert_eq!(LaError::validate_tolerance(value), Tolerance::new(value)); + } + + assert_eq!( + LaError::validate_tolerance(-1.0), + Err(LaError::InvalidTolerance { value: -1.0 }) + ); + assert_matches!( + LaError::validate_tolerance(f64::NAN), + Err(LaError::InvalidTolerance { value }) if value.is_nan() + ); + assert_eq!( + LaError::validate_tolerance(f64::INFINITY), + Err(LaError::InvalidTolerance { + value: f64::INFINITY, + }) + ); + assert_eq!( + LaError::validate_tolerance(f64::NEG_INFINITY), + Err(LaError::InvalidTolerance { + value: f64::NEG_INFINITY, + }) + ); + } + + #[test] + fn laerror_display_formats_asymmetric() { + let err = LaError::Asymmetric { + row: 0, + col: 2, + dim: 3, + }; + assert_eq!( + err.to_string(), + "matrix is not symmetric for dimension 3: asymmetric pair (0, 2)" + ); + } + + #[test] + fn laerror_display_formats_not_positive_semidefinite() { + let err = LaError::NotPositiveSemidefinite { + pivot_col: 1, + value: -3.0, + }; + assert_eq!( + err.to_string(), + "matrix is not positive semidefinite at LDLT pivot column 1: diagonal value -3 < 0" + ); + } + + #[test] + fn laerror_is_std_error_with_no_source() { + let err = LaError::Singular { pivot_col: 0 }; + let e: &dyn std::error::Error = &err; + assert!(e.source().is_none()); + } +} diff --git a/src/exact.rs b/src/exact.rs index 5ac3cec..f2aecdb 100644 --- a/src/exact.rs +++ b/src/exact.rs @@ -6,17 +6,19 @@ //! //! ## Determinants //! -//! All three determinant methods (`det_exact`, `det_exact_f64`, `det_sign_exact`) -//! share the same integer-only Bareiss core (`bareiss_det_int`). Each f64 -//! entry is decomposed via `f64_decompose` into `mantissa × 2^exponent`, -//! all entries are scaled to a common `BigInt` matrix (shifting by -//! `e - e_min`), and Bareiss elimination runs entirely in `BigInt` -//! arithmetic — no `BigRational`, no GCD, no denominator tracking. -//! The result is `(det_int, total_exp)` where `det = det_int × 2^(D × e_min)`. -//! `det_exact` wraps this with `bigint_exp_to_bigrational` to reconstruct a -//! reduced `BigRational`; `det_exact_f64` converts the same pair directly to -//! `f64`; and `det_sign_exact` reads the sign directly from `det_int` (the -//! scale factor is always positive). +//! All determinant methods (`det_exact`, `det_exact_f64`, +//! `det_exact_rounded_f64`, and `det_sign_exact`) share the same integer-scaled +//! determinant core. Each f64 entry is decomposed via `f64_decompose` into +//! `mantissa × 2^exponent`, then all entries are scaled to a common `BigInt` +//! matrix (shifting by `e - e_min`). D≤4 uses direct integer expansions; larger +//! matrices use fraction-free Bareiss elimination entirely in `BigInt` +//! arithmetic — no `BigRational`, no GCD, no denominator tracking. The result +//! is `(det_int, total_exp)` where `det = det_int × 2^(D × e_min)`. `det_exact` +//! wraps this with `bigint_exp_to_bigrational` to reconstruct a reduced +//! `BigRational`; `det_exact_f64` converts the same pair only when the exact +//! value is representable as finite binary64; `det_exact_rounded_f64` rounds +//! the same exact value to finite binary64; and `det_sign_exact` reads the sign +//! directly from `det_int` (the scale factor is always positive). //! //! `det_sign_exact` adds a two-stage adaptive-precision optimisation inspired //! by Shewchuk's robust geometric predicates: @@ -29,9 +31,9 @@ //! //! ## Linear system solve //! -//! `solve_exact` and `solve_exact_f64` solve `A x = b` with a hybrid -//! algorithm that reuses the integer-only Bareiss core used for -//! determinants. Matrix and RHS entries are decomposed via +//! `solve_exact`, `solve_exact_f64`, and `solve_exact_rounded_f64` solve +//! `A x = b` with a hybrid algorithm that reuses the integer-only Bareiss core +//! used for determinants. Matrix and RHS entries are decomposed via //! `f64_decompose` into `mantissa × 2^exponent`, scaled to a shared //! base `2^e_min`, and assembled into a `BigInt` augmented system //! `(A | b)`. Forward elimination runs entirely in `BigInt` with @@ -42,8 +44,10 @@ //! overhead is modest. First-non-zero pivoting is used throughout; //! since all arithmetic is exact, any non-zero pivot gives the correct //! result (no numerical stability concern). Every finite `f64` is -//! exactly representable as a rational, so the result is provably -//! correct. +//! exactly representable as a rational, so the result is provably correct. +//! `solve_exact_f64` returns `Vector` only when every exact component is +//! exactly representable as finite binary64; `solve_exact_rounded_f64` returns +//! the exact components rounded to finite binary64. //! //! ## f64 → integer decomposition //! @@ -75,9 +79,17 @@ use num_bigint::{BigInt, Sign}; use num_rational::BigRational; use num_traits::ToPrimitive; -use crate::LaError; -use crate::matrix::{FiniteMatrix, Matrix}; -use crate::vector::{FiniteVector, Vector}; +use crate::matrix::Matrix; +use crate::vector::Vector; +use crate::{LaError, UnrepresentableReason}; + +const F64_SIGNIFICAND_BITS: i64 = 53; +const F64_FRACTION_BITS: i64 = 52; +const F64_MIN_BINARY_EXPONENT: i64 = -1074; +const F64_MIN_NORMAL_EXPONENT: i64 = -1022; +const F64_MAX_BINARY_EXPONENT: i64 = 1023; +const F64_EXPONENT_BIAS: i64 = 1023; +const F64_FRACTION_MASK: u64 = (1u64 << 52) - 1; /// Decompose a finite `f64` into its IEEE 754 components. /// @@ -160,31 +172,183 @@ fn bigint_exp_to_bigrational(mut value: BigInt, mut exp: i32) -> BigRational { } } -/// Convert a `BigInt × 2^exp` determinant pair to finite `f64` without first -/// reducing a public `BigRational` determinant value. -fn bigint_exp_to_finite_f64(value: BigInt, exp: i32) -> Result { +/// Convert a finite `f64` back to the exact rational value it represents. +fn finite_f64_to_bigrational(value: f64) -> Result { + let Some((mantissa, exponent, is_negative)) = f64_decompose(value)? else { + return Ok(BigRational::from_integer(BigInt::from(0))); + }; + + let value = BigInt::from(mantissa.get()); + let value = if is_negative { -value } else { value }; + Ok(bigint_exp_to_bigrational(value, exponent)) +} + +/// Convert an exact rational result to `f64` only when the conversion is exact. +fn exact_rational_to_finite_f64(exact: &BigRational, index: Option) -> Result { + let Some(value) = exact.to_f64() else { + cold_path(); + return Err(LaError::unrepresentable( + index, + UnrepresentableReason::NotFinite, + )); + }; + if !value.is_finite() { + cold_path(); + return Err(LaError::unrepresentable( + index, + UnrepresentableReason::NotFinite, + )); + } + + if finite_f64_to_bigrational(value)? == *exact { + Ok(value) + } else { + cold_path(); + Err(LaError::unrepresentable( + index, + UnrepresentableReason::RequiresRounding, + )) + } +} + +/// Convert an exact rational result to a rounded finite `f64`. +fn exact_rational_to_rounded_f64( + exact: &BigRational, + index: Option, +) -> Result { + let Some(value) = exact.to_f64() else { + cold_path(); + return Err(LaError::unrepresentable( + index, + UnrepresentableReason::NotFinite, + )); + }; + if value.is_finite() { + Ok(value) + } else { + cold_path(); + Err(LaError::unrepresentable( + index, + UnrepresentableReason::NotFinite, + )) + } +} + +/// Convert a `BigInt × 2^exp` determinant pair to an exactly represented finite +/// `f64` without allocating a `BigRational`. +fn bigint_exp_to_finite_f64(mut value: BigInt, exp: i32) -> Result { if value == BigInt::from(0) { return Ok(0.0); } - let exact = if exp >= 0 { - BigRational::new_raw(value << exp.cast_unsigned(), BigInt::from(1u32)) - } else { - BigRational::new_raw(value, BigInt::from(1u32) << exp.unsigned_abs()) + let is_negative = value.sign() == Sign::Minus; + if is_negative { + value = -value; + } + + let mut exp = i64::from(exp); + if let Some(tz) = value.trailing_zeros() { + value >>= tz; + let Ok(tz) = i64::try_from(tz) else { + cold_path(); + return Err(LaError::unrepresentable( + None, + UnrepresentableReason::NotFinite, + )); + }; + let Some(updated_exp) = exp.checked_add(tz) else { + cold_path(); + return Err(LaError::unrepresentable( + None, + UnrepresentableReason::NotFinite, + )); + }; + exp = updated_exp; + } + + let bit_len = value.bits(); + let Ok(bit_len) = i64::try_from(bit_len) else { + cold_path(); + return Err(LaError::unrepresentable( + None, + UnrepresentableReason::NotFinite, + )); }; + let Some(top_bit_exp) = exp.checked_add(bit_len - 1) else { + cold_path(); + return Err(LaError::unrepresentable( + None, + UnrepresentableReason::NotFinite, + )); + }; + if exp < F64_MIN_BINARY_EXPONENT { + cold_path(); + return Err(LaError::unrepresentable( + None, + UnrepresentableReason::RequiresRounding, + )); + } + if top_bit_exp > F64_MAX_BINARY_EXPONENT { + cold_path(); + return Err(LaError::unrepresentable( + None, + UnrepresentableReason::NotFinite, + )); + } + if bit_len > F64_SIGNIFICAND_BITS { + cold_path(); + return Err(LaError::unrepresentable( + None, + UnrepresentableReason::RequiresRounding, + )); + } - let Some(val) = exact.to_f64() else { + let Some(mantissa) = value.to_u64() else { cold_path(); - return Err(LaError::Overflow { index: None }); + return Err(LaError::unrepresentable( + None, + UnrepresentableReason::NotFinite, + )); }; - if val.is_finite() { - Ok(val) + let sign = if is_negative { 1u64 << 63 } else { 0 }; + + if top_bit_exp < F64_MIN_NORMAL_EXPONENT { + let Ok(shift) = u32::try_from(exp - F64_MIN_BINARY_EXPONENT) else { + cold_path(); + return Err(LaError::unrepresentable( + None, + UnrepresentableReason::RequiresRounding, + )); + }; + Ok(f64::from_bits(sign | (mantissa << shift))) } else { - cold_path(); - Err(LaError::Overflow { index: None }) + let Ok(biased_exp) = u64::try_from(top_bit_exp + F64_EXPONENT_BIAS) else { + cold_path(); + return Err(LaError::unrepresentable( + None, + UnrepresentableReason::NotFinite, + )); + }; + let Ok(shift) = u32::try_from(F64_FRACTION_BITS - (bit_len - 1)) else { + cold_path(); + return Err(LaError::unrepresentable( + None, + UnrepresentableReason::RequiresRounding, + )); + }; + let significand = mantissa << shift; + Ok(f64::from_bits( + sign | (biased_exp << F64_FRACTION_BITS) | (significand & F64_FRACTION_MASK), + )) } } +/// Convert a `BigInt × 2^exp` determinant pair to a rounded finite `f64`. +fn bigint_exp_to_rounded_f64(value: BigInt, exp: i32) -> Result { + let exact = bigint_exp_to_bigrational(value, exp); + exact_rational_to_rounded_f64(&exact, None) +} + // ----------------------------------------------------------------------- // Shared integer-Bareiss primitives // ----------------------------------------------------------------------- @@ -219,12 +383,11 @@ enum Component { /// Returns the per-entry components and the minimum exponent across non-zero /// entries. If every entry is zero, the exponent is `i32::MAX`. fn decompose_finite_matrix( - m: &FiniteMatrix, + m: &Matrix, ) -> Result<([[Component; D]; D], i32), LaError> { let mut components = [[Component::default(); D]; D]; let mut e_min = i32::MAX; - let matrix = m.as_matrix(); - for (r, row) in matrix.rows().iter().enumerate() { + for (r, row) in m.rows().iter().enumerate() { for (c, &entry) in row.iter().enumerate() { if let Some((mantissa, exponent, is_negative)) = f64_decompose(entry).map_err(|_| LaError::non_finite_cell(r, c))? @@ -245,9 +408,7 @@ fn decompose_finite_matrix( /// /// Returns the per-entry components and the minimum exponent across non-zero /// entries. If every entry is zero, the exponent is `i32::MAX`. -fn decompose_finite_vec( - v: &FiniteVector, -) -> Result<([Component; D], i32), LaError> { +fn decompose_finite_vec(v: &Vector) -> Result<([Component; D], i32), LaError> { let mut components = [Component::default(); D]; let mut e_min = i32::MAX; let data = v.as_array(); @@ -298,6 +459,77 @@ fn build_bigint_vec(components: &[Component; D], e_min: i32) -> from_fn(|i| component_to_bigint(components[i], e_min)) } +/// Compute a 2×2 determinant from a scaled integer matrix. +#[inline] +fn det2_bigint(a: &[[BigInt; D]; D]) -> BigInt { + &a[0][0] * &a[1][1] - &a[0][1] * &a[1][0] +} + +/// Compute a 3×3 determinant from scaled integer entries. +#[inline] +#[allow(clippy::too_many_arguments)] +fn det3_bigint_entries( + a00: &BigInt, + a01: &BigInt, + a02: &BigInt, + a10: &BigInt, + a11: &BigInt, + a12: &BigInt, + a20: &BigInt, + a21: &BigInt, + a22: &BigInt, +) -> BigInt { + let m00 = a11 * a22 - a12 * a21; + let m01 = a10 * a22 - a12 * a20; + let m02 = a10 * a21 - a11 * a20; + a00 * m00 - a01 * m01 + a02 * m02 +} + +/// Compute a 3×3 determinant from a scaled integer matrix. +#[inline] +fn det3_bigint(a: &[[BigInt; D]; D]) -> BigInt { + det3_bigint_entries( + &a[0][0], &a[0][1], &a[0][2], &a[1][0], &a[1][1], &a[1][2], &a[2][0], &a[2][1], &a[2][2], + ) +} + +/// Compute a 4×4 determinant from a scaled integer matrix. +#[inline] +fn det4_bigint(a: &[[BigInt; D]; D]) -> BigInt { + let mut det = BigInt::from(0); + + if a[0][0].sign() != Sign::NoSign { + let c00 = det3_bigint_entries( + &a[1][1], &a[1][2], &a[1][3], &a[2][1], &a[2][2], &a[2][3], &a[3][1], &a[3][2], + &a[3][3], + ); + det += &a[0][0] * c00; + } + if a[0][1].sign() != Sign::NoSign { + let c01 = det3_bigint_entries( + &a[1][0], &a[1][2], &a[1][3], &a[2][0], &a[2][2], &a[2][3], &a[3][0], &a[3][2], + &a[3][3], + ); + det -= &a[0][1] * c01; + } + if a[0][2].sign() != Sign::NoSign { + let c02 = det3_bigint_entries( + &a[1][0], &a[1][1], &a[1][3], &a[2][0], &a[2][1], &a[2][3], &a[3][0], &a[3][1], + &a[3][3], + ); + det += &a[0][2] * c02; + } + if a[0][3].sign() != Sign::NoSign { + let c03 = det3_bigint_entries( + &a[1][0], &a[1][1], &a[1][2], &a[2][0], &a[2][1], &a[2][2], &a[3][0], &a[3][1], + &a[3][2], + ); + det -= &a[0][3] * c03; + } + + det +} + /// Outcome of a Bareiss forward-elimination pass. #[derive(Debug)] enum BareissResult { @@ -402,7 +634,7 @@ fn determinant_scale_exp(e_min: i32) -> Result { Ok(total_exp) } -/// Compute the exact determinant using integer-only Bareiss elimination. +/// Compute the exact determinant from integer-scaled entries. /// /// Returns `(det_int, scale_exp)` where the true determinant is /// `det_int × 2^scale_exp`. Since the scale factor `2^scale_exp` is always @@ -410,10 +642,11 @@ fn determinant_scale_exp(e_min: i32) -> Result { /// /// All arithmetic is in `BigInt` — no `BigRational`, no GCD, no denominator /// tracking. Each f64 entry is decomposed into `mantissa × 2^exponent` and -/// scaled to a common base `2^e_min` so every entry becomes an integer. -/// The Bareiss inner-loop division is exact (guaranteed by the algorithm). +/// scaled to a common base `2^e_min` so every entry becomes an integer. D≤4 +/// uses direct determinant expansions; larger matrices use Bareiss elimination +/// whose inner-loop division is exact (guaranteed by the algorithm). /// -fn bareiss_det_int_finite(m: &FiniteMatrix) -> Result<(BigInt, i32), LaError> { +fn bareiss_det_int_finite(m: &Matrix) -> Result<(BigInt, i32), LaError> { // D == 0 has no `a[D-1][D-1]` to read; shortcut to the empty-product // determinant. if D == 0 { @@ -428,18 +661,26 @@ fn bareiss_det_int_finite(m: &FiniteMatrix) -> Result<(BigInt } let mut a = build_bigint_matrix(&components, e_min); - let sign = match bareiss_forward_eliminate(&mut a, None) { - BareissResult::Upper { sign } => sign, - BareissResult::Singular { .. } => { - cold_path(); - return Ok((BigInt::from(0), 0)); - } - }; + let det_int = match D { + 1 => a[0][0].clone(), + 2 => det2_bigint(&a), + 3 => det3_bigint(&a), + 4 => det4_bigint(&a), + _ => { + let sign = match bareiss_forward_eliminate(&mut a, None) { + BareissResult::Upper { sign } => sign, + BareissResult::Singular { .. } => { + cold_path(); + return Ok((BigInt::from(0), 0)); + } + }; - let det_int = if sign < 0 { - -&a[D - 1][D - 1] - } else { - a[D - 1][D - 1].clone() + if sign < 0 { + -&a[D - 1][D - 1] + } else { + a[D - 1][D - 1].clone() + } + } }; let total_exp = determinant_scale_exp::(e_min)?; @@ -448,7 +689,7 @@ fn bareiss_det_int_finite(m: &FiniteMatrix) -> Result<(BigInt /// Compute the exact determinant of a `D×D` matrix using integer-only Bareiss /// elimination and return the result as a `BigRational`. -fn bareiss_det_finite(m: &FiniteMatrix) -> Result { +fn bareiss_det_finite(m: &Matrix) -> Result { let (det_int, total_exp) = bareiss_det_int_finite(m)?; Ok(bigint_exp_to_bigrational(det_int, total_exp)) } @@ -462,8 +703,8 @@ fn bareiss_det_finite(m: &FiniteMatrix) -> Result( - m: &FiniteMatrix, - b: &FiniteVector, + m: &Matrix, + b: &Vector, ) -> Result<[BigRational; D], LaError> { let (m_components, m_e_min) = decompose_finite_matrix(m)?; let (b_components, b_e_min) = decompose_finite_vec(b)?; @@ -499,10 +740,12 @@ fn gauss_solve_components( let mut a = build_bigint_matrix(&m_components, e_min); let mut rhs = build_bigint_vec(&b_components, e_min); - if let BareissResult::Singular { pivot_col } = bareiss_forward_eliminate(&mut a, Some(&mut rhs)) - { - cold_path(); - return Err(LaError::Singular { pivot_col }); + match bareiss_forward_eliminate(&mut a, Some(&mut rhs)) { + BareissResult::Upper { .. } => {} + BareissResult::Singular { pivot_col } => { + cold_path(); + return Err(LaError::Singular { pivot_col }); + } } let mut x: [BigRational; D] = from_fn(|_| BigRational::from_integer(BigInt::from(0))); @@ -519,97 +762,159 @@ fn gauss_solve_components( Ok(x) } -impl FiniteMatrix { - /// Exact determinant for an already finite matrix. - #[inline] - fn det_exact(&self) -> Result { - bareiss_det_finite(self) - } +/// Exact determinant for a finite-by-construction matrix. +/// +/// This is the private implementation target for [`Matrix::det_exact`]. Keeping +/// the helper separate from the public method keeps the exact core focused on +/// the Bareiss computation while relying on the public [`Matrix`] finite-storage +/// invariant. +/// +/// # Errors +/// Returns [`LaError::DeterminantScaleOverflow`] if determinant scaling +/// overflows the internal exponent representation. +#[inline] +fn det_exact_finite(m: &Matrix) -> Result { + bareiss_det_finite(m) +} - /// Exact determinant converted to a finite `f64`. - /// - /// # Errors - /// Returns [`LaError::DeterminantScaleOverflow`] if determinant scaling - /// overflows the internal exponent representation. - /// - /// Returns [`LaError::Overflow`] if the exact determinant cannot be - /// represented as a finite `f64`. - #[inline] - fn det_exact_f64(&self) -> Result { - let (det_int, total_exp) = bareiss_det_int_finite(self)?; - bigint_exp_to_finite_f64(det_int, total_exp) - } +/// Exact determinant converted to finite `f64` without rounding. +/// +/// This preserves the strict contract of [`Matrix::det_exact_f64`]: if the exact +/// determinant is not representable as a finite binary64 value, callers receive +/// a typed [`LaError::Unrepresentable`] instead of a rounded result. +/// +/// # Errors +/// Returns [`LaError::DeterminantScaleOverflow`] if determinant scaling +/// overflows the internal exponent representation. +/// +/// Returns [`LaError::Unrepresentable`] with +/// [`UnrepresentableReason::RequiresRounding`] when the exact determinant is +/// finite but not exactly representable as binary64, or +/// [`UnrepresentableReason::NotFinite`] when no finite `f64` can represent it. +#[inline] +fn det_exact_f64_finite(m: &Matrix) -> Result { + let (det_int, total_exp) = bareiss_det_int_finite(m)?; + bigint_exp_to_finite_f64(det_int, total_exp) +} - /// Exact linear solve for finite inputs. - /// - /// # Errors - /// Returns [`LaError::Singular`] if the matrix is exactly singular. - #[inline] - fn solve_exact(&self, b: FiniteVector) -> Result<[BigRational; D], LaError> { - gauss_solve_finite(self, &b) - } +/// Exact determinant rounded to finite `f64`. +/// +/// This is the intentionally lossy counterpart to [`det_exact_f64_finite`] and +/// the private implementation target for [`Matrix::det_exact_rounded_f64`]. +/// +/// # Errors +/// Returns [`LaError::DeterminantScaleOverflow`] if determinant scaling +/// overflows the internal exponent representation. +/// +/// Returns [`LaError::Unrepresentable`] with +/// [`UnrepresentableReason::NotFinite`] if the rounded result would be NaN or +/// infinity. +#[inline] +fn det_exact_rounded_f64_finite(m: &Matrix) -> Result { + let (det_int, total_exp) = bareiss_det_int_finite(m)?; + bigint_exp_to_rounded_f64(det_int, total_exp) +} - /// Exact linear solve converted to a finite `f64` vector. - /// - /// # Errors - /// Returns [`LaError::Singular`] if the matrix is exactly singular. - /// Returns [`LaError::Overflow`] if any exact solution component cannot be - /// represented as a finite `f64`. - #[inline] - fn solve_exact_f64(&self, b: FiniteVector) -> Result, LaError> { - let exact = self.solve_exact(b)?; - let mut result = [0.0f64; D]; - for (i, val) in exact.iter().enumerate() { - let Some(f) = val.to_f64() else { - cold_path(); - return Err(LaError::Overflow { index: Some(i) }); - }; - if !f.is_finite() { - cold_path(); - return Err(LaError::Overflow { index: Some(i) }); - } - result[i] = f; - } - Ok(FiniteVector::new_unchecked(Vector::new_unchecked(result))) - } +/// Exact linear solve for finite inputs. +/// +/// This is the private implementation target for [`Matrix::solve_exact`]. +/// Public [`Matrix`] and [`Vector`] values are finite by construction, so this +/// helper can focus on the exact Bareiss/rational solve. +/// +/// # Errors +/// Returns [`LaError::Singular`] if the matrix is exactly singular. +#[inline] +fn solve_exact_finite( + m: &Matrix, + b: Vector, +) -> Result<[BigRational; D], LaError> { + gauss_solve_finite(m, &b) +} - /// Exact determinant sign for an already finite matrix. - /// - /// The fast `f64` filter may reject overflowed scalar intermediates as - /// inconclusive, then fall back to integer Bareiss sign computation. - /// - /// # Errors - /// Returns [`LaError::NonFinite`] if a direct determinant or error-bound - /// computation detects a non-finite condition that is not an inconclusive - /// scalar overflow. - /// - /// Returns [`LaError::DeterminantScaleOverflow`] if determinant scaling - /// overflows the internal exponent representation. - #[inline] - fn det_sign_exact(&self) -> Result { - match (self.det_direct(), self.det_errbound()) { - (Ok(Some(det_f64)), Ok(Some(err))) => { - if det_f64 > err { - return Ok(1); - } - if det_f64 < -err { - return Ok(-1); - } +/// Exact linear solve converted to finite `f64` components without rounding. +/// +/// This preserves the strict contract of [`Matrix::solve_exact_f64`]: each exact +/// component must already be representable as finite binary64. +/// +/// # Errors +/// Returns [`LaError::Singular`] if the matrix is exactly singular. +/// +/// Returns [`LaError::Unrepresentable`] with the failing component index when an +/// exact solution component requires rounding or cannot be represented as a +/// finite `f64`. +#[inline] +fn solve_exact_f64_finite( + m: &Matrix, + b: Vector, +) -> Result, LaError> { + let exact = solve_exact_finite(m, b)?; + let mut result = [0.0f64; D]; + for (i, val) in exact.iter().enumerate() { + result[i] = exact_rational_to_finite_f64(val, Some(i))?; + } + Ok(Vector::new_unchecked(result)) +} + +/// Exact linear solve rounded to finite `f64` components. +/// +/// This is the intentionally lossy counterpart to [`solve_exact_f64_finite`] and +/// the private implementation target for [`Matrix::solve_exact_rounded_f64`]. +/// +/// # Errors +/// Returns [`LaError::Singular`] if the matrix is exactly singular. +/// +/// Returns [`LaError::Unrepresentable`] with +/// [`UnrepresentableReason::NotFinite`] if any rounded component would be NaN or +/// infinity. +#[inline] +fn solve_exact_rounded_f64_finite( + m: &Matrix, + b: Vector, +) -> Result, LaError> { + let exact = solve_exact_finite(m, b)?; + let mut result = [0.0f64; D]; + for (i, val) in exact.iter().enumerate() { + result[i] = exact_rational_to_rounded_f64(val, Some(i))?; + } + Ok(Vector::new_unchecked(result)) +} + +/// Exact determinant sign for an already finite matrix. +/// +/// The fast `f64` filter may reject overflowed scalar intermediates as +/// inconclusive, then fall back to integer Bareiss sign computation. +/// +/// # Errors +/// Returns [`LaError::NonFinite`] if a direct determinant or error-bound +/// computation detects a non-finite condition that is not an inconclusive scalar +/// overflow. +/// +/// Returns [`LaError::DeterminantScaleOverflow`] if determinant scaling +/// overflows the internal exponent representation. +#[inline] +fn det_sign_exact_finite(m: &Matrix) -> Result { + match (m.det_direct(), m.det_errbound()) { + (Ok(Some(det_f64)), Ok(Some(err))) => { + if det_f64 > err { + return Ok(1); + } + if det_f64 < -err { + return Ok(-1); } - (Err(LaError::NonFinite { row: None, .. }), _) - | (_, Err(LaError::NonFinite { row: None, .. })) => {} - (Err(err), _) | (_, Err(err)) => return Err(err), - _ => {} } - - cold_path(); - let (det_int, _) = bareiss_det_int_finite(self)?; - Ok(match det_int.sign() { - Sign::Plus => 1, - Sign::Minus => -1, - Sign::NoSign => 0, - }) - } + (Err(LaError::NonFinite { row: None, .. }), _) + | (_, Err(LaError::NonFinite { row: None, .. })) => {} + (Err(err), _) | (_, Err(err)) => return Err(err), + _ => {} + } + + cold_path(); + let (det_int, _) = bareiss_det_int_finite(m)?; + Ok(match det_int.sign() { + Sign::Plus => 1, + Sign::Minus => -1, + Sign::NoSign => 0, + }) } impl Matrix { @@ -646,7 +951,7 @@ impl Matrix { /// overflows the internal exponent representation. #[inline] pub fn det_exact(&self) -> Result { - FiniteMatrix::new_unchecked(*self).det_exact() + det_exact_finite(self) } /// Exact determinant converted to `f64`. @@ -655,10 +960,8 @@ impl Matrix { /// /// Computes the exact determinant with the same integer Bareiss core used by /// [`det_exact`](Self::det_exact), then converts the exact scaled integer - /// result to the nearest `f64` without first materializing the public - /// [`BigRational`] determinant. This is useful when you want the most accurate - /// f64 determinant possible without committing to `BigRational` in your - /// downstream code. + /// result to `f64` only if the result is exactly representable as a finite + /// binary64 value. /// /// # Examples /// ``` @@ -676,11 +979,54 @@ impl Matrix { /// Returns [`LaError::DeterminantScaleOverflow`] if determinant scaling /// overflows the internal exponent representation. /// - /// Returns [`LaError::Overflow`] if the exact determinant is too large to - /// represent as a finite `f64`. + /// Returns [`LaError::Unrepresentable`] if the exact determinant cannot be + /// represented exactly as a finite `f64`. #[inline] pub fn det_exact_f64(&self) -> Result { - FiniteMatrix::new_unchecked(*self).det_exact_f64() + det_exact_f64_finite(self) + } + + /// Exact determinant rounded to `f64`. + /// + /// Requires the `exact` Cargo feature. + /// + /// Computes the exact determinant with the same integer Bareiss core used by + /// [`det_exact`](Self::det_exact), then rounds the exact value to a finite + /// binary64 value. Unlike [`det_exact_f64`](Self::det_exact_f64), this method + /// is intentionally lossy and may round non-dyadic or underflowing nonzero + /// exact determinants. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// # fn main() -> Result<(), LaError> { + /// let m = Matrix::<2>::try_from_rows([ + /// [1.0 + f64::EPSILON, 0.0], + /// [0.0, 1.0 - f64::EPSILON], + /// ])?; + /// + /// assert_eq!( + /// m.det_exact_f64(), + /// Err(LaError::Unrepresentable { + /// index: None, + /// reason: UnrepresentableReason::RequiresRounding, + /// }) + /// ); + /// assert_eq!(m.det_exact_rounded_f64()?.to_bits(), 1.0f64.to_bits()); + /// # Ok(()) + /// # } + /// ``` + /// + /// # Errors + /// Returns [`LaError::DeterminantScaleOverflow`] if determinant scaling + /// overflows the internal exponent representation. + /// + /// Returns [`LaError::Unrepresentable`] if the rounded determinant would be + /// NaN or infinite. + #[inline] + pub fn det_exact_rounded_f64(&self) -> Result { + det_exact_rounded_f64_finite(self) } /// Exact linear system solve using hybrid integer/rational arithmetic. @@ -729,9 +1075,7 @@ impl Matrix { /// Returns [`LaError::Singular`] if the matrix is exactly singular. #[inline] pub fn solve_exact(&self, b: Vector) -> Result<[BigRational; D], LaError> { - let finite_m = FiniteMatrix::new_unchecked(*self); - let finite_b = FiniteVector::new_unchecked(b); - finite_m.solve_exact(finite_b) + solve_exact_finite(self, b) } /// Exact linear system solve converted to `f64`. @@ -739,9 +1083,9 @@ impl Matrix { /// Requires the `exact` Cargo feature. /// /// Computes the exact [`BigRational`] solution via - /// [`solve_exact`](Self::solve_exact) and converts each component to the - /// nearest `f64`. This is useful when you want the most accurate f64 - /// solution possible without committing to `BigRational` downstream. + /// [`solve_exact`](Self::solve_exact) and converts each component to `f64` + /// only if that component is exactly representable as a finite binary64 + /// value. /// /// # Examples /// ``` @@ -759,13 +1103,50 @@ impl Matrix { /// /// # Errors /// Returns [`LaError::Singular`] if the matrix is exactly singular. - /// Returns [`LaError::Overflow`] if any component of the exact solution is - /// too large to represent as a finite `f64`. + /// Returns [`LaError::Unrepresentable`] if any component of the exact solution + /// cannot be represented exactly as a finite `f64`. #[inline] pub fn solve_exact_f64(&self, b: Vector) -> Result, LaError> { - let finite_m = FiniteMatrix::new_unchecked(*self); - let finite_b = FiniteVector::new_unchecked(b); - Ok(finite_m.solve_exact_f64(finite_b)?.into_vector()) + solve_exact_f64_finite(self, b) + } + + /// Exact linear system solve rounded to `f64`. + /// + /// Requires the `exact` Cargo feature. + /// + /// Computes the exact [`BigRational`] solution via + /// [`solve_exact`](Self::solve_exact) and rounds each component to a finite + /// binary64 value. Unlike [`solve_exact_f64`](Self::solve_exact_f64), this + /// method is intentionally lossy and may round non-dyadic or underflowing + /// nonzero exact components. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// # fn main() -> Result<(), LaError> { + /// let a = Matrix::<1>::try_from_rows([[3.0]])?; + /// let b = Vector::<1>::try_new([1.0])?; + /// + /// assert_eq!( + /// a.solve_exact_f64(b), + /// Err(LaError::Unrepresentable { + /// index: Some(0), + /// reason: UnrepresentableReason::RequiresRounding, + /// }) + /// ); + /// assert_eq!(a.solve_exact_rounded_f64(b)?.into_array(), [1.0 / 3.0]); + /// # Ok(()) + /// # } + /// ``` + /// + /// # Errors + /// Returns [`LaError::Singular`] if the matrix is exactly singular. + /// Returns [`LaError::Unrepresentable`] if any rounded component would be + /// NaN or infinite. + #[inline] + pub fn solve_exact_rounded_f64(&self, b: Vector) -> Result, LaError> { + solve_exact_rounded_f64_finite(self, b) } /// Exact determinant sign using adaptive-precision arithmetic. @@ -809,7 +1190,7 @@ impl Matrix { /// overflows the internal exponent representation. #[inline] pub fn det_sign_exact(&self) -> Result { - FiniteMatrix::new_unchecked(*self).det_sign_exact() + det_sign_exact_finite(self) } } @@ -866,22 +1247,22 @@ mod tests { // Macro-generated per-dimension tests (D=2..5) // ----------------------------------------------------------------------- - macro_rules! gen_internal_finite_exact_tests { + macro_rules! gen_internal_matrix_exact_tests { ($d:literal) => { paste! { #[test] - fn []() { - let a = FiniteMatrix::<$d>::new_unchecked(Matrix::<$d>::identity()); - let b = FiniteVector::<$d>::new_unchecked(Vector::<$d>::new([1.0; $d])); + fn []() { + let a = Matrix::<$d>::identity(); + let b = Vector::<$d>::new([1.0; $d]); assert_eq!( - a.det_exact().unwrap(), + det_exact_finite(&a).unwrap(), BigRational::from_integer(BigInt::from(1)) ); - assert!((a.det_exact_f64().unwrap() - 1.0).abs() <= f64::EPSILON); - assert_eq!(a.det_sign_exact().unwrap(), 1); + assert!((det_exact_f64_finite(&a).unwrap() - 1.0).abs() <= f64::EPSILON); + assert_eq!(det_sign_exact_finite(&a).unwrap(), 1); - let exact = a.solve_exact(b).unwrap(); + let exact = solve_exact_finite(&a, b).unwrap(); for value in exact { assert_eq!(value, BigRational::from_integer(BigInt::from(1))); } @@ -895,10 +1276,10 @@ mod tests { }; } - gen_internal_finite_exact_tests!(2); - gen_internal_finite_exact_tests!(3); - gen_internal_finite_exact_tests!(4); - gen_internal_finite_exact_tests!(5); + gen_internal_matrix_exact_tests!(2); + gen_internal_matrix_exact_tests!(3); + gen_internal_matrix_exact_tests!(4); + gen_internal_matrix_exact_tests!(5); macro_rules! gen_det_exact_tests { ($d:literal) => { @@ -934,30 +1315,25 @@ mod tests { gen_det_exact_f64_tests!(4); gen_det_exact_f64_tests!(5); - /// For D ≤ 4, `det_exact_f64` should agree with `det_direct` on - /// well-conditioned matrices. + /// For D ≤ 4, `det_exact_f64` should agree with `det_direct` on matrices + /// whose exact determinant is representable in f64. macro_rules! gen_det_exact_f64_agrees_with_det_direct { ($d:literal) => { paste! { #[test] - #[allow(clippy::cast_precision_loss)] fn []() { - // Diagonally dominant → well-conditioned. + // Power-of-two diagonal entries make the determinant + // exactly representable in binary64. let mut rows = [[0.0f64; $d]; $d]; - for r in 0..$d { - for c in 0..$d { - rows[r][c] = if r == c { - (r as f64) + f64::from($d) + 1.0 - } else { - 0.1 / ((r + c + 1) as f64) - }; - } + let mut value = 2.0; + for (i, row) in rows.iter_mut().enumerate() { + row[i] = value; + value *= 2.0; } - let m = Matrix::<$d>::from_rows(rows); + let m = Matrix::<$d>::try_from_rows(rows).unwrap(); let exact = m.det_exact_f64().unwrap(); let direct = m.det_direct().unwrap().unwrap(); - let eps = direct.abs().mul_add(1e-12, 1e-12); - assert!((exact - direct).abs() <= eps); + assert_eq!(exact.to_bits(), direct.to_bits()); } } }; @@ -974,19 +1350,19 @@ mod tests { #[test] fn det_sign_exact_d1_positive() { - let m = Matrix::<1>::from_rows([[42.0]]); + let m = Matrix::<1>::try_from_rows([[42.0]]).unwrap(); assert_eq!(m.det_sign_exact().unwrap(), 1); } #[test] fn det_sign_exact_d1_negative() { - let m = Matrix::<1>::from_rows([[-3.5]]); + let m = Matrix::<1>::try_from_rows([[-3.5]]).unwrap(); assert_eq!(m.det_sign_exact().unwrap(), -1); } #[test] fn det_sign_exact_d1_zero() { - let m = Matrix::<1>::from_rows([[0.0]]); + let m = Matrix::<1>::try_from_rows([[0.0]]).unwrap(); assert_eq!(m.det_sign_exact().unwrap(), 0); } @@ -1012,39 +1388,43 @@ mod tests { #[test] fn det_sign_exact_singular_duplicate_rows() { - let m = Matrix::<3>::from_rows([ + let m = Matrix::<3>::try_from_rows([ [1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [1.0, 2.0, 3.0], // duplicate of row 0 - ]); + ]) + .unwrap(); assert_eq!(m.det_sign_exact().unwrap(), 0); } #[test] fn det_sign_exact_singular_linear_combination() { // Row 2 = row 0 + row 1 in exact arithmetic. - let m = Matrix::<3>::from_rows([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [5.0, 7.0, 9.0]]); + let m = Matrix::<3>::try_from_rows([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [5.0, 7.0, 9.0]]) + .unwrap(); assert_eq!(m.det_sign_exact().unwrap(), 0); } #[test] fn det_sign_exact_negative_det_row_swap() { // Swapping two rows of the identity negates the determinant. - let m = Matrix::<3>::from_rows([[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]); + let m = Matrix::<3>::try_from_rows([[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) + .unwrap(); assert_eq!(m.det_sign_exact().unwrap(), -1); } #[test] fn det_sign_exact_negative_det_known() { // det([[1,2],[3,4]]) = 1*4 - 2*3 = -2 - let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); + let m = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]]).unwrap(); assert_eq!(m.det_sign_exact().unwrap(), -1); } #[test] fn det_sign_exact_agrees_with_det_for_spd() { // SPD matrix → positive determinant. - let m = Matrix::<3>::from_rows([[4.0, 2.0, 0.0], [2.0, 5.0, 1.0], [0.0, 1.0, 3.0]]); + let m = Matrix::<3>::try_from_rows([[4.0, 2.0, 0.0], [2.0, 5.0, 1.0], [0.0, 1.0, 3.0]]) + .unwrap(); assert_eq!(m.det_sign_exact().unwrap(), 1); assert!(m.det().unwrap() > 0.0); } @@ -1058,11 +1438,12 @@ mod tests { #[test] fn det_sign_exact_near_singular_perturbation() { let perturbation = f64::from_bits(0x3CD0_0000_0000_0000); // 2^-50 - let m = Matrix::<3>::from_rows([ + let m = Matrix::<3>::try_from_rows([ [1.0 + perturbation, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0], - ]); + ]) + .unwrap(); // Exact: det = perturbation × (5×9 − 6×8) = perturbation × (−3) < 0. assert_eq!(m.det_sign_exact().unwrap(), -1); } @@ -1072,12 +1453,13 @@ mod tests { /// but we verify correctness for a range of known signs. #[test] fn det_sign_exact_fast_filter_positive_4x4() { - let m = Matrix::<4>::from_rows([ + let m = Matrix::<4>::try_from_rows([ [2.0, 1.0, 0.0, 0.0], [1.0, 3.0, 1.0, 0.0], [0.0, 1.0, 4.0, 1.0], [0.0, 0.0, 1.0, 5.0], - ]); + ]) + .unwrap(); // SPD tridiagonal → positive det. assert_eq!(m.det_sign_exact().unwrap(), 1); } @@ -1085,12 +1467,13 @@ mod tests { #[test] fn det_sign_exact_fast_filter_negative_4x4() { // Swap rows 0 and 1 of the above → negate det. - let m = Matrix::<4>::from_rows([ + let m = Matrix::<4>::try_from_rows([ [1.0, 3.0, 1.0, 0.0], [2.0, 1.0, 0.0, 0.0], [0.0, 1.0, 4.0, 1.0], [0.0, 0.0, 1.0, 5.0], - ]); + ]) + .unwrap(); assert_eq!(m.det_sign_exact().unwrap(), -1); } @@ -1100,7 +1483,7 @@ mod tests { let tiny = 5e-324_f64; // smallest positive subnormal assert!(tiny.is_subnormal()); - let m = Matrix::<2>::from_rows([[tiny, 0.0], [0.0, tiny]]); + let m = Matrix::<2>::try_from_rows([[tiny, 0.0], [0.0, tiny]]).unwrap(); // det = tiny^2 > 0 assert_eq!(m.det_sign_exact().unwrap(), 1); } @@ -1127,6 +1510,35 @@ mod tests { ); } + #[test] + fn exact_public_methods_reject_unchecked_nonfinite_matrix_before_computation() { + let m = Matrix::<2>::from_rows_unchecked([[1.0, 0.0], [f64::NAN, 1.0]]); + let b = Vector::<2>::new([1.0, 1.0]); + let expected = Err(LaError::NonFinite { + row: Some(1), + col: 0, + }); + + assert_eq!(m.det_exact().map(|_| ()), expected); + assert_eq!(m.det_exact_f64().map(|_| ()), expected); + assert_eq!(m.det_exact_rounded_f64().map(|_| ()), expected); + assert_eq!(m.det_sign_exact().map(|_| ()), expected); + assert_eq!(m.solve_exact(b).map(|_| ()), expected); + assert_eq!(m.solve_exact_f64(b).map(|_| ()), expected); + assert_eq!(m.solve_exact_rounded_f64(b).map(|_| ()), expected); + } + + #[test] + fn exact_solve_public_methods_reject_unchecked_nonfinite_rhs_before_computation() { + let a = Matrix::<2>::identity(); + let b = Vector::<2>::new_unchecked([1.0, f64::INFINITY]); + let expected = Err(LaError::NonFinite { row: None, col: 1 }); + + assert_eq!(a.solve_exact(b).map(|_| ()), expected); + assert_eq!(a.solve_exact_f64(b).map(|_| ()), expected); + assert_eq!(a.solve_exact_rounded_f64(b).map(|_| ()), expected); + } + #[test] fn det_sign_exact_returns_err_on_nan_5x5() { // D ≥ 5 bypasses the fast filter, exercising the bareiss_det path. @@ -1156,26 +1568,28 @@ mod tests { fn det_sign_exact_pivot_needed_5x5() { // D ≥ 5 skips the fast filter → exercises Bareiss pivoting. // Permutation matrix with a single swap (rows 0↔1) → det = −1. - let m = Matrix::<5>::from_rows([ + let m = Matrix::<5>::try_from_rows([ [0.0, 1.0, 0.0, 0.0, 0.0], [1.0, 0.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, 0.0], [0.0, 0.0, 0.0, 0.0, 1.0], - ]); + ]) + .unwrap(); assert_eq!(m.det_sign_exact().unwrap(), -1); } #[test] fn det_sign_exact_5x5_known() { // det of a permutation matrix with two swaps = +1 (even permutation). - let m = Matrix::<5>::from_rows([ + let m = Matrix::<5>::try_from_rows([ [0.0, 1.0, 0.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 1.0], - ]); + ]) + .unwrap(); // Two transpositions → even permutation → det = +1 assert_eq!(m.det_sign_exact().unwrap(), 1); } @@ -1192,14 +1606,14 @@ mod tests { #[test] fn det_errbound_d1_is_zero() { assert_eq!( - Matrix::<1>::from_rows([[42.0]]).det_errbound(), + Matrix::<1>::try_from_rows([[42.0]]).unwrap().det_errbound(), Ok(Some(0.0)) ); } #[test] fn det_errbound_d2_positive() { - let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); + let m = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]]).unwrap(); let bound = m.det_errbound().unwrap().unwrap(); assert!(bound > 0.0); // bound = ERR_COEFF_2 * (|1*4| + |2*3|) = ERR_COEFF_2 * 10 @@ -1216,7 +1630,8 @@ mod tests { #[test] fn det_errbound_d3_non_identity() { // Non-identity matrix to exercise all code paths in D=3 case - let m = Matrix::<3>::from_rows([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 10.0]]); + let m = Matrix::<3>::try_from_rows([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 10.0]]) + .unwrap(); let bound = m.det_errbound().unwrap().unwrap(); assert!(bound > 0.0); } @@ -1231,12 +1646,13 @@ mod tests { #[test] fn det_errbound_d4_non_identity() { // Non-identity matrix to exercise all code paths in D=4 case - let m = Matrix::<4>::from_rows([ + let m = Matrix::<4>::try_from_rows([ [1.0, 0.0, 0.0, 0.0], [0.0, 2.0, 0.0, 0.0], [0.0, 0.0, 3.0, 0.0], [0.0, 0.0, 0.0, 4.0], - ]); + ]) + .unwrap(); let bound = m.det_errbound().unwrap().unwrap(); assert!(bound > 0.0); } @@ -1354,7 +1770,7 @@ mod tests { #[test] fn bareiss_det_int_d0() { - let m = FiniteMatrix::new_unchecked(Matrix::<0>::zero()); + let m = Matrix::<0>::zero(); let (det, exp) = bareiss_det_int_finite(&m).unwrap(); assert_eq!(det, BigInt::from(1)); assert_eq!(exp, 0); @@ -1375,7 +1791,7 @@ mod tests { (0.5, 1, -1), // 0.5 = 1 × 2^(-1) ]; for &(input, expected_det_int, expected_exp) in cases { - let m = FiniteMatrix::new_unchecked(Matrix::<1>::from_rows([[input]])); + let m = Matrix::<1>::try_from_rows([[input]]).unwrap(); let (det, exp) = bareiss_det_int_finite(&m).unwrap(); assert_eq!( det, @@ -1389,7 +1805,7 @@ mod tests { #[test] fn bareiss_det_int_d2_known() { // det([[1,2],[3,4]]) = -2 - let m = FiniteMatrix::new_unchecked(Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]])); + let m = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]]).unwrap(); let (det_int, total_exp) = bareiss_det_int_finite(&m).unwrap(); // Reconstruct and verify. let det = bigint_exp_to_bigrational(det_int, total_exp); @@ -1398,7 +1814,7 @@ mod tests { #[test] fn bareiss_det_int_all_zeros() { - let m = FiniteMatrix::new_unchecked(Matrix::<3>::zero()); + let m = Matrix::<3>::zero(); let (det, _) = bareiss_det_int_finite(&m).unwrap(); assert_eq!(det, BigInt::from(0)); } @@ -1406,11 +1822,8 @@ mod tests { #[test] fn bareiss_det_int_sign_matches_det_sign_exact() { // The sign of det_int should match det_sign_exact for various matrices. - let m = FiniteMatrix::new_unchecked(Matrix::<3>::from_rows([ - [0.0, 1.0, 0.0], - [1.0, 0.0, 0.0], - [0.0, 0.0, 1.0], - ])); + let m = Matrix::<3>::try_from_rows([[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) + .unwrap(); let (det_int, _) = bareiss_det_int_finite(&m).unwrap(); assert_eq!(det_int.sign(), Sign::Minus); // det = -1 } @@ -1419,7 +1832,7 @@ mod tests { fn bareiss_det_int_fractional_entries() { // Entries with negative exponents: 0.5 = 1×2^(-1), 0.25 = 1×2^(-2). // det([[0.5, 0.25], [1.0, 1.0]]) = 0.5×1.0 − 0.25×1.0 = 0.25 - let m = FiniteMatrix::new_unchecked(Matrix::<2>::from_rows([[0.5, 0.25], [1.0, 1.0]])); + let m = Matrix::<2>::try_from_rows([[0.5, 0.25], [1.0, 1.0]]).unwrap(); let (det_int, total_exp) = bareiss_det_int_finite(&m).unwrap(); let det = bigint_exp_to_bigrational(det_int, total_exp); assert_eq!(det, BigRational::new(BigInt::from(1), BigInt::from(4))); @@ -1428,11 +1841,8 @@ mod tests { #[test] fn bareiss_det_int_d3_with_pivoting() { // Zero on diagonal → exercises pivot swap inside bareiss_det_int. - let m = FiniteMatrix::new_unchecked(Matrix::<3>::from_rows([ - [0.0, 1.0, 0.0], - [1.0, 0.0, 0.0], - [0.0, 0.0, 1.0], - ])); + let m = Matrix::<3>::try_from_rows([[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) + .unwrap(); let (det_int, total_exp) = bareiss_det_int_finite(&m).unwrap(); let det = bigint_exp_to_bigrational(det_int, total_exp); assert_eq!(det, BigRational::from_integer(BigInt::from(-1))); @@ -1444,7 +1854,7 @@ mod tests { paste! { #[test] fn []() { - let m = FiniteMatrix::new_unchecked(Matrix::<$d>::identity()); + let m = Matrix::<$d>::identity(); let (det_int, total_exp) = bareiss_det_int_finite(&m).unwrap(); let det = bigint_exp_to_bigrational(det_int, total_exp); assert_eq!(det, BigRational::from_integer(BigInt::from(1))); @@ -1525,14 +1935,18 @@ mod tests { #[test] fn bareiss_det_d1_returns_entry() { - let det = Matrix::<1>::from_rows([[7.0]]).det_exact().unwrap(); + let det = Matrix::<1>::try_from_rows([[7.0]]) + .unwrap() + .det_exact() + .unwrap(); assert_eq!(det, f64_to_bigrational(7.0)); } #[test] fn bareiss_det_d3_with_pivoting() { // First column has zero on diagonal → exercises pivot swap + break. - let m = Matrix::<3>::from_rows([[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]); + let m = Matrix::<3>::try_from_rows([[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) + .unwrap(); let det = m.det_exact().unwrap(); // det of this permutation matrix = -1 assert_eq!(det, BigRational::from_integer(BigInt::from(-1))); @@ -1541,7 +1955,8 @@ mod tests { #[test] fn bareiss_det_singular_all_zeros_in_column() { // Column 1 is all zeros below diagonal after elimination → singular. - let m = Matrix::<3>::from_rows([[1.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]); + let m = Matrix::<3>::try_from_rows([[1.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) + .unwrap(); let det = m.det_exact().unwrap(); assert_eq!(det, BigRational::from_integer(BigInt::from(0))); } @@ -1553,7 +1968,8 @@ mod tests { // compute the correct positive sign. let big = f64::MAX / 2.0; assert!(big.is_finite()); - let m = Matrix::<3>::from_rows([[0.0, 0.0, 1.0], [big, 0.0, 1.0], [0.0, big, 1.0]]); + let m = Matrix::<3>::try_from_rows([[0.0, 0.0, 1.0], [big, 0.0, 1.0], [0.0, big, 1.0]]) + .unwrap(); // det = big^2 > 0 assert_eq!(m.det_sign_exact().unwrap(), 1); } @@ -1571,15 +1987,32 @@ mod tests { #[test] fn det_exact_known_2x2() { // det([[1,2],[3,4]]) = 1*4 - 2*3 = -2 - let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); + let m = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]]).unwrap(); let det = m.det_exact().unwrap(); assert_eq!(det, BigRational::from_integer(BigInt::from(-2))); } + #[test] + fn det_exact_known_dense_4x4() { + let m = Matrix::<4>::try_from_rows([ + [4.0, 1.0, 3.0, 2.0], + [0.0, 5.0, 2.0, 1.0], + [7.0, 2.0, 6.0, 3.0], + [1.0, 8.0, 4.0, 9.0], + ]) + .unwrap(); + + assert_eq!( + m.det_exact(), + Ok(BigRational::from_integer(BigInt::from(92))) + ); + } + #[test] fn det_exact_singular_returns_zero() { // Rows in arithmetic progression → exactly singular. - let m = Matrix::<3>::from_rows([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]); + 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]]) + .unwrap(); let det = m.det_exact().unwrap(); assert_eq!(det, BigRational::from_integer(BigInt::from(0))); } @@ -1588,11 +2021,12 @@ mod tests { fn det_exact_near_singular_perturbation() { // Same 2^-50 perturbation case: exact det = -3 × 2^-50. let perturbation = f64::from_bits(0x3CD0_0000_0000_0000); // 2^-50 - let m = Matrix::<3>::from_rows([ + let m = Matrix::<3>::try_from_rows([ [1.0 + perturbation, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0], - ]); + ]) + .unwrap(); let det = m.det_exact().unwrap(); // det should be exactly -3 × 2^-50. let expected = BigRational::new(BigInt::from(-3), BigInt::from(1u64 << 50)); @@ -1602,13 +2036,14 @@ mod tests { #[test] fn det_exact_5x5_permutation() { // Single swap (rows 0↔1) → det = -1. - let m = Matrix::<5>::from_rows([ + let m = Matrix::<5>::try_from_rows([ [0.0, 1.0, 0.0, 0.0, 0.0], [1.0, 0.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, 0.0], [0.0, 0.0, 0.0, 0.0, 1.0], - ]); + ]) + .unwrap(); let det = m.det_exact().unwrap(); assert_eq!(det, BigRational::from_integer(BigInt::from(-1))); } @@ -1619,7 +2054,7 @@ mod tests { #[test] fn det_exact_f64_known_2x2() { - let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); + let m = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]]).unwrap(); let det = m.det_exact_f64().unwrap(); assert!((det - (-2.0)).abs() <= f64::EPSILON); } @@ -1628,9 +2063,93 @@ mod tests { fn det_exact_f64_overflow_returns_err() { // Entries near f64::MAX produce a determinant too large for f64. let big = f64::MAX / 2.0; - let m = Matrix::<3>::from_rows([[0.0, 0.0, 1.0], [big, 0.0, 1.0], [0.0, big, 1.0]]); + let m = Matrix::<3>::try_from_rows([[0.0, 0.0, 1.0], [big, 0.0, 1.0], [0.0, big, 1.0]]) + .unwrap(); // det = big^2, which overflows f64. - assert_eq!(m.det_exact_f64(), Err(LaError::Overflow { index: None })); + assert_eq!( + m.det_exact_f64(), + Err(LaError::Unrepresentable { + index: None, + reason: UnrepresentableReason::NotFinite, + }) + ); + } + + #[test] + fn det_exact_rounded_f64_overflow_returns_err() { + let big = f64::MAX / 2.0; + let m = Matrix::<3>::try_from_rows([[0.0, 0.0, 1.0], [big, 0.0, 1.0], [0.0, big, 1.0]]) + .unwrap(); + + assert_eq!( + m.det_exact_rounded_f64(), + Err(LaError::Unrepresentable { + index: None, + reason: UnrepresentableReason::NotFinite, + }) + ); + } + + #[test] + fn det_exact_f64_underflow_returns_err_for_nonzero_exact_result() { + let tiny = f64::from_bits(1); + let m = Matrix::<2>::try_from_rows([[tiny, 0.0], [0.0, tiny]]).unwrap(); + + assert!(m.det_exact().unwrap().is_positive()); + assert_eq!( + m.det_exact_f64(), + Err(LaError::Unrepresentable { + index: None, + reason: UnrepresentableReason::RequiresRounding, + }) + ); + } + + #[test] + fn det_exact_f64_rejects_inexact_rounding() { + let m = Matrix::<2>::try_from_rows([[1.0 + f64::EPSILON, 0.0], [0.0, 1.0 - f64::EPSILON]]) + .unwrap(); + + assert_eq!( + m.det_exact(), + Ok(BigRational::new( + (BigInt::from(1_u128) << 104_u32) - BigInt::from(1), + BigInt::from(1_u128 << 104), + )) + ); + assert_eq!( + m.det_exact_f64(), + Err(LaError::Unrepresentable { + index: None, + reason: UnrepresentableReason::RequiresRounding, + }) + ); + } + + #[test] + fn det_exact_f64_accepts_min_positive_subnormal() { + let tiny = f64::from_bits(1); + let m = Matrix::<1>::try_from_rows([[tiny]]).unwrap(); + + assert_eq!(m.det_exact_f64().unwrap().to_bits(), tiny.to_bits()); + } + + #[test] + fn det_exact_f64_accepts_max_finite_binary64() { + let m = Matrix::<1>::try_from_rows([[f64::MAX]]).unwrap(); + + assert_eq!(m.det_exact_f64().unwrap().to_bits(), f64::MAX.to_bits()); + } + + #[test] + fn det_exact_rounded_f64_rounds_inexact_result() { + let m = Matrix::<2>::try_from_rows([[1.0 + f64::EPSILON, 0.0], [0.0, 1.0 - f64::EPSILON]]) + .unwrap(); + + assert_eq!( + m.det_exact_rounded_f64().unwrap().to_bits(), + 1.0f64.to_bits() + ); } // ----------------------------------------------------------------------- @@ -1703,25 +2222,42 @@ mod tests { ($d:literal) => { paste! { #[test] - #[allow(clippy::cast_precision_loss)] fn []() { - // Diagonally dominant → well-conditioned. + // Diagonally dominant integer matrix with an exactly + // representable target solution. The exact result can + // therefore be parsed into f64 without rounding. let mut rows = [[0.0f64; $d]; $d]; for r in 0..$d { for c in 0..$d { rows[r][c] = if r == c { - (r as f64) + f64::from($d) + 1.0 + f64::from($d) + 1.0 } else { - 0.1 / ((r + c + 1) as f64) + 1.0 }; } } - let a = Matrix::<$d>::from_rows(rows); - let b = arbitrary_rhs::<$d>(); + let a = Matrix::<$d>::try_from_rows(rows).unwrap(); + let x_true = { + let mut arr = [0.0f64; $d]; + for (dst, src) in arr.iter_mut().zip([1.0, -2.0, 3.0, -4.0, 5.0]) { + *dst = src; + } + arr + }; + let mut b_arr = [0.0f64; $d]; + for i in 0..$d { + let mut sum = 0.0; + for j in 0..$d { + sum = rows[i][j].mul_add(x_true[j], sum); + } + b_arr[i] = sum; + } + let b = Vector::<$d>::new(b_arr); let exact = a.solve_exact_f64(b).unwrap().into_array(); let lu_sol = a.lu(DEFAULT_SINGULAR_TOL).unwrap() .solve(b).unwrap().into_array(); for i in 0..$d { + assert_eq!(exact[i].to_bits(), x_true[i].to_bits()); let eps = lu_sol[i].abs().mul_add(1e-12, 1e-12); assert!((exact[i] - lu_sol[i]).abs() <= eps); } @@ -1758,7 +2294,7 @@ mod tests { }; } } - let a = Matrix::<$d>::from_rows(rows); + let a = Matrix::<$d>::try_from_rows(rows).unwrap(); // x0 = [1, 2, ..., D]. let mut x0 = [0.0f64; $d]; @@ -1807,7 +2343,7 @@ mod tests { #[test] fn solve_exact_known_2x2() { // [[1,2],[3,4]] x = [5, 11] → x = [1, 2] - let a = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); + let a = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]]).unwrap(); let b = Vector::<2>::new([5.0, 11.0]); let x = a.solve_exact(b).unwrap(); assert_eq!(x[0], BigRational::from_integer(BigInt::from(1))); @@ -1817,7 +2353,8 @@ mod tests { #[test] fn solve_exact_pivoting_needed() { // First column has zero on diagonal → pivot swap required. - let a = Matrix::<3>::from_rows([[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]); + let a = Matrix::<3>::try_from_rows([[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) + .unwrap(); let b = Vector::<3>::new([2.0, 3.0, 4.0]); let x = a.solve_exact(b).unwrap(); // x = [3, 2, 4] @@ -1829,7 +2366,7 @@ mod tests { #[test] fn solve_exact_fractional_result() { // [[2, 1], [1, 3]] x = [1, 1] → x = [2/5, 1/5] - let a = Matrix::<2>::from_rows([[2.0, 1.0], [1.0, 3.0]]); + let a = Matrix::<2>::try_from_rows([[2.0, 1.0], [1.0, 3.0]]).unwrap(); let b = Vector::<2>::new([1.0, 1.0]); let x = a.solve_exact(b).unwrap(); assert_eq!(x[0], BigRational::new(BigInt::from(2), BigInt::from(5))); @@ -1838,7 +2375,8 @@ mod tests { #[test] fn solve_exact_singular_duplicate_rows() { - let a = Matrix::<3>::from_rows([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [1.0, 2.0, 3.0]]); + let a = Matrix::<3>::try_from_rows([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [1.0, 2.0, 3.0]]) + .unwrap(); let b = Vector::<3>::new([1.0, 2.0, 3.0]); assert_matches!(a.solve_exact(b), Err(LaError::Singular { .. })); } @@ -1846,13 +2384,14 @@ mod tests { #[test] fn solve_exact_5x5_permutation() { // Permutation matrix (swap rows 0↔1): P x = b → x = P^T b. - let a = Matrix::<5>::from_rows([ + let a = Matrix::<5>::try_from_rows([ [0.0, 1.0, 0.0, 0.0, 0.0], [1.0, 0.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, 0.0], [0.0, 0.0, 0.0, 0.0, 1.0], - ]); + ]) + .unwrap(); let b = Vector::<5>::new([10.0, 20.0, 30.0, 40.0, 50.0]); let x = a.solve_exact(b).unwrap(); assert_eq!(x[0], f64_to_bigrational(20.0)); @@ -1879,7 +2418,7 @@ mod tests { for i in 0..$d { rows[i][i] = big; } - let a = Matrix::<$d>::from_rows(rows); + let a = Matrix::<$d>::try_from_rows(rows).unwrap(); // RHS = [big, …, big, 0] → x = [1, …, 1, 0]. let mut b_arr = [big; $d]; b_arr[$d - 1] = 0.0; @@ -1919,7 +2458,7 @@ mod tests { rows[i][i] = val; b_arr[i] = val; } - let a = Matrix::<$d>::from_rows(rows); + let a = Matrix::<$d>::try_from_rows(rows).unwrap(); let b = Vector::<$d>::new(b_arr); let x = a.solve_exact(b).unwrap(); for i in 0..$d { @@ -1998,7 +2537,7 @@ mod tests { for i in 2..$d { rows[i][i] = 1.0; } - let a = Matrix::<$d>::from_rows(rows); + let a = Matrix::<$d>::try_from_rows(rows).unwrap(); // b = [3, 4, 12, 13, …]; padded entries are arbitrary // finite integers so the identity block gives x[i] = b[i]. let mut b_arr = [0.0f64; $d]; @@ -2050,7 +2589,7 @@ mod tests { for i in 3..$d { rows[i][i] = 1.0; } - let a = Matrix::<$d>::from_rows(rows); + let a = Matrix::<$d>::try_from_rows(rows).unwrap(); let mut b_arr = [0.0f64; $d]; b_arr[0] = 6.0; b_arr[1] = 7.0; @@ -2094,7 +2633,7 @@ mod tests { rows[$d - 1][i] = 1.0; } // Last column is left all-zero → rank exactly D-1. - let a = Matrix::<$d>::from_rows(rows); + let a = Matrix::<$d>::try_from_rows(rows).unwrap(); let b = Vector::<$d>::new([1.0; $d]); assert_eq!( a.solve_exact(b), @@ -2125,7 +2664,7 @@ mod tests { for i in 0..$d { rows[i][i] = f64::from($d) + 1.0; } - let a = Matrix::<$d>::from_rows(rows); + let a = Matrix::<$d>::try_from_rows(rows).unwrap(); let b = Vector::<$d>::zero(); let x = a.solve_exact(b).unwrap(); for xi in &x { @@ -2173,11 +2712,12 @@ mod tests { #[test] fn solve_exact_near_singular_3x3_integer_x0() { let perturbation = f64::from_bits(0x3CD0_0000_0000_0000); // 2^-50 - let a = Matrix::<3>::from_rows([ + let a = Matrix::<3>::try_from_rows([ [1.0 + perturbation, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0], - ]); + ]) + .unwrap(); let b = Vector::<3>::new([6.0 + perturbation, 15.0, 24.0]); let x = a.solve_exact(b).unwrap(); let one = BigRational::from_integer(BigInt::from(1)); @@ -2186,6 +2726,30 @@ mod tests { assert_eq!(x[2], one); } + #[test] + fn solve_exact_f64_near_singular_benchmark_rhs_is_rejection_path() { + let perturbation = f64::from_bits(0x3CD0_0000_0000_0000); // 2^-50 + let a = Matrix::<3>::try_from_rows([ + [1.0 + perturbation, 2.0, 3.0], + [4.0, 5.0, 6.0], + [7.0, 8.0, 9.0], + ]) + .unwrap(); + let b = Vector::<3>::new([1.0, 2.0, 3.0]); + + assert_eq!( + a.solve_exact_f64(b), + Err(LaError::Unrepresentable { + index: Some(2), + reason: UnrepresentableReason::RequiresRounding, + }) + ); + assert_eq!( + a.solve_exact_rounded_f64(b).unwrap().into_array()[2].to_bits(), + (1.0f64 / 3.0).to_bits() + ); + } + /// Large-entry 3×3 solve (matches the `exact_large_entries_3x3` /// bench). `A = big · I + (1 - I)` with `big = f64::MAX / 2` and /// `b = [big, 1, 1] = A · [1, 0, 0]`. The `BigInt` augmented system @@ -2196,7 +2760,8 @@ mod tests { fn solve_exact_large_entries_3x3_unit_vector() { let big = f64::MAX / 2.0; assert!(big.is_finite()); - let a = Matrix::<3>::from_rows([[big, 1.0, 1.0], [1.0, big, 1.0], [1.0, 1.0, big]]); + let a = Matrix::<3>::try_from_rows([[big, 1.0, 1.0], [1.0, big, 1.0], [1.0, 1.0, big]]) + .unwrap(); let b = Vector::<3>::new([big, 1.0, 1.0]); let x = a.solve_exact(b).unwrap(); let zero = BigRational::from_integer(BigInt::from(0)); @@ -2210,20 +2775,27 @@ mod tests { /// overflows `f64`. `det_direct()` therefore reports a computed /// [`LaError::NonFinite`], the fast filter inside `det_sign_exact` /// treats that as inconclusive, and the Bareiss fallback resolves the - /// positive sign correctly. `det_exact_f64` must report `Overflow`. + /// positive sign correctly. `det_exact_f64` must report `Unrepresentable`. #[test] fn det_sign_exact_large_entries_3x3_positive() { let big = f64::MAX / 2.0; - let a = Matrix::<3>::from_rows([[big, 1.0, 1.0], [1.0, big, 1.0], [1.0, 1.0, big]]); + let a = Matrix::<3>::try_from_rows([[big, 1.0, 1.0], [1.0, big, 1.0], [1.0, 1.0, big]]) + .unwrap(); // Fast filter is inconclusive (big^3 overflows f64 to +∞), so // this exercises the Bareiss cold path. assert_matches!(a.det_direct(), Err(LaError::NonFinite { row: None, .. })); assert_eq!(a.det_sign_exact().unwrap(), 1); // Cross-validate: the exact `BigRational` determinant must agree - // on sign with `det_sign_exact`, and `det_exact_f64` must overflow - // (the value is representable in BigRational but far exceeds f64). + // on sign with `det_sign_exact`, and `det_exact_f64` must reject the + // conversion (the value is representable in BigRational but far exceeds f64). assert!(a.det_exact().unwrap().is_positive()); - assert_eq!(a.det_exact_f64(), Err(LaError::Overflow { index: None })); + assert_eq!( + a.det_exact_f64(), + Err(LaError::Unrepresentable { + index: None, + reason: UnrepresentableReason::NotFinite, + }) + ); } /// Hilbert matrices are symmetric positive-definite, so @@ -2245,7 +2817,7 @@ mod tests { rows[r][c] = 1.0 / ((r + c + 1) as f64); } } - let h = Matrix::<$d>::from_rows(rows); + let h = Matrix::<$d>::try_from_rows(rows).unwrap(); assert_eq!(h.det_sign_exact().unwrap(), 1); } } @@ -2257,6 +2829,37 @@ mod tests { gen_det_sign_exact_hilbert_positive_tests!(4); gen_det_sign_exact_hilbert_positive_tests!(5); + macro_rules! gen_solve_exact_f64_hilbert_benchmark_rhs_tests { + ($d:literal) => { + paste! { + #[test] + #[allow(clippy::cast_precision_loss)] + fn []() { + let mut rows = [[0.0f64; $d]; $d]; + for r in 0..$d { + for c in 0..$d { + rows[r][c] = 1.0 / ((r + c + 1) as f64); + } + } + let h = Matrix::<$d>::try_from_rows(rows).unwrap(); + let b = Vector::<$d>::new([1.0; $d]); + + assert_matches!( + h.solve_exact_f64(b), + Err(LaError::Unrepresentable { + reason: UnrepresentableReason::RequiresRounding, + .. + }) + ); + assert_matches!(h.solve_exact_rounded_f64(b), Ok(_)); + } + } + }; + } + + gen_solve_exact_f64_hilbert_benchmark_rhs_tests!(4); + gen_solve_exact_f64_hilbert_benchmark_rhs_tests!(5); + /// `solve_exact` on a Hilbert matrix must produce a solution whose /// residual `A · x - b` is *exactly* zero in `BigRational` arithmetic. /// Hilbert entries (`1/3`, `1/5`, `1/6`, `1/7`, …) are non-terminating @@ -2275,7 +2878,7 @@ mod tests { rows[r][c] = 1.0 / ((r + c + 1) as f64); } } - let h = Matrix::<$d>::from_rows(rows); + let h = Matrix::<$d>::try_from_rows(rows).unwrap(); // Use a non-trivial RHS with both positive and negative // entries to avoid accidental structural cancellation. let mut b_arr = [0.0f64; $d]; @@ -2305,7 +2908,7 @@ mod tests { #[test] fn solve_exact_f64_known_2x2() { - let a = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); + let a = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]]).unwrap(); let b = Vector::<2>::new([5.0, 11.0]); let x = a.solve_exact_f64(b).unwrap().into_array(); assert!((x[0] - 1.0).abs() <= f64::EPSILON); @@ -2317,11 +2920,69 @@ mod tests { // [[1/big, 0], [0, 1/big]] x = [big, big] → x = [big², big²], // which overflows f64. let big = f64::MAX / 2.0; - let a = Matrix::<2>::from_rows([[1.0 / big, 0.0], [0.0, 1.0 / big]]); + let a = Matrix::<2>::try_from_rows([[1.0 / big, 0.0], [0.0, 1.0 / big]]).unwrap(); let b = Vector::<2>::new([big, big]); assert_eq!( a.solve_exact_f64(b), - Err(LaError::Overflow { index: Some(0) }) + Err(LaError::Unrepresentable { + index: Some(0), + reason: UnrepresentableReason::NotFinite, + }) + ); + } + + #[test] + fn solve_exact_rounded_f64_overflow_returns_err() { + let big = f64::MAX / 2.0; + let a = Matrix::<2>::try_from_rows([[1.0 / big, 0.0], [0.0, 1.0 / big]]).unwrap(); + let b = Vector::<2>::new([big, big]); + + assert_eq!( + a.solve_exact_rounded_f64(b), + Err(LaError::Unrepresentable { + index: Some(0), + reason: UnrepresentableReason::NotFinite, + }) + ); + } + + #[test] + fn solve_exact_f64_underflow_returns_err_for_nonzero_exact_component() { + let tiny = f64::from_bits(1); + let a = Matrix::<1>::try_from_rows([[2.0]]).unwrap(); + let b = Vector::<1>::new([tiny]); + + assert_eq!( + a.solve_exact_f64(b), + Err(LaError::Unrepresentable { + index: Some(0), + reason: UnrepresentableReason::RequiresRounding, + }) + ); + } + + #[test] + fn solve_exact_f64_rejects_non_dyadic_component() { + let a = Matrix::<1>::try_from_rows([[3.0]]).unwrap(); + let b = Vector::<1>::new([1.0]); + + assert_eq!( + a.solve_exact_f64(b), + Err(LaError::Unrepresentable { + index: Some(0), + reason: UnrepresentableReason::RequiresRounding, + }) + ); + } + + #[test] + fn solve_exact_rounded_f64_rounds_non_dyadic_component() { + let a = Matrix::<1>::try_from_rows([[3.0]]).unwrap(); + let b = Vector::<1>::new([1.0]); + + assert_eq!( + a.solve_exact_rounded_f64(b).unwrap().into_array()[0].to_bits(), + (1.0f64 / 3.0).to_bits() ); } @@ -2338,7 +2999,7 @@ mod tests { #[test] fn gauss_solve_d1() { - let a = Matrix::<1>::from_rows([[2.0]]); + let a = Matrix::<1>::try_from_rows([[2.0]]).unwrap(); let b = Vector::<1>::new([6.0]); let x = a.solve_exact(b).unwrap(); assert_eq!(x[0], f64_to_bigrational(3.0)); @@ -2346,7 +3007,8 @@ mod tests { #[test] fn gauss_solve_singular_column_all_zero() { - let a = Matrix::<3>::from_rows([[1.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]); + let a = Matrix::<3>::try_from_rows([[1.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) + .unwrap(); let b = Vector::<3>::new([1.0, 2.0, 3.0]); assert_eq!(a.solve_exact(b), Err(LaError::Singular { pivot_col: 1 })); } diff --git a/src/ldlt.rs b/src/ldlt.rs index 407d584..5f6339c 100644 --- a/src/ldlt.rs +++ b/src/ldlt.rs @@ -13,8 +13,8 @@ use core::hint::cold_path; -use crate::matrix::{FiniteMatrix, Matrix, SymmetricMatrix}; -use crate::vector::{FiniteVector, Vector}; +use crate::matrix::{Matrix, SymmetricMatrix}; +use crate::vector::Vector; use crate::{LaError, Tolerance}; /// LDLT factorization (`A = L D Lᵀ`) for symmetric positive (semi)definite matrices. @@ -128,7 +128,7 @@ impl Ldlt { } } - let f = FiniteMatrix::new(f)?.into_matrix(); + let f = f.validate_finite()?; Ok(Self { factors: LdltFactors::new_unchecked(f), @@ -200,10 +200,7 @@ impl Ldlt { /// overflows to NaN or infinity. #[inline] pub const fn solve(&self, b: Vector) -> Result, LaError> { - match self.solve_finite(FiniteVector::new_unchecked(b)) { - Ok(x) => Ok(x.into_vector()), - Err(err) => Err(err), - } + self.solve_finite(b) } /// Solve `A x = b` using this LDLT factorization and a finite right-hand side. @@ -215,10 +212,7 @@ impl Ldlt { /// Returns [`LaError::NonFinite`] if a computed substitution intermediate /// overflows to NaN or infinity. #[inline] - pub(crate) const fn solve_finite( - &self, - b: FiniteVector, - ) -> Result, LaError> { + pub(crate) const fn solve_finite(&self, b: Vector) -> Result, LaError> { let mut x = b.into_array(); // Forward substitution: L y = b (L has unit diagonal). @@ -271,7 +265,7 @@ impl Ldlt { ii += 1; } - Ok(FiniteVector::new_unchecked(Vector::new_unchecked(x))) + Ok(Vector::new_unchecked(x)) } } @@ -280,8 +274,6 @@ mod tests { use super::*; use crate::DEFAULT_SINGULAR_TOL; - use crate::matrix::FiniteMatrix; - use core::hint::black_box; use approx::assert_abs_diff_eq; @@ -340,7 +332,7 @@ mod tests { rows[i][i] = diag[i]; } - let a = Matrix::<$d>::from_rows(black_box(rows)); + let a = Matrix::<$d>::try_from_rows(black_box(rows)).unwrap(); let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap(); let expected_det = { @@ -379,11 +371,8 @@ mod tests { #[test] fn solve_2x2_known_spd() { - let a = Matrix::<2>::from_rows(black_box([[4.0, 2.0], [2.0, 3.0]])); - let ldlt = FiniteMatrix::new(a) - .unwrap() - .ldlt(DEFAULT_SINGULAR_TOL) - .unwrap(); + let a = Matrix::<2>::try_from_rows(black_box([[4.0, 2.0], [2.0, 3.0]])).unwrap(); + let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<2>::new(black_box([1.0, 2.0])); let x = ldlt.solve(b).unwrap().into_array(); @@ -395,11 +384,12 @@ mod tests { #[test] fn solve_3x3_spd_tridiagonal_smoke() { - let a = Matrix::<3>::from_rows(black_box([ + let a = Matrix::<3>::try_from_rows(black_box([ [2.0, -1.0, 0.0], [-1.0, 2.0, -1.0], [0.0, -1.0, 2.0], - ])); + ])) + .unwrap(); let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap(); // Choose x = 1 so b = A x is simple: [1, 0, 1]. @@ -414,14 +404,14 @@ mod tests { #[test] fn singular_detected_for_degenerate_psd() { // Rank-1 Gram-like matrix. - let a = Matrix::<2>::from_rows(black_box([[1.0, 1.0], [1.0, 1.0]])); + let a = Matrix::<2>::try_from_rows(black_box([[1.0, 1.0], [1.0, 1.0]])).unwrap(); let err = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap_err(); assert_eq!(err, LaError::Singular { pivot_col: 1 }); } #[test] fn negative_initial_diagonal_reports_not_positive_semidefinite() { - let a = Matrix::<2>::from_rows(black_box([[-1.0, 0.0], [0.0, 1.0]])); + let a = Matrix::<2>::try_from_rows(black_box([[-1.0, 0.0], [0.0, 1.0]])).unwrap(); let err = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap_err(); assert_eq!( err, @@ -434,7 +424,7 @@ mod tests { #[test] fn negative_updated_diagonal_reports_not_positive_semidefinite() { - let a = Matrix::<2>::from_rows(black_box([[1.0, 2.0], [2.0, 1.0]])); + let a = Matrix::<2>::try_from_rows(black_box([[1.0, 2.0], [2.0, 1.0]])).unwrap(); let err = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap_err(); assert_eq!( err, @@ -472,7 +462,7 @@ mod tests { #[test] fn nonfinite_l_multiplier_overflow() { // d = 1e-11 > tol, but l = 1e300 / 1e-11 = 1e311 overflows f64. - let a = Matrix::<2>::from_rows([[1e-11, 1e300], [1e300, 1.0]]); + let a = Matrix::<2>::try_from_rows([[1e-11, 1e300], [1e300, 1.0]]).unwrap(); let err = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap_err(); assert_eq!( err, @@ -487,7 +477,7 @@ mod tests { fn nonfinite_trailing_submatrix_overflow() { // L multiplier is finite (1e200), but the rank-1 update // (-1e200 * 1.0) * 1e200 + 1.0 overflows. - let a = Matrix::<2>::from_rows([[1.0, 1e200], [1e200, 1.0]]); + let a = Matrix::<2>::try_from_rows([[1.0, 1e200], [1e200, 1.0]]).unwrap(); let err = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap_err(); assert_eq!( err, @@ -502,11 +492,12 @@ mod tests { fn nonfinite_solve_forward_substitution_overflow() { // SPD matrix with large L multiplier: L[1,0] = 1e153. // Forward substitution overflows: y[1] = 0 - 1e153 * 1e156 = -inf. - let a = Matrix::<3>::from_rows([ + let a = Matrix::<3>::try_from_rows([ [1.0, 1e153, 0.0], [1e153, 1e306 + 1.0, 0.0], [0.0, 0.0, 1.0], - ]); + ]) + .unwrap(); let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<3>::new([1e156, 0.0, 0.0]); @@ -520,7 +511,8 @@ mod tests { // D=[1,1,1], L[2,1]=2. Forward sub and diagonal solve produce // z=[0,0,1e308]. Back-substitution: x[2]=1e308 then // x[1] = 0 - 2*1e308 = -inf (overflows f64). - let a = Matrix::<3>::from_rows([[1.0, 0.0, 0.0], [0.0, 1.0, 2.0], [0.0, 2.0, 5.0]]); + let a = Matrix::<3>::try_from_rows([[1.0, 0.0, 0.0], [0.0, 1.0, 2.0], [0.0, 2.0, 5.0]]) + .unwrap(); let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<3>::new([0.0, 0.0, 1e308]); @@ -535,7 +527,7 @@ mod tests { // large RHS unchanged, then the diagonal solve z[1] = y[1] / D[1] // = 1e300 / 1e-11 = 1e311 overflows f64, exercising the // `!v.is_finite()` branch of the diagonal solve. - let a = Matrix::<2>::from_rows([[1.0, 0.0], [0.0, 1.0e-11]]); + let a = Matrix::<2>::try_from_rows([[1.0, 0.0], [0.0, 1.0e-11]]).unwrap(); let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<2>::new([0.0, 1.0e300]); @@ -545,13 +537,14 @@ mod tests { #[test] fn det_rejects_product_overflow() { - let a = Matrix::<5>::from_rows([ + let a = Matrix::<5>::try_from_rows([ [1.0e100, 0.0, 0.0, 0.0, 0.0], [0.0, 1.0e100, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0e100, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0e100, 0.0], [0.0, 0.0, 0.0, 0.0, 1.0e100], - ]); + ]) + .unwrap(); let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap(); assert_eq!(ldlt.det(), Err(LaError::NonFinite { row: None, col: 3 })); } @@ -559,7 +552,8 @@ mod tests { #[test] fn asymmetric_input_returns_typed_error() { // a[0][1] = 2.0 but a[1][0] = -2.0 → clearly asymmetric. - let a = Matrix::<3>::from_rows([[4.0, 2.0, 0.0], [-2.0, 5.0, 1.0], [0.0, 1.0, 3.0]]); + let a = Matrix::<3>::try_from_rows([[4.0, 2.0, 0.0], [-2.0, 5.0, 1.0], [0.0, 1.0, 3.0]]) + .unwrap(); assert_eq!( a.ldlt(DEFAULT_SINGULAR_TOL), Err(LaError::Asymmetric { diff --git a/src/lib.rs b/src/lib.rs index 63109eb..20576ab 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -101,6 +101,19 @@ mod readme_doctests { /// let det_f64 = m.det_exact_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 rounded_det = match inexact.det_exact_f64() { + /// Ok(det) => det, + /// Err(err) if err.requires_rounding() => inexact.det_exact_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([ @@ -109,7 +122,12 @@ mod readme_doctests { /// [0.0, big, 1.0], /// ])?; /// let huge_det = huge.det_exact()?; - /// assert_eq!(huge.det_exact_f64(), Err(LaError::Overflow { index: None })); + /// assert_eq!( + /// huge.det_exact_f64() + /// .err() + /// .and_then(|err| err.unrepresentable_reason()), + /// Some(UnrepresentableReason::NotFinite) + /// ); /// println!("exact determinant = {huge_det}"); /// /// // Exact linear system solve @@ -124,8 +142,15 @@ mod readme_doctests { fn exact_arithmetic_example() {} } +mod error; #[cfg(feature = "exact")] mod exact; +mod ldlt; +mod lu; +mod matrix; +mod tolerance; +mod vector; + #[cfg(feature = "exact")] pub use num_bigint::BigInt; #[cfg(feature = "exact")] @@ -133,13 +158,6 @@ pub use num_rational::BigRational; #[cfg(feature = "exact")] pub use num_traits::{FromPrimitive, Signed, ToPrimitive}; -mod ldlt; -mod lu; -mod matrix; -mod vector; - -use core::fmt; - // --------------------------------------------------------------------------- // Error-bound constants for `Matrix::det_errbound()`. // @@ -274,443 +292,12 @@ pub const ERR_COEFF_4: f64 = 12.0 * EPS + 128.0 * EPS * EPS; /// dispatch surface explicit. pub const MAX_STACK_MATRIX_DISPATCH_DIM: usize = 7; -/// Finite, non-negative tolerance used by numerical predicates and factorizations. -/// -/// Construct with [`Tolerance::new`] when accepting raw caller input. Once -/// constructed, the stored value is guaranteed to be finite and `>= 0`, so -/// downstream algorithms do not need to revalidate the tolerance. -/// -/// This is the crate-wide tolerance contract: raw negative, NaN, and infinite -/// values are rejected with [`LaError::InvalidTolerance`] at construction time. -#[must_use] -#[derive(Clone, Copy, Debug, PartialEq)] -pub struct Tolerance { - value: f64, -} - -impl Tolerance { - /// Construct a tolerance without checking the raw value. - /// - /// This crate-internal escape hatch is only for constants whose finite, - /// non-negative value is visible at the call site. Public callers should - /// use [`Tolerance::new`] so the returned value carries the validation - /// proof. - pub(crate) const fn new_unchecked(value: f64) -> Self { - Self { value } - } - - /// Construct a finite, non-negative tolerance. - /// - /// # Examples - /// ``` - /// use la_stack::prelude::*; - /// - /// # fn main() -> Result<(), LaError> { - /// let tol = Tolerance::new(1e-12)?; - /// assert_eq!(tol.get(), 1e-12); - /// # Ok(()) - /// # } - /// ``` - /// - /// # Errors - /// Returns [`LaError::InvalidTolerance`] when `value` is NaN, infinite, or - /// negative. - #[inline] - pub const fn new(value: f64) -> Result { - if value >= 0.0 && value.is_finite() { - Ok(Self::new_unchecked(value)) - } else { - Err(LaError::invalid_tolerance(value)) - } - } - - /// Return the raw finite, non-negative tolerance value. - /// - /// # Examples - /// ``` - /// use la_stack::prelude::*; - /// - /// # fn main() -> Result<(), LaError> { - /// let tol = Tolerance::new(0.0)?; - /// assert_eq!(tol.get(), 0.0); - /// # Ok(()) - /// # } - /// ``` - #[inline] - #[must_use] - pub const fn get(self) -> f64 { - self.value - } -} - -/// Default absolute threshold used for singularity/degeneracy detection. -/// -/// This is intentionally conservative for geometric predicates and small systems. -/// -/// Conceptually, this is an absolute bound for deciding when a scalar should be treated -/// as "numerically zero" (e.g. LU pivots, LDLT diagonal entries). -pub const DEFAULT_SINGULAR_TOL: Tolerance = Tolerance::new_unchecked(1e-12); - -/// Relative tolerance used to validate matrices for LDLT factorization. -pub(crate) const LDLT_SYMMETRY_REL_TOL: Tolerance = Tolerance::new_unchecked(1e-12); - -/// Linear algebra errors. -/// -/// This enum is `#[non_exhaustive]` — downstream `match` arms must include a -/// wildcard (`_`) pattern to compile, allowing new variants to be added in -/// future minor releases without breaking existing code. -#[derive(Clone, Copy, Debug, PartialEq)] -#[non_exhaustive] -pub enum LaError { - /// The matrix is (numerically) singular. - Singular { - /// The factorization column/step where a suitable pivot/diagonal could not be found. - pivot_col: usize, - }, - /// A non-finite value (NaN/∞) was encountered. - /// - /// The `(row, col)` coordinate follows a consistent convention across the crate: - /// - /// - `row: Some(r), col: c` — the non-finite value is tied to a matrix/factor - /// cell at `(r, c)`, either because a stored input/factor cell is already - /// non-finite or because factorization computed a non-finite value for - /// that cell before storing it. - /// - `row: None, col: c` — the non-finite value is tied to a vector entry, - /// determinant product, tolerance-scale accumulator, solve accumulator, or - /// other scalar/intermediate that has no matrix row coordinate. - NonFinite { - /// Row of the non-finite entry for a stored matrix cell, or `None` for - /// a vector-input entry or a computed intermediate. See the variant - /// docs for the full convention. - row: Option, - /// Column index (stored cell), vector index, or factorization/solve - /// step where the non-finite value was detected. - col: usize, - }, - /// The exact result overflows the target representation (e.g. `f64`). - /// - /// Returned by `Matrix::det_exact_f64` and `Matrix::solve_exact_f64` - /// (requires `exact` feature) when an exact value is too large to - /// represent as a finite `f64`. - Overflow { - /// For vector results (e.g. `solve_exact_f64`), the index of the - /// component that overflowed. `None` for scalar results. - index: Option, - }, - /// Exact determinant scaling overflowed the internal exponent representation. - DeterminantScaleOverflow { - /// Matrix dimension `D`. - dim: usize, - /// Minimum decomposed f64 exponent among non-zero matrix entries. - min_exponent: i32, - }, - /// A requested runtime matrix dimension has no stack-dispatch arm. - UnsupportedDimension { - /// Runtime dimension requested by the caller. - requested: usize, - /// Largest runtime dimension supported by the dispatch helper. - max: usize, - }, - /// A matrix index is outside the `D×D` storage domain. - IndexOutOfBounds { - /// Requested row index. - row: usize, - /// Requested column index. - col: usize, - /// Matrix dimension `D`; valid row and column indices are `< dim`. - dim: usize, - }, - /// A tolerance value is not finite and non-negative. - InvalidTolerance { - /// Raw tolerance supplied by the caller. - value: f64, - }, - /// A matrix required to be symmetric has an asymmetric off-diagonal pair. - Asymmetric { - /// Row index of the first asymmetric pair. - row: usize, - /// Column index of the first asymmetric pair. - col: usize, - /// Matrix dimension `D`. - dim: usize, - }, - /// A symmetric matrix failed the positive-semidefinite LDLT domain check. - NotPositiveSemidefinite { - /// Factorization column/step where a negative LDLT diagonal was found. - pivot_col: usize, - /// Negative diagonal value observed at that step. - value: f64, - }, -} - -impl LaError { - /// Construct a [`LaError::NonFinite`] pinpointing a stored matrix cell at `(row, col)`. - /// - /// Use this for non-finite values read from a stored `Matrix` entry or - /// factorization cell, and for non-finite factorization updates that would - /// be stored at `(row, col)` if accepted. The resulting error has - /// `row: Some(row), col`, matching the matrix/factor-cell convention - /// documented on [`NonFinite`](Self::NonFinite). For vector-input entries - /// or scalar intermediates without a matrix row coordinate, use - /// [`non_finite_at`](Self::non_finite_at). - /// - /// # Examples - /// ``` - /// use la_stack::prelude::*; - /// - /// assert_eq!( - /// LaError::non_finite_cell(1, 2), - /// LaError::NonFinite { - /// row: Some(1), - /// col: 2, - /// } - /// ); - /// ``` - #[inline] - #[must_use] - pub const fn non_finite_cell(row: usize, col: usize) -> Self { - Self::NonFinite { - row: Some(row), - col, - } - } - - /// Construct a [`LaError::NonFinite`] pinpointing a vector-input entry or - /// computed scalar/intermediate at index `col`. - /// - /// Use this for non-finite values in a `Vector` input, determinant scalar, - /// tolerance-scale accumulator, or solve accumulator that overflowed during - /// forward/back substitution. The resulting error has `row: None, col`, - /// matching the vector/scalar-intermediate convention documented on - /// [`NonFinite`](Self::NonFinite). For stored matrix cells or computed - /// factorization updates tied to a matrix cell, use - /// [`non_finite_cell`](Self::non_finite_cell). - /// - /// # Examples - /// ``` - /// use la_stack::prelude::*; - /// - /// assert_eq!( - /// LaError::non_finite_at(2), - /// LaError::NonFinite { row: None, col: 2 } - /// ); - /// ``` - #[inline] - #[must_use] - pub const fn non_finite_at(col: usize) -> Self { - Self::NonFinite { row: None, col } - } - - /// Construct a [`LaError::DeterminantScaleOverflow`] for exact determinant scaling. - /// - /// # Examples - /// ``` - /// use la_stack::prelude::*; - /// - /// assert_eq!( - /// LaError::determinant_scale_overflow(3, -1074), - /// LaError::DeterminantScaleOverflow { - /// dim: 3, - /// min_exponent: -1074, - /// } - /// ); - /// ``` - #[inline] - #[must_use] - pub const fn determinant_scale_overflow(dim: usize, min_exponent: i32) -> Self { - Self::DeterminantScaleOverflow { dim, min_exponent } - } - - /// Construct a [`LaError::UnsupportedDimension`] for runtime stack dispatch. - /// - /// # Examples - /// ``` - /// use la_stack::prelude::*; - /// - /// assert_eq!( - /// LaError::unsupported_dimension(8, MAX_STACK_MATRIX_DISPATCH_DIM), - /// LaError::UnsupportedDimension { - /// requested: 8, - /// max: MAX_STACK_MATRIX_DISPATCH_DIM, - /// } - /// ); - /// ``` - #[inline] - #[must_use] - pub const fn unsupported_dimension(requested: usize, max: usize) -> Self { - Self::UnsupportedDimension { requested, max } - } - - /// Construct a [`LaError::IndexOutOfBounds`] for a `D×D` matrix index. - /// - /// # Examples - /// ``` - /// use la_stack::prelude::*; - /// - /// assert_eq!( - /// LaError::index_out_of_bounds(2, 0, 2), - /// LaError::IndexOutOfBounds { - /// row: 2, - /// col: 0, - /// dim: 2, - /// } - /// ); - /// ``` - #[inline] - #[must_use] - pub const fn index_out_of_bounds(row: usize, col: usize, dim: usize) -> Self { - Self::IndexOutOfBounds { row, col, dim } - } - - /// Construct a [`LaError::InvalidTolerance`] for a raw tolerance value. - /// - /// # Examples - /// ``` - /// use la_stack::prelude::*; - /// - /// assert_eq!( - /// LaError::invalid_tolerance(-1.0), - /// LaError::InvalidTolerance { value: -1.0 } - /// ); - /// ``` - #[inline] - #[must_use] - pub const fn invalid_tolerance(value: f64) -> Self { - Self::InvalidTolerance { value } - } - - /// Construct a [`LaError::Asymmetric`] for a `D×D` matrix. - /// - /// # Examples - /// ``` - /// use la_stack::prelude::*; - /// - /// assert_eq!( - /// LaError::asymmetric(0, 1, 3), - /// LaError::Asymmetric { - /// row: 0, - /// col: 1, - /// dim: 3, - /// } - /// ); - /// ``` - #[inline] - #[must_use] - pub const fn asymmetric(row: usize, col: usize, dim: usize) -> Self { - Self::Asymmetric { row, col, dim } - } - - /// Construct a [`LaError::NotPositiveSemidefinite`] for LDLT factorization. - /// - /// # Examples - /// ``` - /// use la_stack::prelude::*; - /// - /// assert_eq!( - /// LaError::not_positive_semidefinite(1, -3.0), - /// LaError::NotPositiveSemidefinite { - /// pivot_col: 1, - /// value: -3.0, - /// } - /// ); - /// ``` - #[inline] - #[must_use] - pub const fn not_positive_semidefinite(pivot_col: usize, value: f64) -> Self { - Self::NotPositiveSemidefinite { pivot_col, value } - } - - /// Parse a raw tolerance into a finite, non-negative [`Tolerance`]. - /// - /// # Examples - /// ``` - /// use la_stack::prelude::*; - /// - /// assert_eq!(LaError::validate_tolerance(1e-12)?.get(), 1e-12); - /// - /// let raw = 0.0; - /// let tol = LaError::validate_tolerance(raw)?; - /// let _lu = Matrix::<2>::identity().lu(tol)?; - /// - /// assert_eq!( - /// LaError::validate_tolerance(-1.0), - /// Err(LaError::InvalidTolerance { value: -1.0 }) - /// ); - /// # Ok::<(), LaError>(()) - /// ``` - /// - /// # Errors - /// Returns [`LaError::InvalidTolerance`] when `value` is NaN, infinite, or - /// negative. - #[inline] - pub const fn validate_tolerance(value: f64) -> Result { - Tolerance::new(value) - } -} - -impl fmt::Display for LaError { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - match *self { - Self::Singular { pivot_col } => { - write!(f, "singular matrix at pivot column {pivot_col}") - } - Self::NonFinite { row: Some(r), col } => { - write!(f, "non-finite value at ({r}, {col})") - } - Self::NonFinite { row: None, col } => { - write!(f, "non-finite value at index {col}") - } - Self::Overflow { index: Some(i) } => { - write!( - f, - "exact result overflows the target representation at index {i}" - ) - } - Self::Overflow { index: None } => { - write!(f, "exact result overflows the target representation") - } - Self::DeterminantScaleOverflow { dim, min_exponent } => { - write!( - f, - "exact determinant scale exponent overflows for dimension {dim} with minimum entry exponent {min_exponent}" - ) - } - Self::UnsupportedDimension { requested, max } => { - write!( - f, - "unsupported matrix dimension {requested}; maximum stack-dispatch dimension is {max}" - ) - } - Self::IndexOutOfBounds { row, col, dim } => { - write!( - f, - "matrix index ({row}, {col}) is out of bounds for dimension {dim}" - ) - } - Self::InvalidTolerance { value } => { - write!(f, "invalid tolerance {value}; expected finite value >= 0") - } - Self::Asymmetric { row, col, dim } => { - write!( - f, - "matrix is not symmetric for dimension {dim}: asymmetric pair ({row}, {col})" - ) - } - Self::NotPositiveSemidefinite { pivot_col, value } => { - write!( - f, - "matrix is not positive semidefinite at LDLT pivot column {pivot_col}: diagonal value {value} < 0" - ) - } - } - } -} - -impl std::error::Error for LaError {} - +pub use error::{LaError, UnrepresentableReason}; pub use ldlt::Ldlt; pub use lu::Lu; pub use matrix::Matrix; +pub(crate) use tolerance::LDLT_SYMMETRY_REL_TOL; +pub use tolerance::{DEFAULT_SINGULAR_TOL, Tolerance}; pub use vector::Vector; /// Fallibly dispatch a runtime dimension to a concrete stack-allocated matrix. @@ -798,9 +385,9 @@ macro_rules! try_with_stack_matrix { /// /// This prelude re-exports the primary types and constants: [`Matrix`], /// [`Vector`], [`Lu`], [`Ldlt`], [`Tolerance`], [`LaError`], -/// [`DEFAULT_SINGULAR_TOL`], and the determinant error bound coefficients -/// [`ERR_COEFF_2`], [`ERR_COEFF_3`], and [`ERR_COEFF_4`]. It also re-exports -/// [`MAX_STACK_MATRIX_DISPATCH_DIM`] and +/// [`UnrepresentableReason`], [`DEFAULT_SINGULAR_TOL`], and the determinant +/// error bound coefficients [`ERR_COEFF_2`], [`ERR_COEFF_3`], and +/// [`ERR_COEFF_4`]. It also re-exports [`MAX_STACK_MATRIX_DISPATCH_DIM`] and /// [`try_with_stack_matrix!`] for runtime-to-const matrix dispatch. /// /// When the `exact` feature is enabled, `BigInt` and `BigRational` are also @@ -813,7 +400,8 @@ macro_rules! try_with_stack_matrix { pub mod prelude { pub use crate::{ DEFAULT_SINGULAR_TOL, ERR_COEFF_2, ERR_COEFF_3, ERR_COEFF_4, LaError, Ldlt, Lu, - MAX_STACK_MATRIX_DISPATCH_DIM, Matrix, Tolerance, Vector, try_with_stack_matrix, + MAX_STACK_MATRIX_DISPATCH_DIM, Matrix, Tolerance, UnrepresentableReason, Vector, + try_with_stack_matrix, }; #[cfg(feature = "exact")] @@ -822,215 +410,36 @@ pub mod prelude { #[cfg(test)] mod tests { - use core::assert_matches; - use super::*; use approx::assert_abs_diff_eq; - #[test] - fn default_singular_tol_is_expected() { - assert_abs_diff_eq!(DEFAULT_SINGULAR_TOL.get(), 1e-12, epsilon = 0.0); - } + mod prelude_tests { + use approx::assert_abs_diff_eq; - #[test] - fn tolerance_new_accepts_finite_non_negative_values() { - assert_eq!( - Tolerance::new(0.0).unwrap().get().to_bits(), - 0.0f64.to_bits() - ); - assert_eq!( - Tolerance::new(1e-12).unwrap().get().to_bits(), - 1e-12f64.to_bits() - ); - assert_eq!( - Tolerance::new(f64::MAX).unwrap().get().to_bits(), - f64::MAX.to_bits() - ); - } - - #[test] - fn tolerance_new_rejects_negative_nan_and_infinity() { - assert_eq!( - Tolerance::new(-1.0), - Err(LaError::InvalidTolerance { value: -1.0 }) - ); - assert_matches!( - Tolerance::new(f64::NAN), - Err(LaError::InvalidTolerance { value }) if value.is_nan() - ); - assert_eq!( - Tolerance::new(f64::INFINITY), - Err(LaError::InvalidTolerance { - value: f64::INFINITY, - }) - ); - assert_eq!( - Tolerance::new(f64::NEG_INFINITY), - Err(LaError::InvalidTolerance { - value: f64::NEG_INFINITY, - }) - ); - } - - #[test] - fn validate_tolerance_matches_tolerance_new() { - for value in [0.0, 1e-12, f64::MAX] { - assert_eq!(LaError::validate_tolerance(value), Tolerance::new(value)); - } - - assert_eq!( - LaError::validate_tolerance(-1.0), - Err(LaError::InvalidTolerance { value: -1.0 }) - ); - assert_matches!( - LaError::validate_tolerance(f64::NAN), - Err(LaError::InvalidTolerance { value }) if value.is_nan() - ); - assert_eq!( - LaError::validate_tolerance(f64::INFINITY), - Err(LaError::InvalidTolerance { - value: f64::INFINITY, - }) - ); - assert_eq!( - LaError::validate_tolerance(f64::NEG_INFINITY), - Err(LaError::InvalidTolerance { - value: f64::NEG_INFINITY, - }) - ); - } - - #[test] - fn laerror_display_formats_singular() { - let err = LaError::Singular { pivot_col: 3 }; - assert_eq!(err.to_string(), "singular matrix at pivot column 3"); - } - - #[test] - fn laerror_display_formats_nonfinite_with_row() { - let err = LaError::NonFinite { - row: Some(1), - col: 2, - }; - assert_eq!(err.to_string(), "non-finite value at (1, 2)"); - } - - #[test] - fn laerror_display_formats_nonfinite_without_row() { - let err = LaError::NonFinite { row: None, col: 3 }; - assert_eq!(err.to_string(), "non-finite value at index 3"); - } - - #[test] - fn laerror_display_formats_overflow() { - let err = LaError::Overflow { index: None }; - assert_eq!( - err.to_string(), - "exact result overflows the target representation" - ); - } - - #[test] - fn laerror_display_formats_overflow_with_index() { - let err = LaError::Overflow { index: Some(2) }; - assert_eq!( - err.to_string(), - "exact result overflows the target representation at index 2" - ); - } - - #[test] - fn laerror_display_formats_determinant_scale_overflow() { - let err = LaError::DeterminantScaleOverflow { - dim: 3, - min_exponent: -1074, - }; - assert_eq!( - err.to_string(), - "exact determinant scale exponent overflows for dimension 3 with minimum entry exponent -1074" - ); - } - - #[test] - fn laerror_display_formats_unsupported_dimension() { - let err = LaError::UnsupportedDimension { - requested: 8, - max: MAX_STACK_MATRIX_DISPATCH_DIM, - }; - assert_eq!( - err.to_string(), - "unsupported matrix dimension 8; maximum stack-dispatch dimension is 7" - ); - } - - #[test] - fn laerror_display_formats_index_out_of_bounds() { - let err = LaError::IndexOutOfBounds { - row: 3, - col: 0, - dim: 3, - }; - assert_eq!( - err.to_string(), - "matrix index (3, 0) is out of bounds for dimension 3" - ); - } - - #[test] - fn laerror_display_formats_invalid_tolerance() { - let err = LaError::InvalidTolerance { value: -1.0 }; - assert_eq!( - err.to_string(), - "invalid tolerance -1; expected finite value >= 0" - ); - } - - #[test] - fn laerror_display_formats_asymmetric() { - let err = LaError::Asymmetric { - row: 0, - col: 2, - dim: 3, - }; - assert_eq!( - err.to_string(), - "matrix is not symmetric for dimension 3: asymmetric pair (0, 2)" - ); - } - - #[test] - fn laerror_display_formats_not_positive_semidefinite() { - let err = LaError::NotPositiveSemidefinite { - pivot_col: 1, - value: -3.0, - }; - assert_eq!( - err.to_string(), - "matrix is not positive semidefinite at LDLT pivot column 1: diagonal value -3 < 0" - ); - } - - #[test] - fn laerror_is_std_error_with_no_source() { - let err = LaError::Singular { pivot_col: 0 }; - let e: &dyn std::error::Error = &err; - assert!(e.source().is_none()); - } - - #[test] - fn prelude_reexports_compile_and_work() -> Result<(), LaError> { use crate::prelude::*; - // Use the items so we know they are in scope and usable. - let m = Matrix::<2>::identity(); - let v = Vector::<2>::try_new([1.0, 2.0])?; - assert_abs_diff_eq!(m.inf_norm()?, 1.0, epsilon = 0.0); - assert_abs_diff_eq!(v.norm2_sq()?, 5.0, epsilon = 0.0); - let _ = m.lu(DEFAULT_SINGULAR_TOL)?.solve(v)?; - let _ = m.ldlt(DEFAULT_SINGULAR_TOL)?.solve(v)?; - assert_eq!(MAX_STACK_MATRIX_DISPATCH_DIM, 7); - Ok(()) + #[test] + fn prelude_reexports_compile_and_work() -> Result<(), LaError> { + // Use the items so we know they are in scope and usable. + let m = Matrix::<2>::identity(); + let v = Vector::<2>::try_new([1.0, 2.0])?; + let tol = Tolerance::new(0.0)?; + assert_abs_diff_eq!(tol.get(), 0.0, epsilon = 0.0); + assert_abs_diff_eq!(m.inf_norm()?, 1.0, epsilon = 0.0); + assert_abs_diff_eq!(v.norm2_sq()?, 5.0, epsilon = 0.0); + let _ = m.lu(DEFAULT_SINGULAR_TOL)?.solve(v)?; + let _ = m.ldlt(DEFAULT_SINGULAR_TOL)?.solve(v)?; + assert_eq!( + LaError::unrepresentable(None, UnrepresentableReason::RequiresRounding), + LaError::Unrepresentable { + index: None, + reason: UnrepresentableReason::RequiresRounding, + } + ); + assert_eq!(MAX_STACK_MATRIX_DISPATCH_DIM, 7); + Ok(()) + } } macro_rules! gen_stack_matrix_dispatch_tests { diff --git a/src/lu.rs b/src/lu.rs index 5520ba6..11c1976 100644 --- a/src/lu.rs +++ b/src/lu.rs @@ -2,8 +2,8 @@ use core::hint::cold_path; -use crate::matrix::{FiniteMatrix, Matrix}; -use crate::vector::{FiniteVector, Vector}; +use crate::matrix::Matrix; +use crate::vector::Vector; use crate::{LaError, Tolerance}; /// LU decomposition (PA = LU) with partial pivoting. @@ -48,7 +48,7 @@ impl LuFactors { impl Lu { /// Factor a finite square matrix into in-place LU storage for - /// `FiniteMatrix::lu`. + /// [`Matrix::lu`]. /// /// The input has already proven finite entries, so LU construction rejects /// numerically singular pivots and non-finite elimination intermediates @@ -57,8 +57,8 @@ impl Lu { /// value produced during elimination. #[inline] #[allow(clippy::needless_range_loop)] - pub(crate) fn factor_finite(a: FiniteMatrix, tol: Tolerance) -> Result { - let mut lu = a.into_matrix(); + pub(crate) fn factor_finite(a: Matrix, tol: Tolerance) -> Result { + let mut lu = a; let tol = tol.get(); let mut piv = [0usize; D]; @@ -110,7 +110,7 @@ impl Lu { } } - let lu = FiniteMatrix::new(lu)?.into_matrix(); + let lu = lu.validate_finite()?; Ok(Self { factors: LuFactors::new_unchecked(lu), @@ -148,10 +148,7 @@ impl Lu { /// overflows to NaN or infinity. #[inline] pub const fn solve(&self, b: Vector) -> Result, LaError> { - match self.solve_finite(FiniteVector::new_unchecked(b)) { - Ok(x) => Ok(x.into_vector()), - Err(err) => Err(err), - } + self.solve_finite(b) } /// Solve `A x = b` using this LU factorization and a finite right-hand side. @@ -163,10 +160,7 @@ impl Lu { /// Returns [`LaError::NonFinite`] if a computed substitution intermediate /// overflows to NaN or infinity. #[inline] - pub(crate) const fn solve_finite( - &self, - b: FiniteVector, - ) -> Result, LaError> { + pub(crate) const fn solve_finite(&self, b: Vector) -> Result, LaError> { let mut x = [0.0; D]; let mut i = 0; let b = b.as_array(); @@ -220,7 +214,7 @@ impl Lu { ii += 1; } - Ok(FiniteVector::new_unchecked(Vector::new_unchecked(x))) + Ok(Vector::new_unchecked(x)) } /// Determinant of the original matrix. @@ -284,7 +278,7 @@ mod tests { } rows.swap(0, 1); - let a = Matrix::<$d>::from_rows(black_box(rows)); + let a = Matrix::<$d>::try_from_rows(black_box(rows)).unwrap(); let lu_fn: fn(Matrix<$d>, Tolerance) -> Result, LaError> = black_box(Matrix::<$d>::lu); let lu = lu_fn(a, DEFAULT_SINGULAR_TOL).unwrap(); @@ -324,7 +318,7 @@ mod tests { } rows.swap(0, 1); - let a = Matrix::<$d>::from_rows(black_box(rows)); + let a = Matrix::<$d>::try_from_rows(black_box(rows)).unwrap(); let lu_fn: fn(Matrix<$d>, Tolerance) -> Result, LaError> = black_box(Matrix::<$d>::lu); let lu = lu_fn(a, DEFAULT_SINGULAR_TOL).unwrap(); @@ -364,7 +358,7 @@ mod tests { } } - let a = Matrix::<$d>::from_rows(black_box(rows)); + let a = Matrix::<$d>::try_from_rows(black_box(rows)).unwrap(); let lu_fn: fn(Matrix<$d>, Tolerance) -> Result, LaError> = black_box(Matrix::<$d>::lu); let lu = lu_fn(a, DEFAULT_SINGULAR_TOL).unwrap(); @@ -403,7 +397,7 @@ mod tests { } } - let a = Matrix::<$d>::from_rows(black_box(rows)); + let a = Matrix::<$d>::try_from_rows(black_box(rows)).unwrap(); let lu_fn: fn(Matrix<$d>, Tolerance) -> Result, LaError> = black_box(Matrix::<$d>::lu); let lu = lu_fn(a, DEFAULT_SINGULAR_TOL).unwrap(); @@ -422,11 +416,8 @@ mod tests { #[test] fn solve_1x1() { - let a = Matrix::<1>::from_rows(black_box([[2.0]])); - let lu = FiniteMatrix::new(a) - .unwrap() - .lu(DEFAULT_SINGULAR_TOL) - .unwrap(); + let a = Matrix::<1>::try_from_rows(black_box([[2.0]])).unwrap(); + let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<1>::new(black_box([6.0])); let solve_fn: fn(&Lu<1>, Vector<1>) -> Result, LaError> = @@ -440,11 +431,8 @@ mod tests { #[test] fn solve_2x2_basic() { - let a = Matrix::<2>::from_rows(black_box([[1.0, 2.0], [3.0, 4.0]])); - let lu = FiniteMatrix::new(a) - .unwrap() - .lu(DEFAULT_SINGULAR_TOL) - .unwrap(); + let a = Matrix::<2>::try_from_rows(black_box([[1.0, 2.0], [3.0, 4.0]])).unwrap(); + let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<2>::new(black_box([5.0, 11.0])); let solve_fn: fn(&Lu<2>, Vector<2>) -> Result, LaError> = @@ -457,7 +445,7 @@ mod tests { #[test] fn det_2x2_basic() { - let a = Matrix::<2>::from_rows(black_box([[1.0, 2.0], [3.0, 4.0]])); + let a = Matrix::<2>::try_from_rows(black_box([[1.0, 2.0], [3.0, 4.0]])).unwrap(); let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); let det_fn: fn(&Lu<2>) -> Result = black_box(Lu::<2>::det); @@ -467,7 +455,7 @@ mod tests { #[test] fn det_requires_pivot_sign() { // Row swap ⇒ determinant sign flip. - let a = Matrix::<2>::from_rows(black_box([[0.0, 1.0], [1.0, 0.0]])); + let a = Matrix::<2>::try_from_rows(black_box([[0.0, 1.0], [1.0, 0.0]])).unwrap(); let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); let det_fn: fn(&Lu<2>) -> Result = black_box(Lu::<2>::det); @@ -476,7 +464,7 @@ mod tests { #[test] fn solve_requires_pivoting() { - let a = Matrix::<2>::from_rows(black_box([[0.0, 1.0], [1.0, 0.0]])); + let a = Matrix::<2>::try_from_rows(black_box([[0.0, 1.0], [1.0, 0.0]])).unwrap(); let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<2>::new(black_box([1.0, 2.0])); @@ -490,7 +478,7 @@ mod tests { #[test] fn singular_detected() { - let a = Matrix::<2>::from_rows(black_box([[1.0, 2.0], [2.0, 4.0]])); + let a = Matrix::<2>::try_from_rows(black_box([[1.0, 2.0], [2.0, 4.0]])).unwrap(); let err = a.lu(DEFAULT_SINGULAR_TOL).unwrap_err(); assert_eq!(err, LaError::Singular { pivot_col: 1 }); } @@ -498,7 +486,7 @@ mod tests { #[test] fn singular_due_to_tolerance_at_first_pivot() { // Not exactly singular, but below DEFAULT_SINGULAR_TOL. - let a = Matrix::<2>::from_rows(black_box([[1e-13, 0.0], [0.0, 1.0]])); + let a = Matrix::<2>::try_from_rows(black_box([[1e-13, 0.0], [0.0, 1.0]])).unwrap(); let err = a.lu(DEFAULT_SINGULAR_TOL).unwrap_err(); assert_eq!(err, LaError::Singular { pivot_col: 0 }); } @@ -529,8 +517,12 @@ mod tests { #[test] fn nonfinite_detected_in_trailing_update() { - let a = - Matrix::<3>::from_rows([[1.0, f64::MAX, 0.0], [-1.0, f64::MAX, 0.0], [0.0, 0.0, 1.0]]); + let a = Matrix::<3>::try_from_rows([ + [1.0, f64::MAX, 0.0], + [-1.0, f64::MAX, 0.0], + [0.0, 0.0, 1.0], + ]) + .unwrap(); let err = a.lu(DEFAULT_SINGULAR_TOL).unwrap_err(); assert_eq!( @@ -545,7 +537,8 @@ mod tests { #[test] fn solve_nonfinite_forward_substitution_overflow() { // L has a -1 multiplier, and a large RHS makes forward substitution overflow. - let a = Matrix::<3>::from_rows([[1.0, 0.0, 0.0], [-1.0, 1.0, 0.0], [0.0, 0.0, 1.0]]); + let a = Matrix::<3>::try_from_rows([[1.0, 0.0, 0.0], [-1.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) + .unwrap(); let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<3>::new([1.0e308, 1.0e308, 0.0]); @@ -556,7 +549,7 @@ mod tests { #[test] fn solve_nonfinite_back_substitution_overflow() { // Make x[1] overflow during back substitution, then ensure it is detected on the next row. - let a = Matrix::<2>::from_rows([[1.0, 1.0], [0.0, 2.0e-12]]); + let a = Matrix::<2>::try_from_rows([[1.0, 1.0], [0.0, 2.0e-12]]).unwrap(); let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<2>::new([0.0, 1.0e300]); @@ -572,7 +565,8 @@ mod tests { // reducing row 1, so the failure is detected via the `!sum.is_finite()` // branch of the combined diag/sum check (distinct from the // `q = sum / diag` overflow path covered above). - let a = Matrix::<3>::from_rows([[1.0, 0.0, 0.0], [0.0, 1.0, 1.0e200], [0.0, 0.0, 1.0]]); + let a = Matrix::<3>::try_from_rows([[1.0, 0.0, 0.0], [0.0, 1.0, 1.0e200], [0.0, 0.0, 1.0]]) + .unwrap(); let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<3>::new([0.0, 0.0, 1.0e200]); @@ -582,13 +576,14 @@ mod tests { #[test] fn det_rejects_product_overflow() { - let a = Matrix::<5>::from_rows([ + let a = Matrix::<5>::try_from_rows([ [1.0e100, 0.0, 0.0, 0.0, 0.0], [0.0, 1.0e100, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0e100, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0e100, 0.0], [0.0, 0.0, 0.0, 0.0, 1.0e100], - ]); + ]) + .unwrap(); let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); assert_eq!(lu.det(), Err(LaError::NonFinite { row: None, col: 3 })); } diff --git a/src/matrix.rs b/src/matrix.rs index 5995b20..10a5c9b 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -18,9 +18,10 @@ use crate::{ERR_COEFF_2, ERR_COEFF_3, ERR_COEFF_4, LDLT_SYMMETRY_REL_TOL, LaErro /// [`try_from_rows`](Self::try_from_rows), [`set`](Self::set), and /// [`set_checked`](Self::set_checked). The storage field is private, so a /// `Matrix` value carries the invariant that every stored entry is finite. -/// Algorithms therefore do not re-scan stored entries before using a `Matrix`; -/// they only report non-finite values computed during arithmetic, such as -/// overflowed elimination or determinant intermediates. +/// Algorithms therefore do not re-scan stored entries at every use; user-visible +/// non-finite errors come from construction/mutation boundaries or from values +/// computed during arithmetic, such as overflowed elimination or determinant +/// intermediates. /// /// Direct field construction is intentionally unavailable to downstream callers: /// @@ -37,313 +38,6 @@ pub struct Matrix { rows: [[f64; D]; D], } -/// Fixed-size square matrix whose stored entries are all finite. -/// -/// This internal proof wrapper is used where carrying the invariant explicitly -/// makes algorithm boundaries clearer. Public callers use [`Matrix`], which is -/// already finite by construction. -#[must_use] -#[derive(Clone, Copy, Debug, PartialEq)] -pub struct FiniteMatrix { - matrix: Matrix, -} - -impl FiniteMatrix { - /// Construct a finite matrix without checking the invariant. - /// - /// This crate-internal escape hatch is only for paths with a local proof - /// that every stored entry is finite. - #[inline] - pub(crate) const fn new_unchecked(matrix: Matrix) -> Self { - Self { matrix } - } - - /// Validate a matrix and wrap it for algorithms that carry the finite - /// invariant explicitly. - /// - /// # Errors - /// Returns [`LaError::NonFinite`] with matrix coordinates for the first - /// stored NaN or infinity in row-major order. - #[inline] - pub(crate) const fn new(matrix: Matrix) -> Result { - if let Some((row, col)) = Matrix::::first_non_finite_cell_in(&matrix.rows) { - Err(LaError::non_finite_cell(row, col)) - } else { - Ok(Self::new_unchecked(matrix)) - } - } - - /// Consume the wrapper and return the underlying raw matrix. - #[inline] - pub(crate) const fn into_matrix(self) -> Matrix { - self.matrix - } - - /// Borrow the underlying raw matrix without revalidating stored entries. - #[cfg(feature = "exact")] - #[inline] - pub(crate) const fn as_matrix(&self) -> &Matrix { - &self.matrix - } - - /// Infinity norm (maximum absolute row sum) for a finite matrix. - /// - /// Stored entries are known finite, so this path only checks whether a row - /// sum overflows to NaN or infinity. - /// - /// # Errors - /// Returns [`LaError::NonFinite`] when a row sum overflows to NaN or infinity. - #[inline] - pub(crate) const fn inf_norm(&self) -> Result { - let mut max_row_sum: f64 = 0.0; - - let mut r = 0; - while r < D { - let row = &self.matrix.rows[r]; - let mut row_sum: f64 = 0.0; - let mut c = 0; - while c < D { - row_sum += row[c].abs(); - c += 1; - } - if !row_sum.is_finite() { - cold_path(); - return Err(LaError::non_finite_at(r)); - } - if row_sum > max_row_sum { - max_row_sum = row_sum; - } - r += 1; - } - - Ok(max_row_sum) - } - - /// Returns `true` if the finite matrix is symmetric within a relative tolerance. - /// - /// # Errors - /// Returns [`LaError::NonFinite`] when computing the scaled symmetry tolerance - /// overflows to NaN or infinity. - #[inline] - pub(crate) fn is_symmetric(&self, rel_tol: Tolerance) -> Result { - Ok(self.first_asymmetry(rel_tol)?.is_none()) - } - - /// Returns the first asymmetric off-diagonal pair, if any. - /// - /// # Errors - /// Returns [`LaError::NonFinite`] when computing the scaled symmetry tolerance - /// overflows to NaN or infinity. - #[inline] - pub(crate) fn first_asymmetry( - &self, - rel_tol: Tolerance, - ) -> Result, LaError> { - let eps = self.symmetry_epsilon(rel_tol)?; - for r in 0..D { - for c in (r + 1)..D { - let upper = self.matrix.rows[r][c]; - let lower = self.matrix.rows[c][r]; - - let diff = (upper - lower).abs(); - if !diff.is_finite() || diff > eps { - cold_path(); - return Ok(Some((r, c))); - } - } - } - Ok(None) - } - - /// Compute the symmetry tolerance scale for a finite matrix. - fn symmetry_epsilon(&self, rel_tol: Tolerance) -> Result { - let rel_tol = rel_tol.get(); - let mut eps = rel_tol; - - for r in 0..D { - let mut row_eps = 0.0; - for c in 0..D { - row_eps = rel_tol.mul_add(self.matrix.rows[r][c].abs(), row_eps); - if !row_eps.is_finite() { - cold_path(); - return Err(LaError::non_finite_at(c)); - } - } - if row_eps > eps { - eps = row_eps; - } - } - - Ok(eps) - } - - /// Compute an LU decomposition with partial pivoting. - /// - /// # Errors - /// Returns [`LaError::Singular`] for an unusable pivot, or - /// [`LaError::NonFinite`] if elimination computes a non-finite intermediate. - #[inline] - pub(crate) fn lu(self, tol: Tolerance) -> Result, LaError> { - Lu::factor_finite(self, tol) - } - - /// Compute an LDLT factorization (`A = L D Lᵀ`) without pivoting. - /// - /// # Errors - /// Returns [`LaError::Asymmetric`] if the matrix is not symmetric, - /// [`LaError::NotPositiveSemidefinite`] for a negative LDLT diagonal, - /// [`LaError::Singular`] for a zero or too-small non-negative diagonal, or - /// [`LaError::NonFinite`] if factorization computes a non-finite intermediate. - #[inline] - pub(crate) fn ldlt(self, tol: Tolerance) -> Result, LaError> { - Ldlt::factor_symmetric(SymmetricMatrix::try_new(self)?, tol) - } - - /// Closed-form determinant for dimensions 0–4, bypassing LU factorization. - /// - /// # Errors - /// Returns [`LaError::NonFinite`] when the closed-form determinant overflows - /// to NaN or infinity. - #[inline] - pub(crate) const fn det_direct(&self) -> Result, LaError> { - match D { - 0 => Ok(Some(1.0)), - 1 => Self::computed_scalar_result(Some(self.matrix.rows[0][0])), - 2 => Self::computed_scalar_result(Some(self.matrix.rows[0][0].mul_add( - self.matrix.rows[1][1], - -(self.matrix.rows[0][1] * self.matrix.rows[1][0]), - ))), - 3 => { - let m00 = self.matrix.rows[1][1].mul_add( - self.matrix.rows[2][2], - -(self.matrix.rows[1][2] * self.matrix.rows[2][1]), - ); - let m01 = self.matrix.rows[1][0].mul_add( - self.matrix.rows[2][2], - -(self.matrix.rows[1][2] * self.matrix.rows[2][0]), - ); - let m02 = self.matrix.rows[1][0].mul_add( - self.matrix.rows[2][1], - -(self.matrix.rows[1][1] * self.matrix.rows[2][0]), - ); - Self::computed_scalar_result(Some(self.matrix.rows[0][0].mul_add( - m00, - (-self.matrix.rows[0][1]).mul_add(m01, self.matrix.rows[0][2] * m02), - ))) - } - 4 => { - let r = &self.matrix.rows; - - let s23 = r[2][2].mul_add(r[3][3], -(r[2][3] * r[3][2])); - let s13 = r[2][1].mul_add(r[3][3], -(r[2][3] * r[3][1])); - let s12 = r[2][1].mul_add(r[3][2], -(r[2][2] * r[3][1])); - let s03 = r[2][0].mul_add(r[3][3], -(r[2][3] * r[3][0])); - let s02 = r[2][0].mul_add(r[3][2], -(r[2][2] * r[3][0])); - let s01 = r[2][0].mul_add(r[3][1], -(r[2][1] * r[3][0])); - - let c00 = r[1][1].mul_add(s23, (-r[1][2]).mul_add(s13, r[1][3] * s12)); - let c01 = r[1][0].mul_add(s23, (-r[1][2]).mul_add(s03, r[1][3] * s02)); - let c02 = r[1][0].mul_add(s13, (-r[1][1]).mul_add(s03, r[1][3] * s01)); - let c03 = r[1][0].mul_add(s12, (-r[1][1]).mul_add(s02, r[1][2] * s01)); - - Self::computed_scalar_result(Some(r[0][0].mul_add( - c00, - (-r[0][1]).mul_add(c01, r[0][2].mul_add(c02, -(r[0][3] * c03))), - ))) - } - _ => { - cold_path(); - Ok(None) - } - } - } - - /// Floating-point determinant for a finite matrix. - /// - /// # Errors - /// Returns [`LaError::NonFinite`] if a computed determinant or LU fallback - /// intermediate overflows to NaN or infinity. - #[inline] - pub(crate) fn det(self) -> Result { - if let Some(d) = self.det_direct()? { - return Ok(d); - } - match self.lu(Tolerance::new_unchecked(0.0)) { - Ok(lu) => lu.det(), - Err(LaError::Singular { .. }) => Ok(0.0), - Err(err) => Err(err), - } - } - - /// Conservative absolute error bound for [`det_direct`](Self::det_direct). - /// - /// # Errors - /// Returns [`LaError::NonFinite`] when the bound computation overflows to NaN - /// or infinity. - #[inline] - pub(crate) const fn det_errbound(&self) -> Result, LaError> { - match D { - 0 | 1 => Self::computed_scalar_result(Some(0.0)), - 2 => { - let r = &self.matrix.rows; - let permanent = (r[0][0] * r[1][1]).abs() + (r[0][1] * r[1][0]).abs(); - Self::computed_scalar_result(Some(ERR_COEFF_2 * permanent)) - } - 3 => { - let r = &self.matrix.rows; - let pm00 = (r[1][1] * r[2][2]).abs() + (r[1][2] * r[2][1]).abs(); - let pm01 = (r[1][0] * r[2][2]).abs() + (r[1][2] * r[2][0]).abs(); - let pm02 = (r[1][0] * r[2][1]).abs() + (r[1][1] * r[2][0]).abs(); - let permanent = r[0][2] - .abs() - .mul_add(pm02, r[0][1].abs().mul_add(pm01, r[0][0].abs() * pm00)); - Self::computed_scalar_result(Some(ERR_COEFF_3 * permanent)) - } - 4 => { - let r = &self.matrix.rows; - let sp23 = (r[2][2] * r[3][3]).abs() + (r[2][3] * r[3][2]).abs(); - let sp13 = (r[2][1] * r[3][3]).abs() + (r[2][3] * r[3][1]).abs(); - let sp12 = (r[2][1] * r[3][2]).abs() + (r[2][2] * r[3][1]).abs(); - let sp03 = (r[2][0] * r[3][3]).abs() + (r[2][3] * r[3][0]).abs(); - let sp02 = (r[2][0] * r[3][2]).abs() + (r[2][2] * r[3][0]).abs(); - let sp01 = (r[2][0] * r[3][1]).abs() + (r[2][1] * r[3][0]).abs(); - let pc0 = r[1][3] - .abs() - .mul_add(sp12, r[1][2].abs().mul_add(sp13, r[1][1].abs() * sp23)); - let pc1 = r[1][3] - .abs() - .mul_add(sp02, r[1][2].abs().mul_add(sp03, r[1][0].abs() * sp23)); - let pc2 = r[1][3] - .abs() - .mul_add(sp01, r[1][1].abs().mul_add(sp03, r[1][0].abs() * sp13)); - let pc3 = r[1][2] - .abs() - .mul_add(sp01, r[1][1].abs().mul_add(sp02, r[1][0].abs() * sp12)); - let permanent = r[0][3].abs().mul_add( - pc3, - r[0][2] - .abs() - .mul_add(pc2, r[0][1].abs().mul_add(pc1, r[0][0].abs() * pc0)), - ); - Self::computed_scalar_result(Some(ERR_COEFF_4 * permanent)) - } - _ => { - cold_path(); - Ok(None) - } - } - } - - /// Return a computed scalar result for a matrix with finite stored entries. - const fn computed_scalar_result(value: Option) -> Result, LaError> { - match value { - Some(value) if value.is_finite() => Ok(Some(value)), - Some(_) => Err(LaError::non_finite_at(0)), - None => Ok(None), - } - } -} - /// Matrix proven finite and symmetric under the crate's LDLT symmetry tolerance. #[must_use] #[derive(Clone, Copy, Debug, PartialEq)] @@ -363,7 +57,7 @@ impl SymmetricMatrix { Self { matrix } } - /// Validate that a finite matrix is symmetric under the LDLT symmetry tolerance. + /// Validate that a matrix is symmetric under the LDLT symmetry tolerance. /// /// The predicate is the same one used by [`Matrix::ldlt`]: /// `|A[i][j] - A[j][i]| <= 1e-12 * max(1, inf_norm(A))`, with scaling that @@ -376,12 +70,12 @@ impl SymmetricMatrix { /// Returns [`LaError::NonFinite`] when computing the scaled symmetry /// tolerance overflows to NaN or infinity. #[inline] - pub(crate) fn try_new(matrix: FiniteMatrix) -> Result { + pub(crate) fn try_new(matrix: Matrix) -> Result { if let Some((row, col)) = matrix.first_asymmetry(LDLT_SYMMETRY_REL_TOL)? { cold_path(); Err(LaError::asymmetric(row, col, D)) } else { - Ok(Self::new_unchecked(matrix.into_matrix())) + Ok(Self::new_unchecked(matrix)) } } @@ -393,16 +87,6 @@ impl SymmetricMatrix { } impl Matrix { - /// Test-only infallible constructor for finite literal fixtures. - #[cfg(test)] - #[inline] - pub(crate) const fn from_rows(rows: [[f64; D]; D]) -> Self { - match Self::try_from_rows(rows) { - Ok(matrix) => matrix, - Err(_) => panic!("Matrix::from_rows requires finite entries"), - } - } - /// Try to create a finite matrix from row-major storage. /// /// This is the public raw-storage boundary for matrices. Successful @@ -633,7 +317,8 @@ impl Matrix { /// # Non-finite handling /// Public constructors and setters reject raw non-finite entries, but /// `Matrix` values are finite by construction. `inf_norm` returns - /// [`LaError::NonFinite`] if a row sum overflows to a non-finite value. + /// [`LaError::NonFinite`] with the matrix cell whose addition first makes a + /// row sum non-finite. /// /// Row sums are accumulated in `f64` with ordinary addition. This method /// checks for overflowed accumulators, but it does not provide a certified @@ -660,10 +345,32 @@ impl Matrix { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] when a row sum overflows to NaN or infinity. + /// Returns [`LaError::NonFinite`] with matrix coordinates when a row sum + /// overflows to NaN or infinity. #[inline] pub const fn inf_norm(&self) -> Result { - FiniteMatrix::new_unchecked(*self).inf_norm() + let mut max_row_sum: f64 = 0.0; + + let mut r = 0; + while r < D { + let row = &self.rows[r]; + let mut row_sum: f64 = 0.0; + let mut c = 0; + while c < D { + row_sum += row[c].abs(); + if !row_sum.is_finite() { + cold_path(); + return Err(LaError::non_finite_cell(r, c)); + } + c += 1; + } + if row_sum > max_row_sum { + max_row_sum = row_sum; + } + r += 1; + } + + Ok(max_row_sum) } /// Returns `true` if the matrix is symmetric within a relative tolerance. @@ -684,10 +391,10 @@ impl Matrix { /// [`LaError::InvalidTolerance`]. /// /// # Overflow handling - /// A finite matrix can return [`LaError::NonFinite`] if computing the scaled - /// symmetry tolerance overflows to NaN or infinity. If both stored entries - /// are finite but their difference overflows to ±∞, the pair is reported as - /// asymmetric. + /// A finite matrix can return [`LaError::NonFinite`] with matrix coordinates + /// if computing the scaled symmetry tolerance overflows to NaN or infinity. + /// If both stored entries are finite but their difference overflows to ±∞, + /// the pair is reported as asymmetric. /// /// # Examples /// ``` @@ -705,11 +412,11 @@ impl Matrix { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] when computing the scaled symmetry - /// tolerance overflows to NaN or infinity. + /// Returns [`LaError::NonFinite`] with matrix coordinates when computing the + /// scaled symmetry tolerance overflows to NaN or infinity. #[inline] pub fn is_symmetric(&self, rel_tol: Tolerance) -> Result { - FiniteMatrix::new_unchecked(*self).is_symmetric(rel_tol) + Ok(self.first_asymmetry(rel_tol)?.is_none()) } /// Returns the indices `(r, c)` (with `r < c`) of the first off-diagonal @@ -721,10 +428,10 @@ impl Matrix { /// predicate is the same as [`is_symmetric`](Self::is_symmetric): /// `|self[r][c] - self[c][r]| <= rel_tol * max(1.0, inf_norm(self))`. /// - /// A finite matrix can return [`LaError::NonFinite`] if computing the scaled - /// symmetry tolerance overflows to NaN or infinity. If both stored entries - /// are finite but their difference overflows to ±∞, the pair is reported as - /// asymmetric. + /// A finite matrix can return [`LaError::NonFinite`] with matrix coordinates + /// if computing the scaled symmetry tolerance overflows to NaN or infinity. + /// If both stored entries are finite but their difference overflows to ±∞, + /// the pair is reported as asymmetric. /// /// The `rel_tol` argument is a [`Tolerance`], so raw caller input must be /// finite and non-negative before it can reach this predicate. Use @@ -750,11 +457,24 @@ impl Matrix { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] when computing the scaled symmetry - /// tolerance overflows to NaN or infinity. + /// Returns [`LaError::NonFinite`] with matrix coordinates when computing the + /// scaled symmetry tolerance overflows to NaN or infinity. #[inline] pub fn first_asymmetry(&self, rel_tol: Tolerance) -> Result, LaError> { - FiniteMatrix::new_unchecked(*self).first_asymmetry(rel_tol) + let eps = self.symmetry_epsilon(rel_tol)?; + for r in 0..D { + for c in (r + 1)..D { + let upper = self.rows[r][c]; + let lower = self.rows[c][r]; + + let diff = (upper - lower).abs(); + if !diff.is_finite() || diff > eps { + cold_path(); + return Ok(Some((r, c))); + } + } + } + Ok(None) } /// Compute an LU decomposition with partial pivoting. @@ -789,7 +509,7 @@ impl Matrix { /// to NaN/∞ before it can be stored in the returned [`Lu`]. #[inline] pub fn lu(self, tol: Tolerance) -> Result, LaError> { - FiniteMatrix::new_unchecked(self).lu(tol) + Lu::factor_finite(self, tol) } /// Compute an LDLT factorization (`A = L D Lᵀ`) without pivoting. @@ -843,7 +563,7 @@ impl Matrix { /// Returns [`LaError::Asymmetric`] if the input matrix is not symmetric. #[inline] pub fn ldlt(self, tol: Tolerance) -> Result, LaError> { - FiniteMatrix::new_unchecked(self).ldlt(tol) + Ldlt::factor_symmetric(SymmetricMatrix::try_new(self)?, tol) } /// Return the first non-finite stored cell in row-major order. @@ -862,6 +582,47 @@ impl Matrix { None } + /// Validate storage after unchecked internal construction or mutation. + /// + /// Public constructors and setters already maintain this invariant. This + /// helper is reserved for internal factorization temporaries and test + /// fixtures that intentionally bypass those boundaries. + #[inline] + pub(crate) const fn validate_finite(self) -> Result { + if let Some((row, col)) = Self::first_non_finite_cell_in(&self.rows) { + Err(LaError::non_finite_cell(row, col)) + } else { + Ok(self) + } + } + + /// Compute the symmetry tolerance scale for a finite matrix. + /// + /// This helper protects the public [`is_symmetric`](Self::is_symmetric), + /// [`first_asymmetry`](Self::first_asymmetry), and [`ldlt`](Self::ldlt) + /// error contracts: an overflowed row-scale accumulator is reported with + /// the matrix cell whose contribution made it non-finite. + fn symmetry_epsilon(&self, rel_tol: Tolerance) -> Result { + let rel_tol = rel_tol.get(); + let mut eps = rel_tol; + + for r in 0..D { + let mut row_eps = 0.0; + for c in 0..D { + row_eps = rel_tol.mul_add(self.rows[r][c].abs(), row_eps); + if !row_eps.is_finite() { + cold_path(); + return Err(LaError::non_finite_cell(r, c)); + } + } + if row_eps > eps { + eps = row_eps; + } + } + + Ok(eps) + } + /// Closed-form determinant for dimensions 0–4, bypassing LU factorization. /// /// Returns `Ok(Some(det))` for `D` ∈ {0, 1, 2, 3, 4}, `Ok(None)` for D ≥ 5. @@ -894,7 +655,70 @@ impl Matrix { /// to NaN or infinity. #[inline] pub const fn det_direct(&self) -> Result, LaError> { - FiniteMatrix::new_unchecked(*self).det_direct() + match D { + 0 => Ok(Some(1.0)), + 1 => Self::computed_scalar_result(self.rows[0][0]), + 2 => { + let det = if self.rows[0][1] == 0.0 { + self.rows[0][0] * self.rows[1][1] + } else { + self.rows[0][0].mul_add(self.rows[1][1], -(self.rows[0][1] * self.rows[1][0])) + }; + Self::computed_scalar_result(det) + } + 3 => { + let det = Self::det3_elements( + [self.rows[0][0], self.rows[0][1], self.rows[0][2]], + [self.rows[1][0], self.rows[1][1], self.rows[1][2]], + [self.rows[2][0], self.rows[2][1], self.rows[2][2]], + ); + Self::computed_scalar_result(det) + } + 4 => { + let r = &self.rows; + + let mut det = if r[0][3] == 0.0 { + 0.0 + } else { + let c03 = Self::det3_elements( + [r[1][0], r[1][1], r[1][2]], + [r[2][0], r[2][1], r[2][2]], + [r[3][0], r[3][1], r[3][2]], + ); + -(r[0][3] * c03) + }; + if r[0][2] != 0.0 { + let c02 = Self::det3_elements( + [r[1][0], r[1][1], r[1][3]], + [r[2][0], r[2][1], r[2][3]], + [r[3][0], r[3][1], r[3][3]], + ); + det = r[0][2].mul_add(c02, det); + } + if r[0][1] != 0.0 { + let c01 = Self::det3_elements( + [r[1][0], r[1][2], r[1][3]], + [r[2][0], r[2][2], r[2][3]], + [r[3][0], r[3][2], r[3][3]], + ); + det = (-r[0][1]).mul_add(c01, det); + } + if r[0][0] != 0.0 { + let c00 = Self::det3_elements( + [r[1][1], r[1][2], r[1][3]], + [r[2][1], r[2][2], r[2][3]], + [r[3][1], r[3][2], r[3][3]], + ); + det = r[0][0].mul_add(c00, det); + } + + Self::computed_scalar_result(det) + } + _ => { + cold_path(); + Ok(None) + } + } } /// Floating-point determinant, using closed-form formulas for D ≤ 4 and @@ -929,7 +753,14 @@ impl Matrix { /// factorization cell, or the determinant product overflows to NaN or infinity. #[inline] pub fn det(self) -> Result { - FiniteMatrix::new_unchecked(self).det() + if let Some(d) = self.det_direct()? { + return Ok(d); + } + match self.lu(Tolerance::new_unchecked(0.0)) { + Ok(lu) => lu.det(), + Err(LaError::Singular { .. }) => Ok(0.0), + Err(err) => Err(err), + } } /// Conservative absolute error bound for `det_direct()`. @@ -997,7 +828,121 @@ impl Matrix { /// NaN or infinity. #[inline] pub const fn det_errbound(&self) -> Result, LaError> { - FiniteMatrix::new_unchecked(*self).det_errbound() + match D { + 0 | 1 => Self::computed_scalar_result(0.0), + 2 => { + let r = &self.rows; + let permanent = (r[0][0] * r[1][1]).abs() + (r[0][1] * r[1][0]).abs(); + Self::computed_scalar_result(ERR_COEFF_2 * permanent) + } + 3 => { + let r = &self.rows; + let permanent = Self::det3_abs_permanent_elements( + [r[0][0], r[0][1], r[0][2]], + [r[1][0], r[1][1], r[1][2]], + [r[2][0], r[2][1], r[2][2]], + ); + Self::computed_scalar_result(ERR_COEFF_3 * permanent) + } + 4 => { + let r = &self.rows; + let mut permanent = if r[0][3] == 0.0 { + 0.0 + } else { + let pc3 = Self::det3_abs_permanent_elements( + [r[1][0], r[1][1], r[1][2]], + [r[2][0], r[2][1], r[2][2]], + [r[3][0], r[3][1], r[3][2]], + ); + r[0][3].abs() * pc3 + }; + if r[0][2] != 0.0 { + let pc2 = Self::det3_abs_permanent_elements( + [r[1][0], r[1][1], r[1][3]], + [r[2][0], r[2][1], r[2][3]], + [r[3][0], r[3][1], r[3][3]], + ); + permanent = r[0][2].abs().mul_add(pc2, permanent); + } + if r[0][1] != 0.0 { + let pc1 = Self::det3_abs_permanent_elements( + [r[1][0], r[1][2], r[1][3]], + [r[2][0], r[2][2], r[2][3]], + [r[3][0], r[3][2], r[3][3]], + ); + permanent = r[0][1].abs().mul_add(pc1, permanent); + } + if r[0][0] != 0.0 { + let pc0 = Self::det3_abs_permanent_elements( + [r[1][1], r[1][2], r[1][3]], + [r[2][1], r[2][2], r[2][3]], + [r[3][1], r[3][2], r[3][3]], + ); + permanent = r[0][0].abs().mul_add(pc0, permanent); + } + Self::computed_scalar_result(ERR_COEFF_4 * permanent) + } + _ => { + cold_path(); + Ok(None) + } + } + } + + /// Evaluate a 3×3 determinant expansion while skipping zero coefficients. + /// + /// This helper protects the public [`det_direct`](Self::det_direct) + /// contract: a mathematically absent term must not compute an overflowing + /// minor and poison the determinant with `0.0 * inf == NaN`. Nonzero terms + /// keep the same fused multiply-add ordering as the closed-form expansion. + const fn det3_elements(r0: [f64; 3], r1: [f64; 3], r2: [f64; 3]) -> f64 { + let mut det = if r0[2] == 0.0 { + 0.0 + } else { + let m02 = r1[0].mul_add(r2[1], -(r1[1] * r2[0])); + r0[2] * m02 + }; + if r0[1] != 0.0 { + let m01 = r1[0].mul_add(r2[2], -(r1[2] * r2[0])); + det = (-r0[1]).mul_add(m01, det); + } + if r0[0] != 0.0 { + let m00 = r1[1].mul_add(r2[2], -(r1[2] * r2[1])); + det = r0[0].mul_add(m00, det); + } + det + } + + /// Evaluate a 3×3 absolute permanent while skipping zero coefficients. + /// + /// This mirrors [`det3_elements`](Self::det3_elements) for error-bound + /// computation: absent determinant terms should not force evaluation of an + /// overflowing absolute minor. + const fn det3_abs_permanent_elements(r0: [f64; 3], r1: [f64; 3], r2: [f64; 3]) -> f64 { + let mut permanent = if r0[2] == 0.0 { + 0.0 + } else { + let pm02 = (r1[0] * r2[1]).abs() + (r1[1] * r2[0]).abs(); + r0[2].abs() * pm02 + }; + if r0[1] != 0.0 { + let pm01 = (r1[0] * r2[2]).abs() + (r1[2] * r2[0]).abs(); + permanent = r0[1].abs().mul_add(pm01, permanent); + } + if r0[0] != 0.0 { + let pm00 = (r1[1] * r2[2]).abs() + (r1[2] * r2[1]).abs(); + permanent = r0[0].abs().mul_add(pm00, permanent); + } + permanent + } + + /// Return a computed scalar result for a matrix with finite stored entries. + const fn computed_scalar_result(value: f64) -> Result, LaError> { + if value.is_finite() { + Ok(Some(value)) + } else { + Err(LaError::non_finite_at(0)) + } } } @@ -1012,41 +957,23 @@ impl Default for Matrix { mod tests { use super::*; use crate::DEFAULT_SINGULAR_TOL; - use crate::vector::{FiniteVector, Vector}; + use crate::vector::Vector; use approx::assert_abs_diff_eq; use core::assert_matches; use pastey::paste; use std::hint::black_box; - fn assert_matrix_abs_eq(actual: &Matrix, expected: &Matrix) { - for r in 0..D { - for c in 0..D { - assert_abs_diff_eq!( - actual.get(r, c).unwrap(), - expected.get(r, c).unwrap(), - epsilon = 0.0 - ); - } - } - } - - fn assert_array_abs_eq(actual: &[f64; D], expected: &[f64; D]) { - for i in 0..D { - assert_abs_diff_eq!(actual[i], expected[i], epsilon = 0.0); - } - } - macro_rules! gen_public_api_matrix_tests { ($d:literal) => { paste! { #[test] - fn []() { + fn []() { let mut rows = [[0.0f64; $d]; $d]; rows[0][0] = 1.0; rows[$d - 1][$d - 1] = -2.0; - let mut m = Matrix::<$d>::from_rows(rows); + let mut m = Matrix::<$d>::try_from_rows(rows).unwrap(); assert_eq!(m.get(0, 0), Some(1.0)); assert_eq!(m.get($d - 1, $d - 1), Some(-2.0)); @@ -1125,7 +1052,7 @@ mod tests { rows[1][c] = 0.5; } - let m = Matrix::<$d>::from_rows(rows); + let m = Matrix::<$d>::try_from_rows(rows).unwrap(); assert_abs_diff_eq!(m.inf_norm().unwrap(), f64::from($d), epsilon = 0.0); } @@ -1166,91 +1093,19 @@ mod tests { } #[test] - fn []() { - let mut rows = [[0.0f64; $d]; $d]; - for r in 0..$d { - rows[r][r] = 1.0; - } - - let finite = FiniteMatrix::<$d>::new(Matrix::<$d>::from_rows(rows)).unwrap(); - - assert_matrix_abs_eq(&finite.into_matrix(), &Matrix::<$d>::from_rows(rows)); - assert_eq!(finite.into_matrix().get(0, 0), Some(1.0)); - assert_eq!(finite.into_matrix().get($d, 0), None); - assert_eq!( - finite.into_matrix().get_checked($d, 0), - Err(LaError::IndexOutOfBounds { - row: $d, - col: 0, - dim: $d, - }) - ); - } - - #[test] - fn []() { + fn []() { let mut rows = [[1.0f64; $d]; $d]; rows[$d - 1][0] = f64::NAN; let raw = Matrix::<$d>::from_rows_unchecked(rows); assert_eq!( - FiniteMatrix::<$d>::new(raw), + raw.validate_finite(), Err(LaError::NonFinite { row: Some($d - 1), col: 0, }) ); } - - #[test] - fn []() { - let mut rows = [[0.0f64; $d]; $d]; - let diag_values = [2.0f64, 3.0, 5.0, 7.0, 11.0]; - for i in 0..$d { - rows[i][i] = diag_values[i]; - } - - let raw = Matrix::<$d>::from_rows(rows); - let finite = FiniteMatrix::<$d>::new(raw).unwrap(); - - assert_abs_diff_eq!(finite.inf_norm().unwrap(), raw.inf_norm().unwrap(), epsilon = 0.0); - assert_eq!(finite.det_direct(), raw.det_direct()); - assert_abs_diff_eq!(finite.det().unwrap(), raw.det().unwrap(), epsilon = 0.0); - assert_eq!(finite.det_errbound(), raw.det_errbound()); - assert_eq!( - finite.first_asymmetry(Tolerance::new(1e-12).unwrap()).unwrap(), - raw.first_asymmetry(Tolerance::new(1e-12).unwrap()).unwrap() - ); - assert_eq!( - finite.is_symmetric(Tolerance::new(1e-12).unwrap()).unwrap(), - raw.is_symmetric(Tolerance::new(1e-12).unwrap()).unwrap() - ); - } - - #[test] - fn []() { - let finite = FiniteMatrix::<$d>::new(Matrix::<$d>::identity()).unwrap(); - let rhs = { - let mut arr = [0.0f64; $d]; - let values = [1.0f64, 2.0, 3.0, 4.0, 5.0]; - for (dst, src) in arr.iter_mut().zip(values.iter()) { - *dst = *src; - } - FiniteVector::<$d>::new_unchecked(Vector::<$d>::new(arr)) - }; - - let lu = finite.lu(DEFAULT_SINGULAR_TOL).unwrap(); - assert_array_abs_eq( - &lu.solve_finite(rhs).unwrap().into_array(), - &rhs.into_array() - ); - - let ldlt = finite.ldlt(DEFAULT_SINGULAR_TOL).unwrap(); - assert_array_abs_eq( - &ldlt.solve_finite(rhs).unwrap().into_array(), - &rhs.into_array() - ); - } } }; } @@ -1270,7 +1125,7 @@ mod tests { #[test] fn det_direct_d1_returns_element() { - let m = Matrix::<1>::from_rows([[42.0]]); + let m = Matrix::<1>::try_from_rows([[42.0]]).unwrap(); assert_eq!(m.det_direct(), Ok(Some(42.0))); } @@ -1278,32 +1133,43 @@ mod tests { fn det_direct_d2_known_value() { // [[1,2],[3,4]] → det = 1*4 - 2*3 = -2 // black_box prevents compile-time constant folding of the const fn. - let m = black_box(Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]])); + let m = black_box(Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]]).unwrap()); assert_abs_diff_eq!(m.det_direct().unwrap().unwrap(), -2.0, epsilon = 1e-15); } #[test] fn det_direct_d3_known_value() { // Classic 3×3: det = 0 - let m = black_box(Matrix::<3>::from_rows([ - [1.0, 2.0, 3.0], - [4.0, 5.0, 6.0], - [7.0, 8.0, 9.0], - ])); + let m = black_box( + Matrix::<3>::try_from_rows([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]) + .unwrap(), + ); assert_abs_diff_eq!(m.det_direct().unwrap().unwrap(), 0.0, epsilon = 1e-12); } #[test] fn det_direct_d3_nonsingular() { // [[2,1,0],[0,3,1],[1,0,2]] → det = 2*(6-0) - 1*(0-1) + 0 = 13 - let m = black_box(Matrix::<3>::from_rows([ - [2.0, 1.0, 0.0], - [0.0, 3.0, 1.0], - [1.0, 0.0, 2.0], - ])); + let m = black_box( + Matrix::<3>::try_from_rows([[2.0, 1.0, 0.0], [0.0, 3.0, 1.0], [1.0, 0.0, 2.0]]) + .unwrap(), + ); assert_abs_diff_eq!(m.det_direct().unwrap().unwrap(), 13.0, epsilon = 1e-12); } + #[test] + fn det_direct_d3_skips_zero_coefficient_minor_that_would_overflow() { + let m = black_box( + Matrix::<3>::try_from_rows([ + [1.0, 0.0, 0.0], + [1.0e300, 1.0, 1.0e300], + [1.0e300, 0.0, 1.0e300], + ]) + .unwrap(), + ); + assert_eq!(m.det_direct(), Ok(Some(1.0e300))); + } + #[test] fn det_direct_d4_identity() { let m = black_box(Matrix::<4>::identity()); @@ -1318,10 +1184,38 @@ mod tests { rows[1][1] = 3.0; rows[2][2] = 5.0; rows[3][3] = 7.0; - let m = black_box(Matrix::<4>::from_rows(rows)); + let m = black_box(Matrix::<4>::try_from_rows(rows).unwrap()); assert_abs_diff_eq!(m.det_direct().unwrap().unwrap(), 210.0, epsilon = 1e-12); } + #[test] + fn det_direct_d4_dense_known_value() { + let m = black_box( + Matrix::<4>::try_from_rows([ + [4.0, 1.0, 3.0, 2.0], + [0.0, 5.0, 2.0, 1.0], + [7.0, 2.0, 6.0, 3.0], + [1.0, 8.0, 4.0, 9.0], + ]) + .unwrap(), + ); + assert_abs_diff_eq!(m.det_direct().unwrap().unwrap(), 92.0, epsilon = 1e-12); + } + + #[test] + fn det_direct_d4_skips_zero_coefficient_cofactors_that_would_overflow() { + let m = black_box( + Matrix::<4>::try_from_rows([ + [0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0], + [1.0e300, 0.0, 1.0e300, 1.0e300], + [1.0e300, 0.0, 1.0e300, -1.0e300], + ]) + .unwrap(), + ); + assert_eq!(m.det_direct(), Ok(Some(0.0))); + } + #[test] fn det_direct_d5_returns_none() { assert_eq!(Matrix::<5>::identity().det_direct(), Ok(None)); @@ -1357,7 +1251,7 @@ mod tests { #[test] fn det_direct_rejects_computed_overflow() { - let m = Matrix::<2>::from_rows([[1e300, 0.0], [0.0, 1e300]]); + let m = Matrix::<2>::try_from_rows([[1e300, 0.0], [0.0, 1e300]]).unwrap(); assert_eq!( m.det_direct(), Err(LaError::NonFinite { row: None, col: 0 }) @@ -1366,25 +1260,27 @@ mod tests { #[test] fn det_d5_rejects_lu_product_overflow() { - let m = Matrix::<5>::from_rows([ + let m = Matrix::<5>::try_from_rows([ [1.0e100, 0.0, 0.0, 0.0, 0.0], [0.0, 1.0e100, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0e100, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0e100, 0.0], [0.0, 0.0, 0.0, 0.0, 1.0e100], - ]); + ]) + .unwrap(); assert_eq!(m.det(), Err(LaError::NonFinite { row: None, col: 3 })); } #[test] fn det_d5_rejects_lu_trailing_update_overflow() { - let m = Matrix::<5>::from_rows([ + let m = Matrix::<5>::try_from_rows([ [1.0, f64::MAX, 0.0, 0.0, 0.0], [-1.0, f64::MAX, 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], - ]); + ]) + .unwrap(); assert_eq!( m.det(), @@ -1412,7 +1308,7 @@ mod tests { }; } } - let m = Matrix::<$d>::from_rows(rows); + let m = Matrix::<$d>::try_from_rows(rows).unwrap(); let direct = m.det_direct().unwrap().unwrap(); let lu_det = m.lu(DEFAULT_SINGULAR_TOL).unwrap().det().unwrap(); let eps = lu_det.abs().mul_add(1e-12, 1e-12); @@ -1495,13 +1391,14 @@ mod tests { // A small nonzero determinant is still a determinant. `det` must not // flatten the value to zero merely because the default LU tolerance // would reject a pivot this small. - let m = Matrix::<5>::from_rows([ + let m = Matrix::<5>::try_from_rows([ [1e-13, 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, 0.0], [0.0, 0.0, 0.0, 0.0, 1.0], - ]); + ]) + .unwrap(); assert_abs_diff_eq!(m.det().unwrap(), 1e-13, epsilon = 0.0); assert_eq!( @@ -1542,7 +1439,7 @@ mod tests { // though every matrix entry is finite. The entry scan in `det` // falls through and returns NonFinite { row: None, col: 0 } to signal // a computed overflow rather than a NaN/∞ input. - let m = Matrix::<2>::from_rows([[1e300, 0.0], [0.0, 1e300]]); + let m = Matrix::<2>::try_from_rows([[1e300, 0.0], [0.0, 1e300]]).unwrap(); assert_eq!(m.det(), Err(LaError::NonFinite { row: None, col: 0 })); } @@ -1591,7 +1488,7 @@ mod tests { #[test] fn det_errbound_matches_documented_coefficient_scale() { - let m2 = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); + let m2 = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]]).unwrap(); let expected_2 = ERR_COEFF_2 * ((1.0_f64 * 4.0).abs() + (2.0_f64 * 3.0).abs()); assert_abs_diff_eq!( m2.det_errbound().unwrap().unwrap(), @@ -1611,6 +1508,31 @@ mod tests { ); } + #[test] + fn det_errbound_d3_skips_zero_coefficient_minor_that_would_overflow() { + let m = Matrix::<3>::try_from_rows([ + [1.0, 0.0, 0.0], + [1.0e300, 1.0, 1.0e300], + [1.0e300, 0.0, 1.0e300], + ]) + .unwrap(); + + assert_eq!(m.det_errbound(), Ok(Some(ERR_COEFF_3 * 1.0e300))); + } + + #[test] + fn det_errbound_d4_skips_zero_coefficient_cofactors_that_would_overflow() { + let m = Matrix::<4>::try_from_rows([ + [0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0], + [1.0e300, 0.0, 1.0e300, 1.0e300], + [1.0e300, 0.0, 1.0e300, -1.0e300], + ]) + .unwrap(); + + assert_eq!(m.det_errbound(), Ok(Some(0.0))); + } + #[test] fn det_errbound_d5_returns_none() { // D=5 has no fast filter @@ -1653,7 +1575,7 @@ mod tests { #[test] fn det_errbound_rejects_computed_overflow() { - let m = Matrix::<2>::from_rows([[1e300, 0.0], [0.0, 1e300]]); + let m = Matrix::<2>::try_from_rows([[1e300, 0.0], [0.0, 1e300]]).unwrap(); assert_eq!( m.det_errbound(), Err(LaError::NonFinite { row: None, col: 0 }) @@ -1803,7 +1725,7 @@ mod tests { sym[r][c] = m[r][c] + m[c][r]; } } - let a = Matrix::<$d>::from_rows(sym); + let a = Matrix::<$d>::try_from_rows(sym).unwrap(); assert!(a.is_symmetric(Tolerance::new(1e-12).unwrap()).unwrap()); assert_eq!(a.first_asymmetry(Tolerance::new(1e-12).unwrap()).unwrap(), None); } @@ -1817,7 +1739,7 @@ mod tests { } rows[0][$d - 1] = 1.0; rows[$d - 1][0] = -1.0; // breaks symmetry - let a = Matrix::<$d>::from_rows(rows); + let a = Matrix::<$d>::try_from_rows(rows).unwrap(); assert!(!a.is_symmetric(Tolerance::new(1e-12).unwrap()).unwrap()); assert_eq!( a.first_asymmetry(Tolerance::new(1e-12).unwrap()).unwrap(), @@ -1871,9 +1793,7 @@ mod tests { rows[$d - 1][0] = -1.0; assert_eq!( - Matrix::<$d>::try_from_rows(rows) - .and_then(FiniteMatrix::new) - .and_then(SymmetricMatrix::try_new), + Matrix::<$d>::try_from_rows(rows).and_then(SymmetricMatrix::try_new), Err(LaError::Asymmetric { row: 0, col: $d - 1, @@ -1903,8 +1823,8 @@ mod tests { #[test] fn symmetric_matrix_into_matrix_roundtrips_storage_internally() { - let a = Matrix::<2>::from_rows([[2.0, 1.0], [1.0, 3.0]]); - let symmetric = SymmetricMatrix::try_new(FiniteMatrix::new(a).unwrap()).unwrap(); + let a = Matrix::<2>::try_from_rows([[2.0, 1.0], [1.0, 3.0]]).unwrap(); + let symmetric = SymmetricMatrix::try_new(a).unwrap(); assert_eq!(symmetric.into_matrix(), a); } @@ -1914,7 +1834,7 @@ mod tests { // Off-diagonal entries differ by 1e-6. With inf_norm ≈ 2e6, the // relative tolerance 1e-12 yields eps ≈ 2e-6, which accepts the gap; // a stricter tol of 1e-15 rejects it. - let a = Matrix::<2>::from_rows([[1.0e6, 1.0e6 + 1.0e-6], [1.0e6, 1.0e6]]); + let a = Matrix::<2>::try_from_rows([[1.0e6, 1.0e6 + 1.0e-6], [1.0e6, 1.0e6]]).unwrap(); assert!(a.is_symmetric(Tolerance::new(1e-12).unwrap()).unwrap()); assert!(!a.is_symmetric(Tolerance::new(1e-15).unwrap()).unwrap()); } @@ -1922,7 +1842,8 @@ mod tests { #[test] fn first_asymmetry_returns_lexicographically_first_pair() { // Two asymmetric pairs: (0, 2) and (1, 2). We must get (0, 2) first. - let a = Matrix::<3>::from_rows([[1.0, 0.0, 2.0], [0.0, 1.0, 3.0], [-2.0, -3.0, 1.0]]); + let a = Matrix::<3>::try_from_rows([[1.0, 0.0, 2.0], [0.0, 1.0, 3.0], [-2.0, -3.0, 1.0]]) + .unwrap(); assert_eq!( a.first_asymmetry(Tolerance::new(1e-12).unwrap()).unwrap(), Some((0, 2)) @@ -1953,13 +1874,20 @@ mod tests { #[test] fn first_asymmetry_strict_tol_survives_row_sum_overflow() { - let a = Matrix::<3>::from_rows([ - [f64::MAX, 1.0, f64::MAX], - [2.0, f64::MAX, 0.0], - [f64::MAX, 0.0, f64::MAX], - ]); + let a = Matrix::<3>::try_from_rows([ + [1.0, 1.0, 0.0], + [2.0, f64::MAX, f64::MAX], + [0.0, 0.0, 1.0], + ]) + .unwrap(); - assert_eq!(a.inf_norm(), Err(LaError::NonFinite { row: None, col: 0 })); + assert_eq!( + a.inf_norm(), + Err(LaError::NonFinite { + row: Some(1), + col: 2 + }) + ); assert_eq!( a.first_asymmetry(Tolerance::new(0.0).unwrap()).unwrap(), Some((0, 1)) @@ -1969,22 +1897,28 @@ mod tests { #[test] fn first_asymmetry_rejects_scaled_epsilon_overflow() { - let a = Matrix::<2>::from_rows([[2.0, 0.0], [0.0, 1.0]]); + let a = Matrix::<2>::try_from_rows([[0.0, 0.0], [2.0, 1.0]]).unwrap(); let tol = Tolerance::new(f64::MAX).unwrap(); assert_eq!( a.first_asymmetry(tol), - Err(LaError::NonFinite { row: None, col: 0 }) + Err(LaError::NonFinite { + row: Some(1), + col: 0 + }) ); assert_eq!( a.is_symmetric(tol), - Err(LaError::NonFinite { row: None, col: 0 }) + Err(LaError::NonFinite { + row: Some(1), + col: 0 + }) ); } #[test] fn first_asymmetry_flags_overflowed_finite_difference() { - let a = Matrix::<2>::from_rows([[1.0, f64::MAX], [-f64::MAX, 1.0]]); + let a = Matrix::<2>::try_from_rows([[1.0, f64::MAX], [-f64::MAX, 1.0]]).unwrap(); assert_eq!( a.first_asymmetry(Tolerance::new(1e-12).unwrap()).unwrap(), Some((0, 1)) diff --git a/src/tolerance.rs b/src/tolerance.rs new file mode 100644 index 0000000..cd0a48f --- /dev/null +++ b/src/tolerance.rs @@ -0,0 +1,153 @@ +//! Finite tolerance values used by numerical predicates and factorizations. + +use crate::LaError; + +/// Finite, non-negative tolerance used by numerical predicates and factorizations. +/// +/// Construct with [`Tolerance::new`] when accepting raw caller input. Once +/// constructed, the stored value is guaranteed to be finite and `>= 0`, so +/// downstream algorithms do not need to revalidate the tolerance. +/// +/// This is the crate-wide tolerance contract: raw negative, NaN, and infinite +/// values are rejected with [`LaError::InvalidTolerance`] at construction time. +#[must_use] +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Tolerance { + value: f64, +} + +impl Tolerance { + /// Construct a tolerance without checking the raw value. + /// + /// This crate-internal escape hatch is only for constants whose finite, + /// non-negative value is visible at the call site. Public callers should + /// use [`Tolerance::new`] so the returned value carries the validation + /// proof. + pub(crate) const fn new_unchecked(value: f64) -> Self { + Self { value } + } + + /// Construct a finite, non-negative tolerance. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// # fn main() -> Result<(), LaError> { + /// let tol = Tolerance::new(1e-12)?; + /// assert_eq!(tol.get(), 1e-12); + /// # Ok(()) + /// # } + /// ``` + /// + /// # Errors + /// Returns [`LaError::InvalidTolerance`] when `value` is NaN, infinite, or + /// negative. + #[inline] + pub const fn new(value: f64) -> Result { + if value >= 0.0 && value.is_finite() { + Ok(Self::new_unchecked(value)) + } else { + Err(LaError::invalid_tolerance(value)) + } + } + + /// Return the raw finite, non-negative tolerance value. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// # fn main() -> Result<(), LaError> { + /// let tol = Tolerance::new(0.0)?; + /// assert_eq!(tol.get(), 0.0); + /// # Ok(()) + /// # } + /// ``` + #[inline] + #[must_use] + pub const fn get(self) -> f64 { + self.value + } +} + +/// Default absolute threshold used for singularity/degeneracy detection. +/// +/// This is intentionally conservative for geometric predicates and small systems. +/// +/// Conceptually, this is an absolute bound for deciding when a scalar should be treated +/// as "numerically zero" (e.g. LU pivots, LDLT diagonal entries). +/// +/// # Examples +/// ``` +/// use la_stack::prelude::*; +/// +/// # fn main() -> Result<(), LaError> { +/// let lu = Matrix::<2>::identity().lu(DEFAULT_SINGULAR_TOL)?; +/// assert_eq!(lu.det()?, 1.0); +/// # Ok(()) +/// # } +/// ``` +pub const DEFAULT_SINGULAR_TOL: Tolerance = Tolerance::new_unchecked(1e-12); + +/// Relative tolerance used to validate matrices for LDLT factorization. +/// +/// This is crate-internal because LDLT callers provide the factorization +/// tolerance separately; the symmetry tolerance is a fixed domain check used to +/// parse a public [`Matrix`](crate::Matrix) into the internal symmetric proof +/// type before factorization. +pub const LDLT_SYMMETRY_REL_TOL: Tolerance = Tolerance::new_unchecked(1e-12); + +#[cfg(test)] +mod tests { + use core::assert_matches; + + use approx::assert_abs_diff_eq; + + use super::*; + + #[test] + fn default_singular_tol_is_expected() { + assert_abs_diff_eq!(DEFAULT_SINGULAR_TOL.get(), 1e-12, epsilon = 0.0); + } + + #[test] + fn tolerance_new_accepts_finite_non_negative_values() { + assert_eq!( + Tolerance::new(0.0).unwrap().get().to_bits(), + 0.0f64.to_bits() + ); + assert_eq!( + Tolerance::new(1e-12).unwrap().get().to_bits(), + 1e-12f64.to_bits() + ); + assert_eq!( + Tolerance::new(f64::MAX).unwrap().get().to_bits(), + f64::MAX.to_bits() + ); + } + + #[test] + fn tolerance_new_rejects_negative_nan_and_infinity() { + assert_eq!( + Tolerance::new(-1.0), + Err(LaError::InvalidTolerance { value: -1.0 }) + ); + assert_matches!( + Tolerance::new(f64::NAN), + Err(LaError::InvalidTolerance { value }) if value.is_nan() + ); + assert_eq!( + Tolerance::new(f64::INFINITY), + Err(LaError::InvalidTolerance { + value: f64::INFINITY, + }) + ); + assert_eq!( + Tolerance::new(f64::NEG_INFINITY), + Err(LaError::InvalidTolerance { + value: f64::NEG_INFINITY, + }) + ); + } +} diff --git a/src/vector.rs b/src/vector.rs index e9f0e7b..0f5ad65 100644 --- a/src/vector.rs +++ b/src/vector.rs @@ -9,8 +9,9 @@ use crate::LaError; /// Public construction rejects NaN and infinity through [`try_new`](Self::try_new), /// and the storage field is private, so a `Vector` value carries the invariant /// that every stored entry is finite. Algorithms therefore do not re-scan stored -/// entries before using a `Vector`; they only report non-finite values computed -/// during arithmetic, such as overflowed accumulators. +/// entries at every use; user-visible non-finite errors come from construction +/// boundaries or from values computed during arithmetic, such as overflowed +/// accumulators. /// /// Direct field construction is intentionally unavailable to downstream callers: /// @@ -27,97 +28,6 @@ pub struct Vector { data: [f64; D], } -/// Fixed-size vector whose stored entries are all finite. -/// -/// This internal proof wrapper is used where carrying the invariant explicitly -/// makes algorithm boundaries clearer. Public callers use [`Vector`], which is -/// already finite by construction. -#[must_use] -#[derive(Clone, Copy, Debug, PartialEq)] -pub struct FiniteVector { - vector: Vector, -} - -impl FiniteVector { - /// Construct a finite vector without checking the invariant. - /// - /// This crate-internal escape hatch is only for paths with a local proof - /// that every stored entry is finite. - #[inline] - pub(crate) const fn new_unchecked(vector: Vector) -> Self { - Self { vector } - } - - /// Consume the wrapper and return the underlying raw vector. - #[inline] - pub(crate) const fn into_vector(self) -> Vector { - self.vector - } - - /// Borrow the backing array without revalidating stored entries. - #[inline] - #[must_use] - pub(crate) const fn as_array(&self) -> &[f64; D] { - self.vector.as_array() - } - - /// Consume and return the backing array. - #[inline] - #[must_use] - pub(crate) const fn into_array(self) -> [f64; D] { - self.vector.into_array() - } - - /// Dot product for already finite vectors. - /// - /// Stored entries are known finite, so this path only checks whether the - /// accumulation overflows to NaN or infinity. - /// - /// # Errors - /// Returns [`LaError::NonFinite`] when the accumulated dot product overflows - /// to NaN or infinity. - #[inline] - pub(crate) const fn dot(self, other: Self) -> Result { - let lhs = self.as_array(); - let rhs = other.as_array(); - let mut acc = 0.0; - let mut i = 0; - while i < D { - acc = lhs[i].mul_add(rhs[i], acc); - if !acc.is_finite() { - cold_path(); - return Err(LaError::non_finite_at(i)); - } - i += 1; - } - Ok(acc) - } - - /// Squared Euclidean norm for an already finite vector. - /// - /// Stored entries are known finite, so this path only checks whether the - /// accumulation overflows to NaN or infinity. - /// - /// # Errors - /// Returns [`LaError::NonFinite`] when the accumulated norm overflows to NaN - /// or infinity. - #[inline] - pub(crate) const fn norm2_sq(self) -> Result { - let data = self.as_array(); - let mut acc = 0.0; - let mut i = 0; - while i < D { - acc = data[i].mul_add(data[i], acc); - if !acc.is_finite() { - cold_path(); - return Err(LaError::non_finite_at(i)); - } - i += 1; - } - Ok(acc) - } -} - impl Vector { /// Test-only infallible constructor for finite literal fixtures. #[cfg(test)] @@ -170,8 +80,8 @@ impl Vector { /// Return the first non-finite stored entry in index order. /// - /// Shared by the public raw-storage boundary and the internal finite-vector - /// proof wrapper so both paths report the same first offending index with + /// Shared by the public raw-storage boundary and crate-internal reparsing + /// paths so both report the same first offending index with /// [`LaError::NonFinite`]. const fn first_non_finite_entry_in(data: &[f64; D]) -> Option { let mut i = 0; @@ -260,7 +170,19 @@ impl Vector { /// to NaN or infinity. #[inline] pub const fn dot(self, other: Self) -> Result { - FiniteVector::new_unchecked(self).dot(FiniteVector::new_unchecked(other)) + let lhs = self.as_array(); + let rhs = other.as_array(); + let mut acc = 0.0; + let mut i = 0; + while i < D { + acc = lhs[i].mul_add(rhs[i], acc); + if !acc.is_finite() { + cold_path(); + return Err(LaError::non_finite_at(i)); + } + i += 1; + } + Ok(acc) } /// Squared Euclidean norm. @@ -288,7 +210,7 @@ impl Vector { /// or infinity. #[inline] pub const fn norm2_sq(self) -> Result { - FiniteVector::new_unchecked(self).norm2_sq() + self.dot(self) } } @@ -308,12 +230,6 @@ mod tests { use approx::assert_abs_diff_eq; use pastey::paste; - fn assert_array_abs_eq(actual: &[f64; D], expected: &[f64; D]) { - for i in 0..D { - assert_abs_diff_eq!(actual[i], expected[i], epsilon = 0.0); - } - } - macro_rules! gen_public_api_vector_tests { ($d:literal) => { paste! { @@ -458,22 +374,6 @@ mod tests { assert_eq!(a.norm2_sq(), Err(LaError::NonFinite { row: None, col: 0 })); } - #[test] - fn []() { - let arr = { - let mut arr = [0.0f64; $d]; - let values = [1.0f64, 2.0, 3.0, 4.0, 5.0]; - for (dst, src) in arr.iter_mut().zip(values.iter()) { - *dst = *src; - } - arr - }; - - let finite = FiniteVector::<$d>::new_unchecked(Vector::<$d>::new(arr)); - - assert_array_abs_eq(&finite.into_array(), &arr); - assert_array_abs_eq(&finite.into_vector().into_array(), &arr); - } } }; } diff --git a/tests/proptest_exact.rs b/tests/proptest_exact.rs index 9ad5c76..f076661 100644 --- a/tests/proptest_exact.rs +++ b/tests/proptest_exact.rs @@ -3,6 +3,8 @@ //! //! Covers: //! - `det_sign_exact` on diagonal and full small-integer matrices +//! - `det_exact` on full small-integer matrices against an independent +//! `BigRational` Leibniz-expansion oracle //! - `solve_exact` round-trip with integer inputs (`A · x0` in f64 is //! exact for small integers, so `solve(A, A · x0) == x0`) //! - `solve_exact` residual property (`A · solve(A, b) == b` in @@ -68,6 +70,73 @@ fn bigrational_matvec(a: &[[f64; D]; D], x: &[BigRational; D]) - }) } +/// Compute an exact determinant via the Leibniz permutation expansion. +/// +/// This is intentionally independent from the production Bareiss core. It is +/// factorial-time, but the proptests only use D=2..=5, so it stays tiny while +/// giving `det_exact` a separate dense-matrix oracle. +fn bigrational_det_leibniz(a: &[[f64; D]; D]) -> BigRational { + let mut det = BigRational::from_integer(BigInt::from(0)); + let mut perm: [usize; D] = from_fn(|i| i); + + loop { + let mut term = BigRational::from_integer(BigInt::from(1)); + for (row, &col) in perm.iter().enumerate() { + let entry = BigRational::from_f64(a[row][col]).expect("small int fits in BigRational"); + term *= entry; + } + + if permutation_is_even(&perm) { + det += term; + } else { + det -= term; + } + + if !next_permutation(&mut perm) { + break; + } + } + + det +} + +fn permutation_is_even(perm: &[usize]) -> bool { + let mut inversions = 0usize; + for i in 0..perm.len() { + for j in (i + 1)..perm.len() { + if perm[i] > perm[j] { + inversions += 1; + } + } + } + inversions.is_multiple_of(2) +} + +fn next_permutation(values: &mut [usize]) -> bool { + if values.len() < 2 { + return false; + } + + let mut pivot = values.len() - 2; + loop { + if values[pivot] < values[pivot + 1] { + break; + } + if pivot == 0 { + return false; + } + pivot -= 1; + } + + let mut successor = values.len() - 1; + while values[successor] <= values[pivot] { + successor -= 1; + } + values.swap(pivot, successor); + values[(pivot + 1)..].reverse(); + true +} + /// Build a strictly diagonally-dominant f64 matrix from: /// - an off-diagonal matrix of small integers (entries in `[-10, 10]` /// per `small_int_f64`), and @@ -255,6 +324,36 @@ gen_solve_exact_residual_proptests!(3); gen_solve_exact_residual_proptests!(4); gen_solve_exact_residual_proptests!(5); +/// Dense determinant value oracle: random small-integer matrices should match +/// an independent `BigRational` Leibniz expansion, not just the sign read back +/// from the same Bareiss determinant core. +macro_rules! gen_det_exact_leibniz_oracle_proptests { + ($d:literal) => { + paste! { + proptest! { + #![proptest_config(ProptestConfig::with_cases(32))] + + #[test] + fn []( + entries in proptest::array::[]( + proptest::array::[](small_int_f64()), + ), + ) { + let m = Matrix::<$d>::try_from_rows(entries).unwrap(); + let expected = bigrational_det_leibniz::<$d>(&entries); + + prop_assert_eq!(m.det_exact().unwrap(), expected); + } + } + } + }; +} + +gen_det_exact_leibniz_oracle_proptests!(2); +gen_det_exact_leibniz_oracle_proptests!(3); +gen_det_exact_leibniz_oracle_proptests!(4); +gen_det_exact_leibniz_oracle_proptests!(5); + /// On full (non-diagonal) random small-integer matrices, /// `det_sign_exact()` must agree with `det_exact().signum()`. This /// exercises the adaptive fast-filter / Bareiss-fallback boundary on @@ -383,23 +482,40 @@ macro_rules! gen_mixed_scale_diagonal_exact_det_proptests { }; let expected_f64 = expected.to_f64(); - prop_assert_eq!(m.det_exact().unwrap(), expected); prop_assert_eq!(m.det_sign_exact().unwrap(), expected_sign); match expected_f64 { - Some(expected_f64) if expected_f64.is_finite() => { + Some(expected_f64) + if expected_f64.is_finite() + && BigRational::from_f64(expected_f64).as_ref() + == Some(&expected) => + { prop_assert_eq!( m.det_exact_f64().unwrap().to_bits(), expected_f64.to_bits() ); } + Some(expected_f64) if expected_f64.is_finite() => { + prop_assert_eq!( + m.det_exact_f64(), + Err(LaError::Unrepresentable { + index: None, + reason: UnrepresentableReason::RequiresRounding, + }) + ); + } _ => { prop_assert_eq!( m.det_exact_f64(), - Err(LaError::Overflow { index: None }) + Err(LaError::Unrepresentable { + index: None, + reason: UnrepresentableReason::NotFinite, + }) ); } } + + prop_assert_eq!(m.det_exact().unwrap(), expected); } } } diff --git a/tests/regressions.rs b/tests/regressions.rs new file mode 100644 index 0000000..973b0b4 --- /dev/null +++ b/tests/regressions.rs @@ -0,0 +1,108 @@ +//! Regression tests for bugs caught in public API behavior. + +use la_stack::prelude::*; + +#[test] +#[cfg(feature = "exact")] +fn det_exact_f64_preserves_min_positive_subnormal() -> Result<(), LaError> { + let tiny = f64::from_bits(1); + let m = Matrix::<1>::try_from_rows([[tiny]])?; + + assert_eq!(m.det_exact_f64().unwrap().to_bits(), tiny.to_bits()); + Ok(()) +} + +#[test] +#[cfg(feature = "exact")] +fn det_exact_f64_strict_vs_rounded_inexact_det() -> Result<(), LaError> { + let m = Matrix::<2>::try_from_rows([[1.0 + f64::EPSILON, 0.0], [0.0, 1.0 - f64::EPSILON]])?; + + assert_eq!( + m.det_exact_f64(), + Err(LaError::Unrepresentable { + index: None, + reason: UnrepresentableReason::RequiresRounding, + }) + ); + assert_eq!( + m.det_exact_rounded_f64().unwrap().to_bits(), + 1.0f64.to_bits() + ); + Ok(()) +} + +#[test] +#[cfg(feature = "exact")] +fn solve_exact_f64_strict_vs_rounded_non_dyadic() -> Result<(), LaError> { + let a = Matrix::<1>::try_from_rows([[3.0]])?; + let b = Vector::<1>::try_new([1.0])?; + + assert_eq!( + a.solve_exact_f64(b), + Err(LaError::Unrepresentable { + index: Some(0), + reason: UnrepresentableReason::RequiresRounding, + }) + ); + assert_eq!( + a.solve_exact_rounded_f64(b).unwrap().into_array()[0].to_bits(), + (1.0f64 / 3.0).to_bits() + ); + Ok(()) +} + +#[test] +#[cfg(feature = "exact")] +fn requires_rounding_error_can_fall_back_to_rounded_solve() -> Result<(), LaError> { + let a = Matrix::<1>::try_from_rows([[3.0]])?; + let b = Vector::<1>::try_new([1.0])?; + + let rounded = match a.solve_exact_f64(b) { + Ok(x) => x, + Err(err) if err.requires_rounding() => a.solve_exact_rounded_f64(b)?, + Err(err) => return Err(err), + }; + + assert_eq!(rounded.into_array()[0].to_bits(), (1.0f64 / 3.0).to_bits()); + Ok(()) +} + +#[test] +fn det_direct_skips_zero_coefficient_terms_that_would_overflow() -> Result<(), LaError> { + let d3 = Matrix::<3>::try_from_rows([ + [1.0, 0.0, 0.0], + [1.0e300, 1.0, 1.0e300], + [1.0e300, 0.0, 1.0e300], + ])?; + assert_eq!(d3.det_direct(), Ok(Some(1.0e300))); + + let d4 = Matrix::<4>::try_from_rows([ + [0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0], + [1.0e300, 0.0, 1.0e300, 1.0e300], + [1.0e300, 0.0, 1.0e300, -1.0e300], + ])?; + assert_eq!(d4.det_direct(), Ok(Some(0.0))); + + Ok(()) +} + +#[test] +fn det_errbound_skips_zero_coefficient_terms_that_would_overflow() -> Result<(), LaError> { + let d3 = Matrix::<3>::try_from_rows([ + [1.0, 0.0, 0.0], + [1.0e300, 1.0, 1.0e300], + [1.0e300, 0.0, 1.0e300], + ])?; + assert_eq!(d3.det_errbound(), Ok(Some(ERR_COEFF_3 * 1.0e300))); + + let d4 = Matrix::<4>::try_from_rows([ + [0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0], + [1.0e300, 0.0, 1.0e300, 1.0e300], + [1.0e300, 0.0, 1.0e300, -1.0e300], + ])?; + assert_eq!(d4.det_errbound(), Ok(Some(0.0))); + + Ok(()) +} diff --git a/tests/semgrep/src/project_rules/finite_api_contract.rs b/tests/semgrep/src/project_rules/finite_api_contract.rs index fa8cf04..acb856f 100644 --- a/tests/semgrep/src/project_rules/finite_api_contract.rs +++ b/tests/semgrep/src/project_rules/finite_api_contract.rs @@ -10,23 +10,6 @@ pub mod vector; // ok: la-stack.rust.no-public-raw-linear-algebra-modules mod exact; -pub mod reexports { - // ruleid: la-stack.rust.no-public-finite-proof-reexports - pub use crate::matrix::FiniteMatrix; - - // ruleid: la-stack.rust.no-public-finite-proof-reexports - pub use crate::vector::{FiniteVector, Vector}; - - // ruleid: la-stack.rust.no-public-finite-proof-reexports - pub use crate::matrix::FiniteMatrix as ParsedMatrix; - - // ok: la-stack.rust.no-public-finite-proof-reexports - pub(crate) use crate::vector::FiniteVector as InternalFiniteVector; - - // ok: la-stack.rust.no-public-finite-proof-reexports - pub use crate::matrix::Matrix; -} - #[must_use] pub struct Matrix { // ruleid: la-stack.rust.no-public-matrix-vector-storage-fields diff --git a/tests/semgrep/src/project_rules/raw_f64_constructors.rs b/tests/semgrep/src/project_rules/raw_f64_constructors.rs index da2c7ee..1710cd0 100644 --- a/tests/semgrep/src/project_rules/raw_f64_constructors.rs +++ b/tests/semgrep/src/project_rules/raw_f64_constructors.rs @@ -17,18 +17,13 @@ pub enum LaError { } impl Matrix { - // ruleid: la-stack.rust.no-public-infallible-raw-f64-constructors - pub const fn from_rows(rows: [[f64; D]; D]) -> Self { - Self { rows } - } - // ok: la-stack.rust.no-public-infallible-raw-f64-constructors pub const fn try_from_rows(rows: [[f64; D]; D]) -> Result { Ok(Self { rows }) } // ok: la-stack.rust.no-public-infallible-raw-f64-constructors - pub(crate) const fn from_rows_literal(rows: [[f64; D]; D]) -> Self { + pub(crate) const fn matrix_literal(rows: [[f64; D]; D]) -> Self { Self { rows } } } diff --git a/tests/semgrep/src/project_rules/readme_doctest_mirrors.rs b/tests/semgrep/src/project_rules/readme_doctest_mirrors.rs index 2a02e04..7a5314c 100644 --- a/tests/semgrep/src/project_rules/readme_doctest_mirrors.rs +++ b/tests/semgrep/src/project_rules/readme_doctest_mirrors.rs @@ -2,7 +2,7 @@ #![allow(dead_code)] // ruleid: la-stack.rust.no-unwrap-expect-in-readme-doctest-mirrors -mod readme_doctests { +mod readme_doctests_unwrap { #[test] fn readme_mirror_uses_unwrap() { let _ = Some(1_u32).unwrap(); @@ -10,7 +10,7 @@ mod readme_doctests { } // ruleid: la-stack.rust.no-unwrap-expect-in-readme-doctest-mirrors -mod readme_doctests { +mod readme_doctests_expect { #[test] fn readme_mirror_uses_expect() { let _ = Ok::(1).expect("README mirrors should use ?"); @@ -26,7 +26,7 @@ mod tests { } // ok: la-stack.rust.no-unwrap-expect-in-readme-doctest-mirrors -mod readme_doctests { +mod readme_doctests_result { #[test] fn readme_mirror_uses_result() -> Result<(), &'static str> { let _ = Ok::(1)?;