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.

Requires the exact Cargo feature.

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>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
let det = m.det_exact()?;
// det = 1*4 - 2*3 = -2  (exact)
assert_eq!(det, BigRational::from_integer((-2).into()));

§Errors

Returns LaError::NonFinite if stored matrix entries are NaN or infinity.

Returns LaError::Overflow if determinant scaling overflows the internal exponent representation.

Returns LaError::UnsupportedDimension if D cannot be represented in the internal determinant exponent calculation.

Source

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

Exact determinant converted to f64.

Requires the exact Cargo feature.

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>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
let det = m.det_exact_f64()?;
assert!((det - (-2.0)).abs() <= f64::EPSILON);

§Errors

Returns LaError::NonFinite if stored matrix entries are NaN or infinity.

Returns LaError::Overflow if determinant scaling overflows the internal exponent representation or 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.

Requires the exact Cargo feature.

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>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
let b = Vector::<2>::try_new([5.0, 11.0])?;
let x = a.solve_exact(b)?;
assert_eq!(x[0], BigRational::from_integer(1.into()));
assert_eq!(x[1], BigRational::from_integer(2.into()));

§Errors

Returns LaError::NonFinite if stored matrix or right-hand-side entries are NaN or infinity.

Returns LaError::Singular if the matrix is exactly singular.

Source

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

Exact linear system solve converted to f64.

Requires the exact Cargo feature.

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>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
let b = Vector::<2>::try_new([5.0, 11.0])?;
let x = a.solve_exact_f64(b)?.into_array();
assert!((x[0] - 1.0).abs() <= f64::EPSILON);
assert!((x[1] - 2.0).abs() <= f64::EPSILON);

§Errors

Returns LaError::NonFinite if stored matrix or right-hand-side entries are NaN or infinity.

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

Source

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

Exact determinant sign using adaptive-precision arithmetic.

Requires the exact Cargo feature.

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>::try_from_rows([
    [1.0, 2.0, 3.0],
    [4.0, 5.0, 6.0],
    [7.0, 8.0, 9.0],
])?;
// This matrix is singular (row 3 = row 1 + row 2 in exact arithmetic).
assert_eq!(m.det_sign_exact()?, 0);

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

§Errors

Returns LaError::NonFinite if stored matrix entries are NaN or infinity. This exact sign path has no additional runtime errors for finite matrices.

Source §

impl<const D: usize> Matrix<D>

Source

pub const fn try_from_rows(rows: [[f64; D]; D]) -> Result<Self, LaError>

Try to create a finite matrix from row-major storage.

This is the public raw-storage boundary for matrices. Public compute methods parse stored rows into crate-internal proof-bearing types before arithmetic, including when crate-internal unchecked storage exists.

§Examples

use la_stack::prelude::*;

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

§Errors

Returns LaError::NonFinite with matrix coordinates for the first offending entry in row-major order when rows contains NaN or infinity.

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>::try_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 get_checked(&self, row: usize, col: usize) -> Result<f64, LaError>

Get an element, preserving index context on failure.

Prefer get for const or hot paths that only need Option-style absence. Use this method at public runtime boundaries where row, column, and dimension context should survive in a typed error.

§Examples

use la_stack::prelude::*;

let m = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
assert_eq!(m.get_checked(1, 0)?, 3.0);
assert_eq!(
    m.get_checked(2, 0),
    Err(LaError::IndexOutOfBounds {
        row: 2,
        col: 0,
        dim: 2,
    })
);

§Errors

Returns LaError::IndexOutOfBounds when either index is not < D.

Source

pub const fn set( &mut self, row: usize, col: usize, value: f64, ) -> Result<(), LaError>

Set a finite element with bounds checking.

§Examples

use la_stack::prelude::*;

let mut m = Matrix::<2>::zero();
assert_eq!(m.set(0, 1, 2.5), Ok(()));
assert_eq!(m.get(0, 1), Some(2.5));
assert_eq!(
    m.set(10, 0, 1.0),
    Err(LaError::IndexOutOfBounds {
        row: 10,
        col: 0,
        dim: 2,
    })
);

§Errors

Returns LaError::IndexOutOfBounds when either index is not < D. Returns LaError::NonFinite when value is NaN or infinity.

Source

pub const fn set_checked( &mut self, row: usize, col: usize, value: f64, ) -> Result<(), LaError>

Set an element, preserving index context on failure.

The matrix is mutated only when (row, col) is in bounds and value is finite.

§Examples

use la_stack::prelude::*;

let mut m = Matrix::<2>::zero();
m.set_checked(0, 1, 2.5)?;
assert_eq!(m.get_checked(0, 1)?, 2.5);

assert_eq!(
    m.set_checked(10, 0, 1.0),
    Err(LaError::IndexOutOfBounds {
        row: 10,
        col: 0,
        dim: 2,
    })
);

§Errors

Returns LaError::IndexOutOfBounds when either index is not < D. Returns LaError::NonFinite when value is NaN or infinity.

Source

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

Infinity norm (maximum absolute row sum).

§Non-finite handling

Public constructors and setters reject raw non-finite entries, but crate-internal unchecked storage can still contain NaN or infinity. inf_norm returns LaError::NonFinite if it encounters stored NaN/∞ or if a row sum overflows to a non-finite value.

Row sums are accumulated in f64 with ordinary addition. This method checks for overflowed accumulators, but it does not provide a certified absolute rounding bound for the returned norm.

§Examples

use la_stack::prelude::*;

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

// Raw NaN entries are rejected with coordinates.
assert_eq!(
    Matrix::<2>::try_from_rows([[f64::NAN, 1.0], [2.0, 3.0]]),
    Err(LaError::NonFinite {
        row: Some(0),
        col: 0,
    })
);

§Errors

Returns LaError::NonFinite when stored entries are NaN/infinity or a row sum overflows to NaN or infinity.

Source

pub fn is_symmetric(&self, rel_tol: Tolerance) -> Result<bool, LaError>

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, inf_norm(self)). This mirrors the predicate used internally by ldlt, so callers can pre-validate matrices that may come from untrusted sources.

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

The rel_tol argument is a Tolerance, so raw caller input must be finite and non-negative before it can reach this predicate. Use Tolerance::new or LaError::validate_tolerance when accepting a raw f64; negative, NaN, and infinite tolerances return LaError::InvalidTolerance.

§Overflow handling

A finite matrix can return LaError::NonFinite if computing the scaled symmetry tolerance overflows to NaN or infinity. If both stored entries are finite but their difference overflows to ±∞, the pair is reported as asymmetric.

§Examples

use la_stack::prelude::*;

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

let b = Matrix::<2>::try_from_rows([[4.0, 2.0], [3.0, 3.0]])?;
assert!(!b.is_symmetric(tol)?);

§Errors

Returns LaError::NonFinite when stored entries are NaN/infinity or computing the scaled symmetry tolerance overflows to NaN or infinity.

Source

pub fn first_asymmetry( &self, rel_tol: Tolerance, ) -> Result<Option<(usize, usize)>, LaError>

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, inf_norm(self)).

A finite matrix can return LaError::NonFinite if computing the scaled symmetry tolerance overflows to NaN or infinity. If both stored entries are finite but their difference overflows to ±∞, the pair is reported as asymmetric.

The rel_tol argument is a Tolerance, so raw caller input must be finite and non-negative before it can reach this predicate. Use Tolerance::new or LaError::validate_tolerance when accepting a raw f64; negative, NaN, and infinite tolerances return LaError::InvalidTolerance.

§Examples

use la_stack::prelude::*;

let a = Matrix::<3>::try_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
])?;
let tol = Tolerance::new(1e-12)?;
assert_eq!(a.first_asymmetry(tol)?, Some((1, 2)));
assert_eq!(Matrix::<3>::identity().first_asymmetry(tol)?, None);

§Errors

Returns LaError::NonFinite when stored entries are NaN/infinity or computing the scaled symmetry tolerance overflows to NaN or infinity.

Source

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

Compute an LU decomposition with partial pivoting.

§Examples

use la_stack::prelude::*;

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

let b = Vector::<2>::try_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);

The tol argument is a Tolerance, so raw caller input must be finite and non-negative before it can reach factorization. Use Tolerance::new or LaError::validate_tolerance when accepting a raw f64; negative, NaN, and infinite tolerances return LaError::InvalidTolerance.

§Errors

Returns LaError::Singular if, for some column k, the largest-magnitude candidate pivot in that column satisfies |pivot| <= tol (so no numerically usable pivot exists). Returns LaError::NonFinite if stored entries are NaN/infinity or an elimination intermediate overflows to NaN/∞ before it can be stored in the returned Lu.

Source

pub fn ldlt(self, tol: Tolerance) -> 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.

§Symmetry validation

The input matrix self must be symmetric — that is, self[i][j] == self[j][i] within the crate’s LDLT symmetry tolerance (1e-12, scaled like is_symmetric). This is a correctness invariant, not merely a performance hint, so asymmetric inputs return LaError::Asymmetric before factorization starts. If you need a general-purpose factorization that tolerates non-symmetric inputs, use lu instead.

The tol argument is a Tolerance, so raw caller input must be finite and non-negative before it can reach factorization. Use Tolerance::new or LaError::validate_tolerance when accepting a raw f64; negative, NaN, and infinite tolerances return LaError::InvalidTolerance.

§Examples

use la_stack::prelude::*;

// Note the symmetric layout: a[0][1] == a[1][0] == 2.0.
let a = Matrix::<2>::try_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>::try_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::NotPositiveSemidefinite if, for some step k, the required diagonal entry d = D[k,k] is negative. Returns LaError::Singular if 0 <= d <= tol, treating PSD degeneracy as singular/degenerate. Returns LaError::NonFinite if stored entries are NaN/infinity or factorization computes a non-finite intermediate. Returns LaError::Asymmetric if the input matrix is not symmetric.

Source

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

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

Returns Ok(Some(det)) for D ∈ {0, 1, 2, 3, 4}, Ok(None) for D ≥ 5. D = 0 returns Ok(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>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
assert_eq!(m.det_direct()?, Some(-2.0));

// 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());

§Errors

Returns LaError::NonFinite when stored entries are NaN/infinity or the closed-form determinant overflows to NaN or infinity.

Source

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

Floating-point 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).

Finite inputs return a floating-point determinant estimate in every dimension; this method does not surface LaError::Singular. Because it mixes closed-form paths from det_direct with an LU fallback, the returned value has no certified absolute error bound. Use det_errbound for D ≤ 4 bounds, or the exact determinant APIs when exact singularity classification or certified values matter. For D ≥ 5, the LU fallback only maps an exactly zero pivot to Ok(0.0). Use lu directly when you need tolerance-aware singularity detection or the pivot column.

§Examples

use la_stack::prelude::*;

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

§Errors

Returns LaError::NonFinite if stored entries are NaN/infinity, the LU fallback computes a non-finite factorization cell, or the determinant product overflows to NaN or infinity.

Source

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

Conservative absolute error bound for det_direct().

Returns Ok(Some(bound)) such that |det_direct() - det_exact| ≤ bound, or Ok(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>::try_from_rows([
    [1.0, 2.0, 3.0],
    [4.0, 5.0, 6.0],
    [7.0, 8.0, 9.0],
])?;
if let (Some(bound), Some(det_approx)) = (m.det_errbound()?, m.det_direct()?) {
    // If |det_approx| > bound, the sign is guaranteed correct.
    let sign_is_certified = det_approx.abs() > bound;
    assert!(!sign_is_certified);
}

§Adaptive precision pattern (requires `exact` feature)

use la_stack::prelude::*;

let m = Matrix::<3>::identity();
if let Some(bound) = m.det_errbound()? {
    if let Some(det) = m.det_direct()? {
        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()?;
        }
    }
} else {
    // D ≥ 5: no fast filter, use exact directly
    let sign = m.det_sign_exact()?;
}