Module exact 

Available on crate features exact only.

Exact arithmetic over stored binary64 and rational inputs.

The default build has zero runtime dependencies. Enable the optional exact Cargo feature to add exact arithmetic methods using arbitrary-precision rationals (this pulls in num-bigint, num-rational, and num-traits for BigRational):

See the crate-level installation instructions for Cargo configuration.

The feature exposes two deliberate input domains:

• Matrix<D> / Vector<D> store finite binary64 inputs. Their exact methods treat each stored bit pattern as its exact rational value, so the determinant or solve stage introduces no further roundoff. They cannot recover information already lost before construction.
• RationalMatrix<D> / RationalVector<D> accept coefficients already assembled as BigRational. They preserve derived differences, squared norms, affine coefficients, and other rational expressions without an intermediate f64 conversion.

Determinants:

• det_exact() — returns the exact determinant as a BigRational
• det_exact_f64() — returns the exact determinant as f64 only when it is exactly representable (or LaError::Unrepresentable otherwise)
• det_exact_rounded_f64() — returns the exact determinant rounded to a finite f64 using IEEE 754 round-to-nearest, ties-to-even
• det_sign_exact() — infallibly returns the provably correct DeterminantSign variant (Negative, Zero, or Positive)

Linear system solve:

• solve_exact(b) — solves Ax = b exactly, returning a RationalVector<D>
• solve_exact_f64(b) — solves Ax = b exactly, returning Vector<D> only when every component is exactly representable as f64
• solve_exact_rounded_f64(b) — solves Ax = b exactly, returning each component rounded to finite f64 using IEEE 754 round-to-nearest, ties-to-even
• ExactF64Conversion — converts an existing exact determinant or solution under the strict or rounded contract without repeating exact elimination

Already-exact rational input:

• RationalMatrix::det_sign() — returns the exact sign without constructing a rational determinant
• RationalMatrix::det() — returns the exact BigRational determinant
• RationalMatrix::solve(&rhs) — returns a RationalVector<D> exact solution
• try_with_rational_matrix! — dispatches a runtime-selected dimension through D=8 to a const-generic rational matrix on stable Rust

The Matrix::det_exact* value and conversion methods return LaError::DeterminantScaleOverflow if their aggregate power-of-two scaling exceeds the internal exponent representation. RationalMatrix::det() is infallible because it clears rational row denominators without an exponent-scale conversion. The exact solve methods for both input domains return LaError::Singular with SingularityReason::Exact when the stored matrix is exactly singular.

For exact-to-f64 output, strict conversions use UnrepresentableReason::RequiresRounding when explicit rounding can produce a finite value and UnrepresentableReason::NotFinite otherwise. Rounded conversions opt into nearest-even rounding but still report NotFinite when no finite f64 exists.

Preserving rational inputs

The following 5×5 system has exact determinant 2^-60. Its exact rational inputs therefore produce a unique solution through the general Bareiss path. Supplying the same coefficients as f64 inputs loses the 2^-60 perturbation at 1.0, making the leading rows identical and the binary64 system singular.

use core::assert_matches;

use la_stack::prelude::*;

fn main() -> Result<(), LaError> {
    // This is far below one binary64 ULP at 1.0, so 1.0 + 2^-60 rounds to 1.0.
    let epsilon = BigRational::new(1.into(), (1_u64 << 60).into());
    let one = BigRational::from_integer(1.into());
    let zero = BigRational::from_integer(0.into());

    // The leading block is [[1, 1], [1, 1 + 2^-60]]. The remaining diagonal
    // extends the example to D=5, where the general Bareiss path is used.
    let matrix = RationalMatrix::<5>::try_from_fn(|row, col| match (row, col) {
        (0, 0 | 1) | (1, 0) => one.clone(),
        (1, 1) => &one + &epsilon,
        _ if row == col => one.clone(),
        _ => zero.clone(),
    })?;
    assert_eq!(matrix.det_sign(), DeterminantSign::Positive);
    assert_eq!(matrix.det(), epsilon);

    let rhs = RationalVector::try_new([
        zero,
        -&epsilon,
        BigRational::from_integer(2.into()),
        BigRational::from_integer(3.into()),
        BigRational::from_integer(4.into()),
    ])?;
    let exact_solution = matrix.solve(&rhs)?;
    assert_eq!(
        exact_solution.as_array(),
        &[
            BigRational::from_integer(1.into()),
            BigRational::from_integer((-1).into()),
            BigRational::from_integer(2.into()),
            BigRational::from_integer(3.into()),
            BigRational::from_integer(4.into()),
        ]
    );

    // Supplying the same coefficients as f64 inputs destroys the perturbation
    // and makes the matrix singular, even though the exact solution is integral.
    let epsilon_f64 = epsilon.try_to_f64()?;
    assert_eq!((1.0 + epsilon_f64).to_bits(), 1.0_f64.to_bits());
    let f64_matrix = Matrix::<5>::try_from_rows([
        [1.0, 1.0, 0.0, 0.0, 0.0],
        [1.0, 1.0 + epsilon_f64, 0.0, 0.0, 0.0],
        [0.0, 0.0, 1.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 1.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 1.0],
    ])?;
    let f64_rhs = Vector::<5>::try_new([0.0, -epsilon_f64, 2.0, 3.0, 4.0])?;
    let f64_solve = f64_matrix
        .lu(DEFAULT_SINGULAR_TOL)
        .and_then(|lu| lu.solve(f64_rhs));
    assert_matches!(
        f64_solve,
        Err(LaError::Singular { .. })
    );
    Ok(())
}

Stored binary64 inputs and output conversion

use la_stack::prelude::*;

fn main() -> Result<(), LaError> {
    // Exact determinant
    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],
    ])?;
    assert_eq!(m.det_sign_exact(), DeterminantSign::Zero); // exactly singular

    let det = m.det_exact()?;
    assert_eq!(det, BigRational::from_integer(0.into())); // exact zero
    let det_f64 = det.try_to_f64()?;
    assert_eq!(det_f64, 0.0);

    // If strict exact-to-f64 conversion would require rounding, opt in
    // explicitly with the rounded API.
    let inexact = Matrix::<2>::try_from_rows([
        [1.0 + f64::EPSILON, 0.0],
        [0.0, 1.0 - f64::EPSILON],
    ])?;
    let exact_det = inexact.det_exact()?;
    let rounded_det = match exact_det.try_to_f64() {
        Ok(det) => det,
        Err(err) if err.requires_rounding() => exact_det.to_rounded_f64()?,
        Err(err) => return Err(err),
    };
    assert_eq!(rounded_det.to_bits(), 1.0f64.to_bits());

    // If the exact determinant cannot fit in f64, keep the BigRational value.
    let big = f64::MAX / 2.0;
    let huge = Matrix::<3>::try_from_rows([
        [0.0, 0.0, 1.0],
        [big, 0.0, 1.0],
        [0.0, big, 1.0],
    ])?;
    let huge_det = huge.det_exact()?;
    assert_eq!(
        huge_det
            .try_to_f64()
            .err()
            .and_then(|err| err.unrepresentable_reason()),
        Some(UnrepresentableReason::NotFinite)
    );
    println!("exact determinant = {huge_det}");

    // Exact linear system solve
    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 exact_x = a.solve_exact(b)?;
    let x = exact_x.try_to_f64()?.into_array();
    assert!((x[0] - 1.0).abs() <= f64::EPSILON);
    assert!((x[1] - 2.0).abs() <= f64::EPSILON);

    Ok(())
}

With the exact feature enabled, RationalMatrix, RationalVector, DeterminantSign, ExactF64Conversion, BigInt, and BigRational are re-exported from the crate root and prelude, alongside the most commonly needed num-traits items (FromPrimitive, ToPrimitive, Signed). This lets consumers construct exact values (BigRational::from_f64, from_i64), query sign (is_positive / is_negative), and convert back to f64 (try_to_f64, to_rounded_f64, or the raw to_f64) with a single use la_stack::prelude::*; — no need to add num-bigint, num-rational, or num-traits to their own Cargo.toml. Use DeterminantSign::as_i8() only when numeric −1/0/+1 interoperability is required.

For det_sign_exact(), D ≤ 4 matrices first use a fast f64 filter (error-bounded det_direct_with_errbound()) when its rounded intermediates stay in the normal range or are exact structural zeros. An inconclusive filter falls back to the same direct determinant expansion in BigInt. D ≥ 5 skips the closed-form filter and uses fraction-free Bareiss elimination in BigInt. Because Matrix stores only finite entries, arithmetic range failures in the filter are inconclusive rather than errors and the exact fallback is total.

A five-dimensional rational solve

use core::assert_matches;

use la_stack::prelude::*;

// A tridiagonal exact matrix with determinant 6.
let matrix = RationalMatrix::<5>::try_from_fn(|row, col| {
    BigRational::from_integer(if row == col {
        2.into()
    } else if row.abs_diff(col) == 1 {
        1.into()
    } else {
        0.into()
    })
})?;
let numerators = [4, 8, 12, 16, 14];
let rhs = RationalVector::try_from_fn(|row| {
    BigRational::new(numerators[row].into(), 3.into())
})?;
let solution = matrix.solve(&rhs)?;
let expected = [1, 2, 3, 4, 5].map(|n| BigRational::new(n.into(), 3.into()));
assert_eq!(solution.as_array(), &expected);

// Keep the exact solution until the caller explicitly opts into rounding.
assert_matches!(
    solution.try_to_f64(),
    Err(LaError::Unrepresentable {
        index: Some(0),
        reason: UnrepresentableReason::RequiresRounding,
        ..
    })
);
let rounded = solution.to_rounded_f64()?;
assert_eq!(rounded.as_array(), &[1.0 / 3.0, 2.0 / 3.0, 1.0, 4.0 / 3.0, 5.0 / 3.0]);

Rational dimension dispatch

try_with_rational_matrix! selects a concrete RationalMatrix<N> for dimensions 0 through MAX_RATIONAL_MATRIX_DISPATCH_DIM (8). Larger dimensions return LaError::UnsupportedDimension, converted through From<LaError> into the closure’s declared error type. The macro preserves the const-generic representation.

Rational constructors include try_from_rows / try_new and try_from_fn; as_rows / as_array and get borrow their stored values, while into_rows / into_array return the owned arrays.

Adaptive determinant filtering

This example requires exact and illustrates a custom filter with exact fallback. Use Matrix::det_sign_exact directly when no custom filtering policy is needed.

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> {
    let identity = Matrix::<3>::identity();
    assert_eq!(
        adaptive_det_sign(&identity),
        DeterminantSign::Positive
    );

    // A zero determinant cannot pass the f64 sign filter, so this exercises
    // the exact fallback.
    let singular = Matrix::<3>::try_from_rows([
        [1.0, 2.0, 3.0],
        [4.0, 5.0, 6.0],
        [7.0, 8.0, 9.0],
    ])?;
    assert_eq!(adaptive_det_sign(&singular), DeterminantSign::Zero);

    // The f64 filter overflows for this finite matrix, but the exact fallback
    // still resolves its positive determinant sign.
    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(())
}