acgetchell · acgetchell · Sep 8, 2026 · Sep 8, 2026
@@ -194,9 +194,11 @@ Favor the invariant over a convenient edit or faster implementation.
 
 ## Validation
 
-- Select validators proportionally to the changed surfaces. Compose each
-  relevant focused validator once for mixed changes; core Rust or public
-  behavior changes require final `just ci`.
+- Use `just check` during iterative review and fixes. Reserve `just ci` for
+  final validation once those iterations are complete; core Rust or public
+  behavior changes require that final comprehensive pass.
+- Select additional focused validators proportionally to the changed surfaces.
+  Compose each relevant focused validator once for mixed changes.
 - Use [Contributor validation guidance](CONTRIBUTING.md#validation-workflow)
   for the surface-to-command mapping. The [justfile](justfile) and
   `just --list` own the full command catalog.

@@ -38,7 +38,7 @@ include = [
 # All runtime deps are optional; see [features] below.
 # Must stay in sync with num-rational
 num-bigint = { version = "0.4.8", optional = true }
-num-rational = { version = "0.4.2", features = [ "num-bigint-std" ], optional = true }
+num-rational = { version = "0.4.2", optional = true }
 num-traits = { version = "0.2.19", optional = true }
 
 [dev-dependencies]

@@ -165,13 +165,13 @@ to combine operations.
 
 | Capability | Main entry points |
 |---|---|
-| Vector operations and norms | [`Vector<D>`][api-vector] |
+| Certified dot, affine-difference, and determinant estimates | [`ScalarWithErrorBound`][api-scalar-bound], [`DeterminantWithErrorBound`][api-det-bound] |
+| Exact signs, determinants, solves, and output conversion¹ | [Exact arithmetic examples][api-exact] |
 | Floating-point determinants and solves | [`Matrix<D>`][api-matrix], [`Lu<D>`][api-lu], [`Ldlt<D>`][api-ldlt] |
 | Gram matrix construction | [`gram_matrix`][api-gram] |
-| Certified dot, affine-difference, and determinant estimates | [`ScalarWithErrorBound`][api-scalar-bound], [`DeterminantWithErrorBound`][api-det-bound] |
 | Interval expressions and determinant signs | [`Interval`][api-interval], [`IntervalMatrix<D>`][api-interval-matrix] |
-| Exact signs, determinants, solves, and output conversion¹ | [Exact arithmetic examples][api-exact] |
 | Runtime selection of a const-generic matrix dimension | [Dimension dispatch examples][api-dispatch] |
+| Vector operations and norms | [`Vector<D>`][api-vector] |
 
 [`Tolerance`][api-tolerance] validates numerical rejection thresholds.
 [`LaError`][api-error] and its reason/location enums preserve structured

@@ -68,9 +68,9 @@ intermediates stay normal or are exact zeros, the standard
 `gamma_n = n·u / (1 - n·u)` model, with `u = 2^-53`, `n = D` for the dot
 product, and `n = 2D` for the affine difference, bounds the absolute forward error by
 `gamma_n Σ |a_i b_i|` \[[9], [10], [11]\]. The magnitude sum and final bound are
-rounded upward, while `TwoSum` supplies outward endpoints. Gradual underflow or
-proof-only range exhaustion makes the filter unavailable rather than turning an
-inconclusive result into equality. The affine form evaluates alternating
+rounded upward, while magnitude-ordered `FastTwoSum` supplies outward endpoints
+\[[17]\]. Gradual underflow or proof-only range exhaustion makes the filter
+unavailable rather than turning an inconclusive result into equality. The affine form evaluates alternating
 `axis_i × left_i` and `-axis_i × right_i` FMAs, so its certificate covers the
 original coordinates rather than an already-rounded difference vector.
 
@@ -131,11 +131,11 @@ and significand \[[9]\]. For nonzero `x`, it strips trailing zeros from the
 significand so `|x| = m · 2^e` with `m` odd; signed zeros use a separate zero
 component. The integer matrix is then assembled by shifting each mantissa left by
 `exp − e_min`, giving a GCD-free exact-integer starting point. Solves and D ≥ 5 determinants
-then apply Bareiss elimination; D ≤ 4 determinants use direct expansions. The test-only
-fallible wrapper `decompose_f64` verifies rejection of non-finite raw scalars, while the
-test-only `f64_to_big_rational` helper packages the same decomposition into a single
-`BigRational`. See Goldberg \[[10]\] for background on floating-point representation and
-conversion.
+then apply Bareiss elimination; D ≤ 4 determinants use direct expansions. Tests verify
+non-finite input rejection at the `Matrix` and `Vector` constructors. The test-only
+`f64_to_big_rational` helper uses `BigRational::from_f64` as an independent oracle for
+the production decomposition. See Goldberg \[[10]\] for background on floating-point
+representation and conversion.
 
 ### Gram matrices and geometric measures
 
@@ -183,7 +183,9 @@ algorithmic background.
 `Interval` uses IEEE-754 round-to-nearest binary64 operations plus adjacent
 representable values to enclose exact-real addition, subtraction,
 multiplication, and square results \[[9], [10], [11]\]. Addition and subtraction
-use an error-free `TwoSum` residual \[[8]\]; multiplication independently compares
+use a magnitude-ordered `FastTwoSum` residual \[[17]\]. This ordering prevents
+internal overflow when the rounded sum is finite, including opposite-sign inputs
+at the maximum finite magnitude. Multiplication independently compares
 the exact integer-significand product with the rounded binary64 result,
 including gradual underflow to zero. Results whose exact range cannot fit
 between finite binary64 endpoints return a typed range failure rather than
@@ -292,6 +294,12 @@ alphabetized for navigation without renumbering citations.
 16. <a name="ref-16"></a> Kock, Anders. "Square-densities, and volume forms." Notes, December 10, 2020.
     Introduction and §1.2 (Gram's formula).
     [Author's PDF](https://math.au.dk/~kock/heron4.pdf)
+17. <a name="ref-17"></a> Boldo, Sylvie, Stef Graillat, and Jean-Michel Muller.
+    "On the Robustness of the 2Sum and Fast2Sum Algorithms."
+    *ACM Transactions on Mathematical Software* 44.1 (2017), Article 4: 1–14.
+    Algorithms 1–2 and Theorems 5.1, 6.2 (magnitude ordering and overflow).
+    [DOI](https://doi.org/10.1145/3054947) ·
+    [Authors' PDF](https://perso.lip6.fr/Stef.Graillat/papers/a4-boldo.pdf)
 
 [1]: #ref-1
 [2]: #ref-2
@@ -308,3 +316,4 @@ alphabetized for navigation without renumbering citations.
 [14]: #ref-14
 [15]: #ref-15
 [16]: #ref-16
+[17]: #ref-17
@@ -2,6 +2,9 @@
 
 //! Independently validated rational inputs shared by timing and allocation probes.
 
+use core::array::from_fn;
+use core::cmp::Ordering;
+
 use la_stack::{BigInt, BigRational, DeterminantSign, RationalMatrix, RationalVector};
 
 use super::bench_utils::OrAbort;
@@ -64,8 +67,8 @@ impl<const D: usize> ValidatedRationalInput<D> {
 /// # Panics
 /// Panics if construction or the independent determinant/solve checks fail.
 pub fn rational_input<const D: usize>(kind: RationalInputKind) -> ValidatedRationalInput<D> {
-    let mut rows = std::array::from_fn(|row| {
-        std::array::from_fn(|col| {
+    let mut rows = from_fn(|row| {
+        from_fn(|col| {
             if row == col {
                 let diagonal = 2 * D + row + 1;
                 BigRational::from_integer(BigInt::from(diagonal))
@@ -89,9 +92,8 @@ pub fn rational_input<const D: usize>(kind: RationalInputKind) -> ValidatedRatio
             }
         }
     }
-    let expected_solution = std::array::from_fn(|index| {
-        BigRational::new(BigInt::from(index + 1), BigInt::from(index + 2))
-    });
+    let expected_solution =
+        from_fn(|index| BigRational::new(BigInt::from(index + 1), BigInt::from(index + 2)));
     let rhs_data = rational_matvec(&rows, &expected_solution);
     let matrix = RationalMatrix::try_from_rows(rows.clone())
         .or_abort("rational benchmark matrix construction");
@@ -121,7 +123,7 @@ fn rational_matvec<const D: usize>(
     rows: &[[BigRational; D]; D],
     vector: &[BigRational; D],
 ) -> [BigRational; D] {
-    std::array::from_fn(|row| {
+    from_fn(|row| {
         rows[row]
             .iter()
             .zip(vector.iter())
@@ -202,7 +204,7 @@ pub fn rational_solve_gaussian<const D: usize>(
         }
     }
 
-    let mut solution = std::array::from_fn(|_| zero.clone());
+    let mut solution = from_fn(|_| zero.clone());
     for row in (0..D).rev() {
         let mut value = rhs[row].clone();
         for (coefficient, component) in rows[row].iter().zip(solution.iter()).skip(row + 1) {
@@ -216,8 +218,8 @@ pub fn rational_solve_gaussian<const D: usize>(
 fn determinant_sign(value: &BigRational) -> DeterminantSign {
     let zero = BigRational::from_integer(BigInt::from(0));
     match value.cmp(&zero) {
-        std::cmp::Ordering::Less => DeterminantSign::Negative,
-        std::cmp::Ordering::Equal => DeterminantSign::Zero,
-        std::cmp::Ordering::Greater => DeterminantSign::Positive,
+        Ordering::Less => DeterminantSign::Negative,
+        Ordering::Equal => DeterminantSign::Zero,
+        Ordering::Greater => DeterminantSign::Positive,
     }
 }
@@ -2,11 +2,14 @@
 
 //! Shared helpers for the `vs_linalg` benchmark and its smoke tests.
 
+use core::array::from_fn;
+
 use faer::linalg::solvers::{Ldlt as FaerLdlt, PartialPivLu};
 use faer::perm::PermRef;
-use la_stack::{LaError, Matrix, Tolerance, Vector};
 use nalgebra::SMatrix;
 
+use la_stack::{LaError, Matrix, Tolerance, Vector};
+
 #[cfg(not(la_stack_v0_4_3_api))]
 use crate::bench_utils::OrAbort;
 
@@ -327,8 +330,6 @@ impl<const D: usize> ValidatedLuSolveInput<D> {
 pub fn validated_lu_solve_input<const D: usize>(
     scenario: LuSolveScenario,
 ) -> ValidatedLuSolveInput<D> {
-    use core::array::from_fn;
-
     assert!((2..=64).contains(&D));
     let expected = from_fn(|i| f64::from(u32::try_from(i + 1).or_abort("solution index")));
     let (rows, rhs) = match scenario {
@@ -426,7 +427,7 @@ pub fn make_vector_array<const D: usize>(offset: f64) -> [f64; D] {
 #[inline]
 #[must_use]
 pub fn make_norm_descending_array<const D: usize>() -> [f64; D] {
-    std::array::from_fn(|index| {
+    from_fn(|index| {
         let magnitude = vector_entry(D - index - 1, 0.0);
         if index % 2 == 0 {
             magnitude
@@ -440,7 +441,7 @@ pub fn make_norm_descending_array<const D: usize>() -> [f64; D] {
 #[inline]
 #[must_use]
 pub fn make_norm_repeated_scale_array<const D: usize>() -> [f64; D] {
-    std::array::from_fn(|index| if index % 2 == 0 { 3.0 } else { -3.0 })
+    from_fn(|index| if index % 2 == 0 { 3.0 } else { -3.0 })
 }
 
 /// Build a norm input with one non-zero entry and otherwise skipped zeros.
@@ -469,7 +470,7 @@ pub fn make_norm_wide_dynamic_range_array<const D: usize>() -> [f64; D] {
         -1.0e200,
     ];
 
-    std::array::from_fn(|index| VALUES[index % VALUES.len()])
+    from_fn(|index| VALUES[index % VALUES.len()])
 }
 
 /// Build the named Euclidean-norm scenario corpus in stable benchmark order.

@@ -67,7 +67,7 @@ promotion in one command.
 
 ## Benchmark Suites
 
-`la-stack` has four Criterion benchmark suites.
+`la-stack` has five Criterion benchmark suites.
 
 Newly rendered reports use one table per selected suite. Dimension and
 adversarial-input group appear in a `Case` column instead of creating a separate
@@ -92,6 +92,13 @@ suite compares row-cleared Bareiss operations with direct `BigRational` Gaussian
 operations over already-exact rational inputs across D=2-8. Use it to understand
 exact-arithmetic cost and track optimization progress.
 
+**`gram`** (`benches/gram.rs`) compares `gram_matrix` with checked hand-written
+assembly for square and embedded vector sets with coordinate dimensions 2-8.
+The orthogonal, dependent, near-dependent, and mixed-scale fixtures are checked
+against an independent integer matrix-product oracle before timing. Run it with
+`cargo bench --locked --features bench --bench gram`. This focused construction
+signal is not part of the release-to-release report schema.
+
 **`interval`** (`benches/interval.rs`) measures the default-feature,
 division-free interval determinant sign filter. Its fixtures cover a conclusive
 4×4 relative-coordinate lifted predicate, the corresponding inconclusive

@@ -167,7 +167,8 @@ implementation constructs an upper bound on the magnitude sum: exact
 integer-significand comparison determines whether each rounded product must move
 to its next representable value, and every positive accumulation is rounded
 upward. The division forming `γₙ` and its final multiplication are also rounded
-upward. `TwoSum` then selects finite outward endpoints for `estimate ± bound`.
+upward. Magnitude-ordered `FastTwoSum` then selects finite outward endpoints for
+`estimate ± bound` \[17\].
 
 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
@@ -471,10 +472,12 @@ failure. Both signed-zero inputs represent exact real zero and are canonicalized
 to `+0.0`; finite subnormal endpoints remain valid.
 
 Point construction introduces no width. Exact-real subtraction and interval
-addition use an error-free `TwoSum` residual to determine whether the rounded
-result is exact or which adjacent binary64 value is required for the outward
-endpoint [8]. Multiplication decomposes each nonzero binary64 operand into its exact
-integer significand and power of two, compares the exact 106-bit significand
+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\]. Multiplication decomposes each nonzero
+binary64 operand into its exact integer significand and power of two, compares
+the exact 106-bit significand
 product with the rounded result, and widens only in the required direction.
 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

@@ -23,6 +23,8 @@ const DET: Result<Option<f64>, LaError> = match MAT {
 
 fn main() -> Result<(), LaError> {
     let mat = MAT?;
+    // Integer cofactor expansion gives 72, without runtime factorization.
+    assert_eq!(DET?, Some(72.0));
 
     println!("4×4 matrix:");
     for row in mat.as_rows() {

@@ -2,6 +2,8 @@
 
 //! Compute the determinant of a 5×5 matrix via explicit LU factorization.
 
+use approx::assert_abs_diff_eq;
+
 use la_stack::prelude::*;
 
 fn main() -> Result<(), LaError> {
@@ -18,6 +20,7 @@ fn main() -> Result<(), LaError> {
     // Compute via explicit LU factorization.
     let lu = a.lu(DEFAULT_SINGULAR_TOL)?;
     let det = lu.det()?;
+    assert_abs_diff_eq!(det, 4.0, epsilon = 1.0e-12);
 
     println!("det = {det}");
     Ok(())

@@ -30,6 +30,11 @@ fn main() -> Result<(), LaError> {
     };
     let det_exact = m.det_exact()?;
     let det_exact_as_f64 = det_exact.try_to_f64()?;
+    assert_eq!(
+        det_exact,
+        BigRational::new((-3).into(), (1_u64 << 50).into())
+    );
+    assert_eq!(det_exact_as_f64.to_bits(), (-3.0 * perturbation).to_bits());
 
     println!("Near-singular 3×3 matrix (perturbation = 2^-50 ≈ {perturbation:.2e}):");
     for row in m.as_rows() {

@@ -26,6 +26,7 @@ fn main() -> Result<(), LaError> {
     ])?;
 
     let sign = m.det_sign_exact();
+    assert_eq!(sign, DeterminantSign::Negative);
     let det_f64 = m.det()?;
 
     println!("Near-singular 3×3 matrix (perturbation = 2^-50 ≈ {perturbation:.2e}):");

@@ -34,6 +34,15 @@ fn main() -> Result<(), LaError> {
 
     // Exact solve.
     let exact_x = a.solve_exact(b)?;
+    // Only the third component is nonzero: one third of the third column is b.
+    assert_eq!(
+        exact_x.as_array(),
+        &[
+            BigRational::from_integer(0.into()),
+            BigRational::from_integer(0.into()),
+            BigRational::new(1.into(), 3.into()),
+        ],
+    );
     println!("Near-singular 3×3 system (perturbation = 2^-50 ≈ {perturbation:.2e}):");
     for row in a.as_rows() {
         print!("  [");
@@ -62,7 +71,15 @@ fn main() -> Result<(), LaError> {
         exact_x.as_array()[1],
         exact_x.as_array()[2]
     );
-    match exact_x.try_to_f64() {
+    let strict = exact_x.try_to_f64();
+    assert_eq!(
+        strict,
+        Err(LaError::unrepresentable(
+            Some(2),
+            UnrepresentableReason::RequiresRounding,
+        )),
+    );
+    match strict {
         Ok(x) => {
             let x = x.into_array();
             println!(
@@ -73,6 +90,7 @@ fn main() -> Result<(), LaError> {
         Err(err) if err.requires_rounding() => {
             println!("exact try_to_f64(): {err}");
             let x = exact_x.to_rounded_f64()?.into_array();
+            assert_eq!(x.map(f64::to_bits), [0, 0, 0x3fd5_5555_5555_5555]);
             println!(
                 "exact to_rounded_f64(): x = [{:+.6e}, {:+.6e}, {:+.6e}]",
                 x[0], x[1], x[2]

@@ -8,6 +8,8 @@
 //!
 //! Run with: `cargo run --example ldlt_solve_3x3`
 
+use approx::assert_abs_diff_eq;
+
 use la_stack::prelude::*;
 
 fn main() -> Result<(), LaError> {
@@ -20,6 +22,11 @@ fn main() -> Result<(), LaError> {
     let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL)?;
     let x = ldlt.solve(b)?.into_array();
     let det = ldlt.det()?;
+    for (actual, expected) in x.into_iter().zip([1.0, 2.0, 3.0]) {
+        assert_abs_diff_eq!(actual, expected, epsilon = 1.0e-12);
+    }
+    // Tridiagonal determinant recurrence: d3 = 4*(4*4 - 1) - 4.
+    assert_abs_diff_eq!(det, 56.0, epsilon = 1.0e-12);
 
     println!("A (3×3 SPD tridiagonal):");
     for row in a.as_rows() {