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 arbitrary-precision 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.

§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), the Bareiss algorithm runs in exact BigRational arithmetic.

§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 fn inf_norm(&self) -> f64

Infinity norm (maximum absolute row sum).

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

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.

§Examples

use la_stack::prelude::*;

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.

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 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 matrix permanent using Shewchuk-style error analysis. 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();
}