# Mathematical basis
## Contents
- [Introduction](#introduction)
- [Choosing an API](#choosing-an-api)
- [Geometry relationship and scope](#geometry-relationship-and-scope)
- [Represented values and arithmetic model](#represented-values-and-arithmetic-model)
- [Tolerances and typed errors](#tolerances-and-typed-errors)
- [Linear algebra algorithms](#linear-algebra-algorithms)
- [Certified fixed-vector reductions](#certified-fixed-vector-reductions)
- [Determinants and certified sign filtering](#determinants-and-certified-sign-filtering)
- [Derivation of the returned determinant bound](#derivation-of-the-returned-determinant-bound)
- [Exact arithmetic over binary64 inputs](#exact-arithmetic-over-binary64-inputs)
- [Exact arithmetic over rational inputs](#exact-arithmetic-over-rational-inputs)
- [Exact sparse polynomial determinants](#exact-sparse-polynomial-determinants)
- [Exact-to-binary64 conversion](#exact-to-binary64-conversion)
- [Gram matrices and geometric measures](#gram-matrices-and-geometric-measures)
- [Interval linear solves](#interval-linear-solves)
- [LDLT without pivoting](#ldlt-without-pivoting)
- [LU with partial pivoting](#lu-with-partial-pivoting)
- [Outward-rounded interval expressions](#outward-rounded-interval-expressions)
- [Scalar arithmetic and finite bounds](#scalar-arithmetic-and-finite-bounds)
- [Division and reciprocals](#division-and-reciprocals)
- [Determinant enclosures and sign evidence](#determinant-enclosures-and-sign-evidence)
- [Rectangular column selection](#rectangular-column-selection)
- [Scaled determinant products](#scaled-determinant-products)
- [Scaled Euclidean vector norm](#scaled-euclidean-vector-norm)
- [Unsigned vector angles](#unsigned-vector-angles)
## Introduction
`la-stack` provides fixed-dimension numerical linear algebra over two deliberate
point-value input domains plus one bounded layer. `Matrix<D>` and `Vector<D>`
store finite IEEE 754 binary64 values; their default algorithms operate in
binary64 and are therefore approximate. `Interval` and `IntervalMatrix<D>`
enclose exact-real expression values between outward-rounded finite binary64
bounds when interval operations are used throughout expression assembly.
Lifting an already-rounded `Matrix` encloses only its stored values. The
optional `exact` feature can either lift stored binary64 values to
the exact rational numbers they represent or accept caller-supplied
`BigRational` values through `RationalMatrix<D>` and `RationalVector<D>`.
This document separates three questions that are easy to conflate:
1. Which mathematical factorization or determinant identity is being used?
2. Which part of the computation is rounded in binary64?
3. What does a tolerance or typed error prove about the supplied values?
Reference numbers point to [REFERENCES.md](../REFERENCES.md).
## Choosing an API
| General floating solve | `lu(tol)` then `solve` | Approximate; absolute pivot policy |
| Positive-definite floating solve | `ldlt(tol)` then `solve` | Exact symmetry; computed pivots must exceed tolerance; success is not a certificate |
| Floating determinant, any `D` | `det` | No certified bound; zero is not exact singularity |
| `D ≤ 4` error-bounded determinant/sign test | `det_direct_with_errbound` | Sign is certified when estimate magnitude exceeds bound; otherwise inconclusive |
| Derived-expression determinant sign through `D ≤ 7` | `IntervalMatrix::det_sign` | Outward-rounded proof; overlap with zero is explicitly inconclusive |
| Solution sets with uncertain coefficients or RHS | `IntervalMatrix::solve` | Every admitted exact-real solution is enclosed on success; zero-containing pivots are inconclusive |
| Exact determinant sign | `det_sign_exact` | Exact for stored binary64 entries |
| Exact determinant value or solve | `det_exact`, `solve_exact` | Exact for represented inputs |
| Exact operations over preassembled rationals | `RationalMatrix::det_sign`, `det`, `solve` | No intermediate binary64 reconstruction |
| Exact rational dot products | `RationalVector::dot`, `rational_dot` | Exact canonical result; checked slices have no dimension cap |
| Exact polynomial determinant and first nonzero coefficient | `polynomial_determinant`, `SparsePolynomial::leading_term` | Orders 0–8; increasing `u128` exponents; explicit work limits |
| Rectangular column selection | `column_basis` | Explicit absolute tolerance; provisional numerical rank, requiring independent certification |
| Binary64 output from an exact result | Strict or rounded conversions | Strict conversion forbids rounding |
| Unsigned vector angle in radians | `VectorAngle::angle`, `Vector::angle` | Finite nonzero vectors; rounded, without a certified bound |
## Geometry relationship and scope
Orientation, in-sphere, and related geometric predicates can be reduced to
determinant signs, which is why an adaptive exact sign is useful near degeneracy
\[8\]. `la-stack` supplies both point-matrix and bounded-expression determinant
primitives; callers still own problem-specific matrix assembly and semantic
classification. The crate originated to support
[`delaunay`](https://crates.io/crates/delaunay), but its matrix, factorization,
and exact-arithmetic APIs are general numerical infrastructure.
The borrowed-slice angle API also supports runtime ambient coordinate lengths
without changing the const-generic storage model. Spherical callers retain
radius validation and angle-to-arc-length conversion.
The deliberate anti-goals are general dynamically sized matrix storage,
sparse matrix storage, broad decomposition coverage, alternate floating scalar
families, and GPU-, parallel-, BLAS-, or LAPACK-scale throughput. Those problems
need different storage models and numerical policies rather than extensions to
this small fixed-dimension design.
The rectangular column selector operates on borrowed const-generic arrays;
it does not introduce a general rectangular matrix framework.
The polynomial determinant accepts a borrowed const-generic square array of
canonical sparse polynomials. Its exponent encoding, monomial order and work
budgets are explicit. Simulation of Simplicity variable assignment, base-three
encoding proofs, geometric tie-breaking and consumer budget selection remain
outside this library.
## Represented values and arithmetic model
Public matrix and vector construction rejects NaN and infinity. Because the
backing arrays are private and mutation is validated, every stored entry is
finite. Arithmetic can still overflow, underflow, or round; computed non-finite
values are reported through `LaError::NonFinite` rather than stored in a public
matrix, vector, or factorization.
The floating-point arithmetic model uses ordinary `f64` operators and deliberate
`f64::mul_add` with defined single-rounding semantics. Here, **strict** means
preserving those operation contracts and each algorithm's documented evaluation
order; it does not mean exact real arithmetic or a universal bitwise
reproducibility guarantee across all platforms and library functions.
Rust's algebraic floating-point operations are not used. The
[validation decision](validation/algebraic-floating-point.md) records why their
measured speedups do not justify changed overflow behavior or unsupported
numerical guarantees.
Every finite binary64 value is a dyadic rational with an exact representation
```text
x = (-1)^s m × 2^e,
```
for integer `m` and exponent `e` \[9-11\]. Both signed-zero bit patterns
represent rational zero. Exact methods on `Matrix` and `Vector` exploit this
representation, but their exactness begins only after binary64 construction:
they cannot recover information lost when a decimal value or an earlier
computation was rounded to `f64`. Caller-supplied `RationalMatrix` and
`RationalVector` values instead enter the exact domain directly and do not pass
through binary64.
`Matrix`, `Vector`, `IntervalMatrix`, `Lu`, and `Ldlt` use inline fixed-size storage.
Arbitrary-precision `BigInt` and `BigRational` values allocate when the `exact` feature is
used. Const-generic `D` is not itself a mathematical dimension limit; “small,
fixed dimension” is the intended performance scope. `try_with_stack_matrix!` is
separately limited to `D = 0..=MAX_STACK_MATRIX_DISPATCH_DIM` (currently 7)
because it enumerates concrete stack types.
Except for the fixed-vector reduction and determinant filters described below,
the floating-point APIs do not provide certified forward, backward, or absolute
error bounds. This includes plain `Vector::dot`, Euclidean and squared norms,
matrix norms, factorizations, and solves. Some kernels use FMA to reduce rounding
steps, but that does not make them exact.
The [validation index](validation/README.md#arithmetic-contracts-and-evidence)
connects these arithmetic contracts to independent tests and retained studies.
## Tolerances and typed errors
`Tolerance::try_new` accepts finite values greater than or equal to zero.
`DEFAULT_SINGULAR_TOL` is the absolute value `1e-12`. A tolerance is a rejection
policy for numerical pivots or diagnostics, not an error estimate and not a
condition-number threshold.
The error model keeps distinct mathematical conclusions separate:
- `LaError::Singular` distinguishes numerical rejection from exact singularity.
Numerical context retains the factorization kind, observed pivot magnitude,
and tolerance.
- `LaError::Asymmetric` reports the mirrored values that violate LDLT's exact
symmetry precondition.
- `LaError::NotPositiveSemidefinite` records a negative pivot or a zero pivot
with remaining coupling.
- `LaError::NonFinite` distinguishes invalid input locations from arithmetic
operations that overflowed during computation.
- `LaError::Unrepresentable` distinguishes required rounding from the absence
of any finite binary64 output.
`Matrix::is_symmetric` and `Matrix::first_asymmetry` are tolerance-based,
scale-aware diagnostics. They do not establish the exact mirrored equality
required by `Matrix::ldlt`.
## Linear algebra algorithms
### Certified fixed-vector reductions
`Vector::dot_with_errbound()` binds the ordinary left-to-right FMA estimate to
a certified absolute error bound for the exact-real dot product of the stored
binary64 inputs. Starting from `s₀ = 0`, its arithmetic tree is
```text
sᵢ₊₁ = fma(leftᵢ, rightᵢ, sᵢ).
```
`Vector::dot_difference_with_errbound()` targets the exact-real expression
```text
Σᵢ axisᵢ(leftᵢ - rightᵢ)
```
without rounding coordinate differences first. For each coordinate it applies
`fma(axisᵢ, leftᵢ, s)` followed by `fma(-axisᵢ, rightᵢ, s)`, giving a specified
`2D`-event tree.
Let `u = 2^-53` be binary64 unit roundoff and
`γₙ = n·u / (1 - n·u)`. When every estimate FMA result is normal or an exact
zero, standard floating-point reduction analysis gives [9-11]
```text
with `n = D` for a dot product and `n = 2D` for the affine difference. The
implementation constructs an upper bound on the magnitude sum. A first nonzero
product with a unit (`±1`) multiplier is exact; otherwise the rounded product
moves up one representable value. This conservative step avoids computing a
rounding residual for the magnitude bound. Subsequent products enter the bound
through an FMA followed by one upward step. This encloses the previous upper
bound plus the exact product; for the same incoming bound, monotonicity makes
this step no wider than separately rounding the product upward before adding
and widening. The first product's unconditional step can make the final bound
slightly wider than the former residual-based calculation. Each nonzero product
must still round to a normal value, preserving the relative-error model's
admission rule.
When the only nonzero product is normal and has a unit multiplier, the
implementation instead certifies a zero error bound: its FMA has a zero addend,
and all other FMAs add exact zero. Encountering another nonzero product clears
this special proof. This also permits an exact certificate at `f64::MAX` when
adding a conservative positive error bound would have exhausted the endpoint
range.
The division forming `γₙ` and its final multiplication are also rounded
upward. Magnitude-ordered `FastTwoSum` then selects finite outward endpoints for
`estimate ± bound` \[17\].
The certificate stores the estimate and bound; endpoints are derived when
requested. Construction first proves that both endpoints can remain finite:
the exact value `|estimate| + bound` dominates their magnitudes. If its rounded
value is below `f64::MAX`, even the next representable value is finite. If it
rounds to `f64::MAX`, a nonpositive `FastTwoSum` residual proves the exact sum
does not exceed that limit. An infinite rounded value or positive residual at
the limit makes the certificate unavailable. Endpoint accessors can therefore
remain infallible and use the same tight directed rounding as before.
Zero products leave the magnitude bound unchanged. While the proof is
available, the estimate is normal or positive zero: it starts at positive zero,
and exact cancellation of nonzero terms rounds to positive zero. Adding either
signed zero then preserves every estimate bit, so that FMA can be omitted.
After proof loss, the FMA still executes: underflow may have produced negative
zero, whose sign a subsequent zero product can change. Nonzero FMAs keep their
original order, and the returned estimate matches the specified FMA tree.
The relative-error argument is not used across gradual underflow. A nonzero
product or estimate FMA in the subnormal range, an invalid `γₙ`, or finite-range
exhaustion in proof-only arithmetic returns `Ok(None)`. This is inconclusive
evidence, not equality. A non-finite estimate FMA returns `LaError::NonFinite`
with the failing reduction index and operation. A returned
`ScalarWithErrorBound` can certify a sign or threshold comparison from its
lower/upper endpoints; overlap requires an exact fallback that reconstructs the
same expression over the original inputs. The certificate is a roundoff bound
for this arithmetic tree, not a numerical tolerance chosen by the caller.
### Determinants and certified sign filtering
`det_direct()` evaluates closed forms for `D = 0..=4`, with the empty-product
convention `det(Matrix::<0>) = 1`. `Matrix::det()` uses that path through `D = 4`
and zero-tolerance LU for `D ≥ 5`. The general `det()` result has no certified
roundoff bound, and its LU fallback can report numerical singularity even when
exact arithmetic over the stored entries would find a nonzero determinant.
For `D = 0` and `D = 1`, `det_direct_with_errbound()` returns the exact direct
determinant with a zero bound. For `D = 2..=4`, it can return a determinant
estimate and a conservative absolute bound. Let `ι(A)` denote the exact rational
lift of the stored entries and let
```text
When every rounded intermediate in the implemented determinant and `p(|A|)`
evaluation trees is normal or an exact structural zero, the implementation uses
```text
with `ε = f64::EPSILON` and project-specific coefficients
```text
c_2 = 3ε + 16ε²
c_3 = 8ε + 64ε²
c_4 = 12ε + 128ε².
```
The analysis follows Shewchuk's adaptive-filter framework \[8\], while the
FMA evaluation trees and constants are derived for this crate rather than copied
from that source. IEEE 754 and standard floating-point error analysis provide
the arithmetic model \[9-11\].
If `|determinant| > absolute_error_bound`, the sign is certified. Otherwise the
filter is inconclusive; it does not prove that the matrix is singular. The API
returns `Ok(None)` for `D ≥ 5` or when gradual underflow could invalidate the
relative-error model. A non-finite computed determinant or bound returns
`LaError::NonFinite`.
#### Derivation of the returned determinant bound
The implementation uses a rounded permanent `p_hat`, so its actual returned
bound is `B = fl(c_D × p_hat)`. The following argument includes that rounding.
Assume IEEE 754 round-to-nearest, ties-to-even and gradual underflow, with every
rounded operation finite and normal or an exact zero \[9-11\].
For this proof use `ε = 2^-52`, twice the unit roundoff `u = 2^-53`.
Then each nonzero operation satisfies `fl(t) = t(1 + δ)`, `|δ| ≤ ε`.
This allowance also covers an exact `t` just below the smallest normal that
rounds up to a normal result: its relative error is at most `u / (1 - u) < ε`.
Exact zeros require no error term. Per-operation tracking conservatively accepts
only structural zeros; the grid shortcut below also proves cancellation zeros
exact. Subnormal rounded results are excluded by the filter.
Expanding the implemented arithmetic tree expresses each determinant monomial
as its exact signed value times at most `k_D` factors of the form `(1 + δ)`.
FMA contributes one rounding factor, including when its two terms cancel.
The same count applies to the unsigned monomials of the permanent:
| 2 | One multiply, one FMA | One multiply, one addition | 2 |
| 3 | A 2×2 minor, then three outer operations | An absolute 2×2 sum, then three outer operations | 5 |
| 4 | A 3×3 cofactor, then four outer operations | An absolute 3×3 sum, then four outer operations | 9 |
Sharing the six 2×2 minors in the dense D=4 path changes reuse, not the number
of rounding factors per monomial. Skipping zero coefficients in sparse paths
can only remove terms and operations. The standard product-of-errors bound,
`|Π(1 + δ_j) - 1| ≤ γ_k`, with `γ_k = kε / (1 - kε)`, therefore gives \[11\]
```text
|d_hat - det(ι(A))| ≤ γ_k p
p_hat ≥ (1 - γ_k) p,
```
where `p = p(|A|)` is exact. No independence of the rounding errors is assumed.
Accounting also for rounding the final multiplication yields
```text
B ≥ (1 - ε) c_D p_hat ≥ (1 - ε) c_D (1 - γ_k) p.
```
Consequently, a sufficient condition for `B ≥ |d_hat - det(ι(A))|` is
```text
c_D ≥ γ_k / ((1 - ε)(1 - γ_k))
= kε / ((1 - ε)(1 - 2kε)).
```
All three stored `c_D` values are exactly representable in binary64. For
`ε = 2^-52` and `k ≤ 9`, the denominator `(1 - ε)(1 - 2kε)` is greater than
`3/4`. The required coefficient is thus less than `(4k/3)ε`, which is at most
the respective linear term `3ε`, `8ε`, or `12ε`. The positive quadratic terms
provide additional margin. This proves the returned bound, including downward
rounding of the permanent and bound, for the stated domain. If `p = 0`, all
Leibniz terms vanish and the same inequalities give zero error and bound.
The fast underflow check also has a simple grid argument. When every nonzero
input has magnitude at least `2^-16`, every entry is an integer multiple of
`2^-68`. The homogeneous degree-`j` arithmetic nodes preserve the grid
`2^(-68j)` under rounding. Through degree four, a nonzero result is therefore
at least `2^-272`; since each coefficient is a multiple of `2^-100`, a nonzero
final bound is at least `2^-372`. These are far above the smallest normal
`2^-1022`. Smaller inputs use per-operation underflow tracking instead.
Overflow remains a typed failure in either path.
The independent rational Leibniz checks in `tests/proptest_exact.rs` exercise
the error inequality for D=2–4. They support regression detection; the argument
above supplies the bound for every input satisfying the arithmetic assumptions.
### Exact arithmetic over binary64 inputs
The `exact` feature decomposes each stored entry into an integer mantissa and a
power of two. The entries are scaled to integer matrices without changing their
represented rational values \[9-10\].
For determinants, each row uses its own minimum non-zero binary exponent `eᵢ`:
`A[i, :] = 2^eᵢ A_int[i, :]`. Multilinearity gives
`det(A) = 2^(Σ eᵢ) det(A_int)`. This keeps integer entries small when rows have
widely different magnitudes. A zero row has scale zero and determinant zero;
the empty determinant is one. The positive scale preserves the determinant sign.
The exponent sum is accumulated in `i64` before checking its final `i32` range.
An addressable `D × D` binary64 matrix on Rust's 16/32/64-bit targets has
`D < 2^31`; each row exponent lies in `[-1074, 1023]`. Thus every partial sum
has magnitude below `2^42`, and cancellation cannot cause intermediate overflow.
Exact determinants use direct `BigInt` expansions for `D ≤ 4` and fraction-free
Bareiss elimination for `D ≥ 5` \[7\]. Exact solves apply Bareiss updates to an
integer augmented system, then use integer back-substitution over a common
denominator, reducing once per solved row and canonicalizing final outputs. Matrix
and right-hand-side scales start independently. Writing the selected exponents
as `s_A` and `s_b`, the integer forms satisfy `A = 2ˢᴬ · A_int` and
`b = 2ˢᵇ · b_int` \[9\]. Gaps of at most 64 bits use the lower shared scale,
while larger gaps remain separate; multiplying the integer system's solution by
the exact factor `2^(s_b − s_A)` preserves `A x = b`. First-nonzero pivoting is
sufficient for correctness in exact arithmetic, although pivot choice can still
affect computational cost.
For the upper-triangular system `U x = c` produced by Bareiss elimination,
back-substitution maintains a positive shared denominator `q` for the solved
suffix: `x[j] = n[j]/q` for `j > i`, with integer `n[j]`. Initially `q = 1`
and the suffix is empty. Each row computes
```text
s = q c[i] - Σⱼ₌ᵢ₊₁ᴰ⁻¹ U[i,j] n[j]
p/t = s/U[i,i] in lowest terms, t > 0
q ← q t
n[j] ← n[j] t for j > i
n[i] ← p.
```
Each pivot is nonzero. Before the update, `q x[i] = s/U[i,i] = p/t`, so
`x[i] = p/(q t)`; multiplying the existing suffix's numerators and denominator
by `t` preserves its values. The shared denominator remains the least common
multiple of the solved components' reduced denominators: if the old common
denominator was minimal, adjoining `p/(q t)` requires precisely the extra factor
`t`, since `gcd(p,t) = 1`. A unit `t` needs no suffix rescaling. Pivot signs are
handled by the rational reduction; no row-swap sign correction is needed for
solves. Final `BigRational::new(n[i], q)` calls produce lowest terms with positive
denominators, including `0/1`. Empty systems return the empty solution.
Forward elimination retains the first missing-pivot column for exact singularity
errors.
This takes `O(D²)` integer operations after elimination, with one rational
reduction per solved row and one per output before binary scale restoration.
For integer `A_int` and `b_int`, the adjugate identity
`det(A_int) x = adj(A_int) b_int` shows that each reduced denominator divides
`|det(A_int)|`. Their least common multiple therefore also divides that value,
while avoiding determinant factors not needed by the solution. Integer bit
lengths remain input-dependent. No rounding is introduced, including for
ill-conditioned inputs.
`det_sign_exact()` first attempts the certified binary64 filter for `D ≤ 4` and
falls back to exact integer arithmetic when the filter is inconclusive. It
returns the exact determinant sign for every finite stored matrix. `det_exact()`
and `solve_exact()` return values exact for the stored `A` and `b`, subject to
their documented scale and singularity errors.
Strict `*_exact_f64` conversions succeed only when the exact result already has
an exactly representable finite binary64 value. `RequiresRounding` means a
finite output exists only after rounding; `NotFinite` means even the rounded
result cannot be finite. The explicit rounded conversions use round-to-nearest,
ties-to-even \[9-10\]. A nonzero exact value may consequently round to zero.
### Exact arithmetic over rational inputs
`RationalMatrix<D>` and `RationalVector<D>` are a separate input domain for
coefficients assembled exactly before linear algebra begins. Their constructors
reject zero denominators and store each accepted quotient in lowest terms with
a positive denominator, so equivalent raw representations have identical
storage and denominator-clearing cost. For each matrix row `i`, the
implementation selects a positive least common multiple `sᵢ` of
the rational denominators and forms the integer row `A_int[i, :] = sᵢ A[i, :]`.
Therefore
```text
sign(det(A_int)) = sign(det(A))
det(A) = det(A_int) / ∏ᵢ sᵢ.
```
The sign path reads the `BigInt` determinant sign directly and never constructs
the rational determinant value. General solves include the corresponding right-hand
side denominator in `sᵢ`, so each augmented row is multiplied by the same
positive factor and the solution set is unchanged. The integer determinant and
solve states use the same direct-expansion/Bareiss elimination backend \[7\]
as exact operations over binary64 inputs. General solves then use integer
back-substitution over a common denominator, as derived in
[the binary64-input section](#exact-arithmetic-over-binary64-inputs).
Before clearing denominators, rational solves recognize row permutations of
upper or lower triangular matrices. For an upper triangular ordering, each
row must have a distinct first non-zero column. There are D rows and D columns,
so these columns form a permutation of `0..D`; assigning each row to its first
non-zero column gives a non-zero diagonal and zeros strictly below it. The
same argument using last non-zero columns gives a lower triangular ordering.
Only row indices are stored, and the right-hand side follows the same row map.
Back- or forward-substitution then computes
`x[k] = (b[row(k)] - Σⱼ A[row(k), j] x[j]) / A[row(k), k]`, with the sum over
already-solved columns \[11\]. The residual is accumulated as two arbitrary-precision
integers `N/Q`, initially the canonical RHS numerator and positive denominator.
For a non-zero product `a × x = P/R`, cross multiplication updates
`N ← N R - P Q` and `Q ← Q R`; equal denominators need only `N ← N - P`.
Both `Q` and `R` remain positive. Dividing by the non-zero pivot `p/q` constructs
and reduces `N q / (Q p)` once, so every stored solution component is canonical.
This avoids a GCD reduction after each product and subtraction without creating
unchecked rational values. Zero coefficients or solution components need no
multiplication. Recognition and substitution cost at most
`O(D²)` coefficient operations, with `O(D)` additional storage; bit complexity
still depends on integer component growth. Deferred reduction may increase
intermediate bit lengths relative to reducing every operation. Numerical conditioning introduces
no rounding error in this exact domain.
A zero row or repeated edge column rejects that ordering; if neither ordering
works, the unchanged Bareiss backend handles the general system and owns exact
singularity diagnostics. This includes singular triangular inputs. No Bareiss
update is skipped: even a zero elimination coefficient would still require its
pivot-ratio scaling. Empty systems retain their unique empty solution.
The rational types are const-generic and shape-safe after construction. The
`try_with_rational_matrix!` helper provides explicit runtime dispatch through
D=8 without unstable generic const expressions. Conversion to binary64 remains
a separate strict or explicitly rounded `ExactF64Conversion` operation.
`RationalVector<D>::dot(&rhs)` computes `Σ x_i y_i` with exact rational
multiplication and addition. The type guarantees matching lengths and canonical
operands; the empty sum is canonical zero. Arbitrary-precision intermediates
allocate as needed, but the operand vectors are borrowed without cloning and
there is no binary64 conversion. Canonical reduction changes representation,
never the represented rational number \[12\].
`rational_dot(&left, &right)` admits runtime-length slices after checking equal
lengths and nonzero denominators. Raw unreduced or negative-denominator values
are interpreted as mathematical quotients, and each product is canonicalized
before accumulation. Shape errors precede coefficient errors; left operand
validation precedes right operand validation. The adapter has no dispatch limit:
Delaunay's basis lengths through 8 and full objective lengths 14 and 22 in 6D
are ordinary inputs, as are larger lengths. Both APIs perform a linear number
of rational operations; bit complexity depends on numerator/denominator growth.
### Exact sparse polynomial determinants
`SparsePolynomial` represents `Σ c_e t^e`, with `e ∈ [0, u128::MAX]` and
canonical rational `c_e`. Construction rejects zero denominators before any
arithmetic, reduces quotients, combines duplicate exponents, and removes exactly
zero coefficients. Numeric increasing exponent order is fixed. The leading term
is the least nonzero exponent, as appropriate to `t → 0+`; it is not the
conventional highest-degree term. Arbitrary multivariate orderings are not
supported directly. A caller encoding them into `u128` must ensure injectivity,
order preservation and compatibility with addition throughout the calculation.
For a column subset `S` with `k = |S|`, maintain the determinant `p[S]` of the
first `k` rows restricted to the sorted columns in `S`. Initialize `p[∅] = 1`.
For every `j ∉ S`, add
```text
p[S ∪ {j}] += (-1)^(number of columns in S greater than j) * p[S] * A[k,j].
```
Appending column `j` adds exactly that many inversions to each existing
permutation. Each permutation has exactly one predecessor, so the recurrence
sums the Leibniz terms once with their correct signs \[12\]. All predecessors
have smaller numeric masks, allowing a completed partial to be consumed during
increasing-mask traversal. No division or pivot is needed, even when many
coefficients cancel or entire partial determinants vanish.
Every product and sum uses exact rational arithmetic. No term is discarded
based on its exponent or sign. Only after all contributions are complete can
the first nonzero term be read. A successful empty term map proves the zero
polynomial; `LaError::PolynomialBudgetExceeded` and
`LaError::PolynomialExponentOverflow` provide no sign or zero classification.
Overflowing intermediate exponents are rejected even if subsequent cancellation
could have removed them. There is no rounding error or numerical conditioning
restriction; combinatorial growth and integer bit complexity limit resources.
The caller supplies `PolynomialDeterminantBudget`. For positive order `D`, the
complete structural envelope `D × 2^(D-1)` is checked against
`max_subset_transitions` before allocation, even for zero matrices. Order zero
needs no transitions and returns one. The hard order-eight cap limits the
state table to 256 maps and the structural envelope to 1,024 transitions.
`max_term_products` counts each nonzero coefficient pair before exponent
addition and multiplication. Zero entries/partials cost no products; exact
cancellation never refunds completed work. Traversal order is documented in
the public API, and each call starts a fresh budget.
If `P` products are permitted, at most `P + 1` term insertions can occur,
including the initial one. Consuming states releases their storage, and inputs
are borrowed without deep cloning. Neither budget bounds integer bit lengths,
allocation bytes, elapsed time or the construction of caller-owned polynomials.
Allocator failure retains Rust's normal behavior. A term product may be costly
when coefficients have large numerators or denominators.
For a consumer whose column term maxima are `m_j`, the no-cancellation bound for
a subset is `|S|! × Π_{j∈S} m_j`. Therefore a conservative product envelope is
```text
The Delaunay order-eight lifted shape `[2, 2, 2, 2, 2, 2, 13, 1]` gives
55,631,284 products. This is a consumer policy example, not a library default
or a practical memory guarantee. A fully coincident origin fixture with
base-three coordinate perturbations has 241,920 distinct final monomials at
order eight; even this sparse input can create a large result. The release
benchmark and independent permutation tests exercise that shape through order
eight. Dense small-degree tests separately verify coefficient cancellation,
row swaps, repeated rows and non-dyadic coefficients.
### Exact-to-binary64 conversion
For strict conversion, a canonical nonzero rational must have a power-of-two
denominator. Removing powers of two from its numerator produces
`x = sign × m × 2^e`, with `m` positive and odd. Let `b` be the bit length of
`m`. Exact finite binary64 representation is possible precisely when \[9-10\]
```text
b ≤ 53, e ≥ -1074, e + b - 1 ≤ 1023.
```
These conditions bound significand precision, the least representable bit,
and the highest exponent. Subnormal output stores `m × 2^(e + 1074)` in the
fraction field; normal output aligns the significand to 53 bits and stores the
biased exponent. Zero is handled separately. Raw `BigRational` inputs with
common factors or negative denominators are normalized as needed before this
decision; a zero denominator is rejected. Canonical storage allows
`RationalVector` conversion to skip that repeated normalization.
The integer-and-exponent determinant path rounds directly from `BigInt` digits.
It retains 53 significand bits for normal results, or the bits on the `2^-1074`
grid for subnormals. A discarded guard bit triggers an increment only if some
lower discarded bit is set (the sticky condition) or the retained integer is
odd. This implements nearest-even ties without constructing a rational
denominator. A carry can promote a subnormal to normal, advance the exponent,
or overflow. At `f64::MAX + 2^970`, the tie rounds to infinity; at magnitude
`2^-1075`, the tie rounds to signed zero \[9-10\].
For an already-computed `BigRational` value, rounded conversion delegates to
`num-rational`'s `ToPrimitive::to_f64` and checks for finite output. The strict
path uses that rounded result only when needed to distinguish
`RequiresRounding` from `NotFinite`; it never returns a rounded value as an exact
conversion. Neither conversion reruns determinant or solve elimination.
### Gram matrices and geometric measures
For `M` vectors in `N` coordinates, let `V` contain the vectors as rows. Their
Gram matrix is `G = V Vᵀ`, with `G[i,j] = v_i · v_j`. In exact real arithmetic
`G` is positive semidefinite because `zᵀGz = ||Vᵀz||² ≥ 0`; it is positive
definite exactly when the vectors are linearly independent.
For `M ≤ N`, `det(G)` is the squared M-dimensional volume of the parallelotope
spanned by those vectors. When the vectors are edges from a common simplex
vertex, the simplex volume is `sqrt(det(G)) / M!` \[16\]. These identities
apply to lower-dimensional simplices embedded in a larger coordinate space.
`gram_matrix` evaluates each upper-triangle dot product once using
`Vector::dot`'s left-to-right FMA reduction and copies it to the other triangle.
The returned binary64 matrix is therefore exactly symmetric, but rounding can
destroy positive semidefiniteness or rank. Construction does not certify linear
independence or volume accuracy and provides no absolute rounding-error bound.
A dot-product overflow remains a typed `LaError::NonFinite` failure.
For nonempty `V` with full row rank, `κ₂(G) = κ₂(V)²`: forming a Gram matrix
squares the spectral condition number \[11-12\]. Nearly dependent vectors require care
when interpreting a computed determinant or passing the matrix to `ldlt`.
Its exact-symmetry and computed-pivot checks remain in force. The
[Gram references](../REFERENCES.md#gram-matrices-and-geometric-measures)
provide the volume identity and numerical background.
### Interval linear solves
For finite interval entries `[A]` and `[b]`, `IntervalMatrix<D>::solve` returns
`[x]` with the following contract: **for every** real matrix `A ∈ [A]` and
vector `b ∈ [b]`, `A` is nonsingular and its unique exact-real solution lies
componentwise in `[x]`. Entries vary independently, so correlations inherited
from expression construction may be widened. This encloses the united solution
set; it need not give its tightest bounding box. Interval linear systems and
inclusion arithmetic are described by Moore, Kearfott, and Cloud
\[[20](../REFERENCES.md#ref-20)\]. The
[interval solve references](../REFERENCES.md#interval-linear-solves) connect
that background to this implementation's enclosure proof below.
The algorithm uses Gaussian elimination with row pivoting. At column `k` it
chooses the remaining interval with greatest minimum absolute value, keeping
the first row on ties. The same fixed row permutation applies to every admitted
real system. A strictly positive separation proves that every corresponding
real pivot is nonzero. Outward division encloses each multiplier, and outward
multiply/subtract encloses each remaining coefficient and RHS after the exact
row operation, using the [scalar interval guarantees](#outward-rounded-interval-expressions).
The eliminated entry is set to exact zero: algebraically it is
`a - (a/p)p = 0`, even though an independent interval evaluation would widen it.
Induction therefore preserves inclusion of every transformed real system.
Nonzero pivots in every column prove nonsingularity; outward back substitution
then encloses each exact solution component by reverse induction. All rounding
and input uncertainty is propagated in the returned endpoints; no separate
scalar forward-error estimate is promised.
If every remaining candidate contains zero, the algorithm returns
`LaError::IntervalSolveInconclusive { pivot_col }`, including at the final
column. It never reports `LaError::Singular`: interval dependency can obscure
a nonzero pivot even for exact point inputs. Wide data and poor conditioning
increase enclosure width and the chance of an inconclusive result. Overflow
of an intermediate enclosure returns `IntervalRangeExhausted` attributed to
`IntervalSolve`; it does not imply that the true solution is unbounded.
Neither failure supplies partial solution bounds. Invalid finite ordering and
non-finite endpoints are rejected at construction, before solving.
Matrix and RHS shapes are fixed by `D`; the empty system has the empty solution.
Solves have no D=7 determinant-workspace limit, use O(D³) work and O(D²) inline
scratch, and are intended for small dimensions. Delaunay-style D=1–6 equation
assembly is covered by exact-rational regressions, but circumsphere assembly,
radius bounds, coordinate boxes, and exhaustive fallback remain consumer policy.
Callers needing a certified answer must retain their exact or robust fallback
on either inconclusive pivots or range exhaustion.
<a name="floating-point-factorizations"></a>
### LDLT without pivoting
LDLT factorization targets
```text
A = L D Lᵀ,
```
with unit lower-triangular `L` and diagonal `D` \[4-5, 11-12\]. The input
must be exactly symmetric under binary64 comparison: every mirrored pair must
satisfy `A[i][j] == A[j][i]`. Signed zeros compare equal and are accepted.
A successful `Ldlt` requires every computed diagonal pivot to be positive and
greater than the caller's tolerance. An uncoupled computed zero or a positive
pivot at or below tolerance returns `LaError::Singular`; a negative pivot or
zero pivot with remaining coupling returns `LaError::NotPositiveSemidefinite`
with a typed violation.
This is not a pivoted symmetric-indefinite factorization such as Bunch-Kaufman
\[6, 11-12\]. Pivots are computed in binary64, so a singular represented matrix
can produce a small positive pivot above a low tolerance. Successful
factorization is therefore not an exact positive-definiteness certificate for
the stored matrix, much less for ideal values before binary64 input conversion.
### LU with partial pivoting
LU factorization targets
```text
P A = L U,
```
where `P` is a row permutation, `L` is unit lower triangular, and `U` is upper
triangular. At column `k`, the implementation selects the largest-magnitude
remaining entry in the active column. It rejects the factorization when that
magnitude is less than or equal to the caller's tolerance.
The tolerance is an absolute finite, non-negative threshold. It is not divided
by a matrix norm, so rescaling a system can change whether a pivot is accepted.
Partial pivoting is a practical stability strategy, not an unconditional
accuracy guarantee: worst-case element growth and typical behavior are distinct
questions \[1-3, 11-12\]. `Lu::solve` does not estimate conditioning or refine
the result.
`Lu::det` combines permutation parity with the product of the diagonal of `U`.
Scaled product accumulation avoids some premature overflow and underflow, but
the final binary64 determinant is still rounded. A returned zero or a numerical
`LaError::Singular` is therefore not proof that the represented matrix is
exactly singular.
### Outward-rounded interval expressions
`Interval` owns the invariant `-∞ < lower ≤ upper < +∞`. Its public constructor
rejects non-finite or inverted bounds, its fields are private, and every
arithmetic operation either returns another valid enclosure or a typed range
failure. Both signed-zero inputs represent exact real zero and are canonicalized
to `+0.0`; finite subnormal endpoints remain valid.
The [interval arithmetic references](../REFERENCES.md#interval-arithmetic-division-and-reciprocals)
identify the binary64 model, error-free transformations, and interval inclusion
background. The arguments below explain how this crate combines them.
#### Scalar arithmetic and finite bounds
Point construction introduces no width. Exact-real subtraction and interval
addition use an error-free, magnitude-ordered `FastTwoSum` residual to determine
whether the rounded result is exact or which adjacent binary64 value is required
for the outward endpoint. Ordering the operands prevents intermediate overflow
whenever the rounded sum is finite \[[17](../REFERENCES.md#ref-17)\].
Multiplication uses the exact residual
`left.mul_add(right, -rounded)` of `TwoProductFMA`
\[[18](../REFERENCES.md#ref-18)\] when
`|rounded| ≥ 2^-968`. To justify this threshold, let `E` be the rounded product's
leading binary exponent. The exact product has at most 106 significant bits and
leading exponent at least `E-1`, allowing for a rounding carry. Its least bit
therefore has exponent at least `E-106 ≥ -1074`. The residual fits binary64
precision and cannot lose low bits below the subnormal range. Its sign selects
the same tight outward endpoint as exact integer comparison.
In this range a normal power-of-two operand makes the product exact by exponent
scaling, so no residual evaluation is needed.
Below that conservative threshold, multiplication decomposes each nonzero
operand into an integer significand and power of two and compares the exact
106-bit significand product with the rounded result. A zero FMA residual alone
would not prove exactness here, even for a normal rounded product.
This comparison also handles products that underflow to zero: a positive result
is enclosed by `[0, f64::from_bits(1)]`, and a negative result by the mirrored
interval. Squaring uses multiplication bounds but gives every interval spanning
zero the exact lower bound zero. Singleton squares and sums evaluate the shared
endpoint only once, preserving the same enclosure and error provenance.
These guarantees rely on IEEE-754 binary64
round-to-nearest, ties-to-even, and gradual underflow \[9-11\]. IEEE 1788
provides the broader standardized interval arithmetic model \[14\]; this crate's
deliberately smaller, undecorated surface does not claim conformance.
If the exact result lies outside `[-f64::MAX, f64::MAX]`, no interval with finite
binary64 endpoints can contain it. The operation then returns
`LaError::IntervalRangeExhausted` with the responsible interval
`ArithmeticOperation`; it never stores infinity as an interval bound.
#### Division and reciprocals
Reciprocal and division reject every zero-containing divisor with
`LaError::IntervalZeroDivisor`, even when the numerator is zero. For a zero-free
divisor, extrema of `x/y` over the input rectangle occur at endpoint pairs.
Each endpoint quotient `q` is rounded to nearest, then the sign of its exact
error is found by comparing `|x|` and `|q| |y|` as binary integers with exponents.
The product needs at most 106 significand bits, fitting `u128`. Its comparison
selects the adjacent endpoint without an underflow-prone floating residual;
the quotient's sign reverses the magnitude comparison for negative results.
Exact quotients retain point bounds. Underflow to zero widens toward the least
subnormal, and an outward infinite endpoint reports finite-range exhaustion.
Direct quotient evaluation avoids rejecting a representable ratio solely
because `1/y` is unrepresentable. Reciprocal uses the same kernel with `x=1`.
The interval operations follow the inclusion model in
\[[20](../REFERENCES.md#ref-20)\], Chapters 2 and 4. The integer comparison above
is a local implementation argument based on the binary64 representation
\[[9](../REFERENCES.md#ref-9)\]; it does not assume an error-free floating
division residual or apply `TwoProductFMA` to a quotient.
#### Determinant enclosures and sign evidence
For D≤7, `IntervalMatrix::det()` evaluates the Leibniz determinant
\[[12](../REFERENCES.md#ref-12)\] with subset dynamic programming.
The expansion identity is standard; the interval subset-DP organization is
specific to this crate. A state for each column subset stores the determinant
enclosure for the corresponding leading-row minor, requiring `max(2, 2^D)`
inline states and `D × 2^(D-1)` interval products for positive D. The const
dimension selects the workspace size without changing the expansion order.
The expansion performs no division,
so a pivot interval containing zero cannot make the algorithm unsound.
`det_sign()` classifies a strictly positive or negative enclosure accordingly,
returns `Zero` only for the singleton `[0, 0]`, and returns `Inconclusive` for
every other overlap with zero. Inconclusive evidence is not a singularity
classification.
Lifting a completed `Matrix` creates point intervals for its stored values. It
does not recover rounding from earlier subtraction, dot products, or lifted-norm
construction. Robust callers instead assemble those derived coefficients with
interval operations, use `IntervalMatrix::det_sign()` as a fast proof, and
rebuild the expression in `RationalMatrix` or another exact representation when
the filter is inconclusive or loses range.
### Rectangular column selection
For `A ∈ ℝ^(R×C)`, `column_basis(&A, tolerance)` scans columns in increasing
order. It chooses the largest absolute residual entry in the remaining rows,
breaking ties by the first current row, and accepts it only if `|pivot| > tol`.
Accepted rows are swapped into position; multipliers are bounded by one by
partial pivoting. Each trailing entry uses `fma(-factor, pivot_entry, entry)`;
previously skipped columns are no longer used. This is Gaussian elimination
with partial row pivoting \[11, 12\], extended to skip rejected columns.
The tolerance is absolute and inputs are not normalized. Callers wanting the
Delaunay policy can compute `tol = 1e-12 × max(1, max |A_ij|)` and validate it
with `Tolerance::try_new`. Uniform scaling with a correspondingly scaled
tolerance preserves the real-arithmetic decisions; binary64 rounding and
underflow can still change them. Independent row or column scaling changes
pivot magnitudes and can change the selection. There is no absolute error bound,
and this algorithm is not a rank-revealing QR or SVD. Even at zero tolerance,
neither accepted nor rejected pivots prove exact mathematical rank.
The result retains both dimensions and returns unique increasing input-column
indices. Its numerical rank is their count, at most `min(R,C)`; zero rows,
zero columns, and all-zero matrices produce rank zero. All original cells are
validated before selection. Overflow in elimination returns a typed computation
error locating the row-permuted work cell. Successful partial selection never
returns `SingularityReason::Exact`.
For a shared face with `S` vertices in ambient dimension `D`, the affine
constraint rows are `[point_coordinates, 1]`, giving shape `S × (D+1)`.
Callers can dispatch the small supported pairs of const dimensions, convert each
runtime row to an array with a checked `try_into`, and call `column_basis`.
For 6D this includes `S × 7` through `S=6`; intersection LP systems also fit
`8 × 14`. No square-matrix dispatch cap is used. A full-row-rank numerical
selection defines a square minor which the caller can extract into
`RationalMatrix<S>` and solve or test exactly. Geometry-specific affine
construction, certificates, and classification remain with the consumer.
### Scaled determinant products
Both factorizations compute their determinant from the stored diagonal factors,
with the row-permutation sign included for LU. Direct multiplication is retained
while every running product is finite and normal. Otherwise `ScaledProduct`
replays all factors as a sign, a normalized mantissa, and an integer exponent.
IEEE 754 bit decomposition gives each nonzero factor as `m × 2^e`, including
subnormals, with `1 ≤ m < 2` \[9-10\]. Multiplication of two such mantissas
rounds below 4, so at most one exact division by two restores the invariant.
The separate exponent records the corresponding scale without range loss.
The most recent factor stays pending until finalization. If the final exponent
is in the subnormal range, the two remaining operands are scaled by exact
powers of two before multiplication. This places their product directly on the
destination's subnormal grid, avoiding a normal-mantissa rounding followed by a
second rounding to that coarser grid. The final result can be signed zero;
overflow returns `LaError::NonFinite`.
This is range management for a product of already-rounded factor diagonals.
Earlier mantissa multiplications still round, and factorization error remains.
Scaling guarantees neither an exact product nor correct rounding of the exact
determinant. No certified absolute error bound is provided.
### Scaled Euclidean vector norm
`Vector::norm` computes `sqrt(Σᵢ xᵢ²)` with a left-to-right scaled
sum-of-squares recurrence. Its state represents the accumulated squared norm as
`scale² × scaled_sum`. For each nonzero `|xᵢ|`, either `|xᵢ| / scale` is
squared and accumulated, or a larger `|xᵢ|` becomes the new scale and the old
sum is rescaled. Every squared ratio is therefore at most one, which avoids raw
square overflow and scales all-subnormal vectors into a safe range \[15\].
Division, FMA, square root, and final rescaling still round in binary64. The
method is deterministic for a fixed coordinate order but does not generally
claim correct rounding or publish a certified error bound.
Near the upper range, a fixed-size integer accumulator avoids both false and
hidden overflow from the rounded recurrence. Every coordinate square is an
integer multiple of `2^-2148` and is below `2^2048`, so 4196 bits plus
`usize::BITS` carry bits suffice for every representable vector length. The
fallback sums these integer squares exactly, retaining even the smallest
subnormal square. It compares against squared binary64 rounding midpoints to
return the nearest norm, ties to even, without a floating square root. The
overflow midpoint is `f64::MAX + 2^970`; equality rounds to infinity \[9-10\].
The fallback is selected from the largest coordinate magnitude, independently
of the rounded norm. With `b = bit_length(D)`, a largest magnitude at most
`2^(1023-b)` gives the conservative L1 upper bound `D × scale < 2^1023`, so
the ordinary recurrence has ample overflow margin. Larger magnitudes use the
exact boundary calculation. Its fixed storage stays on the stack and requires
no `exact` feature.
A scalar `LaError::NonFinite` tagged with `ArithmeticOperation::VectorNorm`
therefore means the exact norm rounds to infinity. `Vector::norm_squared`
intentionally remains the direct FMA sum `Σᵢ xᵢ²` and may therefore fail even
when `norm` succeeds.
### Unsigned vector angles
The `VectorAngle` extension trait provides `left.angle(right)` on borrowed,
equal-length coordinate slices;
`Vector<D>::angle(&other)` reuses the kernel with validated finite storage.
Both compute an unsigned angle in radians in `[0, π]`, in quadratic time and
constant auxiliary space, without allocation or optional dependencies.
Every positive ambient length is accepted, including Delaunay's lengths 3–6
and lengths beyond the matrix dispatch limit. Shape mismatch is checked
before coordinates; non-finite coordinates are checked left operand first,
then right, in index order. Two empty vectors return `EmptyVector`; an all-zero
operand returns `ZeroVector`, left before right. Non-finite input metadata
retains `VectorOperandEntry { operand, index }` and `NonFiniteOrigin::Input`.
For nonzero vectors `u, v`, Lagrange's identity gives
```text
w² = Σᵢ<ⱼ (uᵢvⱼ - uⱼvᵢ)² = ‖u‖² ‖v‖² - (u·v)²,
θ = atan2(w, u·v).
```
The exterior-product norm `w` is `‖u‖ ‖v‖ sin(θ)`; the dot is the same
positive norm factor times `cos(θ)`. `atan2` cancels that factor without forming
the norms and remains well-conditioned near both endpoints and `π/2`.
The implementation reduces the minors with `hypot`; it never subtracts rounded
Gram products or squares a tiny minor. For length one the exterior sum is empty
and the sign of the dot selects zero or `π`.
Each input has its own positive maximum magnitude `m`. Extract its binary
exponent `e = floor(log2(m))` and scale its coordinates by `2^(500-e)`, using
up to two power-of-two multiplications. The largest coordinate is then in
`[2^500, 2^501)`. This preserves significands without dividing by a rounded norm
or non-power-of-two maximum. Coordinates below roughly `2^-1574` relative to
their maximum may underflow during scaling; their total angular contribution
stays far below the output range even for addressable slice lengths.
Products of scaled coordinates stay below `2^1002`; differences stay below
`2^1003`. Each minor uses two FMA product residuals \[18\] plus a magnitude-ordered
`FastTwoSum` subtraction residual \[17\]. This retains small differences between
large, nearly equal products. All dot terms and completed minors are multiplied
by `2^-500` before reduction. Dot accumulation also tracks addition residuals.
For 64-bit slice lengths the accumulated dot and exterior norm remain below
roughly `2^567`, with ample overflow headroom. Minors contributing to a
representable tiny angle stay well above underflow at the enlarged product scale.
Kahan's norm-weighted half-angle formula \[19\], used in Delaunay's spherical
implementation, was considered and is retained as a benchmark control on valid
fixtures. General direction normalization can still erase representable tiny
angles: with `n = 2^51`, the planar vectors `[n,n-1]` and `[n+1,n]` have
determinant 1 and dot `2n²`, hence angle `atan(2^-103)`. Dividing each vector by
its maximum can make the directions identical in binary64. Compensated minors
avoid that information loss. They also preserve the least positive angle in
`[1,0]`, `[1,2^-1074]` and the separation of `[2,t]`, `[2,-t]` with
`t = 2^-1074`; separately dividing the latter small coordinates by 2 loses both.
For a positive dot and ratio `r = w / (u·v) ≤ 2^-27`, the implementation
returns the quotient directly. The alternating arctangent series gives
`0 ≤ r - atan(r) ≤ r³/3 ≤ r × 2^-55`; this approximation error is below
half an ulp for normal results and negligible for subnormal results. Division
avoids platform-dependent `atan2` underflow, observed in Windows CI for
`[3,t]`, `[3,-t]` and `[1,1,0]`, `[1,1,t]` with `t = 2^-1074`.
Both true angles round to the least positive subnormal. The bound here concerns
only replacing arctangent by its argument, not the preceding dot and exterior
reductions or their rounding. Larger ratios and nonpositive dots use `atan2`.
The exterior norm is nonnegative, so the result stays in `[0, π]` without clamping
or a half-angle intermediate. Identical inputs yield exactly positive zero;
opposite inputs yield the binary64 constant `π`, including signed-zero changes.
After validation, axes where both raw coordinates are zero (including signed
zero) can be skipped: their dot contribution and every incident minor are zero.
The remaining terms retain their original order and arithmetic. With `K` axes
in the union of the two supports, this requires `K(K-1)/2` compensated
differences and `hypot` steps; scanning still takes up to quadratic time in
the ambient length. Dense inputs select a loop without shared-zero checks.
This is intended for small dimensions, including Delaunay's 3–6, while remaining
allocation-free and without an artificial dimension limit.
This is approximate arithmetic over the stored inputs. There is no certified
absolute error bound, exact parallelism classification, or correct-rounding
promise. Products, compensated differences, reductions, and transcendental
evaluation still round; accumulated error depends on ambient length. Subnormal
outputs have coarse relative spacing, unrepresentable angles may round to
zero, and deviations from `π` smaller than its output spacing are lost.
Coordinate rescaling that itself rounds or underflows changes the input
direction. Rust's `hypot` and `atan2` can also vary across platforms.
The guarantees assume ordinary binary64 arithmetic with gradual underflow.