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 and set. 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 exact for the stored binary64 entries. It cannot recover precision lost before matrix construction.

When to use

Use this when you need the exact determinant value — for example, volume computation over stored coordinates or distinguishing simplices that are exactly degenerate at those coordinates 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-scaled 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. The candidate conversion follows IEEE 754 round-to-nearest, ties-to-even, but is returned only when no rounding is required.

When callers also need the exact value or may recover with explicit rounding, compute det_exact once and use ExactF64Conversion on the returned BigRational.

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-scaled core used by det_exact, then rounds the exact value to a finite binary64 value using IEEE 754 round-to-nearest, ties-to-even. Unlike det_exact_f64, this method is intentionally lossy and may round non-dyadic or underflowing nonzero exact determinants.

Examples

use core::assert_matches;
use la_stack::prelude::*;

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

assert_matches!(
    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 rounding cannot produce a finite f64.

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 exact for the stored binary64 entries. It cannot recover precision lost before matrix or vector construction.

When to use

Use this when you need a solution exact for the stored inputs — for example, circumcenter computation over stored coordinates 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 independently scaled to their own power-of-two bases so both sides of the augmented system (A | b) become integer-valued without needless cross-side shifts. After solving that integer system, the exact power-of-two ratio between the RHS and matrix scales is restored. 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. The candidate conversion follows IEEE 754 round-to-nearest, ties-to-even, but is returned only when no rounding is required.

When callers also need the exact solution or may recover with explicit rounding, compute solve_exact once and use ExactF64Conversion on the returned array.

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 using IEEE 754 round-to-nearest, ties-to-even. Unlike solve_exact_f64, this method is intentionally lossy and may round non-dyadic or underflowing nonzero exact components.

Examples

use core::assert_matches;
use la_stack::prelude::*;

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

assert_matches!(
    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 rounding any component cannot produce a finite f64.

pub fn det_sign_exact(&self) -> DeterminantSign

Exact determinant sign using adaptive-precision arithmetic.

Requires the exact Cargo feature.

Returns DeterminantSign::Positive, DeterminantSign::Negative, or DeterminantSign::Zero according to the determinant that is exact for the stored binary64 entries. This cannot recover precision lost before matrix construction.

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, exact integer arithmetic computes the sign without constructing any BigRational values: direct BigInt expansions for D ≤ 4 and Bareiss elimination for D ≥ 5.

When to use

Use this when the sign of the determinant over the stored entries must be correct regardless of floating-point conditioning (e.g. geometric predicates on near-degenerate stored coordinates). 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(), DeterminantSign::Zero);

assert_eq!(Matrix::<3>::identity().det_sign_exact(), DeterminantSign::Positive);

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 as_rows(&self) -> &[[f64; D]; D]

Borrow the finite row-major backing array.

The returned view is tied to this Matrix, so callers can inspect the canonical storage without copying it or bypassing the finite-value invariant.

Examples

use la_stack::prelude::*;

let matrix = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
assert_eq!(matrix.as_rows(), &[[1.0, 2.0], [3.0, 4.0]]);

A live view keeps the matrix immutably borrowed, so validated mutation cannot occur until the view is no longer used:

use la_stack::Matrix;

let mut matrix = Matrix::<2>::identity();
let rows = matrix.as_rows();
assert!(matrix.set(0, 0, 5.0).is_ok());
assert_eq!(rows[0][0], 1.0);

pub const fn into_rows(self) -> [[f64; D]; D]

Consume this matrix and return its finite row-major backing array.

Examples

use la_stack::prelude::*;

let matrix = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
assert_eq!(matrix.into_rows(), [[1.0, 2.0], [3.0, 4.0]]);

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, row: usize, col: 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 try_get(&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 core::assert_matches;
use la_stack::prelude::*;

let m = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
assert_eq!(m.try_get(1, 0)?, 3.0);
assert_matches!(
    m.try_get(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 core::assert_matches;
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_matches!(
    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 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 core::assert_matches;
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_matches!(
    Matrix::<2>::try_from_rows([[f64::NAN, 1.0], [2.0, 3.0]]),
    Err(LaError::NonFinite {
        location: NonFiniteLocation::MatrixCell { row: 0, col: 0, .. },
        origin: NonFiniteOrigin::Input,
        ..
    })
);

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 approximately 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 is a diagnostic predicate for applications that have an approximation-specific symmetry threshold. It is not the precondition used by ldlt, which requires exact mirrored-entry equality so the returned factors represent the original matrix.

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::try_new 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::try_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 approximate 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)). It is intentionally distinct from the exact equality required by ldlt.

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::try_new 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::try_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. Partial pivoting is a practical finite-precision strategy, not a certified accuracy guarantee; see REFERENCES.md [1-3, 11-12].

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::try_new 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 exactly symmetric positive-definite matrices such as nonsingular Gram matrices. Computed zero and tolerance-small positive pivots are diagnosed as singular rather than returned in a usable factorization. Because pivots are computed in binary64, success is not an exact proof that the stored matrix is positive definite. See REFERENCES.md [4-6, 11-12] for Cholesky/LDLT background and the pivoted symmetric-indefinite alternative.

Symmetry validation

The input matrix self must be exactly symmetric: every mirrored pair must satisfy self[i][j] == self[j][i]. IEEE-754 signed zeros compare equal and are therefore accepted. Exact equality is a correctness invariant, not merely a performance hint: LDLT reads only the lower triangle, so accepting an approximate mismatch would factor a different operator than the matrix supplied by the caller. Asymmetric inputs return LaError::Asymmetric with an allowed absolute difference of 0.0 before factorization starts.

is_symmetric remains available as a tolerance-based diagnostic, but Ok(true) from that method does not establish this exact LDLT precondition. If you need a general-purpose factorization for a non-symmetric matrix, 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::try_new 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 a pivot is negative or a zero pivot retains a non-zero coupling below it. Returns LaError::Singular if a zero pivot has no remaining coupling, or if a positive pivot satisfies d <= tol, treating PSD degeneracy as singular. 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).

Because this method 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 zero-tolerance LU fallback surfaces LaError::Singular when elimination cannot produce a non-zero pivot. Floating-point elimination cannot in general distinguish an exactly singular matrix from a non-singular matrix whose intermediate pivot rounded to zero, so this method never converts that numerical failure into an exact 0.0 result.

Examples

use la_stack::prelude::*;

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

The LU fallback accumulates its diagonal product with power-of-two scaling, so factor order cannot cause premature overflow or underflow in the final product. Elimination intermediates remain subject to binary64 rounding and range limits.

Errors

Returns LaError::Singular if the D ≥ 5 LU fallback cannot produce a non-zero pivot, including when a non-zero mathematical intermediate rounds to zero during elimination. Returns LaError::NonFinite if a D ≤ 4 closed-form result is non-finite, if the LU fallback computes a non-finite factorization cell, or if its final scaled determinant cannot be represented as a finite f64.

pub const fn det_direct_with_errbound( &self, ) -> Result<Option<DeterminantWithErrorBound>, LaError>

Evaluate det_direct() and its absolute error bound together.

Returns Ok(Some(result)) for D ≤ 4 when the relative-error analysis is valid. The result contains the closed-form determinant and a bound such that |result.determinant() - det_exact| ≤ result.absolute_error_bound(). Returns Ok(None) when gradual underflow could invalidate that analysis or for D ≥ 5, where no closed-form bound is available.

This is the preferred API when both values are needed: it evaluates the determinant arithmetic tree once, then computes the matching bound for the same matrix within that call.

Examples

use la_stack::prelude::*;

let matrix = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
if let Some(estimate) = matrix.det_direct_with_errbound()? {
    assert_eq!(estimate.determinant(), -2.0);
    assert!(estimate.absolute_error_bound() >= 0.0);
}

Errors

Returns LaError::NonFinite when the determinant or bound computation overflows to NaN or infinity. Underflow-sensitive finite computations return Ok(None) because they remain valid inputs for an exact fallback.

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 when every rounded intermediate used by the closed-form determinant and bound is normal (or an exact structural zero). Returns Ok(None) when gradual underflow could invalidate the relative-error analysis, or 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 det_direct_with_errbound when the determinant and bound are both needed. This accessor is convenient when only the bound is needed.

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) = m.det_errbound()? {
    assert!(bound >= 0.0);
}

Adaptive precision pattern (requires exact feature)

use la_stack::prelude::*;

fn adaptive_det_sign<const D: usize>(
    matrix: &Matrix<D>,
) -> DeterminantSign {
    if let Ok(Some(estimate)) = matrix.det_direct_with_errbound() {
        if estimate.determinant().abs() > estimate.absolute_error_bound() {
            return if estimate.determinant() > 0.0 {
                DeterminantSign::Positive
            } else {
                DeterminantSign::Negative
            };
        }
    }

    matrix.det_sign_exact()
}

fn main() -> Result<(), LaError> {
    assert_eq!(
        adaptive_det_sign(&Matrix::<3>::identity()),
        DeterminantSign::Positive
    );

    let big = f64::MAX / 2.0;
    let overflowing = Matrix::<3>::try_from_rows([
        [0.0, 0.0, 1.0],
        [big, 0.0, 1.0],
        [0.0, big, 1.0],
    ])?;
    assert_eq!(
        adaptive_det_sign(&overflowing),
        DeterminantSign::Positive
    );
    Ok(())
}

Errors

Returns LaError::NonFinite when the bound computation overflows to NaN or infinity. Underflow-sensitive finite computations return Ok(None) instead because they are valid inputs for an exact fallback.

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

Equality operator ==.

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

Inequality operator !=.

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