diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 7f01214..f799acc 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -232,3 +232,4 @@ jobs: cargo test --lib --verbose cargo test --doc --verbose cargo test --tests --verbose + cargo test --features exact --verbose diff --git a/.github/workflows/codacy.yml b/.github/workflows/codacy.yml index dfac01c..d8a157a 100644 --- a/.github/workflows/codacy.yml +++ b/.github/workflows/codacy.yml @@ -41,6 +41,8 @@ jobs: # only required for a private repository by # github/codeql-action/upload-sarif to get the Action run status actions: read + env: + CODACY_PROJECT_TOKEN: ${{ secrets.CODACY_PROJECT_TOKEN }} name: Codacy Security Scan runs-on: ubuntu-latest timeout-minutes: 30 @@ -77,13 +79,13 @@ jobs: # Execute Codacy Analysis CLI and generate a SARIF output with # the security issues identified during the analysis - name: Run Codacy Analysis CLI + if: ${{ env.CODACY_PROJECT_TOKEN != '' }} + id: codacy_analysis uses: codacy/codacy-analysis-cli-action@562ee3e92b8e92df8b67e0a5ff8aa8e261919c08 with: # Check https://github.com/codacy/codacy-analysis-cli#project-token # to get your project token from your Codacy repository. - # You can also omit the token and run the tools that support - # default configurations - project-token: ${{ secrets.CODACY_PROJECT_TOKEN }} + project-token: ${{ env.CODACY_PROJECT_TOKEN }} verbose: true directory: ${{ env.CODACY_WORKDIR }} output: ${{ env.CODACY_SARIF }} @@ -94,9 +96,24 @@ jobs: # Force 0 exit code to allow SARIF file generation # This will handover control about PR rejection to the GitHub side max-allowed-issues: 2147483647 + # Codacy can fail transiently on PRs (e.g. remote config/tools service outages). + # Keep PR checks non-blocking and continue to SARIF fallback/upload. + continue-on-error: ${{ github.event_name == 'pull_request' }} - # Process SARIF file to split by tool - - name: Split SARIF by tool + - name: Warn when Codacy token is unavailable on PR + if: ${{ github.event_name == 'pull_request' && env.CODACY_PROJECT_TOKEN == '' }} + run: | + echo "::warning::CODACY_PROJECT_TOKEN is unavailable for this pull_request." + echo "::warning::Skipping Codacy Analysis CLI and using SARIF fallback." + + - name: Warn when Codacy analysis fails on PR + if: ${{ always() && github.event_name == 'pull_request' && steps.codacy_analysis.outcome == 'failure' }} + run: | + echo "::warning::Codacy Analysis CLI failed on this pull_request run; continuing with SARIF fallback." + + # Validate SARIF output or create an empty fallback for upload + - name: Validate or create SARIF + if: always() run: | # Fail fast and surface errors clearly set -euo pipefail @@ -124,6 +141,7 @@ jobs: # Select SARIF file for upload - name: Select SARIF file for upload + if: always() run: | set -euo pipefail # Honor preselected SARIF_FILE from earlier steps (e.g., empty SARIF case) diff --git a/AGENTS.md b/AGENTS.md index 8f77f26..971b249 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -24,10 +24,12 @@ When making changes in this repo, prioritize (in order): - Format: `cargo fmt` (or `just fmt`) - Integration tests: `just test-integration` - Lint (Clippy): `cargo clippy --all-targets --all-features -- -D warnings` (or `just clippy`) +- Lint (Clippy, exact feature): `cargo clippy --features exact --all-targets -- -D warnings` (or `just clippy-exact`) - Lint/validate: `just check` - Pre-commit validation: `just ci` - Python tests: `just test-python` - Run a single test (by name filter): `cargo test solve_2x2_basic` (or the full path: `cargo test lu::tests::solve_2x2_basic`) +- Run exact-feature tests: `cargo test --features exact --verbose` (or `just test-exact`) - Run examples: `just examples` (or `cargo run --example det_5x5` / `cargo run --example solve_5x5` / `cargo run --example const_det_4x4`) - Spell check: `just spell-check` (uses `typos.toml` at repo root; add false positives to `[default.extend-words]`) @@ -40,6 +42,8 @@ When making changes in this repo, prioritize (in order): - `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/exact.rs`: `det_sign_exact()` — adaptive-precision determinant sign + (Shewchuk-style f64 filter + Bareiss in `BigRational`); `features = ["exact"]` - A minimal `justfile` exists for common workflows (see `just --list`). - The public API re-exports these items from `src/lib.rs`. - Dev-only benchmarks live in `benches/vs_linalg.rs` (Criterion + nalgebra/faer comparison). diff --git a/Cargo.lock b/Cargo.lock index 9950998..5f9ba00 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -675,12 +675,14 @@ dependencies = [ [[package]] name = "la-stack" -version = "0.1.3" +version = "0.2.0" dependencies = [ "approx", "criterion", "faer", "nalgebra", + "num-bigint", + "num-rational", "pastey", "proptest", ] diff --git a/Cargo.toml b/Cargo.toml index 4c2af4a..0df0e97 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "la-stack" -version = "0.1.3" +version = "0.2.0" edition = "2024" rust-version = "1.94" license = "BSD-3-Clause" @@ -13,10 +13,12 @@ categories = [ "mathematics", "science" ] keywords = [ "linear-algebra", "geometry", "const-generics" ] [dependencies] -# Intentionally empty (bench-only deps are optional features below) +# All runtime deps are optional; see [features] below. criterion = { version = "0.8.2", features = [ "html_reports" ], optional = true } faer = { version = "0.24.0", default-features = false, features = [ "std", "linalg" ], optional = true } nalgebra = { version = "0.34.1", optional = true } +num-bigint = { version = "0.4", optional = true } +num-rational = { version = "0.4", features = [ "num-bigint-std" ], optional = true } [dev-dependencies] approx = "0.5.1" @@ -26,6 +28,11 @@ proptest = "1.10.0" [features] default = [ ] bench = [ "criterion", "faer", "nalgebra" ] +exact = [ "num-bigint", "num-rational" ] + +[[example]] +name = "exact_sign_3x3" +required-features = [ "exact" ] [[bench]] name = "vs_linalg" @@ -36,6 +43,9 @@ required-features = [ "bench" ] lto = "fat" codegen-units = 1 +[package.metadata.docs.rs] +features = [ "exact" ] + [lints.rust] unsafe_code = "forbid" missing_docs = "warn" diff --git a/README.md b/README.md index da85d72..376dfa2 100644 --- a/README.md +++ b/README.md @@ -31,7 +31,8 @@ while keeping the API intentionally small and explicit. - ✅ Const-generic dimensions (no dynamic sizes) - ✅ `const fn` where possible (compile-time evaluation of determinants, dot products, etc.) - ✅ Explicit algorithms (LU, solve, determinant) -- ✅ No runtime dependencies (dev-dependencies are for contributors only) +- ✅ Robust geometric predicates via optional exact arithmetic (`det_sign_exact`) +- ✅ No runtime dependencies by default (optional features may add deps) - ✅ Stack storage only (no heap allocation in core types) - ✅ `unsafe` forbidden @@ -44,8 +45,8 @@ while keeping the API intentionally small and explicit. ## 🔢 Scalar types Today, the core types are implemented for `f64`. The intent is to support `f32` and `f64` -(and `f128` if/when Rust gains a stable primitive for it). Longer term, we may add optional -arbitrary-precision support (e.g. via `rug`) depending on performance. +(and `f128` if/when Rust gains a stable primitive for it). Arbitrary-precision arithmetic +is available via the optional `"exact"` feature (see below). ## 🚀 Quickstart @@ -53,7 +54,7 @@ Add this to your `Cargo.toml`: ```toml [dependencies] -la-stack = "0.1" +la-stack = "0.2" ``` Solve a 5×5 system via LU: @@ -130,27 +131,62 @@ assert_eq!(DET, Some(30.0)); The public `det()` method automatically dispatches through the closed-form path for D ≤ 4 and falls back to LU for D ≥ 5 — no API change needed. +## 🔬 Exact determinant sign (`"exact"` feature) + +The default build has **zero runtime dependencies**. Enable the optional +`exact` Cargo feature to add `det_sign_exact()`, which returns the provably +correct sign (−1, 0, or +1) of the determinant using adaptive-precision +arithmetic (this pulls in `num-bigint` and `num-rational` for `BigRational`): + +```toml +[dependencies] +la-stack = { version = "0.2", features = ["exact"] } +``` + +```rust,ignore +use la_stack::prelude::*; + +let m = Matrix::<3>::from_rows([ + [1.0, 2.0, 3.0], + [4.0, 5.0, 6.0], + [7.0, 8.0, 9.0], +]); +assert_eq!(m.det_sign_exact(), 0); // exactly singular +``` + +For D ≤ 4, a fast f64 filter (error-bounded `det_direct()`) resolves the sign +without allocating. Only near-degenerate or large (D ≥ 5) matrices fall through +to the exact Bareiss algorithm in `BigRational`. + ## 🧩 API at a glance | Type | Storage | Purpose | Key methods | |---|---|---|---| | `Vector` | `[f64; D]` | Fixed-length vector | `new`, `zero`, `dot`, `norm2_sq` | -| `Matrix` | `[[f64; D]; D]` | Fixed-size square matrix | `from_rows`, `zero`, `identity`, `lu`, `ldlt`, `det`, `det_direct` | +| `Matrix` | `[[f64; D]; D]` | Fixed-size square matrix | `from_rows`, `zero`, `identity`, `lu`, `ldlt`, `det`, `det_direct`, `det_sign_exact`¹ | | `Lu` | `Matrix` + pivot array | Factorization for solves/det | `solve_vec`, `det` | | `Ldlt` | `Matrix` | Factorization for symmetric SPD/PSD solves/det | `solve_vec`, `det` | Storage shown above reflects the current `f64` implementation. +¹ Requires `features = ["exact"]`. + ## 📋 Examples The `examples/` directory contains small, runnable programs: +- **`solve_5x5`** — solve a 5×5 system via LU with partial pivoting +- **`det_5x5`** — determinant of a 5×5 matrix via LU +- **`const_det_4x4`** — compile-time 4×4 determinant via `det_direct()` +- **`exact_sign_3x3`** — exact determinant sign of a near-singular 3×3 matrix (requires `exact` feature) + ```bash just examples -# or: +# or individually: cargo run --example solve_5x5 cargo run --example det_5x5 cargo run --example const_det_4x4 +cargo run --features exact --example exact_sign_3x3 ``` ## 🤝 Contributing @@ -174,7 +210,7 @@ If you use this library in academic work, please cite it using [CITATION.cff](CI ## 📚 References -For canonical references to LU / `LDL^T` algorithms used by this crate, see [REFERENCES.md](REFERENCES.md). +For canonical references to the algorithms used by this crate, see [REFERENCES.md](REFERENCES.md). ## 📊 Benchmarks (vs nalgebra/faer) diff --git a/REFERENCES.md b/REFERENCES.md index 89d4a8c..05de553 100644 --- a/REFERENCES.md +++ b/REFERENCES.md @@ -26,6 +26,11 @@ pivoting. For background on the SPD/PSD setting, see [4-5]. For pivoted variants used for symmetric *indefinite* matrices, see [6]. +### Exact determinant sign (adaptive-precision Bareiss) + +`det_sign_exact()` uses a Shewchuk-style f64 error-bound filter [8] backed by exact Bareiss +elimination [7] in `BigRational`. See `src/exact.rs` for the full architecture description. + ## References ### LU / Gaussian elimination @@ -53,3 +58,15 @@ matrices, see [6]. 6. Bunch, J. R., L. Kaufman, and B. N. Parlett. "Decomposition of a Symmetric Matrix." *Numerische Mathematik* 27 (1976/1977): 95–110. [Full text](https://eudml.org/doc/132435) + +### Exact determinant sign (`det_sign_exact`) + +7. Bareiss, Erwin H. "Sylvester's Identity and Multistep Integer-Preserving Gaussian + Elimination." *Mathematics of Computation* 22.103 (1968): 565–578. + [DOI](https://doi.org/10.1090/S0025-5718-1968-0226829-0) · + [PDF](https://www.ams.org/journals/mcom/1968-22-103/S0025-5718-1968-0226829-0/S0025-5718-1968-0226829-0.pdf) +8. Shewchuk, Jonathan Richard. "Adaptive Precision Floating-Point Arithmetic and Fast + Robust Geometric Predicates." *Discrete & Computational Geometry* 18.3 (1997): 305–363. + [DOI](https://doi.org/10.1007/PL00009321) · + [PDF](https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf) + Also: Technical Report CMU-CS-96-140, Carnegie Mellon University, May 1996. diff --git a/examples/exact_sign_3x3.rs b/examples/exact_sign_3x3.rs new file mode 100644 index 0000000..4a9a078 --- /dev/null +++ b/examples/exact_sign_3x3.rs @@ -0,0 +1,45 @@ +//! Exact determinant sign for a near-singular 3×3 matrix. +//! +//! This example demonstrates `det_sign_exact()`, which uses adaptive-precision +//! arithmetic to return the provably correct sign of the determinant — even when +//! the matrix is so close to singular that f64 rounding could give the wrong answer. +//! +//! Run with: `cargo run --features exact --example exact_sign_3x3` + +use la_stack::prelude::*; + +fn main() { + // Base matrix: rows in arithmetic progression → exactly singular (det = 0). + // [[1, 2, 3], + // [4, 5, 6], + // [7, 8, 9]] + // + // Perturb entry (0,0) by 2^-50 ≈ 8.9e-16. + // Exact det = 2^-50 × cofactor(0,0) = 2^-50 × (5×9 − 6×8) = −3 × 2^-50 < 0. + let perturbation = f64::from_bits(0x3CD0_0000_0000_0000); // 2^-50 + let m = Matrix::<3>::from_rows([ + [1.0 + perturbation, 2.0, 3.0], + [4.0, 5.0, 6.0], + [7.0, 8.0, 9.0], + ]); + + let sign = m.det_sign_exact(); + let det_f64 = m.det(DEFAULT_PIVOT_TOL).unwrap(); + + println!("Near-singular 3×3 matrix (perturbation = 2^-50 ≈ {perturbation:.2e}):"); + for r in 0..3 { + print!(" ["); + for c in 0..3 { + if c > 0 { + print!(", "); + } + print!("{:22.18}", m.get(r, c).unwrap()); + } + println!("]"); + } + println!(); + println!("f64 det() = {det_f64:+.6e}"); + println!("det_sign_exact() = {sign}"); + println!(); + println!("The exact sign is −1 (negative), matching the analytical result."); +} diff --git a/justfile b/justfile index b238f7c..4f6e578 100644 --- a/justfile +++ b/justfile @@ -130,6 +130,10 @@ check-fast: ci: check bench-compile test-all examples @echo "🎯 CI checks complete!" +# Clippy for the "exact" feature (catches feature-gated lint issues) +clippy-exact: + cargo clippy --features exact --all-targets -- -D warnings -W clippy::pedantic + # Clean build artifacts clean: cargo clean @@ -145,6 +149,7 @@ clippy: # Common tarpaulin arguments for all coverage runs # Note: -t 300 sets per-test timeout to 5 minutes (needed for slow CI environments) _coverage_base_args := '''--exclude-files 'benches/*' --exclude-files 'examples/*' \ + --features exact \ --workspace --lib --tests \ -t 300 --verbose --implicit-test-threads''' @@ -179,15 +184,16 @@ coverage-ci: default: @just --list -# Documentation build check +# Documentation build check (includes exact feature for full API coverage) doc-check: - RUSTDOCFLAGS='-D warnings' cargo doc --no-deps + RUSTDOCFLAGS='-D warnings' cargo doc --no-deps --features exact # Examples examples: cargo run --quiet --example det_5x5 cargo run --quiet --example solve_5x5 cargo run --quiet --example const_det_4x4 + cargo run --quiet --features exact --example exact_sign_3x3 # Fix (mutating): apply formatters/auto-fixes fix: toml-fmt fmt python-fix shell-fmt markdown-fix yaml-fix @@ -500,9 +506,13 @@ test: cargo test --lib --verbose cargo test --doc --verbose -test-all: test test-integration test-python +test-all: test test-integration test-exact test-python @echo "✅ All tests passed" +# Tests for the "exact" feature (det_sign_exact + BigRational Bareiss) +test-exact: + cargo test --features exact --verbose + test-integration: cargo test --tests --verbose diff --git a/src/exact.rs b/src/exact.rs new file mode 100644 index 0000000..97b53d0 --- /dev/null +++ b/src/exact.rs @@ -0,0 +1,503 @@ +//! Exact determinant sign via arbitrary-precision rational arithmetic. +//! +//! This module is only compiled when the `"exact"` Cargo feature is enabled. +//! +//! # Architecture +//! +//! `det_sign_exact` uses a two-stage adaptive-precision approach inspired by +//! Shewchuk's robust geometric predicates: +//! +//! 1. **Fast filter (D ≤ 4)**: compute `det_direct()` and a conservative error +//! bound. If `|det| > bound`, the f64 sign is provably correct — return +//! immediately without allocating. +//! 2. **Exact fallback**: convert entries to `BigRational` and run the Bareiss +//! algorithm (fraction-free Gaussian elimination) for a guaranteed-correct +//! sign. + +use num_bigint::{BigInt, Sign}; +use num_rational::BigRational; + +use crate::matrix::Matrix; + +// --------------------------------------------------------------------------- +// Error-bound constants for the f64 fast filter. +// +// Each constant bounds the absolute error of `det_direct()` relative to the +// *permanent* (sum of absolute products in the Leibniz expansion). The +// constants are conservative over-estimates following Shewchuk's methodology; +// over-estimating just means we fall through to Bareiss more often. +// --------------------------------------------------------------------------- + +const EPS: f64 = f64::EPSILON; // 2^-52 + +/// D=2: one f64 multiply + one FMA → 2 rounding events. +const ERR_COEFF_2: f64 = 3.0 * EPS + 16.0 * EPS * EPS; + +/// D=3: three 2×2 FMA minors + nested FMA combination. +const ERR_COEFF_3: f64 = 8.0 * EPS + 64.0 * EPS * EPS; + +/// D=4: six hoisted 2×2 minors → four 3×3 cofactors → FMA row combination. +const ERR_COEFF_4: f64 = 12.0 * EPS + 128.0 * EPS * EPS; + +/// Conservative absolute error bound for `det_direct()`, or `None` for D ≥ 5. +/// +/// Returns `Some(bound)` such that `|det_direct() - det_exact| ≤ bound`. +fn det_errbound(m: &Matrix) -> Option { + match D { + 0 | 1 => Some(0.0), // No arithmetic — result is exact. + 2 => { + let r = &m.rows; + let permanent = (r[0][0] * r[1][1]).abs() + (r[0][1] * r[1][0]).abs(); + Some(ERR_COEFF_2 * permanent) + } + 3 => { + let r = &m.rows; + let pm00 = (r[1][1] * r[2][2]).abs() + (r[1][2] * r[2][1]).abs(); + let pm01 = (r[1][0] * r[2][2]).abs() + (r[1][2] * r[2][0]).abs(); + let pm02 = (r[1][0] * r[2][1]).abs() + (r[1][1] * r[2][0]).abs(); + let permanent = r[0][2] + .abs() + .mul_add(pm02, r[0][1].abs().mul_add(pm01, r[0][0].abs() * pm00)); + Some(ERR_COEFF_3 * permanent) + } + 4 => { + let r = &m.rows; + // 2×2 minor permanents from rows 2–3. + let sp23 = (r[2][2] * r[3][3]).abs() + (r[2][3] * r[3][2]).abs(); + let sp13 = (r[2][1] * r[3][3]).abs() + (r[2][3] * r[3][1]).abs(); + let sp12 = (r[2][1] * r[3][2]).abs() + (r[2][2] * r[3][1]).abs(); + let sp03 = (r[2][0] * r[3][3]).abs() + (r[2][3] * r[3][0]).abs(); + let sp02 = (r[2][0] * r[3][2]).abs() + (r[2][2] * r[3][0]).abs(); + let sp01 = (r[2][0] * r[3][1]).abs() + (r[2][1] * r[3][0]).abs(); + // 3×3 cofactor permanents from row 1. + let pc0 = r[1][3] + .abs() + .mul_add(sp12, r[1][2].abs().mul_add(sp13, r[1][1].abs() * sp23)); + let pc1 = r[1][3] + .abs() + .mul_add(sp02, r[1][2].abs().mul_add(sp03, r[1][0].abs() * sp23)); + let pc2 = r[1][3] + .abs() + .mul_add(sp01, r[1][1].abs().mul_add(sp03, r[1][0].abs() * sp13)); + let pc3 = r[1][2] + .abs() + .mul_add(sp01, r[1][1].abs().mul_add(sp02, r[1][0].abs() * sp12)); + let permanent = r[0][3].abs().mul_add( + pc3, + r[0][2] + .abs() + .mul_add(pc2, r[0][1].abs().mul_add(pc1, r[0][0].abs() * pc0)), + ); + Some(ERR_COEFF_4 * permanent) + } + _ => None, + } +} + +/// Convert an `f64` to an exact `BigRational`. +/// +/// Every finite `f64` is exactly representable as a rational number (`m × 2^e`), +/// so this conversion is lossless. +/// +/// # Panics +/// Panics if `x` is NaN or infinite. +fn f64_to_bigrational(x: f64) -> BigRational { + BigRational::from_float(x).expect("non-finite matrix entry in det_sign_exact") +} + +/// Compute the exact determinant of a `D×D` matrix using the Bareiss algorithm +/// (fraction-free Gaussian elimination) in `BigRational` arithmetic. +/// +/// Returns the determinant as an exact `BigRational`. +fn bareiss_det(m: &Matrix) -> BigRational { + if D == 0 { + return BigRational::from_integer(BigInt::from(1)); + } + if D == 1 { + return f64_to_bigrational(m.rows[0][0]); + } + + // Convert f64 entries to exact BigRational. + let mut a: Vec> = Vec::with_capacity(D); + for r in 0..D { + let mut row = Vec::with_capacity(D); + for c in 0..D { + row.push(f64_to_bigrational(m.rows[r][c])); + } + a.push(row); + } + + let zero = BigRational::from_integer(BigInt::from(0)); + let mut prev_pivot = BigRational::from_integer(BigInt::from(1)); + let mut sign: i8 = 1; + + for k in 0..D { + // Partial pivoting: find non-zero entry in column k at or below row k. + if a[k][k] == zero { + let mut found = false; + for i in (k + 1)..D { + if a[i][k] != zero { + a.swap(k, i); + sign = -sign; + found = true; + break; + } + } + if !found { + // Entire column below (and including) diagonal is zero → singular. + return zero; + } + } + + // Bareiss elimination for rows below k. + for i in (k + 1)..D { + for j in (k + 1)..D { + // a[i][j] = (a[k][k] * a[i][j] - a[i][k] * a[k][j]) / prev_pivot + a[i][j] = (&a[k][k] * &a[i][j] - &a[i][k] * &a[k][j]) / &prev_pivot; + } + // Zero out the eliminated column entry (not needed for det, but keeps + // the matrix consistent for potential debugging). + a[i][k] = zero.clone(); + } + + prev_pivot = a[k][k].clone(); + } + + let det = &a[D - 1][D - 1]; + if sign < 0 { -det } else { det.clone() } +} + +impl Matrix { + /// Exact determinant sign using adaptive-precision arithmetic. + /// + /// Returns `1` if `det > 0`, `-1` if `det < 0`, and `0` if `det == 0` (singular). + /// + /// For D ≤ 4, a fast f64 filter is tried first: `det_direct()` is compared + /// against a conservative error bound derived from the matrix permanent. + /// If the f64 result clearly exceeds the bound, the sign is returned + /// immediately without allocating. Otherwise (and always for D ≥ 5), the + /// Bareiss algorithm runs in exact [`BigRational`] arithmetic. + /// + /// # When to use + /// + /// Use this when the sign of the determinant must be correct regardless of + /// floating-point conditioning (e.g. geometric predicates on near-degenerate + /// configurations). For well-conditioned matrices the fast filter resolves + /// the sign without touching `BigRational`, so the overhead is minimal. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// let m = Matrix::<3>::from_rows([ + /// [1.0, 2.0, 3.0], + /// [4.0, 5.0, 6.0], + /// [7.0, 8.0, 9.0], + /// ]); + /// // This matrix is singular (row 3 = row 1 + row 2 in exact arithmetic). + /// assert_eq!(m.det_sign_exact(), 0); + /// + /// assert_eq!(Matrix::<3>::identity().det_sign_exact(), 1); + /// ``` + /// + /// # Panics + /// Panics if any matrix entry is NaN or infinite. + #[must_use] + pub fn det_sign_exact(&self) -> i8 { + // Stage 1: f64 fast filter for D ≤ 4. + if let (Some(det_f64), Some(err)) = (self.det_direct(), det_errbound(self)) { + assert!( + det_f64.is_finite(), + "non-finite matrix entry in det_sign_exact" + ); + if det_f64 > err { + return 1; + } + if det_f64 < -err { + return -1; + } + } + + // Stage 2: exact Bareiss fallback. + let det = bareiss_det(self); + match det.numer().sign() { + Sign::Plus => 1, + Sign::Minus => -1, + Sign::NoSign => 0, + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::DEFAULT_PIVOT_TOL; + + #[test] + fn det_sign_exact_d0_is_positive() { + assert_eq!(Matrix::<0>::zero().det_sign_exact(), 1); + } + + #[test] + fn det_sign_exact_d1_positive() { + let m = Matrix::<1>::from_rows([[42.0]]); + assert_eq!(m.det_sign_exact(), 1); + } + + #[test] + fn det_sign_exact_d1_negative() { + let m = Matrix::<1>::from_rows([[-3.5]]); + assert_eq!(m.det_sign_exact(), -1); + } + + #[test] + fn det_sign_exact_d1_zero() { + let m = Matrix::<1>::from_rows([[0.0]]); + assert_eq!(m.det_sign_exact(), 0); + } + + #[test] + fn det_sign_exact_identity_2d() { + assert_eq!(Matrix::<2>::identity().det_sign_exact(), 1); + } + + #[test] + fn det_sign_exact_identity_3d() { + assert_eq!(Matrix::<3>::identity().det_sign_exact(), 1); + } + + #[test] + fn det_sign_exact_identity_4d() { + assert_eq!(Matrix::<4>::identity().det_sign_exact(), 1); + } + + #[test] + fn det_sign_exact_identity_5d() { + assert_eq!(Matrix::<5>::identity().det_sign_exact(), 1); + } + + #[test] + fn det_sign_exact_singular_duplicate_rows() { + let m = Matrix::<3>::from_rows([ + [1.0, 2.0, 3.0], + [4.0, 5.0, 6.0], + [1.0, 2.0, 3.0], // duplicate of row 0 + ]); + assert_eq!(m.det_sign_exact(), 0); + } + + #[test] + fn det_sign_exact_singular_linear_combination() { + // Row 2 = row 0 + row 1 in exact arithmetic. + let m = Matrix::<3>::from_rows([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [5.0, 7.0, 9.0]]); + assert_eq!(m.det_sign_exact(), 0); + } + + #[test] + fn det_sign_exact_negative_det_row_swap() { + // Swapping two rows of the identity negates the determinant. + let m = Matrix::<3>::from_rows([[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]); + assert_eq!(m.det_sign_exact(), -1); + } + + #[test] + fn det_sign_exact_negative_det_known() { + // det([[1,2],[3,4]]) = 1*4 - 2*3 = -2 + let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); + assert_eq!(m.det_sign_exact(), -1); + } + + #[test] + fn det_sign_exact_agrees_with_det_for_spd() { + // SPD matrix → positive determinant. + let m = Matrix::<3>::from_rows([[4.0, 2.0, 0.0], [2.0, 5.0, 1.0], [0.0, 1.0, 3.0]]); + assert_eq!(m.det_sign_exact(), 1); + assert!(m.det(DEFAULT_PIVOT_TOL).unwrap() > 0.0); + } + + /// Near-singular matrix with an exact perturbation. + /// + /// The base matrix `[[1,2,3],[4,5,6],[7,8,9]]` is exactly singular (rows in + /// arithmetic progression). Adding `2^-50` to entry (0,0) makes + /// `det = 2^-50 × cofactor(0,0) = 2^-50 × (5×9 − 6×8) = −3 × 2^-50 < 0`. + /// Both f64 `det_direct()` and `det_sign_exact()` should agree here. + #[test] + fn det_sign_exact_near_singular_perturbation() { + let perturbation = f64::from_bits(0x3CD0_0000_0000_0000); // 2^-50 + let m = Matrix::<3>::from_rows([ + [1.0 + perturbation, 2.0, 3.0], + [4.0, 5.0, 6.0], + [7.0, 8.0, 9.0], + ]); + // Exact: det = perturbation × (5×9 − 6×8) = perturbation × (−3) < 0. + assert_eq!(m.det_sign_exact(), -1); + } + + /// For D ≤ 4, well-conditioned matrices should hit the fast filter + /// and never allocate `BigRational`. We can't directly observe this, + /// but we verify correctness for a range of known signs. + #[test] + fn det_sign_exact_fast_filter_positive_4x4() { + let m = Matrix::<4>::from_rows([ + [2.0, 1.0, 0.0, 0.0], + [1.0, 3.0, 1.0, 0.0], + [0.0, 1.0, 4.0, 1.0], + [0.0, 0.0, 1.0, 5.0], + ]); + // SPD tridiagonal → positive det. + assert_eq!(m.det_sign_exact(), 1); + } + + #[test] + fn det_sign_exact_fast_filter_negative_4x4() { + // Swap rows 0 and 1 of the above → negate det. + let m = Matrix::<4>::from_rows([ + [1.0, 3.0, 1.0, 0.0], + [2.0, 1.0, 0.0, 0.0], + [0.0, 1.0, 4.0, 1.0], + [0.0, 0.0, 1.0, 5.0], + ]); + assert_eq!(m.det_sign_exact(), -1); + } + + #[test] + fn det_sign_exact_subnormal_entries() { + // Subnormal f64 values should convert losslessly. + let tiny = 5e-324_f64; // smallest positive subnormal + assert!(tiny.is_subnormal()); + + let m = Matrix::<2>::from_rows([[tiny, 0.0], [0.0, tiny]]); + // det = tiny^2 > 0 + assert_eq!(m.det_sign_exact(), 1); + } + + #[test] + #[should_panic(expected = "non-finite matrix entry")] + fn det_sign_exact_panics_on_nan() { + let m = Matrix::<2>::from_rows([[f64::NAN, 0.0], [0.0, 1.0]]); + let _ = m.det_sign_exact(); + } + + #[test] + #[should_panic(expected = "non-finite matrix entry")] + fn det_sign_exact_panics_on_infinity() { + let m = Matrix::<2>::from_rows([[f64::INFINITY, 0.0], [0.0, 1.0]]); + let _ = m.det_sign_exact(); + } + + #[test] + #[should_panic(expected = "non-finite matrix entry")] + fn det_sign_exact_panics_on_nan_5x5() { + // D ≥ 5 bypasses the fast filter, exercising the bareiss_det path. + let mut m = Matrix::<5>::identity(); + m.set(2, 3, f64::NAN); + let _ = m.det_sign_exact(); + } + + #[test] + #[should_panic(expected = "non-finite matrix entry")] + fn det_sign_exact_panics_on_infinity_5x5() { + let mut m = Matrix::<5>::identity(); + m.set(0, 0, f64::INFINITY); + let _ = m.det_sign_exact(); + } + + #[test] + fn det_sign_exact_pivot_needed_5x5() { + // D ≥ 5 skips the fast filter → exercises Bareiss pivoting. + // Permutation matrix with a single swap (rows 0↔1) → det = −1. + let m = Matrix::<5>::from_rows([ + [0.0, 1.0, 0.0, 0.0, 0.0], + [1.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 1.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 1.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 1.0], + ]); + assert_eq!(m.det_sign_exact(), -1); + } + + #[test] + fn det_sign_exact_5x5_known() { + // det of a permutation matrix with two swaps = +1 (even permutation). + let m = Matrix::<5>::from_rows([ + [0.0, 1.0, 0.0, 0.0, 0.0], + [1.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 1.0, 0.0], + [0.0, 0.0, 1.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 1.0], + ]); + // Two transpositions → even permutation → det = +1 + assert_eq!(m.det_sign_exact(), 1); + } + + // ----------------------------------------------------------------------- + // Direct tests for internal helpers (coverage of private functions) + // ----------------------------------------------------------------------- + + #[test] + fn det_errbound_d0_is_zero() { + assert_eq!(det_errbound(&Matrix::<0>::zero()), Some(0.0)); + } + + #[test] + fn det_errbound_d1_is_zero() { + assert_eq!(det_errbound(&Matrix::<1>::from_rows([[42.0]])), Some(0.0)); + } + + #[test] + fn det_errbound_d2_positive() { + let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); + let bound = det_errbound(&m).unwrap(); + assert!(bound > 0.0); + // bound = ERR_COEFF_2 * (|1*4| + |2*3|) = ERR_COEFF_2 * 10 + assert!(ERR_COEFF_2.mul_add(-10.0, bound).abs() < 1e-30); + } + + #[test] + fn det_errbound_d3_positive() { + let m = Matrix::<3>::identity(); + let bound = det_errbound(&m).unwrap(); + assert!(bound > 0.0); + } + + #[test] + fn det_errbound_d4_positive() { + let m = Matrix::<4>::identity(); + let bound = det_errbound(&m).unwrap(); + assert!(bound > 0.0); + } + + #[test] + fn det_errbound_d5_is_none() { + assert_eq!(det_errbound(&Matrix::<5>::identity()), None); + } + + #[test] + fn bareiss_det_d0_is_one() { + let det = bareiss_det(&Matrix::<0>::zero()); + assert_eq!(det, BigRational::from_integer(BigInt::from(1))); + } + + #[test] + fn bareiss_det_d1_returns_entry() { + let det = bareiss_det(&Matrix::<1>::from_rows([[7.0]])); + assert_eq!(det, f64_to_bigrational(7.0)); + } + + #[test] + fn bareiss_det_d3_with_pivoting() { + // First column has zero on diagonal → exercises pivot swap + break. + let m = Matrix::<3>::from_rows([[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]); + let det = bareiss_det(&m); + // det of this permutation matrix = -1 + assert_eq!(det, BigRational::from_integer(BigInt::from(-1))); + } + + #[test] + fn bareiss_det_singular_all_zeros_in_column() { + // Column 1 is all zeros below diagonal after elimination → singular. + let m = Matrix::<3>::from_rows([[1.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]); + let det = bareiss_det(&m); + assert_eq!(det, BigRational::from_integer(BigInt::from(0))); + } +} diff --git a/src/lib.rs b/src/lib.rs index 3210188..46e72e5 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -48,6 +48,8 @@ mod readme_doctests { fn det_5x5_ldlt_example() {} } +#[cfg(feature = "exact")] +mod exact; mod ldlt; mod lu; mod matrix; diff --git a/tests/proptest_exact.rs b/tests/proptest_exact.rs new file mode 100644 index 0000000..1ef388c --- /dev/null +++ b/tests/proptest_exact.rs @@ -0,0 +1,72 @@ +//! Property-based tests for the `det_sign_exact` API (requires `exact` feature). + +#![cfg(feature = "exact")] + +use pastey::paste; +use proptest::prelude::*; + +use la_stack::prelude::*; + +fn small_nonzero_f64() -> impl Strategy { + prop_oneof![(-1000i16..=-1i16), (1i16..=1000i16)].prop_map(|x| f64::from(x) / 10.0) +} + +macro_rules! gen_det_sign_exact_proptests { + ($d:literal) => { + paste! { + proptest! { + #![proptest_config(ProptestConfig::with_cases(64))] + + #[test] + fn []( + diag in proptest::array::[](small_nonzero_f64()), + ) { + // Diagonal matrix: determinant sign = product of diagonal signs. + let mut rows = [[0.0f64; $d]; $d]; + for i in 0..$d { + rows[i][i] = diag[i]; + } + let m = Matrix::<$d>::from_rows(rows); + + let exact_sign = m.det_sign_exact(); + + // Expected sign from the product of diagonal entries. + let neg_count = diag.iter().filter(|&&x| x < 0.0).count(); + let expected_sign: i8 = if neg_count % 2 == 0 { 1 } else { -1 }; + + prop_assert_eq!(exact_sign, expected_sign); + } + + #[test] + fn []( + diag in proptest::array::[](small_nonzero_f64()), + ) { + // For well-conditioned diagonal matrices, det().signum() + // should agree with det_sign_exact(). + let mut rows = [[0.0f64; $d]; $d]; + for i in 0..$d { + rows[i][i] = diag[i]; + } + let m = Matrix::<$d>::from_rows(rows); + + let exact_sign = m.det_sign_exact(); + let fp_det = m.det(DEFAULT_PIVOT_TOL).unwrap(); + let fp_sign: i8 = if fp_det > 0.0 { + 1 + } else if fp_det < 0.0 { + -1 + } else { + 0 + }; + + prop_assert_eq!(exact_sign, fp_sign); + } + } + } + }; +} + +gen_det_sign_exact_proptests!(2); +gen_det_sign_exact_proptests!(3); +gen_det_sign_exact_proptests!(4); +gen_det_sign_exact_proptests!(5); diff --git a/typos.toml b/typos.toml index 2a83eba..0af273f 100644 --- a/typos.toml +++ b/typos.toml @@ -29,3 +29,4 @@ Trefethen = "Trefethen" faer = "faer" frhs = "frhs" nrhs = "nrhs" +numer = "numer"