Struct Matrix 

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

Finite fixed-size square matrix D×D, stored inline.

Matrix is designed for small, robustness-sensitive systems where stack allocation and const-generic dimensions are useful. For large, dynamic, sparse, or parallel workloads, prefer a broader linear-algebra crate such as nalgebra or faer.

Public construction and mutation reject NaN and infinity through try_from_rows, set, and set_checked. The storage field is private, so a Matrix value carries the invariant that every stored entry is finite. Algorithms therefore do not re-scan stored entries at every use; user-visible non-finite errors come from construction/mutation boundaries or from values computed during arithmetic, such as overflowed elimination or determinant intermediates.

Direct field construction is intentionally unavailable to downstream callers:

use la_stack::Matrix;

let _ = Matrix::<2> {
    rows: [[1.0, f64::NAN], [0.0, 1.0]],
};

Implementations

impl<const D: usize> Matrix<D>

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::DeterminantScaleOverflow if determinant scaling overflows the internal exponent representation.

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

Exact determinant converted to f64.

Requires the exact Cargo feature.

Computes the exact determinant with the same integer Bareiss core used by det_exact, then converts the exact scaled integer result to f64 only if the result is exactly representable as a finite binary64 value.

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::DeterminantScaleOverflow if determinant scaling overflows the internal exponent representation.

Returns LaError::Unrepresentable if the exact determinant cannot be represented exactly as a finite f64.

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

Exact determinant rounded to f64.

Requires the exact Cargo feature.

Computes the exact determinant with the same integer Bareiss core used by det_exact, then rounds the exact value to a finite binary64 value. Unlike det_exact_f64, this method is intentionally lossy and may round non-dyadic or underflowing nonzero exact determinants.

Examples

use la_stack::prelude::*;

let m = Matrix::<2>::try_from_rows([
    [1.0 + f64::EPSILON, 0.0],
    [0.0, 1.0 - f64::EPSILON],
])?;

assert_eq!(
    m.det_exact_f64(),
    Err(LaError::Unrepresentable {
        index: None,
        reason: UnrepresentableReason::RequiresRounding,
    })
);
assert_eq!(m.det_exact_rounded_f64()?.to_bits(), 1.0f64.to_bits());

Errors

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

Returns LaError::Unrepresentable if the rounded determinant would be NaN or infinite.

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::Singular if the matrix is exactly singular.

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 f64 only if that component is exactly representable as a finite binary64 value.

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::Singular if the matrix is exactly singular. Returns LaError::Unrepresentable if any component of the exact solution cannot be represented exactly as a finite f64.

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

Exact linear system solve rounded to f64.

Requires the exact Cargo feature.

Computes the exact BigRational solution via solve_exact and rounds each component to a finite binary64 value. Unlike solve_exact_f64, this method is intentionally lossy and may round non-dyadic or underflowing nonzero exact components.

Examples

use la_stack::prelude::*;

let a = Matrix::<1>::try_from_rows([[3.0]])?;
let b = Vector::<1>::try_new([1.0])?;

assert_eq!(
    a.solve_exact_f64(b),
    Err(LaError::Unrepresentable {
        index: Some(0),
        reason: UnrepresentableReason::RequiresRounding,
    })
);
assert_eq!(a.solve_exact_rounded_f64(b)?.into_array(), [1.0 / 3.0]);

Errors

Returns LaError::Singular if the matrix is exactly singular. Returns LaError::Unrepresentable if any rounded component would be NaN or infinite.

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::DeterminantScaleOverflow if determinant scaling overflows the internal exponent representation.

impl<const D: usize> Matrix<D>

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. Successful construction makes the returned Matrix a finite-storage proof.

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.

pub const fn zero() -> Self

All-zeros finite matrix.

Examples

use la_stack::prelude::*;

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

pub const fn identity() -> Self

Finite 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));

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

Get a finite 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);

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

Get a finite 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.

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.

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.

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 Matrix values are finite by construction. inf_norm returns LaError::NonFinite with the matrix cell whose addition first makes a row sum non-finite.

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 with matrix coordinates when a row sum overflows to NaN or infinity.

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 with matrix coordinates 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 with matrix coordinates when computing the scaled symmetry tolerance overflows to NaN or infinity.

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 with matrix coordinates 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 with matrix coordinates when computing the scaled symmetry tolerance overflows to NaN or infinity.

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

Compute an LU decomposition with partial pivoting.

D = 0 follows the empty-matrix convention: factorization succeeds, Lu::det returns 1.0, and solving a length-zero right-hand side returns a length-zero Vector.

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_SINGULAR_TOL)?;

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

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

Empty matrices use the standard empty-product convention:

use la_stack::prelude::*;

let lu = Matrix::<0>::zero().lu(DEFAULT_SINGULAR_TOL)?;

assert_eq!(lu.det()?, 1.0);
assert!(lu.solve(Vector::<0>::zero())?.into_array().is_empty());

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 an elimination intermediate overflows to NaN/∞ before it can be stored in the returned Lu.

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

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

D = 0 follows the empty-matrix convention: factorization succeeds, Ldlt::det returns 1.0, and solving a length-zero right-hand side returns a length-zero Vector.

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(b)?.into_array();
assert!((x[0] - (-0.125)).abs() <= 1e-12);
assert!((x[1] - 0.75).abs() <= 1e-12);

Empty matrices use the standard empty-product convention:

use la_stack::prelude::*;

let ldlt = Matrix::<0>::zero().ldlt(DEFAULT_SINGULAR_TOL)?;

assert_eq!(ldlt.det()?, 1.0);
assert!(ldlt.solve(Vector::<0>::zero())?.into_array().is_empty());

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 factorization computes a non-finite intermediate. Returns LaError::Asymmetric if the input matrix is not symmetric.

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 the closed-form determinant overflows to NaN or infinity.

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 the LU fallback computes a non-finite factorization cell, or the determinant product overflows to NaN or infinity.

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()?;
}

Errors

Returns LaError::NonFinite when the bound computation overflows to NaN or infinity.

Trait Implementations

impl<const D: usize> Clone for Matrix<D>

fn clone(&self) -> Matrix<D>

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl<const D: usize> Copy for Matrix<D>

impl<const D: usize> Debug for Matrix<D>

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl<const D: usize> Default for Matrix<D>

fn default() -> Self

Returns the “default value” for a type.

impl<const D: usize> PartialEq for Matrix<D>

fn eq(&self, other: &Matrix<D>) -> bool

Tests for self and other values to be equal, and is used by ==.

fn ne(&self, other: &Rhs) -> bool

Tests for !=. The default implementation is almost always sufficient, and should not be overridden without very good reason.

impl<const D: usize> StructuralPartialEq for Matrix<D>