diff --git a/.github/workflows/benchmarks.yml b/.github/workflows/benchmarks.yml index d28363a..5469388 100644 --- a/.github/workflows/benchmarks.yml +++ b/.github/workflows/benchmarks.yml @@ -114,8 +114,11 @@ jobs: if [ -d target/criterion/exact_d2/det/main ]; then echo "::notice::Baseline found — comparing against main" + # --baseline-lenient rather than --baseline: benches added on the + # PR branch that don't yet exist in the main baseline get a + # "no baseline data" notice instead of aborting the whole run. cargo bench --features bench,exact --bench exact \ - -- --baseline main 2>&1 | tee bench-output.txt + -- --baseline-lenient main 2>&1 | tee bench-output.txt else echo "::notice::No baseline found — running without comparison" cargo bench --features bench,exact --bench exact \ diff --git a/AGENTS.md b/AGENTS.md index 26cb01f..80184d9 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -221,7 +221,10 @@ When creating or updating issues: - `exact` — enables exact arithmetic methods via `BigRational`: `det_exact()`, `det_exact_f64()`, `det_sign_exact()`, `solve_exact()`, and `solve_exact_f64()`. - Also re-exports `BigRational` from the crate root and prelude. + Re-exports `BigInt`, `BigRational`, and the commonly needed `num-traits` + items (`FromPrimitive`, `ToPrimitive`, and `Signed`) from the crate root and prelude + (so consumers get usable `from_f64` / `to_f64` / `is_positive` etc. without adding + `num-bigint` / `num-rational` / `num-traits` as their own deps). Gates `src/exact.rs`, additional tests, and the `exact_det_3x3`/`exact_sign_3x3`/`exact_solve_3x3` examples. Clippy, doc builds, and test commands have dedicated `--features exact` variants. @@ -245,7 +248,8 @@ When creating or updating issues: - The public API re-exports these items from `src/lib.rs`. - The `justfile` defines all dev workflows (see `just --list`). - Dev-only benchmarks live in `benches/vs_linalg.rs` (Criterion + nalgebra/faer comparison) - and `benches/exact.rs` (exact arithmetic across D=2–5). + and `benches/exact.rs` (exact arithmetic across D=2–5, plus adversarial-input groups + `exact_near_singular_3x3`, `exact_large_entries_3x3`, `exact_hilbert_4x4`, `exact_hilbert_5x5`). - Python scripts under `scripts/`: - `bench_compare.py`: exact-arithmetic benchmark comparison across releases (generates `docs/PERFORMANCE.md`) - `criterion_dim_plot.py`: benchmark plotting (CSV + SVG + README table update) diff --git a/README.md b/README.md index b4a1507..3e37f39 100644 --- a/README.md +++ b/README.md @@ -190,8 +190,13 @@ assert!((x[0] - 1.0).abs() <= f64::EPSILON); assert!((x[1] - 2.0).abs() <= f64::EPSILON); ``` -`BigRational` is re-exported from the crate root and prelude when the `exact` -feature is enabled, so consumers don't need to depend on `num-rational` directly. +With the `exact` feature enabled, `BigInt` and `BigRational` are re-exported +from the crate root and prelude, alongside the most commonly needed +`num-traits` items (`FromPrimitive`, `ToPrimitive`, `Signed`). This lets +consumers construct exact values (`BigRational::from_f64`, `from_i64`), query +sign (`is_positive` / `is_negative`), and convert back to `f64` (`to_f64`) +with a single `use la_stack::prelude::*;` — no need to add `num-bigint`, +`num-rational`, or `num-traits` to their own `Cargo.toml`. For `det_sign_exact()`, D ≤ 4 matrices use a fast f64 filter (error-bounded `det_direct()`) that resolves the sign without allocating. Only near-degenerate diff --git a/benches/exact.rs b/benches/exact.rs index 3fd7846..228cb45 100644 --- a/benches/exact.rs +++ b/benches/exact.rs @@ -1,10 +1,22 @@ //! Benchmarks for exact arithmetic operations. //! //! These benchmarks measure the performance of the `exact` feature's -//! arbitrary-precision methods across dimensions D=2..5 (the primary -//! target for geometric predicates). +//! arbitrary-precision methods. They are organised into two classes: +//! +//! 1. **General-case benches** (`exact_d{2..5}`) — a single +//! well-conditioned diagonally-dominant matrix per dimension. These +//! measure typical-case performance and track regressions against a +//! reproducible input. +//! 2. **Adversarial / extreme-input benches** — matrices chosen to +//! stress specific corners of the exact-arithmetic pipeline: +//! near-singularity (forces the Bareiss fallback), large f64 entries +//! (stresses intermediate `BigInt` growth), and Hilbert-style +//! ill-conditioning (wide range of `(mantissa, exponent)` pairs in +//! the `f64_decompose → BigInt` path). These measure tail behaviour +//! that fixed well-conditioned inputs miss and provide stronger +//! empirical evidence for `docs/PERFORMANCE.md`. -use criterion::Criterion; +use criterion::{BenchmarkGroup, Criterion, measurement::WallTime}; use la_stack::{Matrix, Vector}; use pastey::paste; use std::hint::black_box; @@ -47,8 +59,13 @@ fn make_vector_array() -> [f64; D] { } /// Near-singular matrix: base singular matrix + tiny perturbation. -/// This forces the exact Bareiss fallback in `det_sign_exact` (the fast -/// f64 filter cannot resolve the sign). +/// +/// The base `[[1,2,3],[4,5,6],[7,8,9]]` is exactly singular; adding +/// `2^-50` to entry (0,0) makes `det = -3 × 2^-50 ≠ 0`. The f64 filter +/// in `det_sign_exact` cannot resolve this sign, so Bareiss is forced; +/// `solve_exact` is the primary use case for near-degenerate inputs +/// (exact circumcenter etc.) and exercises the largest intermediate +/// `BigInt` values in the hybrid solve. #[inline] fn near_singular_3x3() -> Matrix<3> { let perturbation = f64::from_bits(0x3CD0_0000_0000_0000); // 2^-50 @@ -59,6 +76,97 @@ fn near_singular_3x3() -> Matrix<3> { ]) } +/// Large-entry 3×3: strictly diagonally-dominant matrix with diagonal +/// entries near `f64::MAX / 2` and ones elsewhere. +/// +/// Each big entry decomposes into a 53-bit mantissa with exponent `~970`; +/// the unit off-diagonals have exponent `0`, so the shared `e_min = 0` +/// shift in `component_to_bigint` produces `BigInt`s of `~1023` bits for +/// the diagonal and small integers elsewhere. Bareiss fraction-free +/// updates then multiply these together, stressing the big-integer +/// multiply and allocator along the full `O(D³)` elimination phase. The +/// matrix is non-singular (det ≈ `big³`) so both `det_*` and `solve_*` +/// paths complete. +#[inline] +fn large_entries_3x3() -> Matrix<3> { + let big = f64::MAX / 2.0; + Matrix::<3>::from_rows([[big, 1.0, 1.0], [1.0, big, 1.0], [1.0, 1.0, big]]) +} + +/// Hilbert matrix `H[i][j] = 1 / (i + j + 1)`. +/// +/// Most entries (`1/3`, `1/5`, `1/6`, `1/7`, …) are non-terminating in +/// binary, so every cell has a distinct 53-bit mantissa and a small +/// negative exponent. `f64_decompose` therefore produces a wide mix of +/// `(mantissa, exponent)` pairs with no shared power-of-two factors, +/// and the scaling shift to the common `e_min` yields `BigInt` values +/// of varied bit-lengths — a different kind of adversarial input from +/// the large-entries case. Hilbert matrices are also classically +/// ill-conditioned (condition number grows exponentially with D), so +/// they are a realistic stand-in for the near-degenerate geometric +/// predicate inputs that motivate exact arithmetic. +#[inline] +#[allow(clippy::cast_precision_loss)] +fn hilbert() -> Matrix { + let mut rows = [[0.0; D]; D]; + let mut r = 0; + while r < D { + let mut c = 0; + while c < D { + rows[r][c] = 1.0 / ((r + c + 1) as f64); + c += 1; + } + r += 1; + } + Matrix::::from_rows(rows) +} + +/// Populate a Criterion group with the four headline exact-arithmetic +/// benches on a single `(matrix, rhs)` pair: `det_sign_exact`, +/// `det_exact`, `solve_exact`, `solve_exact_f64`. +/// +/// Used by every adversarial-input group so each one measures the same +/// operations, making the resulting tables directly comparable. +fn bench_extreme_group( + group: &mut BenchmarkGroup<'_, WallTime>, + m: Matrix, + rhs: Vector, +) { + group.bench_function("det_sign_exact", |bencher| { + bencher.iter(|| { + let sign = black_box(m) + .det_sign_exact() + .expect("finite matrix entries"); + black_box(sign); + }); + }); + + group.bench_function("det_exact", |bencher| { + bencher.iter(|| { + let det = black_box(m).det_exact().expect("finite matrix entries"); + black_box(det); + }); + }); + + group.bench_function("solve_exact", |bencher| { + bencher.iter(|| { + let x = black_box(m) + .solve_exact(black_box(rhs)) + .expect("non-singular matrix with finite entries"); + let _ = black_box(x); + }); + }); + + group.bench_function("solve_exact_f64", |bencher| { + bencher.iter(|| { + let x = black_box(m) + .solve_exact_f64(black_box(rhs)) + .expect("solution representable in f64"); + let _ = black_box(x); + }); + }); +} + macro_rules! gen_exact_benches_for_dim { ($c:expr, $d:literal) => { paste! {{ @@ -72,7 +180,7 @@ macro_rules! gen_exact_benches_for_dim { bencher.iter(|| { let det = black_box(a) .det(la_stack::DEFAULT_PIVOT_TOL) - .expect("should not fail"); + .expect("diagonally dominant matrix is non-singular"); black_box(det); }); }); @@ -87,7 +195,7 @@ macro_rules! gen_exact_benches_for_dim { // === det_exact (BigRational result) === [].bench_function("det_exact", |bencher| { bencher.iter(|| { - let det = black_box(a).det_exact().expect("should not fail"); + let det = black_box(a).det_exact().expect("finite matrix entries"); black_box(det); }); }); @@ -95,7 +203,9 @@ macro_rules! gen_exact_benches_for_dim { // === det_exact_f64 (exact → f64) === [].bench_function("det_exact_f64", |bencher| { bencher.iter(|| { - let det = black_box(a).det_exact_f64().expect("should not fail"); + let det = black_box(a) + .det_exact_f64() + .expect("det representable in f64"); black_box(det); }); }); @@ -103,7 +213,9 @@ macro_rules! gen_exact_benches_for_dim { // === det_sign_exact (adaptive: fast filter + exact fallback) === [].bench_function("det_sign_exact", |bencher| { bencher.iter(|| { - let sign = black_box(a).det_sign_exact().expect("should not fail"); + let sign = black_box(a) + .det_sign_exact() + .expect("finite matrix entries"); black_box(sign); }); }); @@ -111,7 +223,9 @@ macro_rules! gen_exact_benches_for_dim { // === solve_exact (BigRational result) === [].bench_function("solve_exact", |bencher| { bencher.iter(|| { - let x = black_box(a).solve_exact(black_box(rhs)).expect("should not fail"); + let x = black_box(a) + .solve_exact(black_box(rhs)) + .expect("diagonally dominant matrix is non-singular"); black_box(x); }); }); @@ -119,7 +233,9 @@ macro_rules! gen_exact_benches_for_dim { // === solve_exact_f64 (exact → f64) === [].bench_function("solve_exact_f64", |bencher| { bencher.iter(|| { - let x = black_box(a).solve_exact_f64(black_box(rhs)).expect("should not fail"); + let x = black_box(a) + .solve_exact_f64(black_box(rhs)) + .expect("solution representable in f64"); black_box(x); }); }); @@ -140,25 +256,50 @@ fn main() { gen_exact_benches_for_dim!(&mut c, 5); } - // Near-singular 3×3: forces Bareiss fallback in det_sign_exact. + // === Adversarial / extreme-input groups === + // + // Each group runs the same four exact-arithmetic benches + // (`det_sign_exact`, `det_exact`, `solve_exact`, `solve_exact_f64`) + // via `bench_extreme_group`, so the resulting tables are directly + // comparable across input classes. + + // Near-singular 3×3: forces Bareiss fallback in det_sign_exact and + // exercises the largest intermediate BigInt values in solve_exact + // (the primary motivating use case for exact solve). { - let m = near_singular_3x3(); let mut group = c.benchmark_group("exact_near_singular_3x3"); + bench_extreme_group( + &mut group, + near_singular_3x3(), + Vector::<3>::new([1.0, 2.0, 3.0]), + ); + group.finish(); + } - group.bench_function("det_sign_exact", |bencher| { - bencher.iter(|| { - let sign = black_box(m).det_sign_exact().expect("should not fail"); - black_box(sign); - }); - }); + // Large-entry 3×3: diagonal entries near `f64::MAX / 2` stress + // BigInt growth during Bareiss forward elimination. + { + let mut group = c.benchmark_group("exact_large_entries_3x3"); + bench_extreme_group( + &mut group, + large_entries_3x3(), + Vector::<3>::new([1.0, 1.0, 1.0]), + ); + group.finish(); + } - group.bench_function("det_exact", |bencher| { - bencher.iter(|| { - let det = black_box(m).det_exact().expect("should not fail"); - black_box(det); - }); - }); + // Hilbert 4×4 and 5×5: classically ill-conditioned matrices whose + // entries span many orders of magnitude in `(mantissa, exponent)` + // space, exercising the f64 → BigInt scaling path. + { + let mut group = c.benchmark_group("exact_hilbert_4x4"); + bench_extreme_group(&mut group, hilbert::<4>(), Vector::<4>::new([1.0; 4])); + group.finish(); + } + { + let mut group = c.benchmark_group("exact_hilbert_5x5"); + bench_extreme_group(&mut group, hilbert::<5>(), Vector::<5>::new([1.0; 5])); group.finish(); } diff --git a/docs/BENCHMARKING.md b/docs/BENCHMARKING.md index def50a0..94bc1ea 100644 --- a/docs/BENCHMARKING.md +++ b/docs/BENCHMARKING.md @@ -14,6 +14,20 @@ la-stack has two Criterion benchmark suites: (`det_exact`, `solve_exact`, `det_sign_exact`, etc.) alongside f64 baselines (`det`, `det_direct`) across D=2–5. Use this to understand the cost of exact arithmetic and track optimization progress. + In addition to the per-dimension groups (`exact_d{2..5}`), the suite + includes four adversarial-input groups designed to stress specific + corners of the pipeline: + - `exact_near_singular_3x3` — a 2^-50 perturbation of a singular base + matrix; forces the Bareiss fallback in `det_sign_exact` and + exercises the largest intermediate `BigInt` values in `solve_exact`. + - `exact_large_entries_3x3` — diagonal entries near `f64::MAX / 2` + stress `BigInt` growth during Bareiss forward elimination. + - `exact_hilbert_4x4` / `exact_hilbert_5x5` — classically + ill-conditioned matrices whose non-terminating-in-binary entries + stress the `f64_decompose → BigInt` scaling path. + Each adversarial group runs the same four benches (`det_sign_exact`, + `det_exact`, `solve_exact`, `solve_exact_f64`) so the resulting tables + are directly comparable across input classes. ## Quick reference diff --git a/scripts/bench_compare.py b/scripts/bench_compare.py index 0892009..4feeed8 100644 --- a/scripts/bench_compare.py +++ b/scripts/bench_compare.py @@ -38,12 +38,22 @@ # --------------------------------------------------------------------------- # Groups and the benchmarks within each group that we track. +# +# Mirrors the structure of `benches/exact.rs`: general-case per-dimension +# groups (`exact_d{2..5}`) plus adversarial/extreme-input groups that +# share a fixed four-bench layout (`det_sign_exact`, `det_exact`, +# `solve_exact`, `solve_exact_f64`). +_EXTREME_BENCHES: list[str] = ["det_sign_exact", "det_exact", "solve_exact", "solve_exact_f64"] + EXACT_GROUPS: dict[str, list[str]] = { "exact_d2": ["det", "det_direct", "det_exact", "det_exact_f64", "det_sign_exact", "solve_exact", "solve_exact_f64"], "exact_d3": ["det", "det_direct", "det_exact", "det_exact_f64", "det_sign_exact", "solve_exact", "solve_exact_f64"], "exact_d4": ["det", "det_direct", "det_exact", "det_exact_f64", "det_sign_exact", "solve_exact", "solve_exact_f64"], "exact_d5": ["det", "det_direct", "det_exact", "det_exact_f64", "det_sign_exact", "solve_exact", "solve_exact_f64"], - "exact_near_singular_3x3": ["det_sign_exact", "det_exact"], + "exact_near_singular_3x3": _EXTREME_BENCHES, + "exact_large_entries_3x3": _EXTREME_BENCHES, + "exact_hilbert_4x4": _EXTREME_BENCHES, + "exact_hilbert_5x5": _EXTREME_BENCHES, } @@ -197,11 +207,16 @@ def _group_by_group(items: list[_T]) -> dict[str, list[_T]]: def _group_heading(group: str) -> str: """Turn a Criterion group name into a readable heading.""" - # exact_d3 -> "D=3", exact_near_singular_3x3 -> "Near-singular 3x3" + # exact_d3 -> "D=3", exact_near_singular_3x3 -> "Near-singular 3x3", + # exact_hilbert_4x4 -> "Hilbert 4x4", etc. if group.startswith("exact_d"): return f"D={group.removeprefix('exact_d')}" if group == "exact_near_singular_3x3": return "Near-singular 3x3" + if group == "exact_large_entries_3x3": + return "Large entries 3x3" + if group.startswith("exact_hilbert_"): + return f"Hilbert {group.removeprefix('exact_hilbert_')}" return group diff --git a/scripts/tests/test_bench_compare.py b/scripts/tests/test_bench_compare.py index 9fd1c4d..f4a9944 100644 --- a/scripts/tests/test_bench_compare.py +++ b/scripts/tests/test_bench_compare.py @@ -77,6 +77,15 @@ def test_dimension_group(self) -> None: def test_near_singular(self) -> None: assert bench_compare._group_heading("exact_near_singular_3x3") == "Near-singular 3x3" + def test_large_entries(self) -> None: + assert bench_compare._group_heading("exact_large_entries_3x3") == "Large entries 3x3" + + def test_hilbert_4x4(self) -> None: + assert bench_compare._group_heading("exact_hilbert_4x4") == "Hilbert 4x4" + + def test_hilbert_5x5(self) -> None: + assert bench_compare._group_heading("exact_hilbert_5x5") == "Hilbert 5x5" + def test_unknown_passthrough(self) -> None: assert bench_compare._group_heading("something_else") == "something_else" diff --git a/src/exact.rs b/src/exact.rs index 356f70a..1645301 100644 --- a/src/exact.rs +++ b/src/exact.rs @@ -699,7 +699,9 @@ mod tests { use super::*; use crate::DEFAULT_PIVOT_TOL; + use num_traits::Signed; use pastey::paste; + use std::array::from_fn; // ----------------------------------------------------------------------- // Test helpers @@ -2006,6 +2008,164 @@ mod tests { gen_solve_exact_zero_rhs_tests!(4); gen_solve_exact_zero_rhs_tests!(5); + // ----------------------------------------------------------------------- + // Adversarial-input coverage mirroring `benches/exact.rs` + // ----------------------------------------------------------------------- + // + // These tests pin the behaviour of the extreme-input benchmark groups + // (`exact_near_singular_3x3`, `exact_large_entries_3x3`, + // `exact_hilbert_{4x4,5x5}`) so a regression would be caught even + // when benchmarks are not running. + + /// Multiply `A · x` entirely in `BigRational`, using `f64_to_bigrational` + /// to lift each matrix entry. Used by residual assertions for inputs + /// whose exact solution has no closed form we can easily type out. + fn bigrational_matvec(a: &Matrix, 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()) { + sum += f64_to_bigrational(*aij) * xj; + } + sum + }) + } + + /// Near-singular 3×3 solve (matches the `exact_near_singular_3x3` + /// bench). With `A = [[1+2^-50, 2, 3], [4, 5, 6], [7, 8, 9]]` and + /// `x0 = [1, 1, 1]`, `A · x0 = [6 + 2^-50, 15, 24]`; every component is + /// exactly representable in `f64` (`6` has ulp `2^-50` at its exponent). + /// `solve_exact` must recover `x0` exactly — the fractional denominator + /// introduced by `det(A) = -3 × 2^-50` cancels cleanly against the + /// augmented RHS. + #[test] + fn solve_exact_near_singular_3x3_integer_x0() { + let perturbation = f64::from_bits(0x3CD0_0000_0000_0000); // 2^-50 + let a = Matrix::<3>::from_rows([ + [1.0 + perturbation, 2.0, 3.0], + [4.0, 5.0, 6.0], + [7.0, 8.0, 9.0], + ]); + let b = Vector::<3>::new([6.0 + perturbation, 15.0, 24.0]); + let x = a.solve_exact(b).unwrap(); + let one = BigRational::from_integer(BigInt::from(1)); + assert_eq!(x[0], one); + assert_eq!(x[1], one); + assert_eq!(x[2], one); + } + + /// Large-entry 3×3 solve (matches the `exact_large_entries_3x3` + /// bench). `A = big · I + (1 - I)` with `big = f64::MAX / 2` and + /// `b = [big, 1, 1] = A · [1, 0, 0]`. The `BigInt` augmented system + /// sees entries of ~1023 bits on the diagonal and unit entries + /// elsewhere; Bareiss elimination still produces the exact integer + /// solution `[1, 0, 0]`. + #[test] + fn solve_exact_large_entries_3x3_unit_vector() { + let big = f64::MAX / 2.0; + assert!(big.is_finite()); + let a = Matrix::<3>::from_rows([[big, 1.0, 1.0], [1.0, big, 1.0], [1.0, 1.0, big]]); + let b = Vector::<3>::new([big, 1.0, 1.0]); + let x = a.solve_exact(b).unwrap(); + let zero = BigRational::from_integer(BigInt::from(0)); + let one = BigRational::from_integer(BigInt::from(1)); + assert_eq!(x[0], one); + assert_eq!(x[1], zero); + assert_eq!(x[2], zero); + } + + /// Determinant of the large-entry 3×3 is roughly `big^3`, which + /// overflows `f64`. `det_direct()` therefore returns `±∞`, the fast + /// filter inside `det_sign_exact` falls through on the `is_finite()` + /// guard, and the Bareiss fallback resolves the positive sign + /// correctly. `det_exact_f64` must report `Overflow`. + #[test] + fn det_sign_exact_large_entries_3x3_positive() { + let big = f64::MAX / 2.0; + let a = Matrix::<3>::from_rows([[big, 1.0, 1.0], [1.0, big, 1.0], [1.0, 1.0, big]]); + // Fast filter is inconclusive (big^3 overflows f64 to +∞), so + // this exercises the Bareiss cold path. + assert!(!a.det_direct().is_some_and(f64::is_finite)); + assert_eq!(a.det_sign_exact().unwrap(), 1); + // Cross-validate: the exact `BigRational` determinant must agree + // on sign with `det_sign_exact`, and `det_exact_f64` must overflow + // (the value is representable in BigRational but far exceeds f64). + assert!(a.det_exact().unwrap().is_positive()); + assert_eq!(a.det_exact_f64(), Err(LaError::Overflow { index: None })); + } + + /// Hilbert matrices are symmetric positive-definite, so + /// `det_sign_exact` must return `1` for every D. For D=2..=4 the + /// fast f64 filter resolves the positive sign without falling + /// through (Hilbert's determinant is tiny but still well above the + /// `det_errbound` cushion); for D=5 the filter is skipped entirely + /// and the Bareiss path handles inputs whose `(mantissa, exponent)` + /// pairs all differ. + macro_rules! gen_det_sign_exact_hilbert_positive_tests { + ($d:literal) => { + paste! { + #[test] + #[allow(clippy::cast_precision_loss)] + fn []() { + let mut rows = [[0.0f64; $d]; $d]; + for r in 0..$d { + for c in 0..$d { + rows[r][c] = 1.0 / ((r + c + 1) as f64); + } + } + let h = Matrix::<$d>::from_rows(rows); + assert_eq!(h.det_sign_exact().unwrap(), 1); + } + } + }; + } + + gen_det_sign_exact_hilbert_positive_tests!(2); + gen_det_sign_exact_hilbert_positive_tests!(3); + gen_det_sign_exact_hilbert_positive_tests!(4); + gen_det_sign_exact_hilbert_positive_tests!(5); + + /// `solve_exact` on a Hilbert matrix must produce a solution whose + /// residual `A · x - b` is *exactly* zero in `BigRational` arithmetic. + /// Hilbert entries (`1/3`, `1/5`, `1/6`, `1/7`, …) are non-terminating + /// in binary, so this is a stronger test than the + /// `gen_solve_exact_roundtrip_tests` construction (which requires the + /// RHS to be representable as an exact `f64` product). + macro_rules! gen_solve_exact_hilbert_residual_tests { + ($d:literal) => { + paste! { + #[test] + #[allow(clippy::cast_precision_loss)] + fn []() { + let mut rows = [[0.0f64; $d]; $d]; + for r in 0..$d { + for c in 0..$d { + rows[r][c] = 1.0 / ((r + c + 1) as f64); + } + } + let h = Matrix::<$d>::from_rows(rows); + // Use a non-trivial RHS with both positive and negative + // entries to avoid accidental structural cancellation. + let mut b_arr = [0.0f64; $d]; + for i in 0usize..$d { + let sign = if i.is_multiple_of(2) { 1.0 } else { -1.0 }; + b_arr[i] = sign * ((i + 1) as f64); + } + let b = Vector::<$d>::new(b_arr); + let x = h.solve_exact(b).unwrap(); + let ax = bigrational_matvec(&h, &x); + for i in 0..$d { + assert_eq!(ax[i], f64_to_bigrational(b_arr[i])); + } + } + } + }; + } + + gen_solve_exact_hilbert_residual_tests!(2); + gen_solve_exact_hilbert_residual_tests!(3); + gen_solve_exact_hilbert_residual_tests!(4); + gen_solve_exact_hilbert_residual_tests!(5); + // ----------------------------------------------------------------------- // solve_exact_f64: dimension-specific tests // ----------------------------------------------------------------------- diff --git a/src/lib.rs b/src/lib.rs index 8b9eb38..c1f9a5c 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -51,7 +51,11 @@ mod readme_doctests { #[cfg(feature = "exact")] mod exact; #[cfg(feature = "exact")] +pub use num_bigint::BigInt; +#[cfg(feature = "exact")] pub use num_rational::BigRational; +#[cfg(feature = "exact")] +pub use num_traits::{FromPrimitive, Signed, ToPrimitive}; mod ldlt; mod lu; @@ -297,8 +301,14 @@ pub use vector::Vector; /// [`Ldlt`], [`LaError`], [`DEFAULT_PIVOT_TOL`], [`DEFAULT_SINGULAR_TOL`], and the determinant /// error bound coefficients [`ERR_COEFF_2`], [`ERR_COEFF_3`], and [`ERR_COEFF_4`]. /// -/// When the `exact` feature is enabled, `BigRational` is also -/// re-exported for use with `Matrix::det_exact`. +/// When the `exact` feature is enabled, [`BigInt`] and [`BigRational`] +/// are also re-exported so callers can construct exact values (e.g. as +/// the expected result of `Matrix::det_exact`) without adding +/// `num-bigint` / `num-rational` to their own dependencies. The most +/// commonly needed `num-traits` items are re-exported alongside them: +/// [`FromPrimitive`] for `BigRational::from_f64` / `from_i64`, +/// [`ToPrimitive`] for `BigRational::to_f64` / `to_i64`, and [`Signed`] +/// 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, @@ -306,7 +316,7 @@ pub mod prelude { }; #[cfg(feature = "exact")] - pub use crate::BigRational; + pub use crate::{BigInt, BigRational, FromPrimitive, Signed, ToPrimitive}; } #[cfg(test)] @@ -377,4 +387,37 @@ mod tests { let _ = m.lu(DEFAULT_PIVOT_TOL).unwrap().solve_vec(v).unwrap(); let _ = m.ldlt(DEFAULT_SINGULAR_TOL).unwrap().solve_vec(v).unwrap(); } + + /// Exercise every exact-feature re-export via the prelude so a future + /// refactor that drops one (e.g. removing `Signed` from the prelude + /// list) fails to compile rather than silently breaking downstream. + #[cfg(feature = "exact")] + #[test] + fn prelude_exact_reexports_compile_and_work() { + use crate::prelude::*; + + // `BigInt` and `BigRational` constructors. + let n = BigInt::from(7); + let r = BigRational::from_integer(n.clone()); + 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(); + assert_eq!( + half.clone() + half.clone(), + BigRational::from_integer(BigInt::from(1)) + ); + + // `Signed::is_positive` / `is_negative` / `abs`. + assert!(half.is_positive()); + assert!(!half.is_negative()); + let neg = -half.clone(); + assert!(neg.is_negative()); + 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); + } } diff --git a/tests/proptest_exact.rs b/tests/proptest_exact.rs index 2998b8e..764ec02 100644 --- a/tests/proptest_exact.rs +++ b/tests/proptest_exact.rs @@ -1,7 +1,17 @@ -//! Property-based tests for the `det_sign_exact` API (requires `exact` feature). +//! Property-based tests for the exact-arithmetic APIs +//! (requires `exact` feature). +//! +//! Covers: +//! - `det_sign_exact` on diagonal and full small-integer matrices +//! - `solve_exact` round-trip with integer inputs (`A · x0` in f64 is +//! exact for small integers, so `solve(A, A · x0) == x0`) +//! - `solve_exact` residual property (`A · solve(A, b) == b` in +//! `BigRational` arithmetic) on random RHS vectors #![cfg(feature = "exact")] +use std::array::from_fn; + use pastey::paste; use proptest::prelude::*; @@ -11,6 +21,67 @@ fn small_nonzero_f64() -> impl Strategy { prop_oneof![(-1000i16..=-1i16), (1i16..=1000i16)].prop_map(|x| f64::from(x) / 10.0) } +/// Small non-zero integers in `[-5, 5] \ {0}`, returned as `f64`. +/// +/// Used to populate the off-diagonal of diagonally-dominant integer +/// matrices in the `solve_exact` proptests. Every value is an exact +/// `f64` integer, so `A · x` never loses precision for small `x`. +fn small_nonzero_int_f64() -> impl Strategy { + prop_oneof![(-5i32..=-1i32), (1i32..=5i32)].prop_map(f64::from) +} + +/// Small signed integers in `[-10, 10]` as `f64` (may be zero). +fn small_int_f64() -> impl Strategy { + (-10i32..=10i32).prop_map(f64::from) +} + +/// Multiply `A · x` entirely in `BigRational`, lifting each f64 matrix +/// entry via `BigRational::from_f64`. Used by residual assertions. +/// +/// All f64 inputs in the proptests are small exact integers, so +/// `from_f64` always succeeds with an exact rational reconstruction. +fn bigrational_matvec(a: &[[f64; D]; D], x: &[BigRational; D]) -> [BigRational; D] { + from_fn(|i| { + let mut sum = BigRational::from_integer(BigInt::from(0)); + for (aij, xj) in a[i].iter().zip(x.iter()) { + let entry = BigRational::from_f64(*aij).expect("small int fits in BigRational"); + sum += entry * xj; + } + sum + }) +} + +/// Build a strictly diagonally-dominant f64 matrix from: +/// - an off-diagonal matrix of small integers (entries in `[-10, 10]` +/// per `small_int_f64`), and +/// - diagonal entries shifted by `D · 10 + 1` so every row satisfies +/// `|A[i][i]| > Σ_{j≠i} |A[i][j]|`, which guarantees invertibility +/// (Levy–Desplanques). +/// +/// The shift must match the off-diagonal strategy's maximum magnitude +/// (`max_off_diag = 10`): with `D - 1` off-diagonals of magnitude ≤ 10 +/// the row sum is at most `10 (D - 1) < 10 D + 1`, so the shifted +/// diagonal strictly dominates. The shift keeps every entry a small +/// exact `f64` integer, so matrix × small-integer-vector products are +/// exact in f64. +fn make_diagonally_dominant( + offdiag: [[f64; D]; D], + diag: [f64; D], +) -> [[f64; D]; D] { + let mut rows = offdiag; + // Must track `small_int_f64`'s `max_off_diag = 10`: `D · 10 + 1` + // strictly dominates the worst-case row sum of `10 (D - 1)`. + let shift = f64::from(u8::try_from(D).unwrap_or(u8::MAX)).mul_add(10.0, 1.0); + for i in 0..D { + rows[i][i] = if diag[i] >= 0.0 { + diag[i] + shift + } else { + diag[i] - shift + }; + } + rows +} + macro_rules! gen_det_sign_exact_proptests { ($d:literal) => { paste! { @@ -70,3 +141,189 @@ gen_det_sign_exact_proptests!(2); gen_det_sign_exact_proptests!(3); gen_det_sign_exact_proptests!(4); gen_det_sign_exact_proptests!(5); + +/// Round-trip property: for random small-integer `x0` and a random +/// diagonally-dominant integer matrix `A`, `A · x0` is exactly +/// representable in `f64` (small integer products stay well under the +/// 53-bit mantissa) and `solve_exact(A, A · x0)` must return `x0` +/// exactly as `BigRational`. This exercises the full Bareiss forward +/// elimination + rational back-substitution pipeline on a different +/// input for every case. +macro_rules! gen_solve_exact_roundtrip_proptests { + ($d:literal) => { + paste! { + proptest! { + #![proptest_config(ProptestConfig::with_cases(64))] + + #[test] + fn []( + offdiag in proptest::array::[]( + proptest::array::[](small_int_f64()), + ), + diag in proptest::array::[](small_nonzero_int_f64()), + x0 in proptest::array::[](small_int_f64()), + ) { + let rows = make_diagonally_dominant::<$d>(offdiag, diag); + let a = Matrix::<$d>::from_rows(rows); + + // b = A · x0, computed in f64. Small integers keep + // every partial sum exact. + let mut b_arr = [0.0f64; $d]; + for i in 0..$d { + let mut sum = 0.0f64; + for j in 0..$d { + sum += rows[i][j] * x0[j]; + } + b_arr[i] = sum; + } + let b = Vector::<$d>::new(b_arr); + let x = a.solve_exact(b).expect("diagonally-dominant A is non-singular"); + + let expected: [BigRational; $d] = from_fn(|i| { + BigRational::from_f64(x0[i]).expect("small int fits in BigRational") + }); + for i in 0..$d { + prop_assert_eq!(&x[i], &expected[i]); + } + } + } + } + }; +} + +gen_solve_exact_roundtrip_proptests!(2); +gen_solve_exact_roundtrip_proptests!(3); +gen_solve_exact_roundtrip_proptests!(4); +gen_solve_exact_roundtrip_proptests!(5); + +/// Residual property: for a random diagonally-dominant integer matrix +/// `A` and a random integer RHS `b`, `solve_exact` must return an `x` +/// such that `A · x` equals `b` *exactly* in `BigRational` +/// arithmetic. Unlike the round-trip test above, the exact solution +/// is generally fractional — this catches back-substitution bugs that +/// preserve integer inputs but mishandle denominators. +macro_rules! gen_solve_exact_residual_proptests { + ($d:literal) => { + paste! { + proptest! { + #![proptest_config(ProptestConfig::with_cases(32))] + + #[test] + fn []( + offdiag in proptest::array::[]( + proptest::array::[](small_int_f64()), + ), + diag in proptest::array::[](small_nonzero_int_f64()), + b_arr in proptest::array::[](small_int_f64()), + ) { + let rows = make_diagonally_dominant::<$d>(offdiag, diag); + let a = Matrix::<$d>::from_rows(rows); + let b = Vector::<$d>::new(b_arr); + let x = a.solve_exact(b).expect("diagonally-dominant A is non-singular"); + + let ax = bigrational_matvec::<$d>(&rows, &x); + for i in 0..$d { + let b_rat = BigRational::from_f64(b_arr[i]) + .expect("small int fits in BigRational"); + prop_assert_eq!(&ax[i], &b_rat); + } + } + } + } + }; +} + +gen_solve_exact_residual_proptests!(2); +gen_solve_exact_residual_proptests!(3); +gen_solve_exact_residual_proptests!(4); +gen_solve_exact_residual_proptests!(5); + +/// On full (non-diagonal) random small-integer matrices, +/// `det_sign_exact()` must agree with `det_exact().signum()`. This +/// exercises the adaptive fast-filter / Bareiss-fallback boundary on +/// inputs the existing diagonal-only proptests don't touch (e.g. +/// matrices where the f64 det is near its error bound and the filter +/// must defer to Bareiss). +macro_rules! gen_det_sign_agrees_with_det_exact_proptests { + ($d:literal) => { + paste! { + proptest! { + #![proptest_config(ProptestConfig::with_cases(64))] + + #[test] + fn []( + entries in proptest::array::[]( + proptest::array::[](small_int_f64()), + ), + ) { + let m = Matrix::<$d>::from_rows(entries); + let sign = m.det_sign_exact().unwrap(); + let det = m.det_exact().unwrap(); + let expected: i8 = if det.is_positive() { + 1 + } else if det.is_negative() { + -1 + } else { + 0 + }; + prop_assert_eq!(sign, expected); + } + } + } + }; +} + +gen_det_sign_agrees_with_det_exact_proptests!(2); +gen_det_sign_agrees_with_det_exact_proptests!(3); +gen_det_sign_agrees_with_det_exact_proptests!(4); +gen_det_sign_agrees_with_det_exact_proptests!(5); + +/// Fast-filter invariant: whenever `|det_direct()| > det_errbound()`, +/// the f64 sign is provably correct — so +/// `det_direct().signum() == det_sign_exact()`. This is the +/// correctness guarantee the Shewchuk-style filter inside +/// `det_sign_exact` relies on. The proptest cross-checks that the +/// fast-filter boundary itself is honoured, independent of whether +/// `det_sign_exact` ended up using the filter or the Bareiss fallback +/// on any particular input. Only D=2..=4 have a closed-form +/// `det_direct` / `det_errbound` pair. +macro_rules! gen_det_sign_fast_filter_boundary_proptests { + ($d:literal) => { + paste! { + proptest! { + #![proptest_config(ProptestConfig::with_cases(64))] + + #[test] + fn []( + entries in proptest::array::[]( + proptest::array::[](small_int_f64()), + ), + ) { + let m = Matrix::<$d>::from_rows(entries); + let det = m.det_direct().expect("D<=4 has closed-form det_direct"); + let bound = m.det_errbound().expect("D<=4 has a det_errbound"); + let sign = m.det_sign_exact().unwrap(); + + // Only assert when the filter is conclusive. When + // `|det| <= bound` the f64 sign may disagree with the + // exact sign; that case is covered by the other + // proptests via the Bareiss fallback. + if det.abs() > bound { + let direct_sign: i8 = if det > 0.0 { + 1 + } else if det < 0.0 { + -1 + } else { + 0 + }; + prop_assert_eq!(direct_sign, sign); + } + } + } + } + }; +} + +gen_det_sign_fast_filter_boundary_proptests!(2); +gen_det_sign_fast_filter_boundary_proptests!(3); +gen_det_sign_fast_filter_boundary_proptests!(4);