Struct Matrix

Source

pub struct Matrix<const D: usize> { /* private fields */ }

Expand description

Fixed-size square matrix D×D, stored inline.

Implementations§

Source §

impl<const D: usize> Matrix<D>

Source

pub fn det_exact(&self) -> Result<BigRational, LaError>

Exact determinant using arbitrary-precision rational arithmetic.

Returns the determinant as an exact BigRational value. Every finite f64 is exactly representable as a rational, so the conversion is lossless and the result is provably correct.

§When to use

Use this when you need the exact determinant value — for example, exact volume computation or distinguishing truly-degenerate simplices from near-degenerate ones. If you only need the sign, prefer det_sign_exact which has a fast f64 filter.

§Examples

use la_stack::prelude::*;

let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
let det = m.det_exact().unwrap();
// det = 1*4 - 2*3 = -2  (exact)
assert_eq!(det, BigRational::from_integer((-2).into()));

§Errors

Returns LaError::NonFinite if any matrix entry is NaN or infinite.

Source

pub fn det_exact_f64(&self) -> Result<f64, LaError>

Exact determinant converted to f64.

Computes the exact BigRational determinant via det_exact and converts it to the nearest f64. This is useful when you want the most accurate f64 determinant possible without committing to BigRational in your downstream code.

§Examples

use la_stack::prelude::*;

let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
let det = m.det_exact_f64().unwrap();
assert!((det - (-2.0)).abs() <= f64::EPSILON);

§Errors

Returns LaError::NonFinite if any matrix entry is NaN or infinite. Returns LaError::Overflow if the exact determinant is too large to represent as a finite f64.

Source

pub fn solve_exact(&self, b: Vector<D>) -> Result<[BigRational; D], LaError>

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 is exactly representable as a rational, so the conversion is lossless and the result is provably correct.

§When to use

Use this when you need a provably correct solution — for example, 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::*;

// A x = b  where A = [[1,2],[3,4]], b = [5, 11]  →  x = [1, 2]
let a = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
let b = Vector::<2>::new([5.0, 11.0]);
let x = a.solve_exact(b).unwrap();
assert_eq!(x[0], BigRational::from_integer(1.into()));
assert_eq!(x[1], BigRational::from_integer(2.into()));

§Errors

Returns LaError::NonFinite if any matrix or vector entry is NaN or infinite. Returns LaError::Singular if the matrix is exactly singular.

Source

pub fn solve_exact_f64(&self, b: Vector<D>) -> Result<Vector<D>, LaError>

Exact linear system solve converted to f64.

Computes the exact BigRational solution via solve_exact and converts each component to the nearest f64. This is useful when you want the most accurate f64 solution possible without committing to BigRational downstream.

§Examples

use la_stack::prelude::*;

let a = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
let b = Vector::<2>::new([5.0, 11.0]);
let x = a.solve_exact_f64(b).unwrap().into_array();
assert!((x[0] - 1.0).abs() <= f64::EPSILON);
assert!((x[1] - 2.0).abs() <= f64::EPSILON);

§Errors

Returns LaError::NonFinite if any matrix or vector entry is NaN or infinite. Returns LaError::Singular if the matrix is exactly singular. Returns LaError::Overflow if any component of the exact solution is too large to represent as a finite f64.

Source

pub fn det_sign_exact(&self) -> Result<i8, LaError>

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), integer-only Bareiss elimination (bareiss_det_int) computes the exact sign without constructing any BigRational values.

§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().unwrap(), 0);

assert_eq!(Matrix::<3>::identity().det_sign_exact().unwrap(), 1);

§Errors

Returns LaError::NonFinite if any matrix entry is NaN or infinite.

Source §

impl<const D: usize> Matrix<D>

Source

pub const fn from_rows(rows: [[f64; D]; D]) -> Self

Construct from row-major storage.

§Examples

use la_stack::prelude::*;

let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
assert_eq!(m.get(0, 1), Some(2.0));

Source

pub const fn zero() -> Self

All-zeros matrix.

§Examples

use la_stack::prelude::*;

let z = Matrix::<2>::zero();
assert_eq!(z.get(1, 1), Some(0.0));

Source

pub const fn identity() -> Self

Identity matrix.

§Examples

use la_stack::prelude::*;

let i = Matrix::<3>::identity();
assert_eq!(i.get(0, 0), Some(1.0));
assert_eq!(i.get(0, 1), Some(0.0));
assert_eq!(i.get(2, 2), Some(1.0));

Source

pub const fn get(&self, r: usize, c: usize) -> Option<f64>

Get an element with bounds checking.

§Examples

use la_stack::prelude::*;

let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
assert_eq!(m.get(1, 0), Some(3.0));
assert_eq!(m.get(2, 0), None);

Source

pub const fn set(&mut self, r: usize, c: usize, value: f64) -> bool

Set an element with bounds checking.

Returns true if the index was in-bounds.

§Examples

use la_stack::prelude::*;

let mut m = Matrix::<2>::zero();
assert!(m.set(0, 1, 2.5));
assert_eq!(m.get(0, 1), Some(2.5));
assert!(!m.set(10, 0, 1.0));

Source

pub const fn inf_norm(&self) -> f64

Infinity norm (maximum absolute row sum).

§Non-finite handling

If any entry is NaN, the result is NaN. NaN is detected explicitly because a naive row_sum > max_row_sum comparison silently skips NaN rows (every ordered comparison against NaN is false). If any entry is infinite (and no entry is NaN), the result is +∞.

§Examples

use la_stack::prelude::*;

let m = Matrix::<2>::from_rows([[1.0, -2.0], [3.0, 4.0]]);
assert!((m.inf_norm() - 7.0).abs() <= 1e-12);

// NaN entries propagate to the norm.
let nan = Matrix::<2>::from_rows([[f64::NAN, 1.0], [2.0, 3.0]]);
assert!(nan.inf_norm().is_nan());

Source

pub fn is_symmetric(&self, rel_tol: f64) -> bool

Returns true if the matrix is symmetric within a relative tolerance.

Two entries self[r][c] and self[c][r] are considered equal (for the purposes of symmetry) when |self[r][c] - self[c][r]| <= rel_tol * max(1.0, self.inf_norm()). This mirrors the predicate used internally by the debug-build symmetry check inside ldlt, so callers can pre-validate matrices that may come from untrusted sources without relying on a debug-only panic.

Use first_asymmetry to locate the first offending pair when this returns false.

§NaN / infinity handling

Any non-finite |self[r][c] - self[c][r]| (NaN or ±∞) causes this predicate to return false. This catches both NaN off-diagonals and asymmetric pairs where one side is infinite and the other is finite (which would otherwise slip through when inf_norm() blows eps up to +∞ and makes diff > eps trivially false). A matrix whose inf_norm is +∞ can still tolerate finite asymmetries under an infinite eps — callers who need strict equality on large-magnitude finite entries should validate finiteness separately.

§Panics

In debug builds, panics if rel_tol is negative or NaN; in release builds these are silently treated as garbage-in garbage-out, matching the convention of lu and ldlt.

§Examples

use la_stack::prelude::*;

let a = Matrix::<2>::from_rows([[4.0, 2.0], [2.0, 3.0]]);
assert!(a.is_symmetric(1e-12));

let b = Matrix::<2>::from_rows([[4.0, 2.0], [3.0, 3.0]]);
assert!(!b.is_symmetric(1e-12));

Source

pub fn first_asymmetry(&self, rel_tol: f64) -> Option<(usize, usize)>

Returns the indices (r, c) (with r < c) of the first off-diagonal pair that violates symmetry, or None if the matrix is symmetric within rel_tol.

Iteration order is row-major over the strict upper triangle, so the returned indices are the lexicographically smallest such pair. The predicate is the same as is_symmetric: |self[r][c] - self[c][r]| <= rel_tol * max(1.0, self.inf_norm()).

§Panics

In debug builds, panics if rel_tol is negative or NaN.

§Examples

use la_stack::prelude::*;

let a = Matrix::<3>::from_rows([
    [1.0, 2.0, 0.0],
    [2.0, 4.0, 5.0],
    [0.0, 6.0, 9.0], // 6.0 breaks symmetry with a[1][2] = 5.0
]);
assert_eq!(a.first_asymmetry(1e-12), Some((1, 2)));
assert_eq!(Matrix::<3>::identity().first_asymmetry(1e-12), None);

Source

pub fn lu(self, tol: f64) -> Result<Lu<D>, LaError>

Compute an LU decomposition with partial pivoting.

§Examples

use la_stack::prelude::*;

let a = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
let lu = a.lu(DEFAULT_PIVOT_TOL)?;

let b = Vector::<2>::new([5.0, 11.0]);
let x = lu.solve_vec(b)?.into_array();

assert!((x[0] - 1.0).abs() <= 1e-12);
assert!((x[1] - 2.0).abs() <= 1e-12);

§Errors

Returns LaError::Singular if, for some column k, the largest-magnitude candidate pivot in that column satisfies |pivot| <= tol (so no numerically usable pivot exists). Returns LaError::NonFinite if NaN/∞ is detected during factorization.

Source

pub fn ldlt(self, tol: f64) -> Result<Ldlt<D>, LaError>

Compute an LDLT factorization (A = L D Lᵀ) without pivoting.

This is intended for symmetric positive definite (SPD) and positive semi-definite (PSD) matrices such as Gram matrices.

§Preconditions

The input matrix self must be symmetric — that is, self[i][j] == self[j][i] (within rounding) for all i, j. This is a correctness precondition, not merely a performance hint.

In debug builds a debug_assert! verifies symmetry via is_symmetric (relative tolerance scaled by the matrix’s infinity norm) and panics if it fails.
In release builds the check is compiled out for performance. An asymmetric input will be accepted silently and produce a mathematically meaningless factorization — subsequent calls to Ldlt::det and Ldlt::solve_vec will return wrong results with no error.

Callers who cannot statically guarantee symmetry should pre-validate with is_symmetric (or locate the offending pair with first_asymmetry) before calling ldlt. If you need a general-purpose factorization that tolerates non-symmetric inputs, use lu instead.

§Examples

use la_stack::prelude::*;

// Note the symmetric layout: a[0][1] == a[1][0] == 2.0.
let a = Matrix::<2>::from_rows([[4.0, 2.0], [2.0, 3.0]]);
let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL)?;

// det(A) = 8
assert!((ldlt.det() - 8.0).abs() <= 1e-12);

// Solve A x = b
let b = Vector::<2>::new([1.0, 2.0]);
let x = ldlt.solve_vec(b)?.into_array();
assert!((x[0] - (-0.125)).abs() <= 1e-12);
assert!((x[1] - 0.75).abs() <= 1e-12);

§Errors

Returns LaError::Singular if, for some step k, the required diagonal entry d = D[k,k] is <= tol (non-positive or too small). This treats PSD degeneracy (and indefinite inputs) as singular/degenerate. Returns LaError::NonFinite if NaN/∞ is detected during factorization.

Note that an asymmetric input is not reported as an error in release builds — see the Preconditions section above.

Source

pub const fn det_direct(&self) -> Option<f64>

Closed-form determinant for dimensions 0–4, bypassing LU factorization.

Returns Some(det) for D ∈ {0, 1, 2, 3, 4}, None for D ≥ 5. D = 0 returns Some(1.0) (empty product). This is a const fn (Rust 1.94+) and uses fused multiply-add (mul_add) for improved accuracy and performance.

For a determinant that works for any dimension (falling back to LU for D ≥ 5), use det.

§Examples

use la_stack::prelude::*;

let m = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
assert!((m.det_direct().unwrap() - (-2.0)).abs() <= 1e-12);

// D = 0 is the empty product.
assert_eq!(Matrix::<0>::zero().det_direct(), Some(1.0));

// D ≥ 5 returns None.
assert!(Matrix::<5>::identity().det_direct().is_none());

Source

pub fn det(self, tol: f64) -> Result<f64, LaError>

Determinant, using closed-form formulas for D ≤ 4 and LU decomposition for D ≥ 5.

For D ∈ {1, 2, 3, 4}, this bypasses LU factorization entirely for a significant speedup (see det_direct). The tol parameter is only used by the LU fallback path for D ≥ 5.

§Examples

use la_stack::prelude::*;

let det = Matrix::<3>::identity().det(DEFAULT_PIVOT_TOL)?;
assert!((det - 1.0).abs() <= 1e-12);

§Errors

Returns LaError::NonFinite if the result contains NaN or infinity. For D ≥ 5, propagates LU factorization errors (e.g. LaError::Singular).

Source

pub const fn det_errbound(&self) -> Option<f64>

Conservative absolute error bound for det_direct().

Returns Some(bound) such that |det_direct() - det_exact| ≤ bound, or None for D ≥ 5 where no fast bound is available.

For D ≤ 4, the bound is derived from the absolute Leibniz sum using Shewchuk-style error analysis (see REFERENCES.md [8] and the per-constant docs on ERR_COEFF_2, ERR_COEFF_3, and ERR_COEFF_4). For D = 0 or 1, returns Some(0.0) since the determinant computation is exact (no arithmetic).

This method does NOT require the exact feature — the bounds use pure f64 arithmetic and are useful for custom adaptive-precision logic.

§When to use

Use this to build adaptive-precision logic: if |det_direct()| > bound, the f64 sign is provably correct. Otherwise fall back to exact arithmetic.

§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],
]);
let bound = m.det_errbound().unwrap();
let det_approx = m.det_direct().unwrap();
// If |det_approx| > bound, the sign is guaranteed correct.

§Adaptive precision pattern (requires `exact` feature)

use la_stack::prelude::*;

let m = Matrix::<3>::identity();
if let Some(bound) = m.det_errbound() {
    let det = m.det_direct().unwrap();
    if det.abs() > bound {
        // f64 sign is guaranteed correct
        let sign = det.signum() as i8;
    } else {
        // Fall back to exact arithmetic (requires `exact` feature)
        let sign = m.det_sign_exact().unwrap();
    }
} else {
    // D ≥ 5: no fast filter, use exact directly
    let sign = m.det_sign_exact().unwrap();
}