diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 48e8f85..02e77e4 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -27,9 +27,9 @@ env: CARGO_NEXTEST_VERSION: "0.9.137" DPRINT_VERSION: "0.54.0" JUST_VERSION: "1.51.0" - RUMDL_VERSION: "0.2.6" + RUMDL_VERSION: "0.2.9" TAPLO_VERSION: "0.10.0" - TYPOS_VERSION: "1.47.1" + TYPOS_VERSION: "1.47.2" UV_VERSION: "0.11.19" ZIZMOR_VERSION: "1.25.2" diff --git a/AGENTS.md b/AGENTS.md index 7b6ba22..ed213ac 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -36,11 +36,12 @@ invariant over the convenient edit. - Error enums are `#[non_exhaustive]`; public wrapper types are `#[must_use]`. -- New functionality is additive: use the prelude for ergonomic re-exports; - never silently rename or remove a public item. -- Pre-1.0 semver: `0.x.Y` is a patch-level additive bump, `0.X.y` is a - minor bump that may include breaking changes. Conventional-commit - types (`feat`, `fix`, `refactor`, …) mirror this convention. +- New functionality is additive by default: use the prelude for ergonomic + re-exports, and avoid churn for its own sake. +- Pre-1.0 semver: any `0.x.y` release may include breaking API changes when + they materially improve correctness, orthogonality, or long-term API clarity. + Do not keep compatibility aliases that weaken the public model; document + intentional breaks clearly in release notes and commit messages. - Do **not** automatically update the library version in `Cargo.toml`, `Cargo.lock`, README dependency snippets, or related docs during ordinary feature, fix, review, or hygiene work. Version bumps are maintainer-driven @@ -68,7 +69,7 @@ invariant over the convenient edit. `Option` instead of relying on `panic!`, `assert!`, `unwrap`, or `expect`. - Borrow by default (`&T`, `&[T]`); return borrowed views when possible. - Type and function names match textbook vocabulary (`Matrix`, `Vector`, - `Lu`, `Ldlt`, `solve_vec`, `det`, `inf_norm`). Avoid Rust-ecosystem + `Lu`, `Ldlt`, `solve`, `det`, `inf_norm`). Avoid Rust-ecosystem abstractions that obscure the math. ### Scientific notation in docs @@ -210,7 +211,7 @@ testable. #### Reference examples - `src/matrix.rs` — `gen_public_api_matrix_tests!` -- `src/lu.rs` — `gen_public_api_pivoting_solve_vec_and_det_tests!`, `gen_public_api_tridiagonal_smoke_solve_vec_and_det_tests!` +- `src/lu.rs` — `gen_public_api_pivoting_solve_and_det_tests!`, `gen_public_api_tridiagonal_smoke_solve_and_det_tests!` - `src/ldlt.rs` — `gen_public_api_ldlt_identity_tests!`, `gen_public_api_ldlt_diagonal_tests!` - `src/exact.rs` — `gen_det_exact_tests!`, `gen_det_exact_f64_tests!`, `gen_solve_exact_tests!`, `gen_solve_exact_f64_tests!` @@ -354,11 +355,11 @@ When creating or updating issues: - This is a single Rust *library crate* (no `src/main.rs`). The crate root is `src/lib.rs`. - The linear algebra implementation is split across: - - `src/lib.rs`: crate root + shared items (`LaError`, `DEFAULT_SINGULAR_TOL`, `DEFAULT_PIVOT_TOL`) + re-exports + - `src/lib.rs`: crate root + shared items (`LaError`, `DEFAULT_SINGULAR_TOL`) + re-exports - `src/vector.rs`: `Vector` (`[f64; D]`) - `src/matrix.rs`: `Matrix` (`[[f64; D]; D]`) + helpers (`get`, `set`, `inf_norm`, `det`, `det_direct`) - - `src/lu.rs`: `Lu` factorization with partial pivoting (`solve_vec`, `det`) - - `src/ldlt.rs`: `Ldlt` factorization without pivoting for symmetric SPD/PSD matrices (`solve_vec`, `det`) + - `src/lu.rs`: `Lu` factorization with partial pivoting (`solve`, `det`) + - `src/ldlt.rs`: `Ldlt` factorization without pivoting for symmetric SPD/PSD matrices (`solve`, `det`) - `src/exact.rs`: exact arithmetic behind `features = ["exact"]`: - Determinants: `det_exact()`, `det_exact_f64()`, `det_sign_exact()` via integer-only Bareiss in `BigInt` (`bareiss_det_int`); `det_sign_exact()` adds a Shewchuk-style diff --git a/CHANGELOG.md b/CHANGELOG.md index 12aa440..34d851e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,11 +13,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Add a generic docs-version sync check that compares Markdown dependency snippets against the Cargo package name and version - - Run the docs-version check from the repository Semgrep policy lane - Refresh README determinant examples with explicit fallible handling and hidden doctest mirrors - - Update CI uv pins to 0.11.19 ### Documentation @@ -52,22 +50,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 Remove the broad restore-keys pattern to ensure only exact version matches are used from the cache, preventing potential version mismatches during CI runs. - - Revert "ci: modernize tooling checks and example execution" [`a4b9f64`](https://github.com/acgetchell/la-stack/commit/a4b9f64235f53274270ef9272634a61f653dad87) - - Reapply "ci: modernize tooling checks and example execution" [`758321a`](https://github.com/acgetchell/la-stack/commit/758321acf872b1f17286ff3bb7bee6a807e4b440) - - Encode nonzero mantissas in exact decomposition [`7a664ed`](https://github.com/acgetchell/la-stack/commit/7a664ede2f4add168c5813f8d24e16732fa03b30) - Replace the exact-arithmetic zero mantissa sentinel with `Option<NonZeroU64>`. - Carry nonzero mantissa proof through matrix/vector decomposition and BigInt scaling. - Clarify determinant documentation around uncertified `det()` bounds. - Keep SPD determinant proptests on the tolerance-aware LU path. - - Simplify finite proof wrappers [`54b603c`](https://github.com/acgetchell/la-stack/commit/54b603c21f0eb5b4d63ff334fac7f8cc325ebc2e) - - Use the checked proof-wrapper constructors as the single internal path for finite matrices and vectors. - Remove exact-arithmetic tests that duplicated the matrix and vector non-finite boundary checks. @@ -93,12 +86,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Add a roadmap covering the v0.4.x stable-Rust issue sequence and the v0.5.0 generic_const_exprs anchor. - Refresh generated changelog entries and archived changelog grouping. - Document finite RHS solve validation [`075aed7`](https://github.com/acgetchell/la-stack/commit/075aed78cf8264fc920258f1f1d977ddd589ffd7) - - Document that LU and LDLT solve_vec reject non-finite RHS entries with LaError::NonFinite metadata. - Cite the Bareiss reference in the exact solve helper docs and describe exact-arithmetic growth and complexity. - Cover finite proof defaults and non-finite RHS solve boundaries in unit tests. - Clarify finite solve and norm guarantees [`fb71485`](https://github.com/acgetchell/la-stack/commit/fb71485cac0b464b2fa9ee949140a558b7738781) - - State that LU and LDLT solve_vec use floating-point substitution without a certified absolute rounding-error bound. - Clarify that inf_norm reports NonFinite for unchecked stored NaN/∞ as well as row-sum overflow. - Exercise the unchecked finite-proof fixture path directly in exact tests. @@ -110,28 +101,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Enforce the tolerance contract around symmetry checks by surfacing scaled tolerance overflow as a typed non-finite intermediate error. - - Document finite, non-negative tolerance requirements across tolerance-taking matrix APIs. - - Add regression coverage for invalid tolerance construction and symmetry tolerance overflow. - - Update exact examples to propagate typed crate errors instead of unwrapping. - Harden Semgrep fixture parsing [`ac44c07`](https://github.com/acgetchell/la-stack/commit/ac44c078cc4435d5beca27f1890fbb4046cf5952) - - Ignore non-canonical todoruleid annotations when counting expected rule hits. - Reject malformed Semgrep JSON results with clear stderr diagnostics instead of propagating KeyError. - Revalidate finite proof conversions [`419a90f`](https://github.com/acgetchell/la-stack/commit/419a90f7267608051736498154ac5e6faf0909c5) Ensure internal finite proof conversions cannot accept raw Matrix or Vector storage without checking the invariant. - - Revalidate TryFrom<Matrix<D>> and TryFrom<Vector<D>> before constructing finite wrappers. - Measure exact random percentile benchmarks over repeated corpus timings and cumulative input sets. - Tighten Codecov status thresholds and extend benchmark workflow timeout. - Keep Semgrep constructor fixtures aligned with public API guardrails. - Revalidate public compute inputs [`ffca00e`](https://github.com/acgetchell/la-stack/commit/ffca00e9dd6fde5e57c8064f69807a76e45a469e) - - Parse Matrix and Vector storage into private finite proof-bearing types at public compute boundaries. - Reject unchecked non-finite storage before LU, LDLT, determinant, norm, dot, and exact-arithmetic paths can proceed. - Keep unchecked proof-wrapper constructors crate-private for internal paths with local finiteness proofs. @@ -148,7 +133,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Add archive-aware changelog generation, post-processing, and tag-release tooling. - Preserve exact arithmetic overflow reporting without non-finite sentinel defaults. - Modernize tooling checks and example execution [`da626bc`](https://github.com/acgetchell/la-stack/commit/da626bcca899aa91d58f728db433a53a46179e92) - - Run examples from prebuilt binaries instead of invoking cargo run for each example. - Add check/fix recipe aliases and guard documented command ordering with Semgrep. - Document the Rust-native Markdown, YAML, TOML, spelling, and workflow action policy. @@ -181,7 +165,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 calculations, correctly detect and return `LaError::NonFinite` when intermediate calculations overflow to infinity despite having finite inputs. - - Defensive-path test coverage for LU and LDLT solve_vec [`87d426f`](https://github.com/acgetchell/la-stack/commit/87d426fca1627445b804fd26b62fc7d9d4f0ae48) Add unit tests to exercise internal safety nets in the LU and LDLT @@ -189,7 +172,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 invalid diagonals (NaN or sub-tolerance) to verify that solve_vec correctly surfaces NonFinite and Singular errors, even though these states are unreachable via the standard factorization APIs. - - Const-ify Lu/Ldlt det + solve_vec and Matrix inf_norm + det_errbound [`81ecb35`](https://github.com/acgetchell/la-stack/commit/81ecb35bdaf159f1f44d1eb24274ecf82c6567d5) @@ -207,11 +189,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 rewritten as `while` loops since they are not const-stable. Fix error-variant correctness in both solve_vec paths: - - A corrupt stored `U` / `D` diagonal at `(i, i)` now surfaces as `LaError::NonFinite { row: Some(i), col: i }`, matching the convention used by `Matrix::det`, `Lu::factor`, and `Ldlt::factor`. - - A computed-intermediate overflow keeps `row: None, col: i`. - Previously both were conflated into `row: None, col: i`, defeating debuggability for callers who construct factorizations directly. @@ -219,7 +199,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 Sharpen `LaError::NonFinite` variant docs and the `# Errors` sections of `Lu::solve_vec` / `Ldlt::solve_vec` to spell out the `(row, col)` contract for each failure mode. - - Fast-filter boundary proptests for exact determinant sign [`6357db3`](https://github.com/acgetchell/la-stack/commit/6357db35c70bca1b93e6bbf9a4fd231913631950) Introduce proptests to verify that the floating-point determinant sign @@ -240,14 +219,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 determinants, adding a dedicated section for AI agent governance, and refining internal error helpers to use a more consistent coordinate convention for non-finite values. - - Consolidate and expand const-evaluability tests via macros [`f8d80a0`](https://github.com/acgetchell/la-stack/commit/f8d80a01e9913f87e1b19b2ad5ffbc0994e2bfdb) Refactor manual const-evaluation tests in `Ldlt` and `Matrix` into macros to standardize coverage across matrix dimensions 2 through 5. This ensures all determinant and norm variants are fully exercised at compile time. - - Refactor solve_exact to use hybrid Bareiss forward elimination [`ecbbe8a`](https://github.com/acgetchell/la-stack/commit/ecbbe8a571ccaeb9cfedbf0269b8db44d43a5773) @@ -257,7 +234,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 This eliminates per-entry GCD overhead during elimination. Internal f64 decomposition is unified via f64_decompose, and f64_to_bigrational is moved to test scope. - - Polish exact module (Component struct, errors, perf) [`53a5be6`](https://github.com/acgetchell/la-stack/commit/53a5be6abecc0af332398236ed6803ed75564b03) Major refactor of `src/exact.rs`: @@ -266,28 +242,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 `component_to_bigint`, `build_bigint_matrix`, `build_bigint_vec`, `bareiss_forward_eliminate`) so `bareiss_det_int` and `gauss_solve` share a single integer-Bareiss core. - - Hybrid BigInt/BigRational solve: forward elimination runs entirely in `BigInt` with fraction-free Bareiss updates on `(A | b)`; only the `O(D²)` back-substitution phase lifts into `BigRational`. - - `mem::take` in back-substitution eliminates per-entry `clone()` calls on `rhs[i]`, `a[i][j]`, `a[i][i]`. Structured errors replace preconditions: - - `decompose_matrix` / `decompose_vec` fold `is_finite()` into the same pass that decomposes each entry and return `Err(LaError::NonFinite { row, col })` on the first non-finite cell. - - `bareiss_det_int`, `bareiss_det`, `gauss_solve` now return `Result<_, LaError>`; error propagates to every public entry point via `?`. - - `validate_finite` and `validate_finite_vec` removed (dead after the refactor); `det_sign_exact` relies on IEEE 754 NaN/∞ propagation through `det_direct()` plus `bareiss_det_int` for Stage-2 validation. - - Add adversarial-input coverage for exact arithmetic [#80](https://github.com/acgetchell/la-stack/pull/80) [`5bf5815`](https://github.com/acgetchell/la-stack/commit/5bf5815cb165c3b6145c5592420a58085a66efaa) @@ -295,56 +265,43 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 APIs to catch tail cases that fixed well-conditioned inputs miss. Benchmarks (benches/exact.rs): - - Factor out `bench_extreme_group` helper running the same four benches (`det_sign_exact`, `det_exact`, `solve_exact`, `solve_exact_f64`) so adversarial groups are directly comparable. - - Extend `exact_near_singular_3x3` with the two solve benches (the primary motivating use case for exact solve was previously unmeasured). - - Add `exact_large_entries_3x3` (entries near `f64::MAX / 2`) to stress intermediate BigInt growth during Bareiss forward elimination. - - Add `exact_hilbert_4x4` / `exact_hilbert_5x5` to stress the `f64_decompose → BigInt` scaling path on ill-conditioned inputs. - - Tighten bench `expect(...)` messages to name the invariant each call relies on (e.g. "non-singular matrix with finite entries") so panics identify both where (Criterion bench name) and why. Unit tests (src/exact.rs): - - `solve_exact_near_singular_3x3_integer_x0` — integer-x0 round-trip through the 2^-50-perturbed matrix. - - `solve_exact_large_entries_3x3_unit_vector` — `A·[1,0,0] = [big,1,1]` round-trip with `f64::MAX/2` diagonal. - - `det_sign_exact_large_entries_3x3_positive` — asserts the fast filter falls through (`det_direct` is non-finite) and `det_exact_f64` returns `Overflow { index: None }`. - - `det_sign_exact_hilbert_positive_{3,4,5}d` — Hilbert is SPD, sign = 1. - `solve_exact_hilbert_residual_{3,4,5}d` — residual `A·x - b` is exactly zero in `BigRational`, stronger than integer round-trips since Hilbert entries are non-terminating in binary. Proptests (tests/proptest_exact.rs): - - `solve_exact_integer_roundtrip_{2..5}d` — random diagonally-dominant integer `A` and small-integer `x0`, verify `solve(A, A·x0) == x0` exactly. - - `solve_exact_residual_{2..5}d` — random `A` + small-integer `b`, verify `A · solve(A, b) == b` exactly (catches back-sub bugs on fractional solutions). - - `det_sign_exact_agrees_with_det_exact_{2..5}d` — on full (non-diagonal) small-integer matrices, asserts `det_sign_exact() == det_exact().sign()` (exercises the filter/fallback boundary previously only diagonal-tested). Prelude (src/lib.rs): - - Re-export `BigInt` alongside `BigRational` (crate root + prelude). - Re-export `FromPrimitive`, `Signed`, `ToPrimitive` from `num-traits` in the prelude so the re-exported `BigRational` is actually usable for @@ -353,23 +310,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 their own Cargo.toml. Additive; no public API breaks. Tooling (scripts/bench_compare.py + test_bench_compare.py): - - Register the three new adversarial groups in `EXACT_GROUPS`. - Extend `_group_heading` with human-readable titles ("Large entries 3x3", "Hilbert 4x4", "Hilbert 5x5"). - - Add group-heading unit tests for the new cases. Test results (`just ci`): - - cargo test --features exact: 368 lib + 20 proptest_exact + 40 other proptests + 34 doc-tests — all pass - - cargo test (no features): 175 lib + 40 proptests + 29 doc-tests — all pass - Python: 104 tests pass (ty, mypy, ruff clean) - Clippy (pedantic + nursery + cargo, `-D warnings`): clean - fmt, taplo, yamllint, shellcheck, spell-check, bench-compile, examples: clean - - Expand exact-arithmetic re-exports and adversarial benchmarks [`b1a491d`](https://github.com/acgetchell/la-stack/commit/b1a491d902ccdaba6f9cd2e6f8e05514b6dfa3de) @@ -382,7 +334,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 large-entry, and Hilbert matrices to better track performance on ill-conditioned data. CI benchmark comparisons are now more resilient to newly added benchmarks via `--baseline-lenient`. - - Update AGENTS.md [`1e0648d`](https://github.com/acgetchell/la-stack/commit/1e0648dad3147ef127c775d8969b7cd214a2a6ed) ### Dependencies @@ -412,7 +363,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Propagate NaN in Matrix::inf_norm [#85](https://github.com/acgetchell/la-stack/pull/85) [`16ffa45`](https://github.com/acgetchell/la-stack/commit/16ffa45ade11b21a179cad8fcecc51d802086a1d) - - Report infinite vs finite off-diagonal pairs as asymmetric [`1805779`](https://github.com/acgetchell/la-stack/commit/1805779dbca49183fbfa95c68ec00984966aa551) Update Matrix::first_asymmetry to flag any non-finite difference between @@ -428,28 +378,24 @@ 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 - - src/lu.rs: Lu::factor and Lu::solve_vec NonFinite / Singular returns - src/ldlt.rs: Ldlt::factor and Ldlt::solve_vec NonFinite / Singular returns - - src/matrix.rs: det_direct D >= 5 fallback arm (legal because cold_path is const fn in 1.95) and det NonFinite / overflow scan - - Refactor det_sign_exact Stage 1 fast filter to use match + if let guard with let-chain, replacing the tuple destructure; semantics unchanged Test results (local `just ci`): - - cargo fmt --all -- --check: clean - cargo clippy --workspace --all-targets --all-features -D warnings -W clippy::pedantic -W clippy::nursery -W clippy::cargo: clean - - cargo test --features exact --lib: 258 passed, 0 failed - cargo test --features exact (doctests + examples): 31 passed, 0 failed - RUSTDOCFLAGS='-D warnings' cargo doc --no-deps --features exact: clean @@ -461,17 +407,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 [`a1d3bdb`](https://github.com/acgetchell/la-stack/commit/a1d3bdbdb6fcb778a78a2a3d0cb66b79484e1472) Replace BigRational-only gauss_solve with a hybrid that runs - Bareiss fraction-free forward elimination in BigInt on the - augmented (A | b) system, then back-substitutes in BigRational. - Eliminates GCD normalization from the O(D^3) phase while keeping - rational overhead limited to the cheaper O(D^2) back-sub. + Bareiss fraction-free forward elimination in BigInt on the + augmented (A | b) system, then back-substitutes in BigRational. + Eliminates GCD normalization from the O(D^3) phase while keeping + rational overhead limited to the cheaper O(D^2) back-sub. - Scope f64_to_bigrational to cfg(test); production code paths now - use f64_decompose directly (shared with bareiss_det_int). + Scope f64_to_bigrational to cfg(test); production code paths now + use f64_decompose directly (shared with bareiss_det_int). - Closes #72 + Closes #72 - Co-Authored-By: Oz <oz-agent@warp.dev> + Co-Authored-By: Oz ## [0.4.0] - 2026-04-11 @@ -497,23 +443,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 Update CITATION.cff to reflect the v0.3.0 release. Keywords are updated to include exact arithmetic and robust predicates, with existing terms standardized to kebab-case. - - Refine benchmark comparison reporting and documentation [`e1b5955`](https://github.com/acgetchell/la-stack/commit/e1b5955fb5024232e34e9df9701bb24fb98efa15) Refactor the `bench_compare.py` script to group results by dimension and add a new test suite. Consolidate all benchmarking workflows into a new dedicated guide and move machine-specific performance snapshots to `.gitignore`. - - Expand test coverage for benchmark comparison edge cases [`bced7d9`](https://github.com/acgetchell/la-stack/commit/bced7d988bd2f42cb6bb5af9c54dabdcf787a5fc) - - Update documentation and tests for integer-only Bareiss [`2ee3f05`](https://github.com/acgetchell/la-stack/commit/2ee3f05caecfdf1a23b61257a7465b3bb6d63614) Update architecture descriptions in AGENTS.md and REFERENCES.md to reflect the move from BigRational to integer-only BigInt arithmetic. Add unit tests for bareiss_det_int covering negative inputs and pivot swapping. Refine gh CLI patterns for automated issue listing. - - Restrict benchmark baselines to main and improve reporting [`9a7caa2`](https://github.com/acgetchell/la-stack/commit/9a7caa241b5f476c2659772f3189468b10fcba2e) Update internal CI logic to ensure Criterion baselines are only saved @@ -554,14 +496,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - On push to main: run exact benchmarks, save Criterion baseline as artifact (30-day retention) - - On PRs: find latest main baseline artifact, download, run benchmarks with --baseline comparison, report regressions in job summary - - Regressions are warning-only — the workflow never blocks PRs - Uses artifact-based baselines (not cache) for robustness: explicit provenance, no silent eviction, fork-safe - - Validated locally: --save-baseline and --baseline flags confirmed working with Criterion @@ -578,23 +517,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - bareiss_det: Vec<Vec<BigRational>> → [[BigRational; D]; D] via std::array::from_fn - - gauss_solve: Vec<Vec<BigRational>> augmented matrix → separate [[BigRational; D]; D] + [BigRational; D] (avoids unstable const-generic [T; D+1]) - - gauss_solve return type: Vec<BigRational> → [BigRational; D] - solve_exact: drop intermediate clone, forward array directly Error enrichment (src/lib.rs, src/lu.rs, src/ldlt.rs, src/matrix.rs, src/exact.rs): - - LaError::NonFinite gains row: Option<usize> for full (row, col) location on matrix inputs; None for vectors/intermediates - - LaError::Overflow becomes a struct variant with index: Option<usize> to identify which solution component overflowed - - Display messages updated to show positional context - Custom f64 → BigRational via IEEE 754 bit decomposition [#63](https://github.com/acgetchell/la-stack/pull/63) [`0a8ce5b`](https://github.com/acgetchell/la-stack/commit/0a8ce5b45e21fa88e4e9832bc6ada38b3ba68eeb) @@ -602,15 +536,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 Replace `BigRational::from_float(x)` in `f64_to_bigrational` with manual IEEE 754 binary64 bit decomposition and `BigRational::new_raw`, bypassing the unnecessary GCD normalization that `from_float` performs internally. - - Decompose f64 into sign, biased exponent, and significand fields - Strip trailing zeros from the significand so the fraction is already in lowest terms (odd numerator over power-of-two denominator) - - Handle zero (±0.0), subnormals, and normal values; panic on NaN/Inf - Add 15 dedicated unit tests: ±zero, integers, fractions, powers of two, subnormals, round-trip fidelity, lowest-terms, and panic cases - - Add IEEE 754-2019 [9] and Goldberg [10] references to REFERENCES.md - Add module-level and function-level doc citations - Integer-only Bareiss determinant via BigInt [#64](https://github.com/acgetchell/la-stack/pull/64) @@ -620,16 +551,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 all f64 entries are decomposed into mantissa × 2^exponent, scaled to a common integer base, and eliminated without any rational arithmetic. The result is reconstructed as BigRational only at the end. - - Add f64_decompose helper (extracted from f64_to_bigrational) - Add bareiss_det_int: integer-only Bareiss returning (BigInt, i32) - Add bigint_exp_to_bigrational: reconstruction with trailing-zero reduction (no full GCD needed) - - Refactor bareiss_det as thin wrapper over bareiss_det_int - Optimize det_sign_exact to read sign from BigInt directly (skip BigRational reconstruction entirely) - - Import std::array::from_fn for shorter call sites - Replace clippy cast suppression with exact try_from conversions - Update AGENTS.md with non-interactive gh CLI patterns @@ -638,10 +566,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 negative exp, reduction, negative values) Performance (vs pre-bigint baseline on Apple M4 Max): - det_exact: 16x (D=2) → 39x (D=5) faster - det_exact_f64: 10x (D=2) → 38x (D=5) faster - det_sign_exact: 40x at D=5 (bypasses BigRational entirely) - Near-singular: 18x faster + det_exact: 16x (D=2) → 39x (D=5) faster + det_exact_f64: 10x (D=2) → 38x (D=5) faster + det_sign_exact: 40x at D=5 (bypasses BigRational entirely) + Near-singular: 18x faster ## Archives diff --git a/README.md b/README.md index a233e02..ed3931b 100644 --- a/README.md +++ b/README.md @@ -105,8 +105,8 @@ fn main() -> Result<(), LaError> { let b = Vector::<5>::try_new([14.0, 13.0, 12.0, 11.0, 10.0])?; - let lu = a.lu(DEFAULT_PIVOT_TOL)?; - let x = lu.solve_vec(b)?.into_array(); + let lu = a.lu(DEFAULT_SINGULAR_TOL)?; + let x = lu.solve(b)?.into_array(); // Floating-point rounding is expected; compare with a tolerance. let expected = [1.0, 2.0, 3.0, 4.0, 5.0]; @@ -330,18 +330,18 @@ compose the same bound themselves. | Type | Storage | Purpose | Key methods | |---|---|---|---| -| `Vector` | `[f64; D]` | Fixed-length vector for input and computation | `try_new`, `zero`, `dot`, `norm2_sq` | -| `Matrix` | `[[f64; D]; D]` | Square matrix for input and computation | See below | -| `Lu` | `Matrix` + pivot array | Factorization for solves/det | `solve_vec`, `det` | -| `Ldlt` | `Matrix` | Factorization for symmetric SPD/PSD solves/det | `solve_vec`, `det` | +| `Vector` | `[f64; D]` | Finite fixed-length vector for input and computation | `try_new`, `zero`, `dot`, `norm2_sq` | +| `Matrix` | `[[f64; D]; D]` | Finite square matrix for input and computation | See below | +| `Lu` | `Matrix` + pivot array | Factorization for solves/det | `solve`, `det` | +| `Ldlt` | `Matrix` | Factorization for symmetric SPD/PSD solves/det | `solve`, `det` | 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`¹. -Matrix and vector methods validate non-finite inputs at public API boundaries and -carry private proof-bearing types through computation so successful factors do -not store NaN or infinity. +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. ¹ Requires `features = ["exact"]`. diff --git a/benches/vs_linalg.rs b/benches/vs_linalg.rs index 3f816e0..8727c90 100644 --- a/benches/vs_linalg.rs +++ b/benches/vs_linalg.rs @@ -16,7 +16,7 @@ use faer::{Mat, Side}; use nalgebra::{SMatrix, SVector}; use pastey::paste; -use la_stack::{DEFAULT_PIVOT_TOL, DEFAULT_SINGULAR_TOL, Matrix, Vector}; +use la_stack::{DEFAULT_SINGULAR_TOL, Matrix, Vector}; mod common { pub mod vs_linalg; @@ -63,7 +63,6 @@ macro_rules! define_vs_linalg_benches_for_dim { Vector::<$d>::try_new(make_vector_array::<$d>(1.0)), "la_stack vector construction", ); - let na = SMatrix::::from_fn(|r, c| matrix_entry::<$d>(r, c)); let nrhs = SVector::::from_fn(|i, _| vector_entry(i, 0.0)); let nv1 = SVector::::from_fn(|i, _| vector_entry(i, 0.0)); @@ -75,7 +74,7 @@ macro_rules! define_vs_linalg_benches_for_dim { let fv2 = Mat::::from_fn($d, 1, |i, _| vector_entry(i, 1.0)); // Precompute LU once for solve-only / det-only benchmarks. - let a_lu = require_ok(a.lu(DEFAULT_PIVOT_TOL), "precomputed la_stack LU"); + let a_lu = require_ok(a.lu(DEFAULT_SINGULAR_TOL), "precomputed la_stack LU"); let a_ldlt = require_ok(a.ldlt(DEFAULT_SINGULAR_TOL), "precomputed la_stack LDLT"); let na_lu = na.clone().lu(); let na_cholesky = require_some(na.clone().cholesky(), "precomputed nalgebra Cholesky"); @@ -88,7 +87,7 @@ macro_rules! define_vs_linalg_benches_for_dim { [].bench_function("la_stack_det_via_lu", |bencher| { bencher.iter(|| { let lu = require_ok( - black_box(a).lu(DEFAULT_PIVOT_TOL), + black_box(a).lu(DEFAULT_SINGULAR_TOL), "la_stack LU factorization", ); let det = require_ok(lu.det(), "la_stack LU determinant"); @@ -124,7 +123,7 @@ macro_rules! define_vs_linalg_benches_for_dim { [].bench_function("la_stack_lu", |bencher| { bencher.iter(|| { let lu = require_ok( - black_box(a).lu(DEFAULT_PIVOT_TOL), + black_box(a).lu(DEFAULT_SINGULAR_TOL), "la_stack LU factorization", ); let _ = black_box(lu); @@ -180,11 +179,11 @@ macro_rules! define_vs_linalg_benches_for_dim { [].bench_function("la_stack_lu_solve", |bencher| { bencher.iter(|| { let lu = require_ok( - black_box(a).lu(DEFAULT_PIVOT_TOL), + black_box(a).lu(DEFAULT_SINGULAR_TOL), "la_stack LU factorization", ); let x = require_ok( - lu.solve_vec(black_box(rhs)), + lu.solve(black_box(rhs)), "la_stack LU solve", ); let _ = black_box(x); @@ -215,7 +214,7 @@ macro_rules! define_vs_linalg_benches_for_dim { "la_stack LDLT factorization", ); let x = require_ok( - ldlt.solve_vec(black_box(rhs)), + ldlt.solve(black_box(rhs)), "la_stack LDLT solve", ); let _ = black_box(x); @@ -248,7 +247,7 @@ macro_rules! define_vs_linalg_benches_for_dim { [].bench_function("la_stack_solve_from_lu", |bencher| { bencher.iter(|| { let x = require_ok( - a_lu.solve_vec(black_box(rhs)), + a_lu.solve(black_box(rhs)), "precomputed la_stack LU solve", ); let _ = black_box(x); @@ -276,7 +275,7 @@ macro_rules! define_vs_linalg_benches_for_dim { [].bench_function("la_stack_solve_from_ldlt", |bencher| { bencher.iter(|| { let x = require_ok( - a_ldlt.solve_vec(black_box(rhs)), + a_ldlt.solve(black_box(rhs)), "precomputed la_stack LDLT solve", ); let _ = black_box(x); diff --git a/docs/roadmap.md b/docs/roadmap.md index 5200b45..ae3d49f 100644 --- a/docs/roadmap.md +++ b/docs/roadmap.md @@ -58,12 +58,9 @@ API-invariant cleanup: - [#126](https://github.com/acgetchell/la-stack/issues/126) is resolved as an internal parse-don't-validate design rather than a public proof-bearing API. - `Matrix` and `Vector` remain the raw boundary types so callers can pass - deserialized, fixture, or otherwise unchecked `f64` storage and receive - structured diagnostics. Crate-private `FiniteMatrix` and - `FiniteVector` are validated domain types that carry the finite-entry proof - behind the scenes for LU, LDLT, determinant, error-bound, and exact-arithmetic - paths. + `Matrix` and `Vector` parse raw `f64` storage at construction and then + carry the finite-entry proof directly. Crate-private finite wrappers remain + only as implementation helpers at algorithm boundaries. - The public prelude stays focused on downstream composition: raw boundary types, factorization handles, tolerances, crate errors, dispatch helpers, and documented constants. Proof-bearing wrappers remain crate-private, and diff --git a/examples/det_5x5.rs b/examples/det_5x5.rs index 5829d67..73b5da7 100644 --- a/examples/det_5x5.rs +++ b/examples/det_5x5.rs @@ -14,7 +14,7 @@ fn main() -> Result<(), LaError> { ])?; // Compute via explicit LU factorization. - let lu = a.lu(DEFAULT_PIVOT_TOL)?; + let lu = a.lu(DEFAULT_SINGULAR_TOL)?; let det = lu.det()?; println!("det = {det}"); diff --git a/examples/exact_solve_3x3.rs b/examples/exact_solve_3x3.rs index 3ac505a..83c3c6b 100644 --- a/examples/exact_solve_3x3.rs +++ b/examples/exact_solve_3x3.rs @@ -25,8 +25,8 @@ fn main() -> Result<(), LaError> { let b = Vector::<3>::try_new([1.0, 2.0, 3.0])?; // f64 LU solve (using zero pivot tolerance since the matrix is nearly singular - // and would be rejected by DEFAULT_PIVOT_TOL). - let lu_x = a.lu(Tolerance::new(0.0)?)?.solve_vec(b)?.into_array(); + // and would be rejected by DEFAULT_SINGULAR_TOL). + let lu_x = a.lu(Tolerance::new(0.0)?)?.solve(b)?.into_array(); // Exact solve. let exact_x = a.solve_exact(b)?; diff --git a/examples/ldlt_solve_3x3.rs b/examples/ldlt_solve_3x3.rs index 2a3b4c7..cbaca7d 100644 --- a/examples/ldlt_solve_3x3.rs +++ b/examples/ldlt_solve_3x3.rs @@ -16,7 +16,7 @@ fn main() -> Result<(), LaError> { let b = Vector::<3>::try_new([2.0, 4.0, 10.0])?; let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL)?; - let x = ldlt.solve_vec(b)?.into_array(); + let x = ldlt.solve(b)?.into_array(); let det = ldlt.det()?; println!("A (3×3 SPD tridiagonal):"); diff --git a/examples/solve_5x5.rs b/examples/solve_5x5.rs index 29870f4..d8ae5c0 100644 --- a/examples/solve_5x5.rs +++ b/examples/solve_5x5.rs @@ -16,8 +16,8 @@ fn main() -> Result<(), LaError> { // Choose x = [1, 2, 3, 4, 5]. Then b = A x = [14, 13, 12, 11, 10]. let b = Vector::<5>::try_new([14.0, 13.0, 12.0, 11.0, 10.0])?; - let lu = a.lu(DEFAULT_PIVOT_TOL)?; - let x = lu.solve_vec(b)?.into_array(); + let lu = a.lu(DEFAULT_SINGULAR_TOL)?; + let x = lu.solve(b)?.into_array(); println!("x = {x:?}"); Ok(()) diff --git a/justfile b/justfile index 6fe02d3..2d392c1 100644 --- a/justfile +++ b/justfile @@ -10,9 +10,9 @@ cargo_nextest_version := "0.9.137" cargo_llvm_cov_version := "0.8.7" dprint_version := "0.54.0" git_cliff_version := "2.13.1" -rumdl_version := "0.2.6" +rumdl_version := "0.2.9" taplo_version := "0.10.0" -typos_version := "1.47.1" +typos_version := "1.47.2" zizmor_version := "1.25.2" # Internal helpers: ensure external tooling is installed diff --git a/scripts/postprocess_changelog.py b/scripts/postprocess_changelog.py index da856c8..a511907 100644 --- a/scripts/postprocess_changelog.py +++ b/scripts/postprocess_changelog.py @@ -9,7 +9,9 @@ and code spans as atomic tokens (MD013). 3. Tag bare fenced code blocks with a language (MD040). 4. Normalize indented commit-body headings (MD023). - 5. Strip trailing blank lines (MD012). + 5. Normalize list continuation indentation (MD077). + 6. Normalize escaped email autolinks (MD034). + 7. Strip trailing blank lines (MD012). Usage: postprocess-changelog # default: CHANGELOG.md @@ -87,6 +89,10 @@ } _ENTRY_HEADING_RE = re.compile(r"^(?P#{2,6})\s+(?P.*?)(?:\s+#+\s*)?$") +# git-cliff HTML-escapes co-author email angle brackets. Markdown wants real +# angle brackets for email autolinks, otherwise rumdl treats the address as bare. +_ESCAPED_EMAIL_RE = re.compile(r"<(?P<email>[A-Z0-9._%+-]+@[A-Z0-9.-]+\.[A-Z]{2,})>", re.IGNORECASE) + # This label set is intentionally broad, including release labels such as # "added", "fixed", "changed", "removed", and "deprecated". Rewriting is # only allowed when _is_isolated_body_heading accepts the line; do not relax @@ -415,21 +421,115 @@ def _deindent_orphan(line: str, lines: list[str], idx: int) -> str: return line[2:] if nearest_parent_indent is not None else stripped -def _needs_blank_before(stripped: str, result: list[str]) -> bool: +def _normalize_list_continuation_indent(line: str, lines: list[str], idx: int) -> str: + """ + Normalize generated list continuation indentation. + + git-cliff commit bodies can contain pre-indented prose under a list item. + Markdown treats that prose as list-continuation content, where rumdl's MD077 + expects exactly two spaces past the parent bullet indentation. + """ + stripped = line.lstrip() + if not line.startswith(" ") or not stripped or stripped.startswith(("- ", "* ")): + return line + + our_indent = len(line) - len(stripped) + + for j in range(idx - 1, -1, -1): + prev = lines[j] + if not prev.strip(): + continue + + prev_stripped = prev.lstrip() + if prev_stripped.startswith(("- ", "* ")): + parent_indent = len(prev) - len(prev_stripped) + expected_indent = parent_indent + 2 + if our_indent > expected_indent: + return " " * expected_indent + stripped + return line + + if not prev.startswith(" "): + return line + + return line + + +def _list_item_indent(line: str) -> int | None: + """Return the indentation of a Markdown list item, if *line* is one.""" + stripped = line.lstrip() + if not stripped.startswith(("- ", "* ")): + return None + return len(line) - len(stripped) + + +def _has_previous_peer_list_item(lines: list[str], idx: int, peer_indent: int) -> bool: + """Return true if a prior list item exists at *peer_indent* before *idx*.""" + for j in range(idx - 1, -1, -1): + prev = lines[j] + if not prev.strip(): + continue + + prev_indent = _list_item_indent(prev) + if prev_indent == peer_indent: + return True + if prev_indent is not None or prev.startswith(" "): + continue + return False + + return False + + +def _next_list_item_indent(lines: list[str], idx: int) -> int | None: + """Return the next nonblank line's list-item indentation, if it is a list item.""" + for next_line in lines[idx + 1 :]: + if not next_line.strip(): + continue + return _list_item_indent(next_line) + return None + + +def _is_blank_between_peer_list_items(lines: list[str], idx: int) -> bool: + """Return true when a blank line separates adjacent items in the same list.""" + if lines[idx].strip(): + return False + + next_indent = _next_list_item_indent(lines, idx) + if next_indent is None: + return False + + return _has_previous_peer_list_item(lines, idx, next_indent) + + +def _normalize_email_autolinks(line: str) -> str: + """Convert git-cliff's escaped email autolinks into Markdown autolinks.""" + return _ESCAPED_EMAIL_RE.sub(r"<\g<email>>", line) + + +def _needs_blank_before(line: str, lines: list[str], idx: int, result: list[str]) -> bool: """ Determine whether a blank line is required before a list item to satisfy Markdown rule MD032. Parameters: - stripped (str): The current line with leading whitespace removed. + line (str): The current line. + lines (list[str]): The source lines being post-processed. + idx (int): The index of the current line in ``lines``. result (list[str]): The lines already emitted immediately before the current line. Returns: bool: `True` if a blank line should be inserted before the list item, `False` otherwise. """ + stripped = line.lstrip() if not stripped.startswith("- ") or not result or not result[-1].strip(): return False + prev = result[-1].lstrip() - return not prev.startswith(("-", "#")) + if prev.startswith("#"): + return False + if prev.startswith("- "): + return False + + current_indent = len(line) - len(stripped) + return not _has_previous_peer_list_item(result, len(result), current_indent) def _normalize_indented_heading(line: str) -> str: @@ -534,13 +634,14 @@ def _normalize_body_line(line: str, lines: list[str], idx: int, result: list[str """Apply markdown hygiene transforms to a non-code line.""" is_isolated_body_heading = _is_isolated_body_heading(lines, idx) line = _deindent_orphan(line, lines, idx) + line = _normalize_list_continuation_indent(line, lines, idx) line = _normalize_indented_heading(line) line = _normalize_entry_heading(line) if is_isolated_body_heading: line = _normalize_squash_heading(line, nested=current_entry_summary is not None) - if _needs_blank_before(line.lstrip(), result): + if _needs_blank_before(line, lines, idx, result): result.append("") return _reflow_line(line) if len(line) > MAX_LINE_WIDTH else line @@ -571,11 +672,15 @@ def postprocess_text(text: str) -> str: result.append(line) continue + if _is_blank_between_peer_list_items(lines, idx): + continue + # --- MD004: normalise ``* `` list markers to ``- `` --- line = _STAR_LIST_RE.sub(r"\1- ", line) # --- MD030: normalise spaces after list marker --- line = _LIST_MARKER_SPACE_RE.sub(r"\1 ", line) + line = _normalize_email_autolinks(line) current_entry_summary = _update_entry_summary(line, current_entry_summary) is_isolated_body_heading = _is_isolated_body_heading(lines, idx) diff --git a/scripts/tests/test_postprocess_changelog.py b/scripts/tests/test_postprocess_changelog.py index 61a7107..18a5edb 100644 --- a/scripts/tests/test_postprocess_changelog.py +++ b/scripts/tests/test_postprocess_changelog.py @@ -10,8 +10,10 @@ _is_duplicate_squash_heading, _is_isolated_body_heading, _max_pr_number, + _normalize_email_autolinks, _normalize_entry_heading, _normalize_indented_heading, + _normalize_list_continuation_indent, _normalize_squash_heading, _plain_summary, _process_code_fence, @@ -83,6 +85,119 @@ def test_empty_file(self, tmp_path: Path) -> None: assert f.read_text(encoding="utf-8") == "\n" +class TestListContinuationIndent: + def test_normalizes_over_indented_top_level_continuation(self) -> None: + lines = [ + "- Replace BigRational-only gauss_solve", + "", + " Bareiss fraction-free forward elimination.", + ] + + result = _normalize_list_continuation_indent(lines[2], lines, 2) + + assert result == " Bareiss fraction-free forward elimination." + + def test_preserves_nested_list_continuation_indent(self) -> None: + lines = [ + "- Parent item", + "", + " - Nested item", + " continued nested prose", + ] + + result = _normalize_list_continuation_indent(lines[3], lines, 3) + + assert result == " continued nested prose" + + def test_full_postprocess_normalizes_git_cliff_continuation(self) -> None: + content = ( + "# Changelog\n\n" + "## [1.0.0] - 2026-01-01\n\n" + "### Performance\n\n" + "- Integer-only Bareiss determinant via BigInt\n" + " [`d422b25`](https://github.com/acgetchell/la-stack/commit/d422b251ca86a914522f80285964d4513bca1817)\n\n" + " det_exact: 16x faster\n" + ) + + result = postprocess_text(content) + + assert "\n det_exact: 16x faster\n" in result + assert "\n det_exact: 16x faster\n" not in result + + +class TestRumdlCompatibility: + def test_removes_blank_between_peer_list_items(self) -> None: + content = ( + "# Changelog\n\n" + "## [1.0.0] - 2026-01-01\n\n" + "### Added\n\n" + "- First item\n" + " [`1111111`](https://github.com/acgetchell/la-stack/commit/1111111111111111111111111111111111111111)\n\n" + "- Second item\n" + ) + + result = postprocess_text(content) + + assert "\n\n- Second item\n" not in result + assert "\n- Second item\n" in result + + def test_removes_blank_after_nested_body_before_top_level_peer(self) -> None: + content = "# Changelog\n\n## [1.0.0] - 2026-01-01\n\n### Added\n\n- First item\n\n - Body point\n - Body point\n\n- Second item\n" + + result = postprocess_text(content) + + assert "\n\n- Second item\n" not in result + assert "\n- Second item\n" in result + + def test_removes_blank_after_mixed_body_before_top_level_peer(self) -> None: + content = ( + "# Changelog\n\n" + "## [1.0.0] - 2026-01-01\n\n" + "### Changed\n\n" + "- First item\n\n" + " Body paragraph.\n\n" + " - Body point\n" + " - Body point\n\n" + " More body prose.\n\n" + "- Second item\n" + ) + + result = postprocess_text(content) + + assert "\n\n- Second item\n" not in result + assert "\n- Second item\n" in result + + def test_preserves_blank_before_list_item_body(self) -> None: + content = "# Changelog\n\n## [1.0.0] - 2026-01-01\n\n### Added\n\n- First item\n\n Body paragraph.\n" + + result = postprocess_text(content) + + assert "- First item\n\n Body paragraph.\n" in result + + def test_keeps_body_item_tight_after_plain_peer(self) -> None: + content = ( + "# Changelog\n\n" + "## [1.0.0] - 2026-01-01\n\n" + "### Added\n\n" + "- Simple item [`1111111`](https://github.com/acgetchell/la-stack/commit/1111111111111111111111111111111111111111)\n" + "- Item with body [`2222222`](https://github.com/acgetchell/la-stack/commit/2222222222222222222222222222222222222222)\n\n" + " Body paragraph.\n" + ) + + result = postprocess_text(content) + + assert "\n- Simple item" in result + assert "\n- Item with body" in result + assert "\n\n- Item with body" not in result + + def test_normalizes_git_cliff_escaped_email_autolink(self) -> None: + line = " Co-Authored-By: Oz <oz-agent@warp.dev>" + + result = _normalize_email_autolinks(line) + + assert result == " Co-Authored-By: Oz <oz-agent@warp.dev>" + + class TestReflowLine: """Unit tests for the _reflow_line helper.""" diff --git a/semgrep.yaml b/semgrep.yaml index 26de90d..793229f 100644 --- a/semgrep.yaml +++ b/semgrep.yaml @@ -76,6 +76,96 @@ rules: - "/tests/semgrep/src/project_rules/raw_f64_constructors.rs" pattern-regex: '(?m)^\s*pub\s+(?:const\s+)?fn\s+(?:new|from_rows)\s*\([^)]*(?:\[\s*f64\s*;\s*D\s*\]|\[\s*\[\s*f64\s*;\s*D\s*\]\s*;\s*D\s*\])[^)]*\)\s*->\s*(?:Self|(?:Matrix|Vector)\s*<)' + - id: la-stack.rust.no-public-unchecked-finite-constructors + languages: + - rust + severity: WARNING + message: "Unchecked Matrix/Vector constructors must stay crate-private; public callers must parse through fallible constructors." + metadata: + category: correctness + rationale: >- + Matrix and Vector are finite proof types. Public unchecked constructors + would let downstream callers bypass parsing and construct invalid finite + proofs. + paths: + include: + - "/src/**/*.rs" + - "/tests/semgrep/src/project_rules/finite_api_contract.rs" + pattern-regex: '(?m)^\s*pub\s+(?:const\s+)?fn\s+(?:new_unchecked|from_rows_unchecked)\s*\(' + + - id: la-stack.rust.no-public-matrix-vector-storage-fields + languages: + - regex + severity: WARNING + message: "Matrix/Vector storage fields must stay private so callers cannot bypass finite parsing." + metadata: + category: correctness + rationale: >- + Matrix and Vector carry the finite-entry invariant. Public or crate-public + raw storage fields allow construction or mutation paths that bypass the + parsing boundary. + paths: + include: + - "/src/**/*.rs" + - "/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 + severity: WARNING + message: "Keep matrix/vector modules private; re-export only the curated public API from crate root/prelude." + 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. + paths: + include: + - "/src/**/*.rs" + - "/tests/semgrep/src/project_rules/finite_api_contract.rs" + pattern-regex: '(?m)^\s*pub\s+mod\s+(?:matrix|vector)\s*;' + + - id: la-stack.rust.no-legacy-solve-vec-api + languages: + - regex + severity: WARNING + message: "Use the clean solve API name; do not reintroduce solve_vec." + metadata: + category: api + rationale: >- + Lu::solve and Ldlt::solve are the clean API because Vector is already + the finite right-hand-side proof type. + paths: + include: + - "/src/**/*.rs" + - "/benches/**/*.rs" + - "/examples/**/*.rs" + - "/tests/**/*.rs" + - "/README.md" + exclude: + - "/tests/semgrep/src/project_rules/raw_f64_constructors.rs" + - "/tests/semgrep/src/project_rules/public_api_panic_paths.rs" + - "/tests/semgrep/src/project_rules/bench_example_usage.rs" + pattern-regex: '\bsolve_vec\b' + - id: la-stack.rust.no-public-api-panic-paths languages: - regex @@ -125,6 +215,23 @@ rules: - "/tests/semgrep/src/**" pattern-regex: '^\s*//[!/]\s*(?:#\s*)?.*(?:\b[\w:]+|[\]\)])\.(unwrap|expect)\s*\(' + - id: la-stack.rust.no-unwrap-expect-in-readme-doctest-mirrors + languages: + - regex + severity: WARNING + message: "Use fallible flow instead of unwrap() or expect() in src/lib.rs README doctest mirrors." + metadata: + category: correctness + rationale: >- + The private readme_doctests module in src/lib.rs keeps README examples + executable for docs.rs-facing documentation. Those mirrors should model + the same fallible Result/? flow as the public examples they verify. + paths: + 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*\()' + - id: la-stack.rust.no-unwrap-expect-in-benches-examples languages: - rust diff --git a/src/exact.rs b/src/exact.rs index 8d978de..5ac3cec 100644 --- a/src/exact.rs +++ b/src/exact.rs @@ -60,11 +60,11 @@ //! //! ## Validation //! -//! Public exact methods parse raw `Matrix` / `Vector` values into internal -//! finite proof-bearing types before reaching the integer-Bareiss core. The -//! decomposition helpers for those validated domain types can then call -//! `f64_decompose` without repeating stored-entry validation; `f64_decompose` -//! itself is therefore never called with non-finite input from the public API. +//! Public `Matrix` / `Vector` values are finite by construction before exact +//! methods reach the integer-Bareiss core. The decomposition helpers for those +//! domain types can then call `f64_decompose` without repeating stored-entry +//! validation; `f64_decompose` itself is therefore never called with non-finite +//! input from the public API. use core::hint::cold_path; use core::mem::take; @@ -224,7 +224,7 @@ fn decompose_finite_matrix<const D: usize>( 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 matrix.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))? @@ -455,9 +455,9 @@ fn bareiss_det_finite<const D: usize>(m: &FiniteMatrix<D>) -> Result<BigRational /// Solve `A x = b` exactly after matrix and RHS finiteness has been proven. /// -/// Public solve APIs parse raw [`Matrix`] / [`Vector`] values before reaching -/// this helper, so decomposition can proceed without rediscovering stored -/// NaN/∞ entries. +/// Public [`Matrix`] / [`Vector`] values are finite by construction before +/// reaching this helper, so decomposition can proceed without rediscovering +/// stored NaN/∞ entries. /// /// # Errors /// Returns [`LaError::Singular`] if the matrix is exactly singular. @@ -642,13 +642,11 @@ impl<const D: usize> Matrix<D> { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] if stored matrix entries are NaN or infinity. - /// /// Returns [`LaError::DeterminantScaleOverflow`] if determinant scaling /// overflows the internal exponent representation. #[inline] pub fn det_exact(&self) -> Result<BigRational, LaError> { - FiniteMatrix::new(*self)?.det_exact() + FiniteMatrix::new_unchecked(*self).det_exact() } /// Exact determinant converted to `f64`. @@ -675,8 +673,6 @@ impl<const D: usize> Matrix<D> { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] if stored matrix entries are NaN or infinity. - /// /// Returns [`LaError::DeterminantScaleOverflow`] if determinant scaling /// overflows the internal exponent representation. /// @@ -684,7 +680,7 @@ impl<const D: usize> Matrix<D> { /// represent as a finite `f64`. #[inline] pub fn det_exact_f64(&self) -> Result<f64, LaError> { - FiniteMatrix::new(*self)?.det_exact_f64() + FiniteMatrix::new_unchecked(*self).det_exact_f64() } /// Exact linear system solve using hybrid integer/rational arithmetic. @@ -730,14 +726,11 @@ impl<const D: usize> Matrix<D> { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] if stored matrix or right-hand-side entries - /// are NaN or infinity. - /// /// Returns [`LaError::Singular`] if the matrix is exactly singular. #[inline] pub fn solve_exact(&self, b: Vector<D>) -> Result<[BigRational; D], LaError> { - let finite_m = FiniteMatrix::new(*self)?; - let finite_b = FiniteVector::new(b)?; + let finite_m = FiniteMatrix::new_unchecked(*self); + let finite_b = FiniteVector::new_unchecked(b); finite_m.solve_exact(finite_b) } @@ -765,16 +758,13 @@ impl<const D: usize> Matrix<D> { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] if stored matrix or right-hand-side entries - /// are NaN or infinity. - /// /// 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`. #[inline] pub fn solve_exact_f64(&self, b: Vector<D>) -> Result<Vector<D>, LaError> { - let finite_m = FiniteMatrix::new(*self)?; - let finite_b = FiniteVector::new(b)?; + let finite_m = FiniteMatrix::new_unchecked(*self); + let finite_b = FiniteVector::new_unchecked(b); Ok(finite_m.solve_exact_f64(finite_b)?.into_vector()) } @@ -815,20 +805,18 @@ impl<const D: usize> Matrix<D> { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] if stored matrix entries are NaN or infinity. - /// /// Returns [`LaError::DeterminantScaleOverflow`] if determinant scaling /// overflows the internal exponent representation. #[inline] pub fn det_sign_exact(&self) -> Result<i8, LaError> { - FiniteMatrix::new(*self)?.det_sign_exact() + FiniteMatrix::new_unchecked(*self).det_sign_exact() } } #[cfg(test)] mod tests { use super::*; - use crate::DEFAULT_PIVOT_TOL; + use crate::DEFAULT_SINGULAR_TOL; use core::assert_matches; use num_traits::Signed; @@ -912,53 +900,6 @@ mod tests { gen_internal_finite_exact_tests!(4); gen_internal_finite_exact_tests!(5); - macro_rules! gen_public_exact_revalidation_tests { - ($d:literal) => { - paste! { - #[test] - fn [<public_exact_methods_revalidate_unchecked_matrix_storage_ $d d>]() { - let mut rows = [[0.0f64; $d]; $d]; - for i in 0..$d { - rows[i][i] = 1.0; - } - rows[0][$d - 1] = f64::NAN; - let raw = Matrix::<$d>::from_rows_unchecked(rows); - let rhs = Vector::<$d>::new([1.0; $d]); - let expected = LaError::NonFinite { - row: Some(0), - col: $d - 1, - }; - - assert_eq!(raw.det_exact(), Err(expected)); - assert_eq!(raw.det_exact_f64(), Err(expected)); - assert_eq!(raw.det_sign_exact(), Err(expected)); - assert_eq!(raw.solve_exact(rhs), Err(expected)); - assert_eq!(raw.solve_exact_f64(rhs), Err(expected)); - } - - #[test] - fn [<public_exact_solve_methods_revalidate_unchecked_rhs_storage_ $d d>]() { - let matrix = Matrix::<$d>::identity(); - let mut rhs = [1.0; $d]; - rhs[$d - 1] = f64::INFINITY; - let rhs = Vector::<$d>::new_unchecked(rhs); - let expected = LaError::NonFinite { - row: None, - col: $d - 1, - }; - - assert_eq!(matrix.solve_exact(rhs), Err(expected)); - assert_eq!(matrix.solve_exact_f64(rhs), Err(expected)); - } - } - }; - } - - gen_public_exact_revalidation_tests!(2); - gen_public_exact_revalidation_tests!(3); - gen_public_exact_revalidation_tests!(4); - gen_public_exact_revalidation_tests!(5); - macro_rules! gen_det_exact_tests { ($d:literal) => { paste! { @@ -1715,7 +1656,7 @@ mod tests { let b = arbitrary_rhs::<$d>(); let x = a.solve_exact(b).unwrap(); for (i, xi) in x.iter().enumerate() { - assert_eq!(*xi, f64_to_bigrational(b.data[i])); + assert_eq!(*xi, f64_to_bigrational(b.as_array()[i])); } } @@ -1744,7 +1685,7 @@ mod tests { let b = arbitrary_rhs::<$d>(); let x = a.solve_exact_f64(b).unwrap().into_array(); for i in 0..$d { - assert!((x[i] - b.data[i]).abs() <= f64::EPSILON); + assert!((x[i] - b.as_array()[i]).abs() <= f64::EPSILON); } } } @@ -1756,7 +1697,7 @@ mod tests { gen_solve_exact_f64_tests!(4); gen_solve_exact_f64_tests!(5); - /// For D ≤ 4, `solve_exact_f64` should agree with `Lu::solve_vec` on + /// For D ≤ 4, `solve_exact_f64` should agree with `Lu::solve` on /// well-conditioned matrices. macro_rules! gen_solve_exact_f64_agrees_with_lu { ($d:literal) => { @@ -1778,8 +1719,8 @@ mod tests { let a = Matrix::<$d>::from_rows(rows); let b = arbitrary_rhs::<$d>(); let exact = a.solve_exact_f64(b).unwrap().into_array(); - let lu_sol = a.lu(DEFAULT_PIVOT_TOL).unwrap() - .solve_vec(b).unwrap().into_array(); + let lu_sol = a.lu(DEFAULT_SINGULAR_TOL).unwrap() + .solve(b).unwrap().into_array(); for i in 0..$d { let eps = lu_sol[i].abs().mul_add(1e-12, 1e-12); assert!((exact[i] - lu_sol[i]).abs() <= eps); @@ -2215,7 +2156,7 @@ mod tests { fn bigrational_matvec<const D: usize>(a: &Matrix<D>, x: &[BigRational; D]) -> [BigRational; D] { from_fn(|i| { let mut sum = BigRational::from_integer(BigInt::from(0)); - for (aij, xj) in a.rows[i].iter().zip(x.iter()) { + for (aij, xj) in a.rows()[i].iter().zip(x.iter()) { sum += f64_to_bigrational(*aij) * xj; } sum diff --git a/src/ldlt.rs b/src/ldlt.rs index 02940d7..407d584 100644 --- a/src/ldlt.rs +++ b/src/ldlt.rs @@ -59,66 +59,71 @@ impl<const D: usize> LdltFactors<D> { #[inline] #[must_use] const fn row(&self, index: usize) -> &[f64; D] { - &self.storage.rows[index] + &self.storage.rows()[index] } /// Return a factor entry. #[inline] #[must_use] const fn entry(&self, row: usize, col: usize) -> f64 { - self.storage.rows[row][col] + self.storage.rows()[row][col] } /// Return a diagonal entry of `D`. #[inline] #[must_use] const fn diag(&self, index: usize) -> f64 { - self.storage.rows[index][index] + self.storage.rows()[index][index] } } impl<const D: usize> Ldlt<D> { /// Factor a matrix that has already passed LDLT symmetry validation. #[inline] + #[allow(clippy::needless_range_loop)] pub(crate) fn factor_symmetric(a: SymmetricMatrix<D>, tol: Tolerance) -> Result<Self, LaError> { let mut f = a.into_matrix(); let tol = tol.get(); - // LDLT via symmetric rank-1 updates, using only the lower triangle. - for j in 0..D { - let d = f.rows[j][j]; - if !d.is_finite() { - cold_path(); - return Err(LaError::non_finite_cell(j, j)); - } - if d < 0.0 { - cold_path(); - return Err(LaError::not_positive_semidefinite(j, d)); - } - if d <= tol { - cold_path(); - return Err(LaError::Singular { pivot_col: j }); - } + { + let rows = f.rows_mut_unchecked(); - // Compute L multipliers below the diagonal in column j. - for i in (j + 1)..D { - let l = f.rows[i][j] / d; - if !l.is_finite() { + // LDLT via symmetric rank-1 updates, using only the lower triangle. + for j in 0..D { + let d = rows[j][j]; + if !d.is_finite() { cold_path(); - return Err(LaError::non_finite_cell(i, j)); + return Err(LaError::non_finite_cell(j, j)); + } + if d < 0.0 { + cold_path(); + return Err(LaError::not_positive_semidefinite(j, d)); + } + if d <= tol { + cold_path(); + return Err(LaError::Singular { pivot_col: j }); + } + + // Compute L multipliers below the diagonal in column j. + for i in (j + 1)..D { + let l = rows[i][j] / d; + if !l.is_finite() { + cold_path(); + return Err(LaError::non_finite_cell(i, j)); + } + rows[i][j] = l; } - f.rows[i][j] = l; - } - // Update the trailing submatrix (lower triangle): A := A - (L_col * d) * L_col^T. - for i in (j + 1)..D { - let l_i = f.rows[i][j]; - let l_i_d = l_i * d; + // Update the trailing submatrix (lower triangle): A := A - (L_col * d) * L_col^T. + for i in (j + 1)..D { + let l_i = rows[i][j]; + let l_i_d = l_i * d; - for k in (j + 1)..=i { - let l_k = f.rows[k][j]; - let new_val = (-l_i_d).mul_add(l_k, f.rows[i][k]); - f.rows[i][k] = new_val; + for k in (j + 1)..=i { + let l_k = rows[k][j]; + let new_val = (-l_i_d).mul_add(l_k, rows[i][k]); + rows[i][k] = new_val; + } } } } @@ -168,8 +173,10 @@ impl<const D: usize> Ldlt<D> { /// Solve `A x = b` using this LDLT factorization. /// - /// `solve_vec` performs floating-point substitution and does not provide a - /// certified absolute rounding-error bound for the returned solution. + /// [`Vector`] is finite by construction, so this method only checks computed + /// substitution overflows. It performs floating-point substitution and does + /// not provide a certified absolute rounding-error bound for the returned + /// solution. /// /// # Examples /// ``` @@ -180,7 +187,7 @@ impl<const D: usize> Ldlt<D> { /// let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL)?; /// /// let b = Vector::<2>::try_new([1.0, 2.0])?; - /// let x = ldlt.solve_vec(b)?.into_array(); + /// let x = ldlt.solve(b)?.into_array(); /// /// assert!((x[0] - (-0.125)).abs() <= 1e-12); /// assert!((x[1] - 0.75).abs() <= 1e-12); @@ -189,18 +196,11 @@ impl<const D: usize> Ldlt<D> { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] if the right-hand side contains - /// NaN/infinity or a computed substitution intermediate overflows to NaN or - /// infinity. The right-hand side is revalidated at this boundary before - /// entering substitution; public callers normally encounter the same - /// rejection earlier through [`Vector::try_new`](crate::Vector::try_new). + /// Returns [`LaError::NonFinite`] if a computed substitution intermediate + /// overflows to NaN or infinity. #[inline] - pub const fn solve_vec(&self, b: Vector<D>) -> Result<Vector<D>, LaError> { - let b = match FiniteVector::new(b) { - Ok(b) => b, - Err(err) => return Err(err), - }; - match self.solve_finite_vec(b) { + pub const fn solve(&self, b: Vector<D>) -> Result<Vector<D>, LaError> { + match self.solve_finite(FiniteVector::new_unchecked(b)) { Ok(x) => Ok(x.into_vector()), Err(err) => Err(err), } @@ -215,7 +215,7 @@ impl<const D: usize> Ldlt<D> { /// Returns [`LaError::NonFinite`] if a computed substitution intermediate /// overflows to NaN or infinity. #[inline] - pub(crate) const fn solve_finite_vec( + pub(crate) const fn solve_finite( &self, b: FiniteVector<D>, ) -> Result<FiniteVector<D>, LaError> { @@ -306,7 +306,7 @@ mod tests { arr }; let b = Vector::<$d>::new(black_box(b_arr)); - let x = ldlt.solve_vec(b).unwrap().into_array(); + let x = ldlt.solve(b).unwrap().into_array(); for i in 0..$d { assert_abs_diff_eq!(x[i], b_arr[i], epsilon = 1e-12); @@ -362,7 +362,7 @@ mod tests { }; let b = Vector::<$d>::new(black_box(b_arr)); - let x = ldlt.solve_vec(b).unwrap().into_array(); + let x = ldlt.solve(b).unwrap().into_array(); for i in 0..$d { assert_abs_diff_eq!(x[i], b_arr[i] / diag[i], epsilon = 1e-12); @@ -386,7 +386,7 @@ mod tests { .unwrap(); let b = Vector::<2>::new(black_box([1.0, 2.0])); - let x = ldlt.solve_vec(b).unwrap().into_array(); + let x = ldlt.solve(b).unwrap().into_array(); assert_abs_diff_eq!(x[0], -0.125, epsilon = 1e-12); assert_abs_diff_eq!(x[1], 0.75, epsilon = 1e-12); @@ -404,7 +404,7 @@ mod tests { // Choose x = 1 so b = A x is simple: [1, 0, 1]. let b = Vector::<3>::new(black_box([1.0, 0.0, 1.0])); - let x = ldlt.solve_vec(b).unwrap().into_array(); + let x = ldlt.solve(b).unwrap().into_array(); for &x_i in &x { assert_abs_diff_eq!(x_i, 1.0, epsilon = 1e-9); @@ -499,7 +499,7 @@ mod tests { } #[test] - fn nonfinite_solve_vec_forward_substitution_overflow() { + 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([ @@ -510,12 +510,12 @@ mod tests { let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<3>::new([1e156, 0.0, 0.0]); - let err = ldlt.solve_vec(b).unwrap_err(); + let err = ldlt.solve(b).unwrap_err(); assert_eq!(err, LaError::NonFinite { row: None, col: 1 }); } #[test] - fn nonfinite_solve_vec_back_substitution_overflow() { + fn nonfinite_solve_back_substitution_overflow() { // SPD matrix: [[1,0,0],[0,1,2],[0,2,5]] has LDLT factors // D=[1,1,1], L[2,1]=2. Forward sub and diagonal solve produce // z=[0,0,1e308]. Back-substitution: x[2]=1e308 then @@ -524,12 +524,12 @@ mod tests { let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<3>::new([0.0, 0.0, 1e308]); - let err = ldlt.solve_vec(b).unwrap_err(); + let err = ldlt.solve(b).unwrap_err(); assert_eq!(err, LaError::NonFinite { row: None, col: 1 }); } #[test] - fn nonfinite_solve_vec_diagonal_solve_overflow() { + fn nonfinite_solve_diagonal_solve_overflow() { // Diagonal SPD matrix with a tiny diagonal entry just above the // singularity tolerance. Forward substitution passes through the // large RHS unchanged, then the diagonal solve z[1] = y[1] / D[1] @@ -539,7 +539,7 @@ mod tests { let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<2>::new([0.0, 1.0e300]); - let err = ldlt.solve_vec(b).unwrap_err(); + let err = ldlt.solve(b).unwrap_err(); assert_eq!(err, LaError::NonFinite { row: None, col: 1 }); } @@ -570,13 +570,13 @@ mod tests { ); } - macro_rules! gen_solve_vec_boundary_tests { + macro_rules! gen_solve_boundary_tests { ($d:literal) => { paste! { /// Raw non-finite right-hand sides are rejected before a /// public caller can construct a `Vector`. #[test] - fn [<solve_vec_rhs_boundary_rejects_non_finite_ $d d>]() { + fn [<solve_rhs_constructor_rejects_non_finite_ $d d>]() { let mut rhs = [1.0; $d]; rhs[$d - 1] = f64::NAN; @@ -588,35 +588,19 @@ mod tests { }) ); } - - #[test] - fn [<solve_vec_revalidates_unchecked_rhs_storage_ $d d>]() { - let ldlt = Matrix::<$d>::identity().ldlt(DEFAULT_SINGULAR_TOL).unwrap(); - let mut rhs = [1.0; $d]; - rhs[$d - 1] = f64::NAN; - let rhs = Vector::<$d>::new_unchecked(rhs); - - assert_eq!( - ldlt.solve_vec(rhs), - Err(LaError::NonFinite { - row: None, - col: $d - 1, - }) - ); - } } }; } - gen_solve_vec_boundary_tests!(2); - gen_solve_vec_boundary_tests!(3); - gen_solve_vec_boundary_tests!(4); - gen_solve_vec_boundary_tests!(5); + gen_solve_boundary_tests!(2); + gen_solve_boundary_tests!(3); + gen_solve_boundary_tests!(4); + gen_solve_boundary_tests!(5); // ----------------------------------------------------------------------- // Const-evaluability tests. // - // These prove that `Ldlt::det` and `Ldlt::solve_vec` are truly `const fn` + // These prove that `Ldlt::det` and `Ldlt::solve` are truly `const fn` // by forcing the compiler to evaluate them inside a `const` initializer. // `Ldlt::factor` is not (yet) `const fn` because the rank-1 update loop // uses array indexing patterns that still require non-const helpers on @@ -633,8 +617,14 @@ mod tests { #[test] fn [<ldlt_det_const_eval_ $d d>]() { const DET: Result<f64, LaError> = { - let mut factors = Matrix::<$d>::identity(); - factors.rows[0][0] = 2.0; + let mut rows = [[0.0f64; $d]; $d]; + let mut i = 0; + while i < $d { + rows[i][i] = 1.0; + i += 1; + } + rows[0][0] = 2.0; + let factors = Matrix::<$d>::from_rows_unchecked(rows); let ldlt = Ldlt::<$d> { factors: LdltFactors::new_unchecked(factors), }; @@ -643,12 +633,12 @@ mod tests { assert_eq!(DET, Ok(2.0)); } - /// `Ldlt::solve_vec` must be fully const-evaluable. Identity + /// `Ldlt::solve` must be fully const-evaluable. Identity /// factors with RHS `b = [1.0, 2.0, …, D]` round-trips `b` /// unchanged, exercising the full forward sub / diagonal solve /// / back sub pipeline inside a `const { … }` initializer. #[test] - fn [<ldlt_solve_vec_const_eval_ $d d>]() { + fn [<ldlt_solve_const_eval_ $d d>]() { #[allow(clippy::cast_precision_loss)] const X: [f64; $d] = { let ldlt = Ldlt::<$d> { @@ -661,7 +651,7 @@ mod tests { i += 1; } let b = Vector::<$d>::new(b_arr); - match ldlt.solve_vec(b) { + match ldlt.solve(b) { Ok(v) => v.into_array(), Err(_) => [0.0f64; $d], } diff --git a/src/lib.rs b/src/lib.rs index d0fee4a..63109eb 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -20,8 +20,8 @@ mod readme_doctests { /// /// let b = Vector::<5>::try_new([14.0, 13.0, 12.0, 11.0, 10.0])?; /// - /// let lu = a.lu(DEFAULT_PIVOT_TOL)?; - /// let x = lu.solve_vec(b)?.into_array(); + /// let lu = a.lu(DEFAULT_SINGULAR_TOL)?; + /// let x = lu.solve(b)?.into_array(); /// /// // Floating-point rounding is expected; compare with a tolerance. /// let expected = [1.0, 2.0, 3.0, 4.0, 5.0]; @@ -351,12 +351,6 @@ impl Tolerance { /// as "numerically zero" (e.g. LU pivots, LDLT diagonal entries). pub const DEFAULT_SINGULAR_TOL: Tolerance = Tolerance::new_unchecked(1e-12); -/// Default absolute pivot magnitude threshold used for LU pivot selection / singularity detection. -/// -/// This name is kept for backwards compatibility; prefer [`DEFAULT_SINGULAR_TOL`] when the -/// tolerance is not specifically about pivot selection. -pub const DEFAULT_PIVOT_TOL: Tolerance = DEFAULT_SINGULAR_TOL; - /// Relative tolerance used to validate matrices for LDLT factorization. pub(crate) const LDLT_SYMMETRY_REL_TOL: Tolerance = Tolerance::new_unchecked(1e-12); @@ -804,9 +798,9 @@ macro_rules! try_with_stack_matrix { /// /// This prelude re-exports the primary types and constants: [`Matrix`], /// [`Vector`], [`Lu`], [`Ldlt`], [`Tolerance`], [`LaError`], -/// [`DEFAULT_PIVOT_TOL`], [`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 +/// [`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 @@ -818,8 +812,8 @@ macro_rules! try_with_stack_matrix { /// for `.is_positive()` / `.is_negative()` / `.abs()`. pub mod prelude { pub use crate::{ - DEFAULT_PIVOT_TOL, 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, + 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, }; #[cfg(feature = "exact")] @@ -837,11 +831,6 @@ mod tests { #[test] fn default_singular_tol_is_expected() { assert_abs_diff_eq!(DEFAULT_SINGULAR_TOL.get(), 1e-12, epsilon = 0.0); - assert_abs_diff_eq!( - DEFAULT_PIVOT_TOL.get(), - DEFAULT_SINGULAR_TOL.get(), - epsilon = 0.0 - ); } #[test] @@ -1030,15 +1019,18 @@ mod tests { } #[test] - fn prelude_reexports_compile_and_work() { + 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>::new([1.0, 2.0]); - let _ = m.lu(DEFAULT_PIVOT_TOL).unwrap().solve_vec(v).unwrap(); - let _ = m.ldlt(DEFAULT_SINGULAR_TOL).unwrap().solve_vec(v).unwrap(); + 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(()) } macro_rules! gen_stack_matrix_dispatch_tests { @@ -1106,7 +1098,7 @@ mod tests { #[test] fn try_with_stack_matrix_converts_unsupported_dimension_error() { let got = try_with_stack_matrix!(9usize, |m| -> Result<usize, DownstreamError> { - assert_abs_diff_eq!(m.inf_norm().unwrap(), 0.0, epsilon = 0.0); + assert_abs_diff_eq!(m.inf_norm()?, 0.0, epsilon = 0.0); Ok(0) }); @@ -1133,8 +1125,10 @@ mod tests { assert_eq!(*r.numer(), n); // `FromPrimitive::from_f64` / `from_i64` on `BigRational`. - let half = BigRational::from_f64(0.5).unwrap(); - let two = BigRational::from_i64(2).unwrap(); + let half = BigRational::new(BigInt::from(1), BigInt::from(2)); + let two = BigRational::from_integer(BigInt::from(2)); + assert_eq!(BigRational::from_f64(0.5), Some(half.clone())); + assert_eq!(BigRational::from_i64(2), Some(two.clone())); assert_eq!( half.clone() + half.clone(), BigRational::from_integer(BigInt::from(1)) @@ -1148,7 +1142,7 @@ mod tests { assert_eq!(neg.abs(), half); // `ToPrimitive::to_f64` / `to_i64`. - assert!((half.to_f64().unwrap() - 0.5).abs() <= f64::EPSILON); - assert_eq!(two.to_i64().unwrap(), 2); + assert_eq!(half.to_f64(), Some(0.5)); + assert_eq!(two.to_i64(), Some(2)); } } diff --git a/src/lu.rs b/src/lu.rs index 21a6eea..5520ba6 100644 --- a/src/lu.rs +++ b/src/lu.rs @@ -35,14 +35,14 @@ impl<const D: usize> LuFactors<D> { #[inline] #[must_use] const fn row(&self, index: usize) -> &[f64; D] { - &self.storage.rows[index] + &self.storage.rows()[index] } /// Return a diagonal entry of `U`. #[inline] #[must_use] const fn diag(&self, index: usize) -> f64 { - self.storage.rows[index][index] + self.storage.rows()[index][index] } } @@ -56,6 +56,7 @@ impl<const D: usize> Lu<D> { /// checked before return so successful factors do not contain a non-finite /// value produced during elimination. #[inline] + #[allow(clippy::needless_range_loop)] pub(crate) fn factor_finite(a: FiniteMatrix<D>, tol: Tolerance) -> Result<Self, LaError> { let mut lu = a.into_matrix(); let tol = tol.get(); @@ -67,40 +68,44 @@ impl<const D: usize> Lu<D> { let mut piv_sign = 1.0; - for k in 0..D { - // Choose pivot row. - let mut pivot_row = k; - let mut pivot_abs = lu.rows[k][k].abs(); + { + let rows = lu.rows_mut_unchecked(); - for r in (k + 1)..D { - let v = lu.rows[r][k].abs(); - if v > pivot_abs { - pivot_abs = v; - pivot_row = r; + for k in 0..D { + // Choose pivot row. + let mut pivot_row = k; + let mut pivot_abs = rows[k][k].abs(); + + for r in (k + 1)..D { + let v = rows[r][k].abs(); + if v > pivot_abs { + pivot_abs = v; + pivot_row = r; + } } - } - if pivot_abs <= tol { - cold_path(); - return Err(LaError::Singular { pivot_col: k }); - } + if pivot_abs <= tol { + cold_path(); + return Err(LaError::Singular { pivot_col: k }); + } - if pivot_row != k { - lu.rows.swap(k, pivot_row); - piv.swap(k, pivot_row); - piv_sign = -piv_sign; - } + if pivot_row != k { + rows.swap(k, pivot_row); + piv.swap(k, pivot_row); + piv_sign = -piv_sign; + } - let pivot = lu.rows[k][k]; + let pivot = rows[k][k]; - // Eliminate below pivot. - for r in (k + 1)..D { - let mult = lu.rows[r][k] / pivot; - lu.rows[r][k] = mult; + // Eliminate below pivot. + for r in (k + 1)..D { + let mult = rows[r][k] / pivot; + rows[r][k] = mult; - for c in (k + 1)..D { - let updated = (-mult).mul_add(lu.rows[k][c], lu.rows[r][c]); - lu.rows[r][c] = updated; + for c in (k + 1)..D { + let updated = (-mult).mul_add(rows[k][c], rows[r][c]); + rows[r][c] = updated; + } } } } @@ -116,11 +121,10 @@ impl<const D: usize> Lu<D> { /// Solve `A x = b` using this LU factorization. /// - /// `solve_vec` performs floating-point forward/back substitution and does - /// not provide a certified absolute rounding-error bound for the returned - /// solution. Callers should not expect an absolute error bound like - /// [`Matrix::det_errbound`](crate::Matrix::det_errbound) or the - /// `ERR_COEFF_*` determinant constants provide. + /// [`Vector`] is finite by construction, so this method only checks computed + /// substitution overflows. It performs floating-point forward/back + /// substitution and does not provide a certified absolute rounding-error + /// bound for the returned solution. /// /// # Examples /// ``` @@ -128,10 +132,10 @@ impl<const D: usize> Lu<D> { /// /// # fn main() -> Result<(), LaError> { /// let a = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?; - /// let lu = a.lu(DEFAULT_PIVOT_TOL)?; + /// let lu = a.lu(DEFAULT_SINGULAR_TOL)?; /// /// let b = Vector::<2>::try_new([5.0, 11.0])?; - /// let x = lu.solve_vec(b)?.into_array(); + /// let x = lu.solve(b)?.into_array(); /// /// assert!((x[0] - 1.0).abs() <= 1e-12); /// assert!((x[1] - 2.0).abs() <= 1e-12); @@ -140,18 +144,11 @@ impl<const D: usize> Lu<D> { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] if the right-hand side contains - /// NaN/infinity or a computed substitution intermediate overflows to NaN or - /// infinity. The right-hand side is revalidated at this boundary before - /// entering substitution; public callers normally encounter the same - /// rejection earlier through [`Vector::try_new`](crate::Vector::try_new). + /// Returns [`LaError::NonFinite`] if a computed substitution intermediate + /// overflows to NaN or infinity. #[inline] - pub const fn solve_vec(&self, b: Vector<D>) -> Result<Vector<D>, LaError> { - let b = match FiniteVector::new(b) { - Ok(b) => b, - Err(err) => return Err(err), - }; - match self.solve_finite_vec(b) { + pub const fn solve(&self, b: Vector<D>) -> Result<Vector<D>, LaError> { + match self.solve_finite(FiniteVector::new_unchecked(b)) { Ok(x) => Ok(x.into_vector()), Err(err) => Err(err), } @@ -166,7 +163,7 @@ impl<const D: usize> Lu<D> { /// Returns [`LaError::NonFinite`] if a computed substitution intermediate /// overflows to NaN or infinity. #[inline] - pub(crate) const fn solve_finite_vec( + pub(crate) const fn solve_finite( &self, b: FiniteVector<D>, ) -> Result<FiniteVector<D>, LaError> { @@ -234,7 +231,7 @@ impl<const D: usize> Lu<D> { /// /// # fn main() -> Result<(), LaError> { /// let a = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?; - /// let lu = a.lu(DEFAULT_PIVOT_TOL)?; + /// let lu = a.lu(DEFAULT_SINGULAR_TOL)?; /// /// let det = lu.det()?; /// assert!((det - (-2.0)).abs() <= 1e-12); @@ -264,20 +261,20 @@ impl<const D: usize> Lu<D> { #[cfg(test)] mod tests { use super::*; - use crate::DEFAULT_PIVOT_TOL; + use crate::DEFAULT_SINGULAR_TOL; use core::hint::black_box; use approx::assert_abs_diff_eq; use pastey::paste; - macro_rules! gen_public_api_pivoting_solve_vec_and_det_tests { + macro_rules! gen_public_api_pivoting_solve_and_det_tests { ($d:literal) => { paste! { #[test] - fn [<public_api_lu_solve_vec_pivoting_ $d d>]() { + fn [<public_api_lu_solve_pivoting_ $d d>]() { // Public API path under test: - // Matrix::lu (pub) -> Lu::solve_vec (pub). + // Matrix::lu (pub) -> Lu::solve (pub). // Permutation matrix that swaps the first two basis vectors. // This forces pivoting in column 0 for any D >= 2. @@ -290,7 +287,7 @@ mod tests { let a = Matrix::<$d>::from_rows(black_box(rows)); let lu_fn: fn(Matrix<$d>, Tolerance) -> Result<Lu<$d>, LaError> = black_box(Matrix::<$d>::lu); - let lu = lu_fn(a, DEFAULT_PIVOT_TOL).unwrap(); + let lu = lu_fn(a, DEFAULT_SINGULAR_TOL).unwrap(); // Pick a simple RHS with unique entries, so the expected swap is obvious. let b_arr = { @@ -307,7 +304,7 @@ mod tests { let b = Vector::<$d>::new(black_box(b_arr)); let solve_fn: fn(&Lu<$d>, Vector<$d>) -> Result<Vector<$d>, LaError> = - black_box(Lu::<$d>::solve_vec); + black_box(Lu::<$d>::solve); let x = solve_fn(&lu, b).unwrap().into_array(); for i in 0..$d { @@ -330,7 +327,7 @@ mod tests { let a = Matrix::<$d>::from_rows(black_box(rows)); let lu_fn: fn(Matrix<$d>, Tolerance) -> Result<Lu<$d>, LaError> = black_box(Matrix::<$d>::lu); - let lu = lu_fn(a, DEFAULT_PIVOT_TOL).unwrap(); + let lu = lu_fn(a, DEFAULT_SINGULAR_TOL).unwrap(); // Row swap ⇒ determinant sign flip. let det_fn: fn(&Lu<$d>) -> Result<f64, LaError> = @@ -341,18 +338,18 @@ mod tests { }; } - gen_public_api_pivoting_solve_vec_and_det_tests!(2); - gen_public_api_pivoting_solve_vec_and_det_tests!(3); - gen_public_api_pivoting_solve_vec_and_det_tests!(4); - gen_public_api_pivoting_solve_vec_and_det_tests!(5); + gen_public_api_pivoting_solve_and_det_tests!(2); + gen_public_api_pivoting_solve_and_det_tests!(3); + gen_public_api_pivoting_solve_and_det_tests!(4); + gen_public_api_pivoting_solve_and_det_tests!(5); - macro_rules! gen_public_api_tridiagonal_smoke_solve_vec_and_det_tests { + macro_rules! gen_public_api_tridiagonal_smoke_solve_and_det_tests { ($d:literal) => { paste! { #[test] - fn [<public_api_lu_solve_vec_tridiagonal_smoke_ $d d>]() { + fn [<public_api_lu_solve_tridiagonal_smoke_ $d d>]() { // Public API path under test: - // Matrix::lu (pub) -> Lu::solve_vec (pub). + // Matrix::lu (pub) -> Lu::solve (pub). // Classic SPD tridiagonal: 2 on diagonal, -1 on sub/super-diagonals. #[allow(clippy::large_stack_arrays)] @@ -370,7 +367,7 @@ mod tests { let a = Matrix::<$d>::from_rows(black_box(rows)); let lu_fn: fn(Matrix<$d>, Tolerance) -> Result<Lu<$d>, LaError> = black_box(Matrix::<$d>::lu); - let lu = lu_fn(a, DEFAULT_PIVOT_TOL).unwrap(); + let lu = lu_fn(a, DEFAULT_SINGULAR_TOL).unwrap(); // Choose x = 1, so b = A x is simple: [1, 0, 0, ..., 0, 1]. let mut b_arr = [0.0f64; $d]; @@ -379,7 +376,7 @@ mod tests { let b = Vector::<$d>::new(black_box(b_arr)); let solve_fn: fn(&Lu<$d>, Vector<$d>) -> Result<Vector<$d>, LaError> = - black_box(Lu::<$d>::solve_vec); + black_box(Lu::<$d>::solve); let x = solve_fn(&lu, b).unwrap().into_array(); for &x_i in &x { @@ -409,7 +406,7 @@ mod tests { let a = Matrix::<$d>::from_rows(black_box(rows)); let lu_fn: fn(Matrix<$d>, Tolerance) -> Result<Lu<$d>, LaError> = black_box(Matrix::<$d>::lu); - let lu = lu_fn(a, DEFAULT_PIVOT_TOL).unwrap(); + let lu = lu_fn(a, DEFAULT_SINGULAR_TOL).unwrap(); let det_fn: fn(&Lu<$d>) -> Result<f64, LaError> = black_box(Lu::<$d>::det); @@ -419,18 +416,21 @@ mod tests { }; } - gen_public_api_tridiagonal_smoke_solve_vec_and_det_tests!(16); - gen_public_api_tridiagonal_smoke_solve_vec_and_det_tests!(32); - gen_public_api_tridiagonal_smoke_solve_vec_and_det_tests!(64); + gen_public_api_tridiagonal_smoke_solve_and_det_tests!(16); + gen_public_api_tridiagonal_smoke_solve_and_det_tests!(32); + gen_public_api_tridiagonal_smoke_solve_and_det_tests!(64); #[test] fn solve_1x1() { let a = Matrix::<1>::from_rows(black_box([[2.0]])); - let lu = FiniteMatrix::new(a).unwrap().lu(DEFAULT_PIVOT_TOL).unwrap(); + let lu = FiniteMatrix::new(a) + .unwrap() + .lu(DEFAULT_SINGULAR_TOL) + .unwrap(); let b = Vector::<1>::new(black_box([6.0])); let solve_fn: fn(&Lu<1>, Vector<1>) -> Result<Vector<1>, LaError> = - black_box(Lu::<1>::solve_vec); + black_box(Lu::<1>::solve); let x = solve_fn(&lu, b).unwrap().into_array(); assert_abs_diff_eq!(x[0], 3.0, epsilon = 1e-12); @@ -441,11 +441,14 @@ 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_PIVOT_TOL).unwrap(); + let lu = FiniteMatrix::new(a) + .unwrap() + .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<Vector<2>, LaError> = - black_box(Lu::<2>::solve_vec); + black_box(Lu::<2>::solve); let x = solve_fn(&lu, b).unwrap().into_array(); assert_abs_diff_eq!(x[0], 1.0, epsilon = 1e-12); @@ -455,7 +458,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 lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); + let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); let det_fn: fn(&Lu<2>) -> Result<f64, LaError> = black_box(Lu::<2>::det); assert_abs_diff_eq!(det_fn(&lu).unwrap(), -2.0, epsilon = 1e-12); @@ -465,7 +468,7 @@ mod tests { 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 lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); + let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); let det_fn: fn(&Lu<2>) -> Result<f64, LaError> = black_box(Lu::<2>::det); assert_abs_diff_eq!(det_fn(&lu).unwrap(), -1.0, epsilon = 0.0); @@ -474,11 +477,11 @@ mod tests { #[test] fn solve_requires_pivoting() { let a = Matrix::<2>::from_rows(black_box([[0.0, 1.0], [1.0, 0.0]])); - let lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); + let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<2>::new(black_box([1.0, 2.0])); let solve_fn: fn(&Lu<2>, Vector<2>) -> Result<Vector<2>, LaError> = - black_box(Lu::<2>::solve_vec); + black_box(Lu::<2>::solve); let x = solve_fn(&lu, b).unwrap().into_array(); assert_abs_diff_eq!(x[0], 2.0, epsilon = 1e-12); @@ -488,15 +491,15 @@ mod tests { #[test] fn singular_detected() { let a = Matrix::<2>::from_rows(black_box([[1.0, 2.0], [2.0, 4.0]])); - let err = a.lu(DEFAULT_PIVOT_TOL).unwrap_err(); + let err = a.lu(DEFAULT_SINGULAR_TOL).unwrap_err(); assert_eq!(err, LaError::Singular { pivot_col: 1 }); } #[test] fn singular_due_to_tolerance_at_first_pivot() { - // Not exactly singular, but below DEFAULT_PIVOT_TOL. + // 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 err = a.lu(DEFAULT_PIVOT_TOL).unwrap_err(); + let err = a.lu(DEFAULT_SINGULAR_TOL).unwrap_err(); assert_eq!(err, LaError::Singular { pivot_col: 0 }); } @@ -529,7 +532,7 @@ mod tests { 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 err = a.lu(DEFAULT_PIVOT_TOL).unwrap_err(); + let err = a.lu(DEFAULT_SINGULAR_TOL).unwrap_err(); assert_eq!( err, LaError::NonFinite { @@ -540,29 +543,29 @@ mod tests { } #[test] - fn solve_vec_nonfinite_forward_substitution_overflow() { + 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 lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); + let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<3>::new([1.0e308, 1.0e308, 0.0]); - let err = lu.solve_vec(b).unwrap_err(); + let err = lu.solve(b).unwrap_err(); assert_eq!(err, LaError::NonFinite { row: None, col: 1 }); } #[test] - fn solve_vec_nonfinite_back_substitution_overflow() { + 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 lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); + let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<2>::new([0.0, 1.0e300]); - let err = lu.solve_vec(b).unwrap_err(); + let err = lu.solve(b).unwrap_err(); assert_eq!(err, LaError::NonFinite { row: None, col: 1 }); } #[test] - fn solve_vec_nonfinite_back_substitution_sum_overflow() { + fn solve_nonfinite_back_substitution_sum_overflow() { // Upper-triangular U with a very large off-diagonal in row 1 and a // very large x[2] produced by the RHS. The back-substitution // accumulator `sum = (-row[j]).mul_add(x[j], sum)` overflows while @@ -570,10 +573,10 @@ mod tests { // 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 lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); + let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<3>::new([0.0, 0.0, 1.0e200]); - let err = lu.solve_vec(b).unwrap_err(); + let err = lu.solve(b).unwrap_err(); assert_eq!(err, LaError::NonFinite { row: None, col: 1 }); } @@ -586,17 +589,17 @@ mod tests { [0.0, 0.0, 0.0, 1.0e100, 0.0], [0.0, 0.0, 0.0, 0.0, 1.0e100], ]); - let lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); + let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); assert_eq!(lu.det(), Err(LaError::NonFinite { row: None, col: 3 })); } - macro_rules! gen_solve_vec_boundary_tests { + macro_rules! gen_solve_boundary_tests { ($d:literal) => { paste! { /// Raw non-finite right-hand sides are rejected before a /// public caller can construct a `Vector`. #[test] - fn [<solve_vec_rhs_boundary_rejects_non_finite_ $d d>]() { + fn [<solve_rhs_constructor_rejects_non_finite_ $d d>]() { let mut rhs = [1.0; $d]; rhs[$d - 1] = f64::NAN; @@ -608,35 +611,19 @@ mod tests { }) ); } - - #[test] - fn [<solve_vec_revalidates_unchecked_rhs_storage_ $d d>]() { - let lu = Matrix::<$d>::identity().lu(DEFAULT_PIVOT_TOL).unwrap(); - let mut rhs = [1.0; $d]; - rhs[$d - 1] = f64::NAN; - let rhs = Vector::<$d>::new_unchecked(rhs); - - assert_eq!( - lu.solve_vec(rhs), - Err(LaError::NonFinite { - row: None, - col: $d - 1, - }) - ); - } } }; } - gen_solve_vec_boundary_tests!(2); - gen_solve_vec_boundary_tests!(3); - gen_solve_vec_boundary_tests!(4); - gen_solve_vec_boundary_tests!(5); + gen_solve_boundary_tests!(2); + gen_solve_boundary_tests!(3); + gen_solve_boundary_tests!(4); + gen_solve_boundary_tests!(5); // ----------------------------------------------------------------------- // Const-evaluability tests. // - // These prove that `Lu::det` and `Lu::solve_vec` are truly `const fn` by + // These prove that `Lu::det` and `Lu::solve` are truly `const fn` by // forcing the compiler to evaluate them inside a `const` initializer. // `Lu::factor` is not (yet) `const fn` because it relies on `<[T]>::swap`, // which is not const-stable; we therefore construct `Lu<D>` directly. @@ -646,9 +633,7 @@ mod tests { fn lu_det_const_eval_d2() { const DET: Result<f64, LaError> = { // Triangular factors with diag [2.0, 3.0] and no row swaps. - let mut factors = Matrix::<2>::identity(); - factors.rows[0][0] = 2.0; - factors.rows[1][1] = 3.0; + let factors = Matrix::<2>::from_rows_unchecked([[2.0, 0.0], [0.0, 3.0]]); let lu = Lu::<2> { factors: LuFactors::new_unchecked(factors), piv: [0, 1], @@ -675,8 +660,8 @@ mod tests { } #[test] - fn lu_solve_vec_const_eval_d2() { - // Identity LU ⇒ solve_vec returns the permuted RHS untouched. + fn lu_solve_const_eval_d2() { + // Identity LU ⇒ solve returns the permuted RHS untouched. const X: [f64; 2] = { let lu = Lu::<2> { factors: LuFactors::new_unchecked(Matrix::<2>::identity()), @@ -684,7 +669,7 @@ mod tests { piv_sign: 1.0, }; let b = Vector::<2>::new([1.0, 2.0]); - match lu.solve_vec(b) { + match lu.solve(b) { Ok(v) => v.into_array(), Err(_) => [0.0, 0.0], } diff --git a/src/matrix.rs b/src/matrix.rs index 32965a9..5995b20 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -6,28 +6,45 @@ use crate::ldlt::Ldlt; use crate::lu::Lu; use crate::{ERR_COEFF_2, ERR_COEFF_3, ERR_COEFF_4, LDLT_SYMMETRY_REL_TOL, LaError, Tolerance}; -/// Fixed-size square matrix `D×D`, stored inline. +/// Finite fixed-size square matrix `D×D`, stored inline. /// /// `Matrix` is designed for small, robustness-sensitive systems where stack /// allocation and const-generic dimensions are useful. For large, dynamic, sparse, /// or parallel workloads, prefer a broader linear-algebra crate such as /// [`nalgebra`](https://crates.io/crates/nalgebra) or /// [`faer`](https://crates.io/crates/faer). +/// +/// Public construction and mutation reject NaN and infinity through +/// [`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. +/// +/// Direct field construction is intentionally unavailable to downstream callers: +/// +/// ```compile_fail +/// use la_stack::Matrix; +/// +/// let _ = Matrix::<2> { +/// rows: [[1.0, f64::NAN], [0.0, 1.0]], +/// }; +/// ``` #[must_use] #[derive(Clone, Copy, Debug, PartialEq)] pub struct Matrix<const D: usize> { - pub(crate) rows: [[f64; D]; D], + rows: [[f64; D]; D], } /// Fixed-size square matrix whose stored entries are all finite. /// -/// This proof-bearing wrapper makes the finite invariant explicit for internal -/// algorithms that should not repeatedly check stored entries for NaN or -/// infinity. +/// 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)] -#[allow(clippy::redundant_pub_crate)] -pub(crate) struct FiniteMatrix<const D: usize> { +pub struct FiniteMatrix<const D: usize> { matrix: Matrix<D>, } @@ -56,24 +73,6 @@ impl<const D: usize> FiniteMatrix<D> { } } - /// Validate raw row-major storage and construct a finite matrix. - /// # Errors - /// Returns [`LaError::NonFinite`] with matrix coordinates for the first - /// offending entry in row-major order when `rows` contains NaN or infinity. - #[inline] - pub(crate) const fn from_rows(rows: [[f64; D]; D]) -> Result<Self, LaError> { - match Matrix::try_from_rows(rows) { - Ok(matrix) => Ok(Self::new_unchecked(matrix)), - Err(err) => Err(err), - } - } - - /// All-zeros finite matrix. - #[inline] - pub(crate) const fn zero() -> Self { - Self::new_unchecked(Matrix::zero()) - } - /// Consume the wrapper and return the underlying raw matrix. #[inline] pub(crate) const fn into_matrix(self) -> Matrix<D> { @@ -345,38 +344,6 @@ impl<const D: usize> FiniteMatrix<D> { } } -impl<const D: usize> Default for FiniteMatrix<D> { - #[inline] - fn default() -> Self { - Self::zero() - } -} - -impl<const D: usize> From<FiniteMatrix<D>> for Matrix<D> { - #[inline] - fn from(value: FiniteMatrix<D>) -> Self { - value.into_matrix() - } -} - -impl<const D: usize> TryFrom<Matrix<D>> for FiniteMatrix<D> { - type Error = LaError; - - #[inline] - fn try_from(value: Matrix<D>) -> Result<Self, Self::Error> { - Self::new(value) - } -} - -impl<const D: usize> TryFrom<[[f64; D]; D]> for FiniteMatrix<D> { - type Error = LaError; - - #[inline] - fn try_from(value: [[f64; D]; D]) -> Result<Self, Self::Error> { - Self::from_rows(value) - } -} - /// Matrix proven finite and symmetric under the crate's LDLT symmetry tolerance. #[must_use] #[derive(Clone, Copy, Debug, PartialEq)] @@ -438,9 +405,8 @@ impl<const D: usize> Matrix<D> { /// Try to create a finite matrix from row-major storage. /// - /// This is the public raw-storage boundary for matrices. Public compute - /// methods parse stored rows into crate-internal proof-bearing types before - /// arithmetic, including when crate-internal unchecked storage exists. + /// This is the public raw-storage boundary for matrices. Successful + /// construction makes the returned [`Matrix`] a finite-storage proof. /// /// # Examples /// ``` @@ -467,14 +433,34 @@ impl<const D: usize> Matrix<D> { /// Construct a matrix without checking that entries are finite. /// - /// This crate-internal escape hatch is reserved for literals and algorithm - /// outputs whose finite invariant is visible at the call site. + /// This crate-internal escape hatch is reserved for finite literals and + /// algorithm outputs whose finite invariant is visible at the call site. + /// Computed outputs must be validated before becoming observable API values. #[inline] pub(crate) const fn from_rows_unchecked(rows: [[f64; D]; D]) -> Self { Self { rows } } - /// All-zeros matrix. + /// Borrow finite row-major storage. + /// + /// This accessor exposes the already validated backing array to internal + /// algorithms without giving them mutable access that could invalidate the + /// [`Matrix`] invariant. + #[inline] + pub(crate) const fn rows(&self) -> &[[f64; D]; D] { + &self.rows + } + + /// Mutably borrow raw row-major storage without preserving the finite invariant. + /// + /// This is reserved for internal factorization temporaries whose results are + /// validated before becoming observable API values. + #[inline] + pub(crate) const fn rows_mut_unchecked(&mut self) -> &mut [[f64; D]; D] { + &mut self.rows + } + + /// All-zeros finite matrix. /// /// # Examples /// ``` @@ -488,7 +474,7 @@ impl<const D: usize> Matrix<D> { Self::from_rows_unchecked([[0.0; D]; D]) } - /// Identity matrix. + /// Finite identity matrix. /// /// # Examples /// ``` @@ -512,7 +498,7 @@ impl<const D: usize> Matrix<D> { m } - /// Get an element with bounds checking. + /// Get a finite element with bounds checking. /// /// # Examples /// ``` @@ -535,7 +521,7 @@ impl<const D: usize> Matrix<D> { } } - /// Get an element, preserving index context on failure. + /// Get a finite element, preserving index context on failure. /// /// Prefer [`get`](Self::get) for const or hot paths that only need /// `Option`-style absence. Use this method at public runtime boundaries @@ -646,9 +632,8 @@ impl<const D: usize> Matrix<D> { /// /// # Non-finite handling /// Public constructors and setters reject raw non-finite entries, but - /// crate-internal unchecked storage can still contain NaN or infinity. - /// `inf_norm` returns [`LaError::NonFinite`] if it encounters stored NaN/∞ - /// or if a row sum overflows to a non-finite value. + /// `Matrix` values are finite by construction. `inf_norm` returns + /// [`LaError::NonFinite`] if a row sum overflows to a non-finite value. /// /// Row sums are accumulated in `f64` with ordinary addition. This method /// checks for overflowed accumulators, but it does not provide a certified @@ -675,14 +660,10 @@ impl<const D: usize> Matrix<D> { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] when stored entries are NaN/infinity or a - /// row sum overflows to NaN or infinity. + /// Returns [`LaError::NonFinite`] when a row sum overflows to NaN or infinity. #[inline] pub const fn inf_norm(&self) -> Result<f64, LaError> { - match FiniteMatrix::new(*self) { - Ok(matrix) => matrix.inf_norm(), - Err(err) => Err(err), - } + FiniteMatrix::new_unchecked(*self).inf_norm() } /// Returns `true` if the matrix is symmetric within a relative tolerance. @@ -724,11 +705,11 @@ impl<const D: usize> Matrix<D> { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] when stored entries are NaN/infinity or - /// computing the scaled symmetry tolerance overflows to NaN or infinity. + /// Returns [`LaError::NonFinite`] when computing the scaled symmetry + /// tolerance overflows to NaN or infinity. #[inline] pub fn is_symmetric(&self, rel_tol: Tolerance) -> Result<bool, LaError> { - FiniteMatrix::new(*self)?.is_symmetric(rel_tol) + FiniteMatrix::new_unchecked(*self).is_symmetric(rel_tol) } /// Returns the indices `(r, c)` (with `r < c`) of the first off-diagonal @@ -769,11 +750,11 @@ impl<const D: usize> Matrix<D> { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] when stored entries are NaN/infinity or - /// computing the scaled symmetry tolerance overflows to NaN or infinity. + /// Returns [`LaError::NonFinite`] when computing the scaled symmetry + /// tolerance overflows to NaN or infinity. #[inline] pub fn first_asymmetry(&self, rel_tol: Tolerance) -> Result<Option<(usize, usize)>, LaError> { - FiniteMatrix::new(*self)?.first_asymmetry(rel_tol) + FiniteMatrix::new_unchecked(*self).first_asymmetry(rel_tol) } /// Compute an LU decomposition with partial pivoting. @@ -784,10 +765,10 @@ impl<const D: usize> Matrix<D> { /// /// # fn main() -> Result<(), LaError> { /// let a = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?; - /// let lu = a.lu(DEFAULT_PIVOT_TOL)?; + /// let lu = a.lu(DEFAULT_SINGULAR_TOL)?; /// /// let b = Vector::<2>::try_new([5.0, 11.0])?; - /// let x = lu.solve_vec(b)?.into_array(); + /// let x = lu.solve(b)?.into_array(); /// /// assert!((x[0] - 1.0).abs() <= 1e-12); /// assert!((x[1] - 2.0).abs() <= 1e-12); @@ -804,12 +785,11 @@ impl<const D: usize> Matrix<D> { /// # Errors /// Returns [`LaError::Singular`] if, for some column `k`, the largest-magnitude candidate pivot /// in that column satisfies `|pivot| <= tol` (so no numerically usable pivot exists). - /// Returns [`LaError::NonFinite`] if stored entries are NaN/infinity or an - /// elimination intermediate overflows to NaN/∞ before it can be stored in - /// the returned [`Lu`]. + /// Returns [`LaError::NonFinite`] if an elimination intermediate overflows + /// to NaN/∞ before it can be stored in the returned [`Lu`]. #[inline] pub fn lu(self, tol: Tolerance) -> Result<Lu<D>, LaError> { - FiniteMatrix::new(self)?.lu(tol) + FiniteMatrix::new_unchecked(self).lu(tol) } /// Compute an LDLT factorization (`A = L D Lᵀ`) without pivoting. @@ -846,7 +826,7 @@ impl<const D: usize> Matrix<D> { /// /// // Solve A x = b /// let b = Vector::<2>::try_new([1.0, 2.0])?; - /// let x = ldlt.solve_vec(b)?.into_array(); + /// let x = ldlt.solve(b)?.into_array(); /// assert!((x[0] - (-0.125)).abs() <= 1e-12); /// assert!((x[1] - 0.75).abs() <= 1e-12); /// # Ok(()) @@ -858,12 +838,12 @@ impl<const D: usize> Matrix<D> { /// diagonal entry `d = D[k,k]` is negative. /// Returns [`LaError::Singular`] if `0 <= d <= tol`, treating PSD degeneracy /// as singular/degenerate. - /// Returns [`LaError::NonFinite`] if stored entries are NaN/infinity or - /// factorization computes a non-finite intermediate. + /// Returns [`LaError::NonFinite`] if factorization computes a non-finite + /// intermediate. /// Returns [`LaError::Asymmetric`] if the input matrix is not symmetric. #[inline] pub fn ldlt(self, tol: Tolerance) -> Result<Ldlt<D>, LaError> { - FiniteMatrix::new(self)?.ldlt(tol) + FiniteMatrix::new_unchecked(self).ldlt(tol) } /// Return the first non-finite stored cell in row-major order. @@ -910,14 +890,11 @@ impl<const D: usize> Matrix<D> { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] when stored entries are NaN/infinity or - /// the closed-form determinant overflows to NaN or infinity. + /// Returns [`LaError::NonFinite`] when the closed-form determinant overflows + /// to NaN or infinity. #[inline] pub const fn det_direct(&self) -> Result<Option<f64>, LaError> { - match FiniteMatrix::new(*self) { - Ok(matrix) => matrix.det_direct(), - Err(err) => Err(err), - } + FiniteMatrix::new_unchecked(*self).det_direct() } /// Floating-point determinant, using closed-form formulas for D ≤ 4 and @@ -948,12 +925,11 @@ impl<const D: usize> Matrix<D> { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] if stored entries are NaN/infinity, the - /// LU fallback computes a non-finite factorization cell, or the determinant - /// product overflows to NaN or infinity. + /// Returns [`LaError::NonFinite`] if the LU fallback computes a non-finite + /// factorization cell, or the determinant product overflows to NaN or infinity. #[inline] pub fn det(self) -> Result<f64, LaError> { - FiniteMatrix::new(self)?.det() + FiniteMatrix::new_unchecked(self).det() } /// Conservative absolute error bound for `det_direct()`. @@ -1017,14 +993,11 @@ impl<const D: usize> Matrix<D> { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] when stored entries are NaN/infinity or - /// the bound computation overflows to NaN or infinity. + /// Returns [`LaError::NonFinite`] when the bound computation overflows to + /// NaN or infinity. #[inline] pub const fn det_errbound(&self) -> Result<Option<f64>, LaError> { - match FiniteMatrix::new(*self) { - Ok(matrix) => matrix.det_errbound(), - Err(err) => Err(err), - } + FiniteMatrix::new_unchecked(*self).det_errbound() } } @@ -1038,7 +1011,7 @@ impl<const D: usize> Default for Matrix<D> { #[cfg(test)] mod tests { use super::*; - use crate::DEFAULT_PIVOT_TOL; + use crate::DEFAULT_SINGULAR_TOL; use crate::vector::{FiniteVector, Vector}; use approx::assert_abs_diff_eq; @@ -1157,7 +1130,7 @@ mod tests { } #[test] - fn [<public_api_matrix_identity_lu_det_solve_vec_ $d d>]() { + fn [<public_api_matrix_identity_lu_det_solve_ $d d>]() { let m = Matrix::<$d>::identity(); // Identity has ones on diag and zeros off diag. @@ -1173,7 +1146,7 @@ mod tests { assert_abs_diff_eq!(det, 1.0, epsilon = 1e-12); // LU solve on identity returns the RHS. - let lu = m.lu(DEFAULT_PIVOT_TOL).unwrap(); + let lu = m.lu(DEFAULT_SINGULAR_TOL).unwrap(); let b_arr = { let mut arr = [0.0f64; $d]; @@ -1185,7 +1158,7 @@ mod tests { }; let b = Vector::<$d>::new(b_arr); - let x = lu.solve_vec(b).unwrap().into_array(); + let x = lu.solve(b).unwrap().into_array(); for (x_i, b_i) in x.iter().zip(b_arr.iter()) { assert_abs_diff_eq!(*x_i, *b_i, epsilon = 1e-12); @@ -1199,11 +1172,9 @@ mod tests { rows[r][r] = 1.0; } - let finite = FiniteMatrix::<$d>::from_rows(rows).unwrap(); + let finite = FiniteMatrix::<$d>::new(Matrix::<$d>::from_rows(rows)).unwrap(); assert_matrix_abs_eq(&finite.into_matrix(), &Matrix::<$d>::from_rows(rows)); - assert_matrix_abs_eq(&FiniteMatrix::<$d>::zero().into_matrix(), &Matrix::<$d>::zero()); - assert_matrix_abs_eq(&FiniteMatrix::<$d>::default().into_matrix(), &Matrix::<$d>::zero()); assert_eq!(finite.into_matrix().get(0, 0), Some(1.0)); assert_eq!(finite.into_matrix().get($d, 0), None); assert_eq!( @@ -1214,44 +1185,19 @@ mod tests { dim: $d, }) ); - assert_matrix_abs_eq(&Matrix::from(finite), &Matrix::<$d>::from_rows(rows)); - assert_matrix_abs_eq( - &FiniteMatrix::<$d>::try_from(Matrix::<$d>::from_rows(rows)) - .unwrap() - .into_matrix(), - &finite.into_matrix() - ); - assert_matrix_abs_eq( - &FiniteMatrix::<$d>::try_from(rows).unwrap().into_matrix(), - &finite.into_matrix() - ); } #[test] fn [<finite_matrix_rejects_nonfinite_with_coordinates_ $d d>]() { let mut rows = [[1.0f64; $d]; $d]; rows[$d - 1][0] = f64::NAN; - - assert_eq!( - FiniteMatrix::<$d>::from_rows(rows), - Err(LaError::NonFinite { - row: Some($d - 1), - col: 0, - }) - ); - } - - #[test] - fn [<finite_matrix_try_from_raw_matrix_revalidates_entries_ $d d>]() { - let mut rows = [[1.0f64; $d]; $d]; - rows[$d - 1][$d - 1] = f64::INFINITY; let raw = Matrix::<$d>::from_rows_unchecked(rows); assert_eq!( - FiniteMatrix::<$d>::try_from(raw), + FiniteMatrix::<$d>::new(raw), Err(LaError::NonFinite { row: Some($d - 1), - col: $d - 1, + col: 0, }) ); } @@ -1265,7 +1211,7 @@ mod tests { } let raw = Matrix::<$d>::from_rows(rows); - let finite = FiniteMatrix::<$d>::from_rows(rows).unwrap(); + 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()); @@ -1290,18 +1236,18 @@ mod tests { for (dst, src) in arr.iter_mut().zip(values.iter()) { *dst = *src; } - FiniteVector::<$d>::from_array(arr).unwrap() + FiniteVector::<$d>::new_unchecked(Vector::<$d>::new(arr)) }; - let lu = finite.lu(DEFAULT_PIVOT_TOL).unwrap(); + let lu = finite.lu(DEFAULT_SINGULAR_TOL).unwrap(); assert_array_abs_eq( - &lu.solve_finite_vec(rhs).unwrap().into_array(), + &lu.solve_finite(rhs).unwrap().into_array(), &rhs.into_array() ); - let ldlt = finite.ldlt(DEFAULT_PIVOT_TOL).unwrap(); + let ldlt = finite.ldlt(DEFAULT_SINGULAR_TOL).unwrap(); assert_array_abs_eq( - &ldlt.solve_finite_vec(rhs).unwrap().into_array(), + &ldlt.solve_finite(rhs).unwrap().into_array(), &rhs.into_array() ); } @@ -1315,61 +1261,6 @@ mod tests { gen_public_api_matrix_tests!(4); gen_public_api_matrix_tests!(5); - macro_rules! gen_public_matrix_revalidation_tests { - ($d:literal) => { - paste! { - #[test] - fn [<public_matrix_compute_methods_revalidate_unchecked_storage_ $d d>]() { - let mut rows = [[0.0f64; $d]; $d]; - for i in 0..$d { - rows[i][i] = 1.0; - } - rows[0][$d - 1] = f64::NAN; - let raw = Matrix::<$d>::from_rows_unchecked(rows); - let expected = LaError::NonFinite { - row: Some(0), - col: $d - 1, - }; - - assert_eq!(raw.inf_norm(), Err(expected)); - assert_eq!( - raw.is_symmetric(Tolerance::new(0.0).unwrap()), - Err(expected) - ); - assert_eq!( - raw.first_asymmetry(Tolerance::new(0.0).unwrap()), - Err(expected) - ); - assert_eq!(raw.lu(DEFAULT_PIVOT_TOL), Err(expected)); - assert_eq!(raw.ldlt(DEFAULT_PIVOT_TOL), Err(expected)); - assert_eq!(raw.det_direct(), Err(expected)); - assert_eq!(raw.det(), Err(expected)); - assert_eq!(raw.det_errbound(), Err(expected)); - } - } - }; - } - - gen_public_matrix_revalidation_tests!(2); - gen_public_matrix_revalidation_tests!(3); - gen_public_matrix_revalidation_tests!(4); - gen_public_matrix_revalidation_tests!(5); - - #[test] - fn public_matrix_fast_none_paths_revalidate_unchecked_storage() { - let mut rows = [[1.0f64; 5]; 5]; - rows[4][1] = f64::INFINITY; - let raw = Matrix::<5>::from_rows_unchecked(rows); - - let expected = Err(LaError::NonFinite { - row: Some(4), - col: 1, - }); - - assert_eq!(raw.det_direct(), expected); - assert_eq!(raw.det_errbound(), expected); - } - // === det_direct tests === #[test] @@ -1523,7 +1414,7 @@ mod tests { } let m = Matrix::<$d>::from_rows(rows); let direct = m.det_direct().unwrap().unwrap(); - let lu_det = m.lu(DEFAULT_PIVOT_TOL).unwrap().det().unwrap(); + let lu_det = m.lu(DEFAULT_SINGULAR_TOL).unwrap().det().unwrap(); let eps = lu_det.abs().mul_add(1e-12, 1e-12); assert_abs_diff_eq!(direct, lu_det, epsilon = eps); } @@ -1614,7 +1505,7 @@ mod tests { assert_abs_diff_eq!(m.det().unwrap(), 1e-13, epsilon = 0.0); assert_eq!( - m.lu(DEFAULT_PIVOT_TOL), + m.lu(DEFAULT_SINGULAR_TOL), Err(LaError::Singular { pivot_col: 0 }) ); } @@ -1966,7 +1857,7 @@ mod tests { paste! { #[test] fn [<matrix_ldlt_accepts_identity_ $d d>]() { - let ldlt = Matrix::<$d>::identity().ldlt(DEFAULT_PIVOT_TOL).unwrap(); + let ldlt = Matrix::<$d>::identity().ldlt(DEFAULT_SINGULAR_TOL).unwrap(); assert_abs_diff_eq!(ldlt.det().unwrap(), 1.0, epsilon = 0.0); } diff --git a/src/vector.rs b/src/vector.rs index 8314b91..e9f0e7b 100644 --- a/src/vector.rs +++ b/src/vector.rs @@ -1,23 +1,40 @@ //! Fixed-size, stack-allocated vectors. +use core::hint::cold_path; + use crate::LaError; -/// Fixed-size vector of length `D`, stored inline. +/// Finite fixed-size vector of length `D`, stored inline. +/// +/// 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. +/// +/// Direct field construction is intentionally unavailable to downstream callers: +/// +/// ```compile_fail +/// use la_stack::Vector; +/// +/// let _ = Vector::<2> { +/// data: [1.0, f64::NAN], +/// }; +/// ``` #[must_use] #[derive(Clone, Copy, Debug, PartialEq)] pub struct Vector<const D: usize> { - pub(crate) data: [f64; D], + data: [f64; D], } /// Fixed-size vector whose stored entries are all finite. /// -/// This proof-bearing wrapper makes the finite invariant explicit for internal -/// algorithms that should not rediscover that no stored entry is NaN or -/// infinite. +/// 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)] -#[allow(clippy::redundant_pub_crate)] -pub(crate) struct FiniteVector<const D: usize> { +pub struct FiniteVector<const D: usize> { vector: Vector<D>, } @@ -31,39 +48,6 @@ impl<const D: usize> FiniteVector<D> { Self { vector } } - /// Validate a vector and wrap it for algorithms that carry the finite - /// invariant explicitly. - /// - /// # Errors - /// Returns [`LaError::NonFinite`] with `row: None` and the first offending - /// entry index when `vector` contains NaN or infinity. - #[inline] - pub(crate) const fn new(vector: Vector<D>) -> Result<Self, LaError> { - if let Some(index) = Vector::<D>::first_non_finite_entry_in(&vector.data) { - Err(LaError::non_finite_at(index)) - } else { - Ok(Self::new_unchecked(vector)) - } - } - - /// Validate raw vector storage and construct a finite vector. - /// # Errors - /// Returns [`LaError::NonFinite`] with `row: None` and the first offending - /// entry index when `data` contains NaN or infinity. - #[inline] - pub(crate) const fn from_array(data: [f64; D]) -> Result<Self, LaError> { - match Vector::try_new(data) { - Ok(vector) => Ok(Self::new_unchecked(vector)), - Err(err) => Err(err), - } - } - - /// All-zeros finite vector. - #[inline] - pub(crate) const fn zero() -> Self { - Self::new_unchecked(Vector::zero()) - } - /// Consume the wrapper and return the underlying raw vector. #[inline] pub(crate) const fn into_vector(self) -> Vector<D> { @@ -72,6 +56,7 @@ impl<const D: usize> FiniteVector<D> { /// Borrow the backing array without revalidating stored entries. #[inline] + #[must_use] pub(crate) const fn as_array(&self) -> &[f64; D] { self.vector.as_array() } @@ -100,6 +85,7 @@ impl<const D: usize> FiniteVector<D> { 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; @@ -123,6 +109,7 @@ impl<const D: usize> FiniteVector<D> { 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; @@ -131,38 +118,6 @@ impl<const D: usize> FiniteVector<D> { } } -impl<const D: usize> Default for FiniteVector<D> { - #[inline] - fn default() -> Self { - Self::zero() - } -} - -impl<const D: usize> From<FiniteVector<D>> for Vector<D> { - #[inline] - fn from(value: FiniteVector<D>) -> Self { - value.into_vector() - } -} - -impl<const D: usize> TryFrom<Vector<D>> for FiniteVector<D> { - type Error = LaError; - - #[inline] - fn try_from(value: Vector<D>) -> Result<Self, Self::Error> { - Self::new(value) - } -} - -impl<const D: usize> TryFrom<[f64; D]> for FiniteVector<D> { - type Error = LaError; - - #[inline] - fn try_from(value: [f64; D]) -> Result<Self, Self::Error> { - Self::from_array(value) - } -} - impl<const D: usize> Vector<D> { /// Test-only infallible constructor for finite literal fixtures. #[cfg(test)] @@ -176,9 +131,8 @@ impl<const D: usize> Vector<D> { /// Try to create a finite vector from a backing array. /// - /// This is the public raw-storage boundary for vectors. Public compute - /// methods parse stored entries into crate-internal proof-bearing types - /// before arithmetic, including when crate-internal unchecked storage exists. + /// This is the public raw-storage boundary for vectors. Successful + /// construction makes the returned [`Vector`] a finite-storage proof. /// /// # Examples /// ``` @@ -205,8 +159,10 @@ impl<const D: usize> Vector<D> { /// Construct a vector without checking that entries are finite. /// - /// This crate-internal escape hatch is reserved for literals and algorithm - /// outputs whose finite invariant is visible at the call site. + /// This crate-internal escape hatch is reserved for finite literals and + /// algorithm outputs whose finite invariant is visible at the call site. + /// Computed outputs must check their intermediates before using this + /// constructor to create an observable [`Vector`]. #[inline] pub(crate) const fn new_unchecked(data: [f64; D]) -> Self { Self { data } @@ -228,7 +184,7 @@ impl<const D: usize> Vector<D> { None } - /// All-zeros vector. + /// All-zeros finite vector. /// /// # Examples /// ``` @@ -242,7 +198,7 @@ impl<const D: usize> Vector<D> { Self::new_unchecked([0.0; D]) } - /// Borrow the backing array. + /// Borrow the finite backing array. /// /// # Examples /// ``` @@ -260,7 +216,7 @@ impl<const D: usize> Vector<D> { &self.data } - /// Consume and return the backing array. + /// Consume and return the finite backing array. /// /// # Examples /// ``` @@ -284,9 +240,8 @@ impl<const D: usize> Vector<D> { /// Terms are accumulated in `f64` using [`f64::mul_add`] at each index. /// Intermediate rounding occurs, and this method does not provide a /// certified absolute rounding bound for the returned dot product. Raw - /// storage is parsed into the crate-internal finite-vector proof before - /// accumulation, so internally unchecked vectors are rejected before - /// arithmetic starts. + /// `Vector` values are finite by construction, so this method only checks + /// whether the accumulation overflows to NaN or infinity. /// /// # Examples /// ``` @@ -301,19 +256,11 @@ impl<const D: usize> Vector<D> { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] when either input contains NaN/infinity or - /// the accumulated dot product overflows to NaN or infinity. + /// Returns [`LaError::NonFinite`] when the accumulated dot product overflows + /// to NaN or infinity. #[inline] pub const fn dot(self, other: Self) -> Result<f64, LaError> { - let lhs = match FiniteVector::new(self) { - Ok(lhs) => lhs, - Err(err) => return Err(err), - }; - let rhs = match FiniteVector::new(other) { - Ok(rhs) => rhs, - Err(err) => return Err(err), - }; - lhs.dot(rhs) + FiniteVector::new_unchecked(self).dot(FiniteVector::new_unchecked(other)) } /// Squared Euclidean norm. @@ -321,9 +268,9 @@ impl<const D: usize> Vector<D> { /// This is computed as `dot(self, self)`, so `norm2_sq` has the same /// `f64` [`mul_add`](f64::mul_add) accumulation behavior as [`dot`](Self::dot). /// Intermediate rounding occurs, and this method does not provide a - /// certified absolute rounding bound for the returned squared norm. Raw - /// storage is parsed into the crate-internal finite-vector proof before - /// accumulation. + /// certified absolute rounding bound for the returned squared norm. + /// `Vector` values are finite by construction, so this method only checks + /// whether the accumulation overflows to NaN or infinity. /// /// # Examples /// ``` @@ -337,14 +284,11 @@ impl<const D: usize> Vector<D> { /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] when stored entries are NaN/infinity or - /// the accumulated norm overflows to NaN or infinity. + /// Returns [`LaError::NonFinite`] when the accumulated norm overflows to NaN + /// or infinity. #[inline] pub const fn norm2_sq(self) -> Result<f64, LaError> { - match FiniteVector::new(self) { - Ok(vector) => vector.norm2_sq(), - Err(err) => Err(err), - } + FiniteVector::new_unchecked(self).norm2_sq() } } @@ -514,38 +458,6 @@ mod tests { assert_eq!(a.norm2_sq(), Err(LaError::NonFinite { row: None, col: 0 })); } - #[test] - fn [<public_api_vector_methods_revalidate_unchecked_storage_ $d d>]() { - let mut lhs_arr = [1.0f64; $d]; - lhs_arr[$d - 1] = f64::NAN; - let lhs = Vector::<$d>::new_unchecked(lhs_arr); - let valid = Vector::<$d>::new([1.0; $d]); - - assert_eq!( - lhs.dot(valid), - Err(LaError::NonFinite { - row: None, - col: $d - 1, - }) - ); - assert_eq!( - lhs.norm2_sq(), - Err(LaError::NonFinite { - row: None, - col: $d - 1, - }) - ); - - let mut rhs_arr = [1.0f64; $d]; - rhs_arr[0] = f64::INFINITY; - let rhs = Vector::<$d>::new_unchecked(rhs_arr); - - assert_eq!( - valid.dot(rhs), - Err(LaError::NonFinite { row: None, col: 0 }) - ); - } - #[test] fn [<finite_vector_accepts_and_roundtrips_ $d d>]() { let arr = { @@ -557,47 +469,10 @@ mod tests { arr }; - let finite = FiniteVector::<$d>::from_array(arr).unwrap(); + 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); - assert_array_abs_eq(&FiniteVector::<$d>::zero().into_array(), &[0.0; $d]); - assert_array_abs_eq(&FiniteVector::<$d>::default().into_array(), &[0.0; $d]); - assert_array_abs_eq(Vector::from(finite).as_array(), &arr); - assert_array_abs_eq( - &FiniteVector::<$d>::try_from(Vector::<$d>::new(arr)).unwrap().into_array(), - &arr - ); - assert_array_abs_eq(&FiniteVector::<$d>::try_from(arr).unwrap().into_array(), &arr); - } - - #[test] - fn [<finite_vector_rejects_nonfinite_with_index_ $d d>]() { - let mut arr = [1.0f64; $d]; - arr[$d - 1] = f64::INFINITY; - - assert_eq!( - FiniteVector::<$d>::from_array(arr), - Err(LaError::NonFinite { - row: None, - col: $d - 1, - }) - ); - } - - #[test] - fn [<finite_vector_try_from_raw_vector_revalidates_entries_ $d d>]() { - let mut arr = [1.0f64; $d]; - arr[$d - 1] = f64::NAN; - let raw = Vector::<$d>::new_unchecked(arr); - - assert_eq!( - FiniteVector::<$d>::try_from(raw), - Err(LaError::NonFinite { - row: None, - col: $d - 1, - }) - ); } } }; diff --git a/tests/proptest_factorizations.rs b/tests/proptest_factorizations.rs index fb8498e..a2d8fc8 100644 --- a/tests/proptest_factorizations.rs +++ b/tests/proptest_factorizations.rs @@ -92,7 +92,7 @@ macro_rules! gen_factorization_proptests { assert_abs_diff_eq!(ldlt.det().unwrap(), expected_det, epsilon = 1e-8); let b = Vector::<$d>::try_new(b_arr).unwrap(); - let x = ldlt.solve_vec(b).unwrap().into_array(); + let x = ldlt.solve(b).unwrap().into_array(); for i in 0..$d { assert_abs_diff_eq!(x[i], x_true[i], epsilon = 1e-8); } @@ -167,12 +167,12 @@ macro_rules! gen_factorization_proptests { } let a = Matrix::<$d>::try_from_rows(a_rows).unwrap(); - let lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); + let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); assert_abs_diff_eq!(lu.det().unwrap(), expected_det, epsilon = 1e-8); let b = Vector::<$d>::try_new(b_arr).unwrap(); - let x = lu.solve_vec(b).unwrap().into_array(); + let x = lu.solve(b).unwrap().into_array(); for i in 0..$d { assert_abs_diff_eq!(x[i], x_true[i], epsilon = 1e-8); } @@ -253,12 +253,12 @@ macro_rules! gen_factorization_proptests { } let a = Matrix::<$d>::try_from_rows(a_rows).unwrap(); - let lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); + let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); assert_abs_diff_eq!(lu.det().unwrap(), expected_det, epsilon = 1e-8); let b = Vector::<$d>::try_new(b_arr).unwrap(); - let x = lu.solve_vec(b).unwrap().into_array(); + let x = lu.solve(b).unwrap().into_array(); for i in 0..$d { assert_abs_diff_eq!(x[i], x_true[i], epsilon = 1e-8); } diff --git a/tests/proptest_matrix.rs b/tests/proptest_matrix.rs index dd685bb..bd79e3a 100644 --- a/tests/proptest_matrix.rs +++ b/tests/proptest_matrix.rs @@ -20,7 +20,7 @@ fn small_ldlt_l_entry() -> impl Strategy<Value = f64> { } fn positive_ldlt_diag() -> impl Strategy<Value = f64> { - // Positive diagonal, comfortably above DEFAULT_PIVOT_TOL. + // Positive diagonal, comfortably above DEFAULT_SINGULAR_TOL. (1i16..=20i16).prop_map(|x| f64::from(x) / 10.0) } @@ -136,7 +136,7 @@ macro_rules! gen_public_api_matrix_proptests { } #[test] - fn [<matrix_det_and_solve_vec_for_diagonal_ $d d>]( + fn [<matrix_det_and_solve_for_diagonal_ $d d>]( diag in proptest::array::[<uniform $d>](small_nonzero_f64()), b_arr in proptest::array::[<uniform $d>](small_f64()), ) { @@ -161,9 +161,9 @@ macro_rules! gen_public_api_matrix_proptests { let eps = expected_det.abs().mul_add(1e-12, 1e-12); assert_abs_diff_eq!(det, expected_det, epsilon = eps); - let lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); + let lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap(); let b = Vector::<$d>::try_new(b_arr).unwrap(); - let x = lu.solve_vec(b).unwrap().into_array(); + let x = lu.solve(b).unwrap().into_array(); for i in 0..$d { let expected_x = b_arr[i] / diag[i]; @@ -219,11 +219,11 @@ macro_rules! gen_public_api_matrix_proptests { let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap(); let det_ldlt = ldlt.det().unwrap(); - let det_lu = a.lu(DEFAULT_PIVOT_TOL).unwrap().det().unwrap(); + let det_lu = a.lu(DEFAULT_SINGULAR_TOL).unwrap().det().unwrap(); assert_abs_diff_eq!(det_ldlt, det_lu, epsilon = 1e-8); let b = Vector::<$d>::try_new(b_arr).unwrap(); - let x = ldlt.solve_vec(b).unwrap().into_array(); + let x = ldlt.solve(b).unwrap().into_array(); for i in 0..$d { assert_abs_diff_eq!(x[i], x_true[i], epsilon = 1e-8); diff --git a/tests/semgrep/src/project_rules/finite_api_contract.rs b/tests/semgrep/src/project_rules/finite_api_contract.rs new file mode 100644 index 0000000..fa8cf04 --- /dev/null +++ b/tests/semgrep/src/project_rules/finite_api_contract.rs @@ -0,0 +1,88 @@ +#![forbid(unsafe_code)] +#![allow(dead_code)] + +// ruleid: la-stack.rust.no-public-raw-linear-algebra-modules +pub mod matrix; + +// ruleid: la-stack.rust.no-public-raw-linear-algebra-modules +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<const D: usize> { + // ruleid: la-stack.rust.no-public-matrix-vector-storage-fields + pub rows: [[f64; D]; D], +} + +#[must_use] +pub struct Vector<const D: usize> { + // ruleid: la-stack.rust.no-public-matrix-vector-storage-fields + pub(crate) data: [f64; D], +} + +pub struct CleanMatrix<const D: usize> { + // ok: la-stack.rust.no-public-matrix-vector-storage-fields + rows: [[f64; D]; D], +} + +impl<const D: usize> Matrix<D> { + // ruleid: la-stack.rust.no-public-unchecked-finite-constructors + pub const fn from_rows_unchecked(rows: [[f64; D]; D]) -> Self { + Self { rows } + } + + // ok: la-stack.rust.no-public-unchecked-finite-constructors + pub(crate) const fn from_rows_unchecked_internal(rows: [[f64; D]; D]) -> Self { + Self { rows } + } +} + +impl<const D: usize> Vector<D> { + // ruleid: la-stack.rust.no-public-unchecked-finite-constructors + pub const fn new_unchecked(data: [f64; D]) -> Self { + Self { data } + } + + // ok: la-stack.rust.no-public-unchecked-finite-constructors + pub(crate) const fn new_unchecked_internal(data: [f64; D]) -> Self { + Self { data } + } +} + +pub struct Lu<const D: usize>; + +impl<const D: usize> Lu<D> { + // ruleid: la-stack.rust.no-legacy-solve-vec-api + pub fn solve_vec(&self, rhs: Vector<D>) -> Vector<D> { + rhs + } + + // ok: la-stack.rust.no-legacy-solve-vec-api + pub fn solve(&self, rhs: Vector<D>) -> Vector<D> { + rhs + } +} + +fn call_sites<const D: usize>(lu: &Lu<D>, rhs: Vector<D>) { + // ruleid: la-stack.rust.no-legacy-solve-vec-api + let _ = lu.solve_vec(rhs); +} diff --git a/tests/semgrep/src/project_rules/readme_doctest_mirrors.rs b/tests/semgrep/src/project_rules/readme_doctest_mirrors.rs new file mode 100644 index 0000000..2a02e04 --- /dev/null +++ b/tests/semgrep/src/project_rules/readme_doctest_mirrors.rs @@ -0,0 +1,35 @@ +#![forbid(unsafe_code)] +#![allow(dead_code)] + +// ruleid: la-stack.rust.no-unwrap-expect-in-readme-doctest-mirrors +mod readme_doctests { + #[test] + fn readme_mirror_uses_unwrap() { + let _ = Some(1_u32).unwrap(); + } +} + +// ruleid: la-stack.rust.no-unwrap-expect-in-readme-doctest-mirrors +mod readme_doctests { + #[test] + fn readme_mirror_uses_expect() { + let _ = Ok::<u32, &'static str>(1).expect("README mirrors should use ?"); + } +} + +// ok: la-stack.rust.no-unwrap-expect-in-readme-doctest-mirrors +mod tests { + #[test] + fn ordinary_internal_tests_may_use_unwrap() { + let _ = Some(1_u32).unwrap(); + } +} + +// ok: la-stack.rust.no-unwrap-expect-in-readme-doctest-mirrors +mod readme_doctests { + #[test] + fn readme_mirror_uses_result() -> Result<(), &'static str> { + let _ = Ok::<u32, &'static str>(1)?; + Ok(()) + } +} diff --git a/tests/vs_linalg_inputs.rs b/tests/vs_linalg_inputs.rs index b2748c5..405215e 100644 --- a/tests/vs_linalg_inputs.rs +++ b/tests/vs_linalg_inputs.rs @@ -6,7 +6,7 @@ use faer::linalg::solvers::Solve; use faer::{Mat, Side}; use nalgebra::{SMatrix, SVector}; -use la_stack::{DEFAULT_PIVOT_TOL, DEFAULT_SINGULAR_TOL, Matrix, Vector}; +use la_stack::{DEFAULT_SINGULAR_TOL, Matrix, Vector}; #[path = "../benches/common/vs_linalg.rs"] mod vs_linalg; @@ -75,7 +75,7 @@ macro_rules! gen_smoke_test { let fv2 = Mat::<f64>::from_fn($d, 1, |i, _| vector_entry(i, 1.0)); let la_lu = a - .lu(DEFAULT_PIVOT_TOL) + .lu(DEFAULT_SINGULAR_TOL) .unwrap_or_else(|err| panic!("la_stack LU factorization failed: {err}")); let na_lu = na.lu(); let fa_lu = fa.partial_piv_lu(); @@ -91,7 +91,7 @@ macro_rules! gen_smoke_test { ); let la_lu_x = la_lu - .solve_vec(rhs) + .solve(rhs) .unwrap_or_else(|err| panic!("la_stack LU solve failed: {err}")); let na_lu_x = na_lu .solve(&nrhs) @@ -134,7 +134,7 @@ macro_rules! gen_smoke_test { ); let la_ldlt_x = la_ldlt - .solve_vec(rhs) + .solve(rhs) .unwrap_or_else(|err| panic!("la_stack LDLT solve failed: {err}")); let na_cholesky_x = na_cholesky.solve(&nrhs); let fa_ldlt_x = fa_ldlt.solve(&frhs);