diff --git a/REFERENCES.md b/REFERENCES.md index bd850c7..02fcf49 100644 --- a/REFERENCES.md +++ b/REFERENCES.md @@ -40,14 +40,27 @@ f64 entry is decomposed into `mantissa × 2^exponent`, scaled to a common intege eliminated without any `BigRational` or GCD overhead. See `src/exact.rs` for the full architecture description. -### f64 → BigRational conversion (`f64_to_bigrational`) - -`f64_to_bigrational` converts an f64 to an exact `BigRational` by decomposing the IEEE 754 -binary64 bit representation [9] into its sign, exponent, and significand fields. Because -every finite f64 is exactly `±m × 2^e` (where `m` is an integer), the rational can be -constructed directly via `BigRational::new_raw` without GCD normalization — trailing zeros -in the significand are stripped first so the fraction is already in lowest terms. See -Goldberg [10] for background on IEEE 754 representation and exact rational reconstruction. +### Exact linear system solve (hybrid Bareiss / BigRational) + +`solve_exact()` and `solve_exact_f64()` share the BigInt core used for determinants. Matrix +and RHS entries are decomposed via IEEE 754 bit extraction [9] and scaled to a shared base +`2^e_min` so the augmented system `(A | b)` becomes a `BigInt` matrix. Forward elimination +runs in `BigInt` using Bareiss fraction-free updates [7] — no `BigRational` and no GCD +normalisation in the `O(D³)` phase. The upper-triangular result is then lifted into +`BigRational` for back-substitution, where fractions are inherent and the cost is only +`O(D²)`. Row swaps from first-non-zero pivoting are applied to both the matrix and the +RHS; because power-of-two scaling is applied uniformly to both sides of `A x = b`, the +solution is unchanged by the scale factor. + +### f64 → integer decomposition (`f64_decompose`) + +Both the determinant and solve paths convert f64 entries via `f64_decompose`, which extracts +the IEEE 754 binary64 sign, unbiased exponent, and significand [9] and strips trailing zeros +from the significand so `|x| = m · 2^e` with `m` odd. The integer matrix is then assembled +by shifting each mantissa left by `exp − e_min`, giving a GCD-free, Bareiss-ready starting +point. A one-shot wrapper `f64_to_bigrational` (used only in tests) packages the same +decomposition into a single `BigRational`. See Goldberg [10] for background on IEEE 754 +representation and exact rational reconstruction. ### LDL^T factorization (symmetric SPD/PSD) diff --git a/src/exact.rs b/src/exact.rs index 83d9562..356f70a 100644 --- a/src/exact.rs +++ b/src/exact.rs @@ -28,22 +28,48 @@ //! //! ## Linear system solve //! -//! `solve_exact` and `solve_exact_f64` solve `A x = b` using Gaussian -//! elimination with first-non-zero pivoting in `BigRational` arithmetic. -//! Since all arithmetic is exact, any non-zero pivot gives the correct result -//! (there is no numerical stability concern). Every finite `f64` is exactly -//! representable as a rational, so the result is provably correct. +//! `solve_exact` and `solve_exact_f64` solve `A x = b` with a hybrid +//! algorithm that reuses the integer-only Bareiss core used for +//! determinants. Matrix and RHS entries are decomposed via +//! `f64_decompose` into `mantissa × 2^exponent`, scaled to a shared +//! base `2^e_min`, and assembled into a `BigInt` augmented system +//! `(A | b)`. Forward elimination runs entirely in `BigInt` with +//! fraction-free Bareiss updates — no `BigRational`, no GCD +//! normalisation in the `O(D³)` phase. Once the system is upper +//! triangular, back-substitution is performed in `BigRational`, where +//! fractions are inherent; this phase is only `O(D²)` so the rational +//! overhead is modest. First-non-zero pivoting is used throughout; +//! since all arithmetic is exact, any non-zero pivot gives the correct +//! result (no numerical stability concern). Every finite `f64` is +//! exactly representable as a rational, so the result is provably +//! correct. //! -//! ## f64 → `BigRational` conversion +//! ## f64 → integer decomposition //! -//! All entry conversions use `f64_to_bigrational`, which decomposes the -//! IEEE 754 binary64 bit representation (\[9\]) into sign, exponent, and -//! significand and constructs a `BigRational` directly — avoiding the GCD -//! normalization that `BigRational::from_float` performs. See Goldberg -//! \[10\] for background on floating-point representation and exact -//! rational reconstruction. Reference numbers refer to `REFERENCES.md`. +//! Both the determinant and solve paths share a single conversion +//! primitive, `f64_decompose`, which extracts `(mantissa, exponent, +//! sign)` from the IEEE 754 binary64 bit representation (\[9\]). The +//! determinant path combines those components into a `BigInt` matrix +//! (for Bareiss) and a `2^(D × e_min)` scale factor, while the solve +//! path builds a `BigInt` augmented system and lifts the +//! upper-triangular result into `BigRational` for back-substitution. +//! See Goldberg \[10\] for background on floating-point representation +//! and exact rational reconstruction. Reference numbers refer to +//! `REFERENCES.md`. +//! +//! ## Validation +//! +//! `decompose_matrix` / `decompose_vec` fold an `is_finite()` check +//! into the same pass that decomposes each entry, returning +//! `Err(LaError::NonFinite { row, col })` on the first NaN or ±∞ +//! encountered. This error is propagated through `bareiss_det_int`, +//! `bareiss_det`, and `gauss_solve` via the `?` operator, so every +//! public entry point that reaches the integer-Bareiss core is +//! automatically validated — `f64_decompose` itself is therefore +//! never called with non-finite input from the public API. use core::hint::cold_path; +use core::mem::take; use std::array::from_fn; use num_bigint::{BigInt, Sign}; @@ -54,36 +80,6 @@ use crate::LaError; use crate::matrix::Matrix; use crate::vector::Vector; -/// Validate that all entries in a `D×D` matrix are finite (not NaN or infinite). -/// -/// Returns `Ok(())` if all entries are finite, or `Err(LaError::NonFinite)` with -/// the column of the first non-finite entry found. -fn validate_finite(m: &Matrix) -> Result<(), LaError> { - for r in 0..D { - for c in 0..D { - if !m.rows[r][c].is_finite() { - cold_path(); - return Err(LaError::non_finite_cell(r, c)); - } - } - } - Ok(()) -} - -/// Validate that all entries in a length-`D` vector are finite. -/// -/// Returns `Ok(())` if all entries are finite, or `Err(LaError::NonFinite)` with -/// the index of the first non-finite entry found. -fn validate_finite_vec(v: &Vector) -> Result<(), LaError> { - for (i, &x) in v.data.iter().enumerate() { - if !x.is_finite() { - cold_path(); - return Err(LaError::non_finite_at(i)); - } - } - Ok(()) -} - /// Decompose a finite `f64` into its IEEE 754 components. /// /// Returns `None` for ±0.0, or `Some((mantissa, exponent, is_negative))` where @@ -124,38 +120,6 @@ fn f64_decompose(x: f64) -> Option<(u64, i32, bool)> { Some((mantissa, exponent, is_negative)) } -/// Convert an `f64` to an exact `BigRational` via IEEE 754 bit decomposition. -/// -/// Every finite `f64` is exactly representable as `±m × 2^e` where `m` is a -/// non-negative integer and `e` is an integer. This function extracts `(m, e)` -/// directly from the IEEE 754 binary64 bit layout \[9\], strips trailing zeros -/// from `m` so the resulting fraction is already in lowest terms, then -/// constructs a `BigRational` via `new_raw` — bypassing the GCD reduction -/// that `BigRational::from_float` performs internally. -/// -/// See `REFERENCES.md` \[9-10\] for the IEEE 754 standard and Goldberg's -/// survey of floating-point representation. -/// -/// # Panics -/// Panics if `x` is NaN or infinite. -fn f64_to_bigrational(x: f64) -> BigRational { - let Some((mantissa, exponent, is_negative)) = f64_decompose(x) else { - return BigRational::from_integer(BigInt::from(0)); - }; - - let numer = if is_negative { - -BigInt::from(mantissa) - } else { - BigInt::from(mantissa) - }; - - if exponent >= 0 { - BigRational::new_raw(numer << exponent.cast_unsigned(), BigInt::from(1u32)) - } else { - BigRational::new_raw(numer, BigInt::from(1u32) << (-exponent).cast_unsigned()) - } -} - /// Convert a `BigInt × 2^exp` pair to a reduced `BigRational`. /// /// When `exp < 0` (denominator is `2^(-exp)`), shared factors of 2 are @@ -182,101 +146,237 @@ fn bigint_exp_to_bigrational(mut value: BigInt, mut exp: i32) -> BigRational { } } -/// Compute the exact determinant using integer-only Bareiss elimination. +// ----------------------------------------------------------------------- +// Shared integer-Bareiss primitives +// ----------------------------------------------------------------------- +// +// Both `bareiss_det_int` (determinants) and `gauss_solve` (linear system +// solve) follow the same pipeline: decompose every f64 entry into +// `(mantissa, exponent, is_negative)`, track the minimum exponent across +// non-zero entries, scale each entry by `2^(exp − e_min)` to build a +// fully-integer `BigInt` matrix, and run Bareiss fraction-free forward +// elimination. The helpers below factor out each stage so the two +// callers differ only in post-processing (± sign for det, back-sub for +// solve) and in whether they carry a RHS through the elimination. + +/// Decomposed finite f64 in the form `(-1)^is_negative · mantissa · 2^exponent`. /// -/// Returns `(det_int, scale_exp)` where the true determinant is -/// `det_int × 2^scale_exp`. Since the scale factor `2^scale_exp` is always -/// positive, `det_int.sign()` gives the sign of the determinant directly. +/// Zero entries have `mantissa == 0`; the other fields are unused in that +/// case. `Default` yields such a zero component, which is what the +/// per-entry initialiser in `decompose_matrix` / `decompose_vec` produces +/// for ±0.0 cells. +#[derive(Clone, Copy, Default)] +struct Component { + mantissa: u64, + exponent: i32, + is_negative: bool, +} + +/// Decompose every entry of a `D×D` matrix via `f64_decompose`, +/// validating finiteness in the same pass. Returns the per-entry +/// components and the minimum exponent across non-zero entries. If +/// every entry is zero, the exponent is `i32::MAX`. /// -/// All arithmetic is in `BigInt` — no `BigRational`, no GCD, no denominator -/// tracking. Each f64 entry is decomposed into `mantissa × 2^exponent` and -/// scaled to a common base `2^e_min` so every entry becomes an integer. -/// The Bareiss inner-loop division is exact (guaranteed by the algorithm). -fn bareiss_det_int(m: &Matrix) -> (BigInt, i32) { - if D == 0 { - return (BigInt::from(1), 0); - } - if D == 1 { - return match f64_decompose(m.rows[0][0]) { - None => (BigInt::from(0), 0), - Some((mant, exp, neg)) => { - let v = if neg { - -BigInt::from(mant) - } else { - BigInt::from(mant) +/// # Errors +/// Returns [`LaError::NonFinite`] with `row: Some(r), col: c` pointing +/// at the first non-finite entry encountered (row-major order). +fn decompose_matrix(m: &Matrix) -> Result<([[Component; D]; D], i32), LaError> { + let mut components = [[Component::default(); D]; D]; + let mut e_min = i32::MAX; + for (r, row) in m.rows.iter().enumerate() { + for (c, &entry) in row.iter().enumerate() { + if !entry.is_finite() { + cold_path(); + return Err(LaError::non_finite_cell(r, c)); + } + if let Some((mantissa, exponent, is_negative)) = f64_decompose(entry) { + components[r][c] = Component { + mantissa, + exponent, + is_negative, }; - (v, exp) + e_min = e_min.min(exponent); } - }; + } } + Ok((components, e_min)) +} - // Decompose all entries and find the minimum exponent. - let mut components = [[(0u64, 0i32, false); D]; D]; +/// Decompose every entry of a length-`D` vector via `f64_decompose`, +/// validating finiteness in the same pass. Returns the per-entry +/// components and the minimum exponent across non-zero entries. If +/// every entry is zero, the exponent is `i32::MAX`. +/// +/// # Errors +/// Returns [`LaError::NonFinite`] with `row: None, col: i` pointing at +/// the first non-finite entry encountered. +fn decompose_vec(v: &Vector) -> Result<([Component; D], i32), LaError> { + let mut components = [Component::default(); D]; let mut e_min = i32::MAX; - - for (r, row) in m.rows.iter().enumerate() { - for (c, &entry) in row.iter().enumerate() { - if let Some((mant, exp, neg)) = f64_decompose(entry) { - components[r][c] = (mant, exp, neg); - e_min = e_min.min(exp); - } - // Zero entries keep the default (0, 0, false); their exponent is - // excluded from e_min. + for (i, &entry) in v.data.iter().enumerate() { + if !entry.is_finite() { + cold_path(); + return Err(LaError::non_finite_at(i)); + } + if let Some((mantissa, exponent, is_negative)) = f64_decompose(entry) { + components[i] = Component { + mantissa, + exponent, + is_negative, + }; + e_min = e_min.min(exponent); } } + Ok((components, e_min)) +} - // All entries are zero → singular. - if e_min == i32::MAX { - return (BigInt::from(0), 0); - } - - // Build the integer matrix: a[r][c] = (±mantissa) × 2^(exp − e_min). - let mut a: [[BigInt; D]; D] = from_fn(|r| { - from_fn(|c| { - let (mant, exp, neg) = components[r][c]; - if mant == 0 { - BigInt::from(0) - } else { - let shift = (exp - e_min).cast_unsigned(); - let v = BigInt::from(mant) << shift; - if neg { -v } else { v } - } - }) - }); +/// Convert a single decomposed component to its scaled `BigInt` +/// representation: `(±mantissa) << (exp − e_min)`. Zero components map +/// to `BigInt::from(0)`. +#[inline] +fn component_to_bigint(c: Component, e_min: i32) -> BigInt { + if c.mantissa == 0 { + BigInt::from(0) + } else { + let v = BigInt::from(c.mantissa) << (c.exponent - e_min).cast_unsigned(); + if c.is_negative { -v } else { v } + } +} - // Bareiss elimination in BigInt. +/// Build a `D×D` integer matrix from a component table, scaled to the +/// shared base `2^e_min`. +fn build_bigint_matrix( + components: &[[Component; D]; D], + e_min: i32, +) -> [[BigInt; D]; D] { + from_fn(|r| from_fn(|c| component_to_bigint(components[r][c], e_min))) +} + +/// Build a length-`D` integer vector from a component array, scaled to +/// the shared base `2^e_min`. +fn build_bigint_vec(components: &[Component; D], e_min: i32) -> [BigInt; D] { + from_fn(|i| component_to_bigint(components[i], e_min)) +} + +/// Outcome of a Bareiss forward-elimination pass. +#[derive(Debug)] +enum BareissResult { + /// Elimination completed; `sign` is `±1` based on the parity of row + /// swaps (relevant for determinants; solves discard it). + Upper { sign: i8 }, + /// Column `pivot_col` has no non-zero pivot at or below its diagonal. + Singular { pivot_col: usize }, +} + +/// Run Bareiss fraction-free forward elimination on the `D×D` integer +/// matrix `a`, optionally augmented with a length-`D` RHS vector. +/// +/// When `rhs` is `Some`, row swaps and the inner-loop Bareiss update are +/// mirrored on the RHS (treating it as column `D+1` of an augmented +/// system). On return, `a` is upper triangular and the last pivot lives +/// in `a[D-1][D-1]`. +/// +/// First-non-zero pivoting is used: since all arithmetic is exact, any +/// non-zero pivot is valid — no tolerance is required. +fn bareiss_forward_eliminate( + a: &mut [[BigInt; D]; D], + mut rhs: Option<&mut [BigInt; D]>, +) -> BareissResult { let zero = BigInt::from(0); let mut prev_pivot = BigInt::from(1); let mut sign: i8 = 1; for k in 0..D { - // Pivot search. + // First-non-zero pivot search. if a[k][k] == zero { let mut found = false; for i in (k + 1)..D { if a[i][k] != zero { a.swap(k, i); + if let Some(r) = &mut rhs { + r.swap(k, i); + } sign = -sign; found = true; break; } } if !found { - return (BigInt::from(0), 0); + cold_path(); + return BareissResult::Singular { pivot_col: k }; } } - // Elimination. + // Elimination. The Bareiss update reads the current `a[i][k]` + // in both the inner `j`-loop and the RHS update, so zero it only + // *after* those reads. 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; } + if let Some(r) = &mut rhs { + r[i] = (&a[k][k] * &r[i] - &a[i][k] * &r[k]) / &prev_pivot; + } a[i][k].clone_from(&zero); } prev_pivot.clone_from(&a[k][k]); } + // Post-conditions (debug builds only): `a` is upper triangular with + // non-zero pivots. These catch future regressions in the inner-loop + // update or pivot-search logic without runtime cost in release. + // Indexed iteration is clearer than iterator chains here because the + // checks read disjoint cells across rows and columns at each step. + #[cfg(debug_assertions)] + #[allow(clippy::needless_range_loop)] + for k in 0..D { + assert_ne!(a[k][k], zero, "pivot at ({k}, {k}) must be non-zero"); + for i in (k + 1)..D { + assert_eq!(a[i][k], zero, "sub-diagonal at ({i}, {k}) must be zero"); + } + } + + BareissResult::Upper { sign } +} + +/// Compute the exact determinant using integer-only Bareiss elimination. +/// +/// Returns `(det_int, scale_exp)` where the true determinant is +/// `det_int × 2^scale_exp`. Since the scale factor `2^scale_exp` is always +/// positive, `det_int.sign()` gives the sign of the determinant directly. +/// +/// All arithmetic is in `BigInt` — no `BigRational`, no GCD, no denominator +/// tracking. Each f64 entry is decomposed into `mantissa × 2^exponent` and +/// scaled to a common base `2^e_min` so every entry becomes an integer. +/// The Bareiss inner-loop division is exact (guaranteed by the algorithm). +/// +/// # Errors +/// Returns [`LaError::NonFinite`] (propagated from `decompose_matrix`) if +/// any matrix entry is NaN or infinite. +fn bareiss_det_int(m: &Matrix) -> Result<(BigInt, i32), LaError> { + // D == 0 has no `a[D-1][D-1]` to read; shortcut to the empty-product + // determinant. + if D == 0 { + return Ok((BigInt::from(1), 0)); + } + + let (components, e_min) = decompose_matrix(m)?; + + // All entries are zero → singular (det = 0). + if e_min == i32::MAX { + return Ok((BigInt::from(0), 0)); + } + + let mut a = build_bigint_matrix(&components, e_min); + let sign = match bareiss_forward_eliminate(&mut a, None) { + BareissResult::Upper { sign } => sign, + BareissResult::Singular { .. } => { + cold_path(); + return Ok((BigInt::from(0), 0)); + } + }; + let det_int = if sign < 0 { -&a[D - 1][D - 1] } else { @@ -289,73 +389,81 @@ fn bareiss_det_int(m: &Matrix) -> (BigInt, i32) { .checked_mul(d_i32) .expect("exponent overflow in bareiss_det_int"); - (det_int, total_exp) + Ok((det_int, total_exp)) } /// Compute the exact determinant of a `D×D` matrix using integer-only Bareiss /// elimination and return the result as a `BigRational`. -fn bareiss_det(m: &Matrix) -> BigRational { - let (det_int, total_exp) = bareiss_det_int(m); - bigint_exp_to_bigrational(det_int, total_exp) +/// +/// # Errors +/// Returns [`LaError::NonFinite`] if any matrix entry is NaN or infinite. +fn bareiss_det(m: &Matrix) -> Result { + let (det_int, total_exp) = bareiss_det_int(m)?; + Ok(bigint_exp_to_bigrational(det_int, total_exp)) } -/// Solve `A x = b` using Gaussian elimination with first-non-zero pivoting -/// in `BigRational` arithmetic. +/// Solve `A x = b` using a hybrid BigInt/BigRational algorithm. /// -/// Since all arithmetic is exact, any non-zero pivot gives the correct result -/// (no numerical stability concern). This matches the pivoting strategy used -/// by `bareiss_det`. +/// Forward elimination runs entirely in `BigInt` using fraction-free +/// Bareiss updates on the augmented system `(A | b)`: every f64 entry +/// is decomposed into `mantissa × 2^exponent` and scaled to a shared +/// base `2^e_min` so both the matrix and the RHS become integer. +/// Because the same power-of-two scaling is applied to both sides of +/// `A x = b`, the solution is unchanged. Row swaps also swap the RHS +/// row; no sign tracking is needed (pivot permutations do not affect +/// the solution of a linear system). +/// +/// After forward elimination, the upper-triangular `BigInt` system and +/// its RHS are lifted into `BigRational` for back-substitution, where +/// fractions are inherent. This keeps the expensive `O(D³)` phase +/// GCD-free and limits `BigRational` work to the cheaper `O(D²)` phase. /// /// Returns the exact solution as `[BigRational; D]`. -/// Returns `Err(LaError::Singular)` if the matrix is exactly singular. +/// +/// # Errors +/// Returns [`LaError::NonFinite`] (propagated from `decompose_matrix` / +/// `decompose_vec`) if any matrix or vector entry is NaN or infinite. +/// The matrix is validated before the vector, matching public-API order. +/// Returns [`LaError::Singular`] if the matrix is exactly singular. fn gauss_solve(m: &Matrix, b: &Vector) -> Result<[BigRational; D], LaError> { - let zero = BigRational::from_integer(BigInt::from(0)); + // Decompose both matrix and RHS (validating finiteness in one pass); + // the shared minimum exponent makes every entry of the augmented + // system an integer after scaling. + let (m_components, m_e_min) = decompose_matrix(m)?; + let (b_components, b_e_min) = decompose_vec(b)?; + let mut e_min = m_e_min.min(b_e_min); + + // All matrix + RHS entries are zero. For `D > 0` this surfaces as + // singular inside forward elimination; for `D == 0` the elimination + // loop body is empty and we return `Ok([])` without touching e_min. + // Pick any finite value so the shift computation is well-defined (the + // resulting BigInts are all zero either way). + if e_min == i32::MAX { + e_min = 0; + } - // Build matrix and RHS separately (cannot use [BigRational; D+1] augmented - // columns because const-generic expressions are unstable). - let mut mat: [[BigRational; D]; D] = from_fn(|r| from_fn(|c| f64_to_bigrational(m.rows[r][c]))); - let mut rhs: [BigRational; D] = from_fn(|r| f64_to_bigrational(b.data[r])); + let mut a = build_bigint_matrix(&m_components, e_min); + let mut rhs = build_bigint_vec(&b_components, e_min); - // Forward elimination with first-non-zero pivoting. - for k in 0..D { - // Find first non-zero pivot in column k at or below row k. - if mat[k][k] == zero { - if let Some(swap_row) = ((k + 1)..D).find(|&i| mat[i][k] != zero) { - mat.swap(k, swap_row); - rhs.swap(k, swap_row); - } else { - cold_path(); - return Err(LaError::Singular { pivot_col: k }); - } - } - - // Eliminate below pivot. - let pivot = mat[k][k].clone(); - for i in (k + 1)..D { - if mat[i][k] != zero { - let factor = &mat[i][k] / &pivot; - // We need index `j` to read mat[k][j] and write mat[i][j] - // (two distinct rows) — iterators can't borrow both. - #[allow(clippy::needless_range_loop)] - for j in (k + 1)..D { - let term = &factor * &mat[k][j]; - mat[i][j] -= term; - } - let rhs_term = &factor * &rhs[k]; - rhs[i] -= rhs_term; - mat[i][k] = zero.clone(); - } - } + if let BareissResult::Singular { pivot_col } = bareiss_forward_eliminate(&mut a, Some(&mut rhs)) + { + cold_path(); + return Err(LaError::Singular { pivot_col }); } - // Back-substitution. - let mut x: [BigRational; D] = from_fn(|_| zero.clone()); + // Back-substitution in `BigRational`. Only the upper triangle of `a` + // and the transformed `rhs` are read, each exactly once — so we + // `mem::take` instead of `clone` to avoid a per-entry allocation. + // `BigInt::default()` is the zero value and does not allocate. + let mut x: [BigRational; D] = from_fn(|_| BigRational::from_integer(BigInt::from(0))); for i in (0..D).rev() { - let mut sum = rhs[i].clone(); + let mut sum = BigRational::from_integer(take(&mut rhs[i])); for j in (i + 1)..D { - sum -= &mat[i][j] * &x[j]; + let a_ij = BigRational::from_integer(take(&mut a[i][j])); + sum -= &a_ij * &x[j]; } - x[i] = sum / &mat[i][i]; + let a_ii = BigRational::from_integer(take(&mut a[i][i])); + x[i] = sum / &a_ii; } Ok(x) @@ -389,8 +497,7 @@ impl Matrix { /// Returns [`LaError::NonFinite`] if any matrix entry is NaN or infinite. #[inline] pub fn det_exact(&self) -> Result { - validate_finite(self)?; - Ok(bareiss_det(self)) + bareiss_det(self) } /// Exact determinant converted to `f64`. @@ -425,7 +532,7 @@ impl Matrix { } } - /// Exact linear system solve using arbitrary-precision rational arithmetic. + /// Exact linear system solve using hybrid integer/rational arithmetic. /// /// Solves `A x = b` where `A` is `self` and `b` is the given vector. /// Returns the exact solution as `[BigRational; D]`. Every finite `f64` @@ -438,6 +545,18 @@ impl Matrix { /// exact circumcenter computation for near-degenerate simplices where /// f64 arithmetic may produce wildly wrong results. /// + /// # Algorithm + /// + /// Matrix and RHS entries are decomposed via IEEE 754 bit extraction and + /// scaled to a shared power-of-two base so the augmented system `(A | b)` + /// becomes integer-valued. Forward elimination runs entirely in `BigInt` + /// with fraction-free Bareiss updates — no `BigRational`, no GCD, no + /// denominator tracking in the `O(D³)` phase. Only the upper-triangular + /// result is lifted into `BigRational` for back-substitution (the `O(D²)` + /// phase where fractions are inherent). First-non-zero pivoting is used + /// throughout; since all arithmetic is exact, any non-zero pivot yields + /// the correct answer (no numerical-stability concerns). + /// /// # Examples /// ``` /// use la_stack::prelude::*; @@ -456,8 +575,6 @@ impl Matrix { /// Returns [`LaError::Singular`] if the matrix is exactly singular. #[inline] pub fn solve_exact(&self, b: Vector) -> Result<[BigRational; D], LaError> { - validate_finite(self)?; - validate_finite_vec(&b)?; gauss_solve(self, &b) } @@ -537,14 +654,15 @@ impl Matrix { /// Returns [`LaError::NonFinite`] if any matrix entry is NaN or infinite. #[inline] pub fn det_sign_exact(&self) -> Result { - validate_finite(self)?; - // Stage 1: f64 fast filter for D ≤ 4. // // When entries are large (e.g. near f64::MAX) the determinant can // overflow to infinity even though every individual entry is finite. // In that case the fast filter is inconclusive; fall through to the - // exact Bareiss path. + // exact Bareiss path. For NaN/±∞ entries IEEE 754 propagates + // non-finite through `det_direct()`, the `det_f64.is_finite()` + // guard fails, and we also fall through — validation then happens + // inside `bareiss_det_int` via `decompose_matrix`. match self.det_direct() { Some(det_f64) if let Some(err) = self.det_errbound() @@ -561,11 +679,13 @@ impl Matrix { } // Stage 2: integer Bareiss fallback — the 2^(D×e_min) scale factor - // is always positive, so det_int.sign() == det(A).sign(). - // This is the cold path: the fast filter resolves the vast majority of - // well-conditioned calls without allocating. + // is always positive, so det_int.sign() == det(A).sign(). This is + // the cold path: the fast filter resolves the vast majority of + // well-conditioned calls without allocating. `bareiss_det_int` + // validates finiteness via `decompose_matrix`, so NaN/±∞ inputs + // surface here as `Err(LaError::NonFinite)`. cold_path(); - let (det_int, _) = bareiss_det_int(self); + let (det_int, _) = bareiss_det_int(self)?; Ok(match det_int.sign() { Sign::Plus => 1, Sign::Minus => -1, @@ -581,6 +701,43 @@ mod tests { use pastey::paste; + // ----------------------------------------------------------------------- + // Test helpers + // ----------------------------------------------------------------------- + + /// Build an exact `BigRational` from an `f64` via IEEE 754 bit decomposition. + /// + /// Thin wrapper over [`f64_decompose`] that packs the mantissa/exponent + /// pair into a fully-formed `BigRational` of the form `±m · 2^e`. The + /// production code paths (`bareiss_det_int`, `gauss_solve`) instead + /// decompose every entry into a shared-scale `BigInt` matrix, which + /// avoids per-entry GCD work in the elimination loops — so this helper + /// is not used by them and lives here to keep test assertions concise + /// (e.g. `assert_eq!(x[0], f64_to_bigrational(3.0))`). + /// + /// See `REFERENCES.md` \[9-10\] for the IEEE 754 standard and Goldberg's + /// survey of floating-point representation. + /// + /// # Panics + /// Panics if `x` is NaN or infinite. + fn f64_to_bigrational(x: f64) -> BigRational { + let Some((mantissa, exponent, is_negative)) = f64_decompose(x) else { + return BigRational::from_integer(BigInt::from(0)); + }; + + let numer = if is_negative { + -BigInt::from(mantissa) + } else { + BigInt::from(mantissa) + }; + + if exponent >= 0 { + BigRational::new_raw(numer << exponent.cast_unsigned(), BigInt::from(1u32)) + } else { + BigRational::new_raw(numer, BigInt::from(1u32) << (-exponent).cast_unsigned()) + } + } + // ----------------------------------------------------------------------- // Macro-generated per-dimension tests (D=2..5) // ----------------------------------------------------------------------- @@ -1010,30 +1167,41 @@ mod tests { #[test] fn bareiss_det_int_d0() { - let (det, exp) = bareiss_det_int(&Matrix::<0>::zero()); + let (det, exp) = bareiss_det_int(&Matrix::<0>::zero()).unwrap(); assert_eq!(det, BigInt::from(1)); assert_eq!(exp, 0); } + /// Table-driven coverage of the D=1 fast-path: each 1×1 matrix + /// decomposes to `(±mant, exp)` directly. Includes an integer, zero, + /// a negative fractional, and a positive fractional case — the + /// combinations that exercise the sign handling, the all-zero early + /// return, trailing-zero stripping, and negative exponent scaling. #[test] - fn bareiss_det_int_d1_value() { - // 7.0 = 7 × 2^0 - let (det, exp) = bareiss_det_int(&Matrix::<1>::from_rows([[7.0]])); - assert_eq!(det, BigInt::from(7)); - assert_eq!(exp, 0); - } - - #[test] - fn bareiss_det_int_d1_zero() { - let (det, _) = bareiss_det_int(&Matrix::<1>::from_rows([[0.0]])); - assert_eq!(det, BigInt::from(0)); + fn bareiss_det_int_d1_cases() { + let cases: &[(f64, i64, i32)] = &[ + // (input, expected_det_int, expected_exp) + (7.0, 7, 0), // integer → (7, 0) + (0.0, 0, 0), // all-zero early return → (0, 0) + (-3.5, -7, -1), // -3.5 = -7 × 2^(-1) + (0.5, 1, -1), // 0.5 = 1 × 2^(-1) + ]; + for &(input, expected_det_int, expected_exp) in cases { + let (det, exp) = bareiss_det_int(&Matrix::<1>::from_rows([[input]])).unwrap(); + assert_eq!( + det, + BigInt::from(expected_det_int), + "det_int for input={input}" + ); + assert_eq!(exp, expected_exp, "exp for input={input}"); + } } #[test] fn bareiss_det_int_d2_known() { // det([[1,2],[3,4]]) = -2 let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); - let (det_int, total_exp) = bareiss_det_int(&m); + let (det_int, total_exp) = bareiss_det_int(&m).unwrap(); // Reconstruct and verify. let det = bigint_exp_to_bigrational(det_int, total_exp); assert_eq!(det, BigRational::from_integer(BigInt::from(-2))); @@ -1041,7 +1209,7 @@ mod tests { #[test] fn bareiss_det_int_all_zeros() { - let (det, _) = bareiss_det_int(&Matrix::<3>::zero()); + let (det, _) = bareiss_det_int(&Matrix::<3>::zero()).unwrap(); assert_eq!(det, BigInt::from(0)); } @@ -1049,7 +1217,7 @@ mod tests { fn bareiss_det_int_sign_matches_det_sign_exact() { // The sign of det_int should match det_sign_exact for various matrices. 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_int, _) = bareiss_det_int(&m); + let (det_int, _) = bareiss_det_int(&m).unwrap(); assert_eq!(det_int.sign(), Sign::Minus); // det = -1 } @@ -1058,34 +1226,46 @@ mod tests { // Entries with negative exponents: 0.5 = 1×2^(-1), 0.25 = 1×2^(-2). // det([[0.5, 0.25], [1.0, 1.0]]) = 0.5×1.0 − 0.25×1.0 = 0.25 let m = Matrix::<2>::from_rows([[0.5, 0.25], [1.0, 1.0]]); - let (det_int, total_exp) = bareiss_det_int(&m); + let (det_int, total_exp) = bareiss_det_int(&m).unwrap(); let det = bigint_exp_to_bigrational(det_int, total_exp); assert_eq!(det, BigRational::new(BigInt::from(1), BigInt::from(4))); } #[test] - fn bareiss_det_int_d1_negative() { - // -3.5 = -7 × 2^(-1) - let (det, exp) = bareiss_det_int(&Matrix::<1>::from_rows([[-3.5]])); - assert_eq!(det, BigInt::from(-7)); - assert_eq!(exp, -1); + fn bareiss_det_int_d3_with_pivoting() { + // Zero on diagonal → exercises pivot swap inside bareiss_det_int. + 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_int, total_exp) = bareiss_det_int(&m).unwrap(); + let det = bigint_exp_to_bigrational(det_int, total_exp); + assert_eq!(det, BigRational::from_integer(BigInt::from(-1))); } + /// Non-finite matrix entries surface as `LaError::NonFinite` with the + /// row/col of the first offending entry. #[test] - fn bareiss_det_int_d1_fractional() { - // 0.5 = 1 × 2^(-1) - let (det, exp) = bareiss_det_int(&Matrix::<1>::from_rows([[0.5]])); - assert_eq!(det, BigInt::from(1)); - assert_eq!(exp, -1); + fn bareiss_det_int_errs_on_nan() { + let mut m = Matrix::<3>::identity(); + m.set(1, 2, f64::NAN); + assert_eq!( + bareiss_det_int(&m), + Err(LaError::NonFinite { + row: Some(1), + col: 2 + }) + ); } #[test] - fn bareiss_det_int_d3_with_pivoting() { - // Zero on diagonal → exercises pivot swap inside bareiss_det_int. - 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_int, total_exp) = bareiss_det_int(&m); - let det = bigint_exp_to_bigrational(det_int, total_exp); - assert_eq!(det, BigRational::from_integer(BigInt::from(-1))); + fn bareiss_det_int_errs_on_inf() { + let mut m = Matrix::<2>::identity(); + m.set(0, 0, f64::INFINITY); + assert_eq!( + bareiss_det_int(&m), + Err(LaError::NonFinite { + row: Some(0), + col: 0 + }) + ); } /// Per AGENTS.md: dimension-generic tests must cover D=2–5. @@ -1094,7 +1274,8 @@ mod tests { paste! { #[test] fn []() { - let (det_int, total_exp) = bareiss_det_int(&Matrix::<$d>::identity()); + let (det_int, total_exp) = + bareiss_det_int(&Matrix::<$d>::identity()).unwrap(); let det = bigint_exp_to_bigrational(det_int, total_exp); assert_eq!(det, BigRational::from_integer(BigInt::from(1))); } @@ -1161,13 +1342,13 @@ mod tests { #[test] fn bareiss_det_d0_is_one() { - let det = bareiss_det(&Matrix::<0>::zero()); + let det = bareiss_det(&Matrix::<0>::zero()).unwrap(); 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]])); + let det = bareiss_det(&Matrix::<1>::from_rows([[7.0]])).unwrap(); assert_eq!(det, f64_to_bigrational(7.0)); } @@ -1175,7 +1356,7 @@ mod tests { 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); + let det = bareiss_det(&m).unwrap(); // det of this permutation matrix = -1 assert_eq!(det, BigRational::from_integer(BigInt::from(-1))); } @@ -1184,7 +1365,7 @@ mod tests { 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); + let det = bareiss_det(&m).unwrap(); assert_eq!(det, BigRational::from_integer(BigInt::from(0))); } @@ -1419,6 +1600,63 @@ mod tests { gen_solve_exact_f64_agrees_with_lu!(4); gen_solve_exact_f64_agrees_with_lu!(5); + /// Round-trip: for a well-conditioned integer matrix `A` and integer + /// target `x0`, solving `A x = A x0` must return `x0` exactly. All + /// intermediate values stay small enough that `A * x0` is exactly + /// representable in `f64`, so the round-trip is a precise equality + /// check on the hybrid BigInt/BigRational path. + macro_rules! gen_solve_exact_roundtrip_tests { + ($d:literal) => { + paste! { + #[test] + #[allow(clippy::cast_precision_loss)] + fn []() { + // A = D * I + J (diag = D+1, off-diag = 1). Invertible + // for any D >= 1 and cheap to multiply by hand. + let mut rows = [[0.0f64; $d]; $d]; + for r in 0..$d { + for c in 0..$d { + rows[r][c] = if r == c { + f64::from($d) + 1.0 + } else { + 1.0 + }; + } + } + let a = Matrix::<$d>::from_rows(rows); + + // x0 = [1, 2, ..., D]. + let mut x0 = [0.0f64; $d]; + for i in 0..$d { + x0[i] = (i + 1) as f64; + } + + // b = A * x0 computed in f64. With small integers the + // multiply-add sequence is exact. + let mut b_arr = [0.0f64; $d]; + for r in 0..$d { + let mut sum = 0.0_f64; + for c in 0..$d { + sum += rows[r][c] * x0[c]; + } + b_arr[r] = sum; + } + let b = Vector::<$d>::new(b_arr); + + let x = a.solve_exact(b).unwrap(); + for i in 0..$d { + assert_eq!(x[i], f64_to_bigrational(x0[i])); + } + } + } + }; + } + + gen_solve_exact_roundtrip_tests!(2); + gen_solve_exact_roundtrip_tests!(3); + gen_solve_exact_roundtrip_tests!(4); + gen_solve_exact_roundtrip_tests!(5); + // ----------------------------------------------------------------------- // solve_exact: dimension-specific tests // ----------------------------------------------------------------------- @@ -1489,6 +1727,285 @@ mod tests { assert_eq!(x[4], f64_to_bigrational(50.0)); } + /// Entries near `f64::MAX / 2` are finite but their product would + /// overflow to ±∞ in pure f64 arithmetic. The `BigInt` augmented-system + /// path computes the correct solution without any overflow. The D×D + /// case uses a diagonal matrix with `big` on every diagonal and a RHS + /// of `[big, …, big, 0]`, giving the known solution `[1, …, 1, 0]`. + macro_rules! gen_solve_exact_large_finite_entries_tests { + ($d:literal) => { + paste! { + #[test] + fn []() { + let big = f64::MAX / 2.0; + assert!(big.is_finite()); + // D×D diagonal matrix with `big` on the diagonal. + let mut rows = [[0.0f64; $d]; $d]; + for i in 0..$d { + rows[i][i] = big; + } + let a = Matrix::<$d>::from_rows(rows); + // RHS = [big, …, big, 0] → x = [1, …, 1, 0]. + let mut b_arr = [big; $d]; + b_arr[$d - 1] = 0.0; + let b = Vector::<$d>::new(b_arr); + let x = a.solve_exact(b).unwrap(); + for i in 0..($d - 1) { + assert_eq!(x[i], BigRational::from_integer(BigInt::from(1))); + } + assert_eq!(x[$d - 1], BigRational::from_integer(BigInt::from(0))); + } + } + }; + } + + gen_solve_exact_large_finite_entries_tests!(2); + gen_solve_exact_large_finite_entries_tests!(3); + gen_solve_exact_large_finite_entries_tests!(4); + gen_solve_exact_large_finite_entries_tests!(5); + + /// Matrix and RHS entries span many orders of magnitude (from + /// `f64::MIN_POSITIVE` up through `1e100`). This exercises the + /// shared `e_min` scaling: even the largest shift keeps every entry a + /// representable `BigInt`. The D×D case alternates `huge`/`tiny` + /// along the diagonal with a matching RHS, giving `x = [1, …, 1]`. + macro_rules! gen_solve_exact_mixed_magnitude_entries_tests { + ($d:literal) => { + paste! { + #[test] + fn []() { + let tiny = f64::MIN_POSITIVE; // 2^-1022, smallest normal + let huge = 1.0e100_f64; + // Alternate huge/tiny along the diagonal. + let mut rows = [[0.0f64; $d]; $d]; + let mut b_arr = [0.0f64; $d]; + for i in 0..$d { + let val = if i % 2 == 0 { huge } else { tiny }; + rows[i][i] = val; + b_arr[i] = val; + } + let a = Matrix::<$d>::from_rows(rows); + let b = Vector::<$d>::new(b_arr); + let x = a.solve_exact(b).unwrap(); + for i in 0..$d { + assert_eq!(x[i], BigRational::from_integer(BigInt::from(1))); + } + } + } + }; + } + + gen_solve_exact_mixed_magnitude_entries_tests!(2); + gen_solve_exact_mixed_magnitude_entries_tests!(3); + gen_solve_exact_mixed_magnitude_entries_tests!(4); + gen_solve_exact_mixed_magnitude_entries_tests!(5); + + /// Subnormal RHS entries must survive the decomposition and + /// back-substitution paths unchanged. The D×D case uses the identity + /// matrix and RHS `[1·tiny, 2·tiny, …, D·tiny]`; each entry remains a + /// valid subnormal f64 (integer multiples of `2^-1074` fit in the + /// 52-bit subnormal mantissa for the small integers used here). + macro_rules! gen_solve_exact_subnormal_rhs_tests { + ($d:literal) => { + paste! { + #[test] + #[allow(clippy::cast_precision_loss)] + fn []() { + let tiny = 5e-324_f64; // smallest positive subnormal + assert!(tiny.is_subnormal()); + let a = Matrix::<$d>::identity(); + // b[i] = (i+1) · tiny — each entry remains a valid subnormal. + let mut b_arr = [0.0f64; $d]; + for i in 0..$d { + b_arr[i] = (i + 1) as f64 * tiny; + assert!(b_arr[i].is_subnormal()); + } + let b = Vector::<$d>::new(b_arr); + let x = a.solve_exact(b).unwrap(); + for i in 0..$d { + assert_eq!(x[i], f64_to_bigrational((i + 1) as f64 * tiny)); + } + } + } + }; + } + + gen_solve_exact_subnormal_rhs_tests!(2); + gen_solve_exact_subnormal_rhs_tests!(3); + gen_solve_exact_subnormal_rhs_tests!(4); + gen_solve_exact_subnormal_rhs_tests!(5); + + /// Pivoting path with a zero top-left entry forces a row swap in the + /// `BigInt` forward-elimination loop and propagates it to the RHS. + /// Combined with a fractional solution, this exercises the + /// `BigRational` back-substitution after integer forward elimination. + /// + /// The 2×2 block `[[0, 1], [2, 1]]` with rhs `[3, 4]` (→ `x = [1/2, 3]`) + /// is embedded into the top-left of a D×D identity matrix. Remaining + /// rows contribute pass-through equalities `x[i] = b[i]`, so the same + /// fractional solution appears at indices 0 and 1 regardless of D. + macro_rules! gen_solve_exact_pivot_swap_fractional_tests { + ($d:literal) => { + paste! { + #[test] + #[allow(clippy::cast_precision_loss)] + // `2..$d` is empty when D=2 (no padded rows); that is the + // intended behaviour of the macro, not a bug. + #[allow(clippy::reversed_empty_ranges)] + fn []() { + // Top-left 2×2: A = [[0, 1], [2, 1]]. After swap: + // [[2, 1], [0, 1]], rhs = [4, 3] → x[1] = 3, x[0] = 1/2. + let mut rows = [[0.0f64; $d]; $d]; + rows[0][1] = 1.0; + rows[1][0] = 2.0; + rows[1][1] = 1.0; + // Identity padding for the remaining rows. + for i in 2..$d { + rows[i][i] = 1.0; + } + let a = Matrix::<$d>::from_rows(rows); + // b = [3, 4, 12, 13, …]; padded entries are arbitrary + // finite integers so the identity block gives x[i] = b[i]. + let mut b_arr = [0.0f64; $d]; + b_arr[0] = 3.0; + b_arr[1] = 4.0; + for i in 2..$d { + b_arr[i] = (i + 10) as f64; + } + let b = Vector::<$d>::new(b_arr); + let x = a.solve_exact(b).unwrap(); + assert_eq!(x[0], BigRational::new(BigInt::from(1), BigInt::from(2))); + assert_eq!(x[1], BigRational::from_integer(BigInt::from(3))); + for i in 2..$d { + assert_eq!(x[i], f64_to_bigrational((i + 10) as f64)); + } + } + } + }; + } + + gen_solve_exact_pivot_swap_fractional_tests!(2); + gen_solve_exact_pivot_swap_fractional_tests!(3); + gen_solve_exact_pivot_swap_fractional_tests!(4); + gen_solve_exact_pivot_swap_fractional_tests!(5); + + /// Mid-elimination pivot swap: the 3×3 block + /// `[[1, 2, 3], [0, 0, 4], [0, 5, 6]]` has a non-zero pivot at k=0 but + /// a zero pivot at k=1, so the swap happens *during* forward + /// elimination rather than at the start. With rhs `[6, 7, 8]` the + /// exact solution is `[7/4, -1/2, 7/4]`. For D > 3 the block is + /// embedded into the top-left of a D×D identity matrix so the same + /// fractional solution appears in `x[0..3]` and `x[i] = b[i]` for + /// `i >= 3`. + macro_rules! gen_solve_exact_mid_pivot_swap_tests { + ($d:literal) => { + paste! { + #[test] + #[allow(clippy::cast_precision_loss)] + // `3..$d` is empty when D=3 (no padded rows); that is the + // intended behaviour of the macro, not a bug. + #[allow(clippy::reversed_empty_ranges)] + fn []() { + let mut rows = [[0.0f64; $d]; $d]; + rows[0][0] = 1.0; rows[0][1] = 2.0; rows[0][2] = 3.0; + // rows[1][0..2] are zero; rows[1][2] = 4. + rows[1][2] = 4.0; + rows[2][1] = 5.0; rows[2][2] = 6.0; + // Identity padding for the remaining rows. + for i in 3..$d { + rows[i][i] = 1.0; + } + let a = Matrix::<$d>::from_rows(rows); + let mut b_arr = [0.0f64; $d]; + b_arr[0] = 6.0; + b_arr[1] = 7.0; + b_arr[2] = 8.0; + for i in 3..$d { + b_arr[i] = (i + 10) as f64; + } + let b = Vector::<$d>::new(b_arr); + let x = a.solve_exact(b).unwrap(); + // x[0..3] = [7/4, -1/2, 7/4]. + assert_eq!(x[0], BigRational::new(BigInt::from(7), BigInt::from(4))); + assert_eq!(x[1], BigRational::new(BigInt::from(-1), BigInt::from(2))); + assert_eq!(x[2], BigRational::new(BigInt::from(7), BigInt::from(4))); + for i in 3..$d { + assert_eq!(x[i], f64_to_bigrational((i + 10) as f64)); + } + } + } + }; + } + + gen_solve_exact_mid_pivot_swap_tests!(3); + gen_solve_exact_mid_pivot_swap_tests!(4); + gen_solve_exact_mid_pivot_swap_tests!(5); + + /// Rank-deficient singular: the last column is identically zero and the + /// leading `(D-1)×(D-1)` block is full rank, so every intermediate + /// pivot is non-zero and the singularity surfaces only at the final + /// column. The matrix is identity in the top-left `(D-1)×(D-1)` with + /// a row of ones as the last row (and an all-zero last column), so the + /// rank is exactly `D-1`. `solve_exact` must return + /// `LaError::Singular { pivot_col: D - 1 }`. + macro_rules! gen_solve_exact_singular_rank_deficient_tests { + ($d:literal) => { + paste! { + #[test] + fn []() { + let mut rows = [[0.0f64; $d]; $d]; + for i in 0..($d - 1) { + rows[i][i] = 1.0; + rows[$d - 1][i] = 1.0; + } + // Last column is left all-zero → rank exactly D-1. + let a = Matrix::<$d>::from_rows(rows); + let b = Vector::<$d>::new([1.0; $d]); + assert_eq!( + a.solve_exact(b), + Err(LaError::Singular { pivot_col: $d - 1 }) + ); + } + } + }; + } + + gen_solve_exact_singular_rank_deficient_tests!(2); + gen_solve_exact_singular_rank_deficient_tests!(3); + gen_solve_exact_singular_rank_deficient_tests!(4); + gen_solve_exact_singular_rank_deficient_tests!(5); + + /// Zero RHS with a non-singular matrix. Every Bareiss update reads + /// `rhs[k]` and `rhs[i]`, both initialised to zero; every update + /// produces zero; back-substitution therefore yields `x = 0` + /// regardless of the matrix entries. This exercises the + /// back-substitution `mem::take` path against an all-zero `rhs`. + macro_rules! gen_solve_exact_zero_rhs_tests { + ($d:literal) => { + paste! { + #[test] + fn []() { + // A = D*I + J (diagonally dominant, invertible). + let mut rows = [[1.0f64; $d]; $d]; + for i in 0..$d { + rows[i][i] = f64::from($d) + 1.0; + } + let a = Matrix::<$d>::from_rows(rows); + let b = Vector::<$d>::zero(); + let x = a.solve_exact(b).unwrap(); + for xi in &x { + assert_eq!(*xi, BigRational::from_integer(BigInt::from(0))); + } + } + } + }; + } + + gen_solve_exact_zero_rhs_tests!(2); + gen_solve_exact_zero_rhs_tests!(3); + gen_solve_exact_zero_rhs_tests!(4); + gen_solve_exact_zero_rhs_tests!(5); + // ----------------------------------------------------------------------- // solve_exact_f64: dimension-specific tests // ----------------------------------------------------------------------- @@ -1666,64 +2183,4 @@ mod tests { fn f64_to_bigrational_panics_on_neg_inf() { f64_to_bigrational(f64::NEG_INFINITY); } - - // ----------------------------------------------------------------------- - // validate_finite_vec tests - // ----------------------------------------------------------------------- - - #[test] - fn validate_finite_vec_ok() { - assert!(validate_finite_vec(&Vector::<3>::new([1.0, 2.0, 3.0])).is_ok()); - } - - #[test] - fn validate_finite_vec_err_on_nan() { - assert_eq!( - validate_finite_vec(&Vector::<2>::new([f64::NAN, 1.0])), - Err(LaError::NonFinite { row: None, col: 0 }) - ); - } - - #[test] - fn validate_finite_vec_err_on_inf() { - assert_eq!( - validate_finite_vec(&Vector::<2>::new([1.0, f64::NEG_INFINITY])), - Err(LaError::NonFinite { row: None, col: 1 }) - ); - } - - // ----------------------------------------------------------------------- - // validate_finite tests - // ----------------------------------------------------------------------- - - #[test] - fn validate_finite_ok_for_finite() { - assert!(validate_finite(&Matrix::<3>::identity()).is_ok()); - } - - #[test] - fn validate_finite_err_on_nan() { - let mut m = Matrix::<2>::identity(); - m.set(1, 0, f64::NAN); - assert_eq!( - validate_finite(&m), - Err(LaError::NonFinite { - row: Some(1), - col: 0 - }) - ); - } - - #[test] - fn validate_finite_err_on_inf() { - let mut m = Matrix::<2>::identity(); - m.set(0, 1, f64::NEG_INFINITY); - assert_eq!( - validate_finite(&m), - Err(LaError::NonFinite { - row: Some(0), - col: 1 - }) - ); - } }