Skip to content

feat: Accept exact rational matrices for determinant and solve operations #216

Description

@acgetchell

Summary

Add an exact-input matrix/vector surface behind the existing exact feature so
callers can supply preassembled BigRational coefficients and obtain exact
determinant signs, determinant values, and linear-system solutions without an
intermediate f64 round trip.

Matrix<D> already treats each stored finite f64 as its exact IEEE-754 value
for det_exact, det_sign_exact, and solve_exact. That is sufficient when
every matrix entry came directly from one binary64 value. It is insufficient
when a caller must form derived coefficients exactly first, such as a
difference, squared norm, affine coefficient, or rational objective. Today
those callers must either round the derived value back through f64 or maintain
their own rational elimination implementation.

This request is independent of the generic_const_exprs redesign in #123 and
can be implemented on stable Rust.

Downstream motivation: delaunay

The delaunay crate needs exact-input matrices in three cold paths:

  1. Relative in-sphere predicates form each row from exact binary64 rationals as
    (p_i - p_0, sum_j (p_ij - p_0j)^2), then need only the determinant sign.
    Rounding the subtraction or squared norm to f64 before exact elimination
    can change a near-degenerate predicate sign.
  2. Circumcenter recovery forms the relative perpendicular-bisector system
    rationally when finite binary64 subtraction or squaring overflows or the
    ordinary LU path is numerically singular, then needs an exact solution that
    is rounded only once at the public coordinate boundary.
  3. Exact simplex-intersection and revised-simplex certification repeatedly
    solve runtime-selected rational basis systems.

For Delaunay's current exact envelope through geometric dimension D = 6, the
required square systems are at most D + 2 = 8. The API therefore needs either
bounded runtime dispatch through at least dimension 8 or an exact workspace
whose validated runtime dimension does not depend on generic_const_exprs.

Requested contract

The exact-input type should:

  • accept BigRational entries without converting through f64;
  • validate square shape once and prevent later shape mismatch;
  • provide determinant sign without constructing an exact determinant value
    when only DeterminantSign is requested;
  • provide the exact determinant as BigRational when requested;
  • solve A x = b from rational matrix and RHS inputs, returning rational
    components;
  • preserve typed singularity, dimension, and conversion diagnostics;
  • remain gated by the existing exact feature;
  • offer explicit strict and rounded conversion through the existing
    ExactF64Conversion vocabulary rather than performing implicit rounding;
  • support Delaunay's stable-Rust D + 1 and runtime-selected basis dimensions.

One possible surface is:

pub struct RationalMatrix<const D: usize> { /* private exact storage */ }
pub struct RationalVector<const D: usize> { /* private exact storage */ }

impl<const D: usize> RationalMatrix<D> {
    pub fn from_rows(rows: [[BigRational; D]; D]) -> Self;
    pub fn from_fn(f: impl FnMut(usize, usize) -> BigRational) -> Self;

    pub fn det_sign(&self) -> DeterminantSign;
    pub fn det(&self) -> BigRational;
    pub fn solve(&self, rhs: &RationalVector<D>)
        -> Result<RationalVector<D>, LaError>;
}

A parallel try_with_rational_matrix! dispatcher through at least dimension 8
would let generic downstream code request D + 1 on stable Rust. A validated
runtime-sized RationalMatrix/borrowed view is also acceptable if it better
fits the exact backend. The important contract is exact rational input, one-time
shape validation, and no binary64 reconstruction.

The existing Matrix<D> should remain the finite binary64 type. Documentation
can then advertise two deliberate input domains:

  • Matrix<D> / Vector<D>: finite f64 input, ordinary floating operations,
    plus exact operations over the stored IEEE-754 values;
  • RationalMatrix / RationalVector: already-exact input and exact output,
    with any conversion to f64 explicit.

This is preferable to making the current Copy-able finite Matrix<D> generic
over scalar storage: BigRational has different ownership, allocation, and
invariant properties, and a separate type keeps those contracts visible.

Suggested implementation

Avoid rational Gaussian elimination in the cubic phase:

  1. Clear denominators with a positive per-row scale. This preserves determinant
    sign; exact determinant values divide by the product of row scales; linear
    solves scale the corresponding RHS entry by the same factor and preserve the
    solution.
  2. Reuse the existing fraction-free BigInt Bareiss elimination and pivoting.
  3. Convert only the upper-triangular solve state to BigRational for
    back-substitution, matching the current Matrix<D>::solve_exact strategy.

For Delaunay-derived coefficients the rationals are dyadic, so row denominator
clearing can use the largest required power of two rather than a general costly
cross-row denominator product.

Acceptance criteria

  • Exact rational matrices expose sign-only determinant, exact determinant, and
    exact solve operations behind exact.
  • Derived rational entries that cannot make an exact BigRational -> f64 -> BigRational round trip are accepted and computed without rounding.
  • Determinant signs agree with an independent BigRational reference for known,
    singular, pivoting, and property-generated matrices through dimension 8.
  • Exact solve satisfies A * x == b rationally for nonsingular generated
    systems and returns LaError::Singular for singular systems.
  • Row-denominator clearing is covered for mixed signs, non-reduced input,
    different per-row denominators, and dyadic exponents.
  • The crate still builds without exact, and existing Matrix<D> behavior is
    unchanged.
  • Documentation and examples explain the distinction between exact operations
    over stored f64 inputs and operations over already-exact rational inputs.
  • Benchmarks compare the rational-input backend with straightforward
    BigRational Gaussian elimination for dimensions used downstream.

Downstream cleanup enabled

Once released, delaunay can remove:

  • its generic rational determinant-sign implementation;
  • two duplicate rational Gaussian solvers;
  • the convert-to-f64, exact-round-trip check used only to reach
    Matrix<D>::solve_exact;
  • associated pivoting and singularity plumbing.

Delaunay will still use BigRational where exact values are themselves the
algorithmic state: assembling derived predicate/circumcenter coefficients,
Simulation-of-Simplicity polynomial coefficients, and exact linear-program
tableaux and witnesses. The proposed API retires most hand-written exact linear
algebra, not all exact arithmetic.

Metadata

Metadata

Assignees

No one assigned

    Labels

    apienhancementNew feature or requestperformancePerformance related issuesrustPull requests that update rust code

    Projects

    No projects

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions