diff --git a/AGENTS.md b/AGENTS.md
index cb0a306..69ed400 100644
--- a/AGENTS.md
+++ b/AGENTS.md
@@ -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.
diff --git a/Cargo.toml b/Cargo.toml
index 882f3c0..540f037 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -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]
diff --git a/README.md b/README.md
index f203612..df47a8d 100644
--- a/README.md
+++ b/README.md
@@ -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
diff --git a/REFERENCES.md b/REFERENCES.md
index 0f0fe86..2700214 100644
--- a/REFERENCES.md
+++ b/REFERENCES.md
@@ -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
diff --git a/benches/common/rational.rs b/benches/common/rational.rs
index e8597d0..2a98cc9 100644
--- a/benches/common/rational.rs
+++ b/benches/common/rational.rs
@@ -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,
     }
 }
diff --git a/benches/common/vs_linalg.rs b/benches/common/vs_linalg.rs
index 5856b17..bfadb5a 100644
--- a/benches/common/vs_linalg.rs
+++ b/benches/common/vs_linalg.rs
@@ -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.
diff --git a/docs/BENCHMARKING.md b/docs/BENCHMARKING.md
index 4753da0..6d745b9 100644
--- a/docs/BENCHMARKING.md
+++ b/docs/BENCHMARKING.md
@@ -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
diff --git a/docs/mathematical_basis.md b/docs/mathematical_basis.md
index 0cc41f2..c04471c 100644
--- a/docs/mathematical_basis.md
+++ b/docs/mathematical_basis.md
@@ -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
diff --git a/examples/const_det_4x4.rs b/examples/const_det_4x4.rs
index 5e192a9..6149019 100644
--- a/examples/const_det_4x4.rs
+++ b/examples/const_det_4x4.rs
@@ -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() {
diff --git a/examples/det_5x5.rs b/examples/det_5x5.rs
index 5acd06a..2051316 100644
--- a/examples/det_5x5.rs
+++ b/examples/det_5x5.rs
@@ -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(())
diff --git a/examples/exact_det_3x3.rs b/examples/exact_det_3x3.rs
index a4bfadc..365a122 100644
--- a/examples/exact_det_3x3.rs
+++ b/examples/exact_det_3x3.rs
@@ -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() {
diff --git a/examples/exact_sign_3x3.rs b/examples/exact_sign_3x3.rs
index 896f94b..e1deb45 100644
--- a/examples/exact_sign_3x3.rs
+++ b/examples/exact_sign_3x3.rs
@@ -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}):");
diff --git a/examples/exact_solve_3x3.rs b/examples/exact_solve_3x3.rs
index 41777a5..979d759 100644
--- a/examples/exact_solve_3x3.rs
+++ b/examples/exact_solve_3x3.rs
@@ -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]
diff --git a/examples/ldlt_solve_3x3.rs b/examples/ldlt_solve_3x3.rs
index abef061..c591c89 100644
--- a/examples/ldlt_solve_3x3.rs
+++ b/examples/ldlt_solve_3x3.rs
@@ -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() {
diff --git a/examples/rational_input_5x5.rs b/examples/rational_input_5x5.rs
index eb01373..a878fd3 100644
--- a/examples/rational_input_5x5.rs
+++ b/examples/rational_input_5x5.rs
@@ -24,6 +24,12 @@ fn main() -> Result<(), LaError> {
     ])?;
 
     let exact_solution = matrix.solve(&rhs)?;
+    assert_eq!(matrix.det(), epsilon);
+    assert_eq!(matrix.det_sign(), DeterminantSign::Positive);
+    assert_eq!(
+        exact_solution.as_array(),
+        &[1, -1, 2, 3, 4].map(|value| BigRational::from_integer(value.into())),
+    );
     println!("exact determinant: {}", matrix.det());
     println!("exact determinant sign: {:?}", matrix.det_sign());
     println!("exact solution: {:?}", exact_solution.as_array());
@@ -41,7 +47,15 @@ fn main() -> Result<(), LaError> {
     let f64_solve = f64_matrix
         .lu(DEFAULT_SINGULAR_TOL)
         .and_then(|lu| lu.solve(f64_rhs));
-    assert!(matches!(&f64_solve, Err(LaError::Singular { .. })));
+    assert_eq!(
+        f64_solve,
+        Err(LaError::singular_numerical(
+            1,
+            FactorizationKind::Lu,
+            0.0,
+            DEFAULT_SINGULAR_TOL.get(),
+        )),
+    );
 
     println!("1.0 + 2^-60 supplied as f64: {}", 1.0 + epsilon_f64);
     println!("solve from f64 inputs: {f64_solve:?}");
diff --git a/examples/solve_5x5.rs b/examples/solve_5x5.rs
index ccc0866..4d5675a 100644
--- a/examples/solve_5x5.rs
+++ b/examples/solve_5x5.rs
@@ -2,6 +2,8 @@
 
 //! Solve a 5×5 linear system via LU factorization (with pivoting).
 
+use approx::assert_abs_diff_eq;
+
 use la_stack::prelude::*;
 
 fn main() -> Result<(), LaError> {
@@ -20,6 +22,9 @@ fn main() -> Result<(), LaError> {
 
     let lu = a.lu(DEFAULT_SINGULAR_TOL)?;
     let x = lu.solve(b)?.into_array();
+    for (actual, expected) in x.into_iter().zip([1.0, 2.0, 3.0, 4.0, 5.0]) {
+        assert_abs_diff_eq!(actual, expected, epsilon = 1.0e-12);
+    }
 
     println!("x = {x:?}");
     Ok(())
diff --git a/src/error.rs b/src/error.rs
index 0c6182c..0d00c20 100644
--- a/src/error.rs
+++ b/src/error.rs
@@ -397,9 +397,9 @@ pub enum LaError {
     /// A matrix algorithm or runtime dispatch helper does not support a dimension.
     #[non_exhaustive]
     UnsupportedDimension {
-        /// Runtime dimension requested by the caller.
+        /// Matrix dimension requested by the caller.
         requested: usize,
-        /// Largest dimension supported by the dispatch helper.
+        /// Largest matrix dimension supported by the operation.
         max: usize,
     },
     /// A matrix index is outside the `D×D` storage domain.
diff --git a/src/exact.rs b/src/exact.rs
index 0bfe65b..6559be8 100644
--- a/src/exact.rs
+++ b/src/exact.rs
@@ -90,8 +90,9 @@
 //!
 //! Public `Matrix` / `Vector` values are finite by construction before exact
 //! methods reach the integer-scaled exact core. The decomposition helpers consume
-//! that proof without repeating stored-entry validation; a fallible raw-f64
-//! decomposition remains only to test rejection at the primitive boundary.
+//! that proof without repeating stored-entry validation. Tests verify rejection
+//! at the owning constructors and check decomposition against independent
+//! `BigRational::from_f64` conversions.
 
 use core::hint::cold_path;
 use core::mem::take;
@@ -259,28 +260,6 @@ const fn decompose_proven_finite_f64(x: f64) -> Component {
     }
 }
 
-/// Parse an arbitrary `f64` into its exact IEEE 754 components.
-///
-/// Returns [`Component::Zero`] for ±0.0, or [`Component::NonZero`] with a
-/// non-zero mantissa where the value is exactly
-/// `(-1)^is_negative × mantissa × 2^exponent` and `mantissa` is odd (trailing
-/// zeros stripped).  See `REFERENCES.md` \[9-10\].
-///
-/// # Errors
-/// Returns [`LaError::NonFinite`] if `x` is NaN or infinite.
-#[cfg(test)]
-const fn decompose_f64(x: f64) -> Result<Component, LaError> {
-    let bits = x.to_bits();
-    let biased_exp = ((bits >> 52) & 0x7FF) as i32;
-
-    if biased_exp == 0x7FF {
-        cold_path();
-        return Err(LaError::non_finite_input_scalar());
-    }
-
-    Ok(decompose_proven_finite_f64(x))
-}
-
 /// Convert a [`BigInt`] × `2^exp` pair to a reduced [`BigRational`].
 ///
 /// When `exp < 0` (denominator is `2^(-exp)`), shared factors of 2 are
@@ -1197,8 +1176,11 @@ enum BareissResult {
 ///
 /// When `rhs` is `Some`, row swaps and the inner-loop Bareiss update are
 /// mirrored on the RHS (treating it as column `D+1` of an augmented
-/// system).  On return, `a` is upper triangular and the last pivot lives
-/// in `a[D-1][D-1]`.
+/// system). On [`BareissResult::Upper`], `a` is upper triangular with non-zero
+/// diagonal entries; for D>0, the last pivot lives in `a[D-1][D-1]`.
+/// [`BareissResult::Singular`] may leave `a` and `rhs` partially eliminated.
+/// Callers own these scratch arrays and discard them on failure; no rollback
+/// is performed.
 ///
 /// First-non-zero pivoting is used: since all arithmetic is exact, any
 /// non-zero pivot is valid — no tolerance is required.
@@ -1818,7 +1800,7 @@ mod tests {
     use core::assert_matches;
     use std::array::from_fn;
 
-    use num_traits::Signed;
+    use num_traits::{FromPrimitive, Signed};
     use pastey::paste;
     use proptest::prelude::*;
 
@@ -1877,57 +1859,15 @@ mod tests {
 
     // Test helpers
 
-    /// Build an exact `BigRational` from an `f64` via IEEE 754 bit decomposition.
-    ///
-    /// Thin wrapper over [`decompose_f64`] that packs the mantissa/exponent
-    /// pair into a fully-formed `BigRational` of the form `±m · 2^e`.  The
-    /// production code paths (`exact_det_int_finite`, `bareiss_solve_finite`) instead
-    /// decompose entries into scaled `BigInt` collections, which avoids
-    /// per-entry GCD work in the elimination loops — so this helper
-    /// is not used by them and lives here to keep test assertions concise
-    /// (e.g. `assert_eq!(x.as_array()[0], f64_to_big_rational(3.0))`).
+    /// Lift binary64 inputs independently of production decomposition and scaling.
     ///
-    /// See `REFERENCES.md` \[9-10\] for the IEEE 754 standard and Goldberg's
-    /// survey of floating-point representation.
+    /// Expected solutions and residuals must not inherit a decomposition bug
+    /// from the implementation they check.
     ///
     /// # Panics
     /// Panics if `x` is NaN or infinite.
     fn f64_to_big_rational(x: f64) -> BigRational {
-        let component = decompose_f64(x).expect("test helper requires finite f64 input");
-        let Component::NonZero {
-            mantissa,
-            exponent,
-            is_negative,
-        } = component
-        else {
-            return BigRational::from_integer(BigInt::from(0));
-        };
-
-        let numer = if is_negative {
-            -BigInt::from(mantissa.get())
-        } else {
-            BigInt::from(mantissa.get())
-        };
-
-        if exponent >= 0 {
-            BigRational::new_raw(numer << exponent.cast_unsigned(), BigInt::from(1u32))
-        } else {
-            BigRational::new_raw(numer, BigInt::from(1u32) << (-exponent).cast_unsigned())
-        }
-    }
-
-    fn assert_non_finite_input_scalar<T>(result: &Result<T, LaError>) {
-        let Err(error) = result else {
-            panic!("expected a non-finite scalar-input error");
-        };
-        assert!(matches!(
-            *error,
-            LaError::NonFinite {
-                location: NonFiniteLocation::Scalar,
-                origin: NonFiniteOrigin::Input,
-                ..
-            }
-        ));
+        BigRational::from_f64(x).expect("test oracle requires finite f64 input")
     }
 
     fn assert_unrepresentable<T>(
@@ -2207,76 +2147,37 @@ mod tests {
     }
 
     // -----------------------------------------------------------------------
-    // Direct tests for internal helpers (coverage of private functions)
-    // -----------------------------------------------------------------------
-
-    #[test]
-    fn det_errbound_d0_is_zero() {
-        assert_eq!(Matrix::<0>::zero().det_errbound(), Ok(Some(0.0)));
-    }
-
-    #[test]
-    fn det_errbound_d1_is_zero() {
-        assert_eq!(
-            Matrix::<1>::try_from_rows([[42.0]]).unwrap().det_errbound(),
-            Ok(Some(0.0))
-        );
-    }
-
-    #[test]
-    fn det_errbound_d3_non_identity() {
-        // Non-identity matrix to exercise all code paths in D=3 case
-        let m = Matrix::<3>::try_from_rows([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 10.0]])
-            .unwrap();
-        let bound = m.det_errbound().unwrap().unwrap();
-        assert!(bound > 0.0);
-    }
-
-    #[test]
-    fn det_errbound_d4_non_identity() {
-        // Non-identity matrix to exercise all code paths in D=4 case
-        let m = Matrix::<4>::try_from_rows([
-            [1.0, 0.0, 0.0, 0.0],
-            [0.0, 2.0, 0.0, 0.0],
-            [0.0, 0.0, 3.0, 0.0],
-            [0.0, 0.0, 0.0, 4.0],
-        ])
-        .unwrap();
-        let bound = m.det_errbound().unwrap().unwrap();
-        assert!(bound > 0.0);
-    }
-
-    // -----------------------------------------------------------------------
-    // decompose_f64 tests
+    // Finite binary64 decomposition tests. Matrix/Vector constructor tests own
+    // rejection of non-finite inputs before they can reach this helper.
     // -----------------------------------------------------------------------
 
     #[test]
     fn decompose_f64_zero() {
-        assert_eq!(decompose_f64(0.0), Ok(Component::Zero));
-        assert_eq!(decompose_f64(-0.0), Ok(Component::Zero));
+        assert_eq!(decompose_proven_finite_f64(0.0), Component::Zero);
+        assert_eq!(decompose_proven_finite_f64(-0.0), Component::Zero);
     }
 
     #[test]
     fn decompose_f64_one() {
         assert_eq!(
-            decompose_f64(1.0),
-            Ok(Component::NonZero {
+            decompose_proven_finite_f64(1.0),
+            Component::NonZero {
                 mantissa: NonZeroU64::new(1).unwrap(),
                 exponent: 0,
                 is_negative: false,
-            })
+            }
         );
     }
 
     #[test]
     fn decompose_f64_negative() {
         assert_eq!(
-            decompose_f64(-3.5),
-            Ok(Component::NonZero {
+            decompose_proven_finite_f64(-3.5),
+            Component::NonZero {
                 mantissa: NonZeroU64::new(7).unwrap(),
                 exponent: -1,
                 is_negative: true,
-            })
+            }
         );
     }
 
@@ -2285,12 +2186,12 @@ mod tests {
         let tiny = f64::from_bits(1);
         assert!(tiny.is_subnormal());
         assert_eq!(
-            decompose_f64(tiny),
-            Ok(Component::NonZero {
+            decompose_proven_finite_f64(tiny),
+            Component::NonZero {
                 mantissa: NonZeroU64::new(1).unwrap(),
                 exponent: -1074,
                 is_negative: false,
-            })
+            }
         );
     }
 
@@ -2299,43 +2200,52 @@ mod tests {
         let value = f64::from_bits(0x000C_0000_0000_0000);
         assert!(value.is_subnormal());
         assert_eq!(
-            decompose_f64(value),
-            Ok(Component::NonZero {
+            decompose_proven_finite_f64(value),
+            Component::NonZero {
                 mantissa: NonZeroU64::new(3).unwrap(),
                 exponent: -1024,
                 is_negative: false,
-            })
+            }
         );
     }
 
     #[test]
     fn decompose_f64_power_of_two() {
         assert_eq!(
-            decompose_f64(1024.0),
-            Ok(Component::NonZero {
+            decompose_proven_finite_f64(1024.0),
+            Component::NonZero {
                 mantissa: NonZeroU64::new(1).unwrap(),
                 exponent: 10,
                 is_negative: false,
-            })
+            }
         );
     }
 
-    #[test]
-    fn decompose_f64_rejects_nan() {
-        assert_non_finite_input_scalar(&decompose_f64(f64::NAN));
-    }
-
     proptest! {
         #[test]
         fn finite_f64_round_trips_through_exact_decomposition(bits in any::<u64>()) {
             let value = f64::from_bits(bits);
             prop_assume!(value.is_finite());
 
-            if let Ok(Component::NonZero { mantissa, .. }) = decompose_f64(value) {
-                prop_assert_eq!(mantissa.get() & 1, 1);
-            }
-
             let exact = f64_to_big_rational(value);
+            let decomposed = match decompose_proven_finite_f64(value) {
+                Component::Zero => BigRational::from_integer(BigInt::from(0)),
+                Component::NonZero { mantissa, exponent, is_negative } => {
+                    prop_assert_eq!(mantissa.get() & 1, 1);
+                    prop_assert!((-1074..=1023).contains(&exponent));
+                    let numerator = if is_negative {
+                        -BigInt::from(mantissa.get())
+                    } else {
+                        BigInt::from(mantissa.get())
+                    };
+                    if exponent >= 0 {
+                        BigRational::from_integer(numerator << exponent.cast_unsigned())
+                    } else {
+                        BigRational::new(numerator, BigInt::from(1) << exponent.unsigned_abs())
+                    }
+                }
+            };
+            prop_assert_eq!(&decomposed, &exact);
             let reconstructed = exact_rational_to_finite_f64(&exact, None);
 
             prop_assert_eq!(reconstructed, Ok(value));
@@ -2373,17 +2283,7 @@ mod tests {
             let exact = big_int_exp_to_big_rational(value, exp);
             let oracle = exact_rational_to_rounded_f64(&exact, None);
 
-            match (direct, oracle) {
-                (Ok(actual), Ok(expected)) => {
-                    prop_assert_eq!(actual.to_bits(), expected.to_bits());
-                }
-                (Err(actual), Err(expected)) => {
-                    prop_assert_eq!(actual, expected);
-                }
-                (actual, expected) => {
-                    prop_assert_eq!(actual, expected);
-                }
-            }
+            prop_assert_eq!(direct.map(f64::to_bits), oracle.map(f64::to_bits));
         }
     }
 
@@ -3732,7 +3632,7 @@ mod tests {
     #[test]
     fn f64_to_big_rational_round_trip() {
         // -0.0 is excluded: it maps to BigRational(0) which round-trips
-        // to +0.0 (correct; tested separately in f64_to_big_rational_negative_zero).
+        // to +0.0 (covered by the negative-zero case in f64_to_big_rational_scalar_cases).
         let values = [
             0.0,
             1.0,
@@ -3758,11 +3658,4 @@ mod tests {
             );
         }
     }
-
-    #[test]
-    fn decompose_f64_rejects_non_finite_inputs() {
-        for value in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
-            assert_non_finite_input_scalar(&decompose_f64(value));
-        }
-    }
 }
diff --git a/src/interval.rs b/src/interval.rs
index d266e5b..19545b0 100644
--- a/src/interval.rs
+++ b/src/interval.rs
@@ -2,7 +2,7 @@
 
 //! Outward-rounded intervals and fixed-size interval determinant signs.
 //!
-//! See `REFERENCES.md` \[8\] for `TwoSum`, \[9-11\] for the binary64 arithmetic
+//! See `REFERENCES.md` \[17\] for `FastTwoSum`, \[9-11\] for the binary64 arithmetic
 //! model, and \[12\] for the Leibniz determinant identity. The column-subset
 //! evaluation is specialized to this crate's small dimensions. Reference
 //! \[14\] describes the broader interval standard; this module does not claim
@@ -34,7 +34,7 @@ pub const MAX_INTERVAL_MATRIX_DIM: usize = 7;
 ///
 /// # Examples
 /// ```
-/// use la_stack::{Interval, LaError};
+/// use la_stack::prelude::*;
 ///
 /// # fn main() -> Result<(), LaError> {
 /// let difference = Interval::try_from_subtraction(1.0, 0.1)?;
@@ -78,7 +78,7 @@ pub enum IntervalDeterminantSign {
 ///
 /// # Examples
 /// ```
-/// use la_stack::{IntervalDeterminantSign, IntervalMatrix, LaError};
+/// use la_stack::prelude::*;
 ///
 /// # fn main() -> Result<(), LaError> {
 /// let matrix = IntervalMatrix::<3>::try_from_point_rows([
@@ -103,7 +103,7 @@ const fn canonical_zero(value: f64) -> f64 {
 }
 
 /// Turn a finite rounded sum into the tight adjacent-float enclosure implied by
-/// its exact `TwoSum` residual.
+/// its exact `FastTwoSum` residual.
 #[inline]
 const fn rounded_add_bounds(
     left: f64,
@@ -182,6 +182,23 @@ impl Interval {
     ///
     /// Signed zero endpoints are canonicalized to `+0.0`.
     ///
+    /// # Examples
+    /// ```
+    /// use core::assert_matches;
+    /// use la_stack::prelude::*;
+    ///
+    /// # fn main() -> Result<(), LaError> {
+    /// let range = Interval::try_new(-2.0, 3.0)?;
+    /// assert!(range.contains(1.0));
+    /// assert!(!range.contains(4.0));
+    /// assert_matches!(
+    ///     Interval::try_new(3.0, -2.0),
+    ///     Err(LaError::InvertedInterval { lower: 3.0, upper: -2.0, .. })
+    /// );
+    /// # Ok(())
+    /// # }
+    /// ```
+    ///
     /// # Errors
     /// Returns [`LaError::NonFinite`] when either endpoint is NaN or infinity.
     /// Returns [`LaError::InvertedInterval`] when `lower > upper`.
@@ -205,6 +222,22 @@ impl Interval {
 
     /// Construct a point interval from a finite binary64 value.
     ///
+    /// This preserves the supplied value, including any earlier rounding.
+    /// Use [`try_from_subtraction`](Self::try_from_subtraction) to enclose a
+    /// subtraction before its rounding uncertainty is lost.
+    ///
+    /// # Examples
+    /// ```
+    /// use la_stack::prelude::*;
+    ///
+    /// # fn main() -> Result<(), LaError> {
+    /// let half = Interval::point(0.5)?;
+    /// assert_eq!((half.lower(), half.upper()), (0.5, 0.5));
+    /// assert_eq!(half.try_add(&half)?, Interval::ONE);
+    /// # Ok(())
+    /// # }
+    /// ```
+    ///
     /// # Errors
     /// Returns [`LaError::NonFinite`] when `value` is NaN or infinity.
     #[inline]
@@ -221,6 +254,23 @@ impl Interval {
     /// Unlike subtracting first and then calling [`point`](Self::point), this
     /// method preserves the rounding uncertainty introduced by the subtraction.
     ///
+    /// # Examples
+    /// ```
+    /// use la_stack::prelude::*;
+    ///
+    /// # fn main() -> Result<(), LaError> {
+    /// // The exact difference 1 - 2^-54 lies between adjacent binary64 values.
+    /// let difference = Interval::try_from_subtraction(1.0, f64::EPSILON / 4.0)?;
+    /// assert_eq!(difference.lower(), 1.0_f64.next_down());
+    /// assert_eq!(difference.upper(), 1.0);
+    ///
+    /// // Subtracting first loses that uncertainty and produces a point at 1.
+    /// let rounded = Interval::point(1.0 - f64::EPSILON / 4.0)?;
+    /// assert_eq!(rounded, Interval::ONE);
+    /// # Ok(())
+    /// # }
+    /// ```
+    ///
     /// # Errors
     /// Returns [`LaError::NonFinite`] for a non-finite input, preserving whether
     /// it was the left or right operand. Returns
@@ -267,6 +317,18 @@ impl Interval {
 
     /// Add two intervals with outward rounding.
     ///
+    /// # Examples
+    /// ```
+    /// use la_stack::prelude::*;
+    ///
+    /// # fn main() -> Result<(), LaError> {
+    /// let left = Interval::try_new(1.0, 2.0)?;
+    /// let right = Interval::try_new(0.5, 1.0)?;
+    /// assert_eq!(left.try_add(&right)?, Interval::try_new(1.5, 3.0)?);
+    /// # Ok(())
+    /// # }
+    /// ```
+    ///
     /// # Errors
     /// Returns [`LaError::IntervalRangeExhausted`] when the exact result range
     /// has no finite binary64 enclosure.
@@ -277,6 +339,21 @@ impl Interval {
 
     /// Multiply two intervals with outward rounding.
     ///
+    /// For a square of the same represented value, prefer
+    /// [`try_square`](Self::try_square), which can give a tighter enclosure.
+    ///
+    /// # Examples
+    /// ```
+    /// use la_stack::prelude::*;
+    ///
+    /// # fn main() -> Result<(), LaError> {
+    /// let left = Interval::try_new(-2.0, 3.0)?;
+    /// let right = Interval::try_new(-4.0, -1.0)?;
+    /// assert_eq!(left.try_mul(&right)?, Interval::try_new(-12.0, 8.0)?);
+    /// # Ok(())
+    /// # }
+    /// ```
+    ///
     /// # Errors
     /// Returns [`LaError::IntervalRangeExhausted`] when the exact result range
     /// has no finite binary64 enclosure.
@@ -286,6 +363,17 @@ impl Interval {
     }
 
     /// Negate an interval exactly by swapping and negating its endpoints.
+    ///
+    /// # Examples
+    /// ```
+    /// use la_stack::prelude::*;
+    ///
+    /// # fn main() -> Result<(), LaError> {
+    /// let range = Interval::try_new(-2.0, 3.0)?;
+    /// assert_eq!(range.negate(), Interval::try_new(-3.0, 2.0)?);
+    /// # Ok(())
+    /// # }
+    /// ```
     #[inline]
     pub const fn negate(&self) -> Self {
         Self::new_unchecked(-self.upper, -self.lower)
@@ -296,6 +384,19 @@ impl Interval {
     /// An interval spanning zero has exact lower bound zero. The upper bound is
     /// the outward-rounded square of the endpoint with greatest magnitude.
     ///
+    /// # Examples
+    /// ```
+    /// use la_stack::prelude::*;
+    ///
+    /// # fn main() -> Result<(), LaError> {
+    /// let range = Interval::try_new(-2.0, 3.0)?;
+    /// assert_eq!(range.try_square()?, Interval::try_new(0.0, 9.0)?);
+    /// // Multiplication treats its two operands independently and is wider.
+    /// assert_eq!(range.try_mul(&range)?, Interval::try_new(-6.0, 9.0)?);
+    /// # Ok(())
+    /// # }
+    /// ```
+    ///
     /// # Errors
     /// Returns [`LaError::IntervalRangeExhausted`] when the exact square range
     /// has no finite binary64 enclosure.
@@ -529,6 +630,8 @@ impl Default for Interval {
 
 impl<const D: usize> IntervalMatrix<D> {
     /// Construct an interval matrix from already-validated interval rows.
+    ///
+    /// See [`det`](Self::det) for an example with non-point entries.
     #[inline]
     pub const fn from_rows(rows: [[Interval; D]; D]) -> Self {
         Self { rows }
@@ -538,6 +641,7 @@ impl<const D: usize> IntervalMatrix<D> {
     ///
     /// This preserves the stored binary64 values exactly; it does not recover
     /// uncertainty from arithmetic performed before this call.
+    /// See [`IntervalMatrix`] for a determinant-sign example using this constructor.
     ///
     /// # Errors
     /// Returns [`LaError::NonFinite`] with matrix coordinates for the first NaN
@@ -566,6 +670,18 @@ impl<const D: usize> IntervalMatrix<D> {
     /// Earlier rounded expression construction is not enclosed; use interval
     /// operations while constructing derived coefficients when that uncertainty
     /// belongs in the proof.
+    ///
+    /// # Examples
+    /// ```
+    /// use la_stack::prelude::*;
+    ///
+    /// # fn main() -> Result<(), LaError> {
+    /// let matrix = Matrix::<2>::try_from_rows([[2.0, 0.0], [0.0, 3.0]])?;
+    /// let intervals = IntervalMatrix::from_matrix(&matrix);
+    /// assert_eq!(intervals.det()?, Interval::point(6.0)?);
+    /// # Ok(())
+    /// # }
+    /// ```
     #[inline]
     pub const fn from_matrix(matrix: &Matrix<D>) -> Self {
         let matrix_rows = matrix.as_rows();
@@ -626,6 +742,8 @@ impl<const D: usize> IntervalMatrix<D> {
 
     /// Get an interval entry while preserving index context on failure.
     ///
+    /// See [`set`](Self::set) for an example of mutation and checked access.
+    ///
     /// # Errors
     /// Returns [`LaError::IndexOutOfBounds`] when either index is not `< D`.
     #[inline]
@@ -640,7 +758,28 @@ impl<const D: usize> IntervalMatrix<D> {
     /// Set an interval entry with bounds checking.
     ///
     /// Validation is unnecessary for the value because [`Interval`] already
-    /// carries the finite ordered-bound proof.
+    /// carries the finite ordered-bound proof. An invalid index leaves the
+    /// matrix unchanged.
+    ///
+    /// # Examples
+    /// ```
+    /// use core::assert_matches;
+    /// use la_stack::prelude::*;
+    ///
+    /// # fn main() -> Result<(), LaError> {
+    /// let mut matrix = IntervalMatrix::<2>::identity();
+    /// let range = Interval::try_new(2.0, 3.0)?;
+    /// matrix.set(0, 0, range)?;
+    /// assert_eq!(matrix.try_get(0, 0)?, range);
+    /// let before = matrix;
+    /// assert_matches!(
+    ///     matrix.set(2, 0, Interval::ZERO),
+    ///     Err(LaError::IndexOutOfBounds { row: 2, col: 0, dim: 2, .. })
+    /// );
+    /// assert_eq!(matrix, before);
+    /// # Ok(())
+    /// # }
+    /// ```
     ///
     /// # Errors
     /// Returns [`LaError::IndexOutOfBounds`] when either index is not `< D`.
@@ -664,6 +803,23 @@ impl<const D: usize> IntervalMatrix<D> {
     ///
     /// The D=0 determinant follows the empty-product convention and is `[1, 1]`.
     ///
+    /// Use [`det_sign`](Self::det_sign) when only sign evidence is needed.
+    ///
+    /// # Examples
+    /// ```
+    /// use la_stack::prelude::*;
+    ///
+    /// # fn main() -> Result<(), LaError> {
+    /// // Every represented diagonal matrix has a determinant in [8, 15].
+    /// let matrix = IntervalMatrix::<2>::from_rows([
+    ///     [Interval::try_new(2.0, 3.0)?, Interval::ZERO],
+    ///     [Interval::ZERO, Interval::try_new(4.0, 5.0)?],
+    /// ]);
+    /// assert_eq!(matrix.det()?, Interval::try_new(8.0, 15.0)?);
+    /// # Ok(())
+    /// # }
+    /// ```
+    ///
     /// # Errors
     /// Returns [`LaError::UnsupportedDimension`] for D>7. Returns
     /// [`LaError::IntervalRangeExhausted`] with interval-determinant provenance
@@ -718,6 +874,25 @@ impl<const D: usize> IntervalMatrix<D> {
     /// singleton interval `[0, 0]` proves `Zero`; every other overlap with zero
     /// is [`IntervalDeterminantSign::Inconclusive`].
     ///
+    /// # Examples
+    /// ```
+    /// use la_stack::prelude::*;
+    ///
+    /// # fn main() -> Result<(), LaError> {
+    /// let mut matrix = IntervalMatrix::<2>::identity();
+    /// assert_eq!(matrix.det_sign()?, IntervalDeterminantSign::Positive);
+    ///
+    /// matrix.set(0, 0, Interval::try_new(-1.0, 1.0)?)?;
+    /// // This range includes nonsingular matrices of both signs; a caller
+    /// // needs tighter or exact input before it can decide singularity.
+    /// assert_eq!(matrix.det_sign()?, IntervalDeterminantSign::Inconclusive);
+    ///
+    /// matrix.set(0, 0, Interval::ZERO)?;
+    /// assert_eq!(matrix.det_sign()?, IntervalDeterminantSign::Zero);
+    /// # Ok(())
+    /// # }
+    /// ```
+    ///
     /// # Errors
     /// Propagates the dimension and arithmetic range failures from
     /// [`det`](Self::det).
@@ -864,13 +1039,18 @@ mod tests {
     #[test]
     fn inexact_operations_expand_only_in_the_required_direction() -> Result<(), LaError> {
         let subtraction = Interval::try_from_subtraction(1.0, 0.1)?;
-        let rounded_subtraction = 1.0 - 0.1;
-        assert!(subtraction.contains(rounded_subtraction));
-        assert!(subtraction.lower() < subtraction.upper());
+        let rounded_subtraction = 1.0_f64 - 0.1;
+        assert_eq!(
+            subtraction,
+            Interval::try_new(rounded_subtraction.next_down(), rounded_subtraction)?
+        );
 
         let product = Interval::point(0.1)?.try_mul(&Interval::point(0.2)?)?;
-        assert!(product.contains(0.1 * 0.2));
-        assert!(product.lower() < product.upper());
+        let rounded_product = 0.1_f64 * 0.2;
+        assert_eq!(
+            product,
+            Interval::try_new(rounded_product.next_down(), rounded_product)?
+        );
 
         let below_one = 1.0 - f64::EPSILON;
         let above_one = 1.0 + f64::EPSILON;
diff --git a/src/lib.rs b/src/lib.rs
index 4d40b30..4eefc0e 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -1012,17 +1012,14 @@ const EPS: f64 = f64::EPSILON; // 2^-52
 ///
 /// # fn main() -> Result<(), LaError> {
 /// let m = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
-/// let Some(det) = m.det_direct()? else {
-///     return Ok(());
-/// };
-/// assert_eq!(det, -2.0);
+/// let det = m.det_direct()?;
+/// assert_eq!(det, Some(-2.0));
 /// // Compute the bound from the raw constant for illustration; most
 /// // callers would match on `m.det_errbound()?` instead.
 /// let p = (1.0_f64 * 4.0).abs() + (2.0_f64 * 3.0).abs();
 /// let bound = ERR_COEFF_2 * p;
-/// if det.abs() > bound {
-///     // The f64 sign is provably correct without exact arithmetic.
-/// }
+/// // The f64 sign is provably correct without exact arithmetic.
+/// assert_eq!(det.map(|value| value.abs() > bound), Some(true));
 /// # Ok(())
 /// # }
 /// ```
@@ -1118,6 +1115,8 @@ pub use vector::{ScalarWithErrorBound, Vector};
 /// The macro creates a zero matrix with type `Matrix<N>` for the selected
 /// runtime dimension `N`, then evaluates the supplied closure body.  Supported
 /// runtime dimensions run from `0` through [`MAX_STACK_MATRIX_DISPATCH_DIM`].
+/// The body may mutate or consume captured values. It is not evaluated for
+/// unsupported dimensions.
 /// Unsupported dimensions return
 /// `Err(LaError::UnsupportedDimension { requested, max })` converted with
 /// `From<LaError>`, so downstream crates can use their own public error type.
@@ -1185,11 +1184,11 @@ macro_rules! try_with_stack_matrix {
         }
     }};
     (@arm $d:literal, $matrix:ident, $ret:ty, $body:block) => {{
-        let __la_stack_body = |$matrix: $crate::Matrix<$d>| -> $ret { $body };
+        let mut __la_stack_body = |$matrix: $crate::Matrix<$d>| -> $ret { $body };
         __la_stack_body($crate::Matrix::<$d>::zero())
     }};
     (@arm_mut $d:literal, $matrix:ident, $ret:ty, $body:block) => {{
-        let __la_stack_body = |mut $matrix: $crate::Matrix<$d>| -> $ret { $body };
+        let mut __la_stack_body = |mut $matrix: $crate::Matrix<$d>| -> $ret { $body };
         __la_stack_body($crate::Matrix::<$d>::zero())
     }};
 }
@@ -1200,6 +1199,8 @@ macro_rules! try_with_stack_matrix {
 /// dimension, then evaluates the closure body. Supported dimensions run from
 /// `0` through [`MAX_INTERVAL_MATRIX_DIM`]. Unsupported dimensions return
 /// [`LaError::UnsupportedDimension`] converted through `From<LaError>`.
+/// The body may mutate or consume captured values. It is not evaluated for
+/// unsupported dimensions.
 ///
 /// # Errors
 /// Returns [`LaError::UnsupportedDimension`] (converted through
@@ -1267,11 +1268,11 @@ macro_rules! try_with_interval_matrix {
         }
     }};
     (@arm $d:literal, $matrix:ident, $ret:ty, $body:block) => {{
-        let __la_stack_body = |$matrix: $crate::IntervalMatrix<$d>| -> $ret { $body };
+        let mut __la_stack_body = |$matrix: $crate::IntervalMatrix<$d>| -> $ret { $body };
         __la_stack_body($crate::IntervalMatrix::<$d>::zero())
     }};
     (@arm_mut $d:literal, $matrix:ident, $ret:ty, $body:block) => {{
-        let __la_stack_body = |mut $matrix: $crate::IntervalMatrix<$d>| -> $ret { $body };
+        let mut __la_stack_body = |mut $matrix: $crate::IntervalMatrix<$d>| -> $ret { $body };
         __la_stack_body($crate::IntervalMatrix::<$d>::zero())
     }};
 }
@@ -1283,6 +1284,8 @@ macro_rules! try_with_interval_matrix {
 /// supported on stable Rust. The closure may fill the matrix through
 /// [`RationalMatrix::set`] or replace it with a value built by
 /// [`RationalMatrix::try_from_fn`].
+/// The body may mutate or consume captured values. It is not evaluated for
+/// unsupported dimensions.
 ///
 /// # Errors
 /// Returns [`LaError::UnsupportedDimension`] (converted through
@@ -1354,11 +1357,11 @@ macro_rules! try_with_rational_matrix {
         }
     }};
     (@arm $d:literal, $matrix:ident, $ret:ty, $body:block) => {{
-        let __la_stack_body = |$matrix: $crate::RationalMatrix<$d>| -> $ret { $body };
+        let mut __la_stack_body = |$matrix: $crate::RationalMatrix<$d>| -> $ret { $body };
         __la_stack_body($crate::RationalMatrix::<$d>::zero())
     }};
     (@arm_mut $d:literal, $matrix:ident, $ret:ty, $body:block) => {{
-        let __la_stack_body = |mut $matrix: $crate::RationalMatrix<$d>| -> $ret { $body };
+        let mut __la_stack_body = |mut $matrix: $crate::RationalMatrix<$d>| -> $ret { $body };
         __la_stack_body($crate::RationalMatrix::<$d>::zero())
     }};
 }
diff --git a/src/lu.rs b/src/lu.rs
index bb7b35d..b504661 100644
--- a/src/lu.rs
+++ b/src/lu.rs
@@ -419,6 +419,7 @@ impl<const D: usize> Lu<D> {
 
 #[cfg(test)]
 mod tests {
+    use core::array::from_fn;
     use core::hint::black_box;
 
     use approx::assert_abs_diff_eq;
@@ -536,24 +537,15 @@ mod tests {
     #[test]
     fn row_permutation_keeps_mapping_and_parity_synchronized() {
         let mut permutation = RowPermutation::<4>::identity();
-        assert_eq!(
-            core::array::from_fn(|row| permutation.source_row(row)),
-            [0, 1, 2, 3]
-        );
+        assert_eq!(from_fn(|row| permutation.source_row(row)), [0, 1, 2, 3]);
         assert!(!permutation.is_odd());
 
         permutation.swap(0, 3);
-        assert_eq!(
-            core::array::from_fn(|row| permutation.source_row(row)),
-            [3, 1, 2, 0]
-        );
+        assert_eq!(from_fn(|row| permutation.source_row(row)), [3, 1, 2, 0]);
         assert!(permutation.is_odd());
 
         permutation.swap(1, 2);
-        assert_eq!(
-            core::array::from_fn(|row| permutation.source_row(row)),
-            [3, 2, 1, 0]
-        );
+        assert_eq!(from_fn(|row| permutation.source_row(row)), [3, 2, 1, 0]);
         assert!(!permutation.is_odd());
     }
 
diff --git a/src/matrix.rs b/src/matrix.rs
index 16abf1c..40cd00c 100644
--- a/src/matrix.rs
+++ b/src/matrix.rs
@@ -1175,10 +1175,12 @@ impl<const D: usize> Matrix<D> {
     ///
     /// # fn main() -> Result<(), LaError> {
     /// let matrix = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
-    /// if let Some(estimate) = matrix.det_direct_with_errbound()? {
-    ///     assert_eq!(estimate.determinant(), -2.0);
-    ///     assert!(estimate.absolute_error_bound() >= 0.0);
-    /// }
+    /// let estimate = matrix.det_direct_with_errbound()?;
+    /// assert_eq!(estimate.map(|value| value.determinant()), Some(-2.0));
+    /// assert_eq!(
+    ///     estimate.map(|value| (0.0..1.0e-12).contains(&value.absolute_error_bound())),
+    ///     Some(true),
+    /// );
     /// # Ok(())
     /// # }
     /// ```
@@ -1240,9 +1242,8 @@ impl<const D: usize> Matrix<D> {
     ///     [4.0, 5.0, 6.0],
     ///     [7.0, 8.0, 9.0],
     /// ])?;
-    /// if let Some(bound) = m.det_errbound()? {
-    ///     assert!(bound >= 0.0);
-    /// }
+    /// let bound = m.det_errbound()?;
+    /// assert_eq!(bound.map(|value| (0.0..1.0e-12).contains(&value)), Some(true));
     /// # Ok(())
     /// # }
     /// ```
@@ -1288,9 +1289,14 @@ impl<const D: usize> Matrix<D> {
     /// ```
     ///
     /// # Errors
-    /// Returns [`LaError::NonFinite`] when the bound computation overflows to
-    /// NaN or infinity. Underflow-sensitive finite computations return
-    /// `Ok(None)` instead because they are valid inputs for an exact fallback.
+    /// Propagates [`LaError::NonFinite`] from
+    /// [`det_direct_with_errbound`](Self::det_direct_with_errbound) when either
+    /// the determinant or bound computation produces NaN or infinity. The error
+    /// retains [`ArithmeticOperation::Determinant`] or
+    /// [`ArithmeticOperation::DeterminantErrorBound`] as its computation origin.
+    /// A non-finite determinant remains an error even if underflow prevents
+    /// computing its bound. Underflow-sensitive finite computations return
+    /// `Ok(None)` because they remain valid inputs for an exact fallback.
     #[inline]
     pub const fn det_errbound(&self) -> Result<Option<f64>, LaError> {
         match self.det_direct_with_errbound() {
@@ -2400,6 +2406,40 @@ mod tests {
 
     // === det_errbound tests (no `exact` feature required) ===
 
+    #[test]
+    fn det_errbound_d0_is_zero() {
+        assert_eq!(Matrix::<0>::zero().det_errbound(), Ok(Some(0.0)));
+    }
+
+    #[test]
+    fn det_errbound_d1_is_zero() {
+        assert_eq!(
+            Matrix::<1>::try_from_rows([[42.0]]).unwrap().det_errbound(),
+            Ok(Some(0.0))
+        );
+    }
+
+    #[test]
+    fn det_errbound_d3_non_identity() {
+        let m = Matrix::<3>::try_from_rows([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 10.0]])
+            .unwrap();
+        let bound = m.det_errbound().unwrap().unwrap();
+        assert!(bound > 0.0);
+    }
+
+    #[test]
+    fn det_errbound_d4_non_identity() {
+        let m = Matrix::<4>::try_from_rows([
+            [1.0, 0.0, 0.0, 0.0],
+            [0.0, 2.0, 0.0, 0.0],
+            [0.0, 0.0, 3.0, 0.0],
+            [0.0, 0.0, 0.0, 4.0],
+        ])
+        .unwrap();
+        let bound = m.det_errbound().unwrap().unwrap();
+        assert!(bound > 0.0);
+    }
+
     #[test]
     fn det_errbound_matches_documented_coefficient_scale() {
         let m2 = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]]).unwrap();
diff --git a/src/rational.rs b/src/rational.rs
index f94c47d..95fbd37 100644
--- a/src/rational.rs
+++ b/src/rational.rs
@@ -17,6 +17,7 @@ use std::array::from_fn;
 
 use num_bigint::{BigInt, Sign};
 use num_rational::BigRational;
+use num_traits::One;
 
 use crate::exact::{det_big_int, solve_big_int};
 use crate::{DeterminantSign, LaError};
@@ -163,6 +164,35 @@ impl<const D: usize> RationalMatrix<D> {
     /// Replace one exact entry while preserving the canonical non-zero-
     /// denominator invariant.
     ///
+    /// As with [`try_from_rows`](Self::try_from_rows), non-reduced values and
+    /// negative denominators are accepted and canonicalized. Rejected indices
+    /// or denominators leave the matrix unchanged.
+    ///
+    /// # Examples
+    /// ```
+    /// use core::assert_matches;
+    /// use la_stack::prelude::*;
+    ///
+    /// # fn main() -> Result<(), LaError> {
+    /// let mut matrix = RationalMatrix::<2>::zero();
+    /// matrix.set(0, 1, BigRational::new_raw((-2).into(), (-4).into()))?;
+    /// let half = BigRational::new(1.into(), 2.into());
+    /// assert_eq!(matrix.get(0, 1), Some(&half));
+    ///
+    /// let before = matrix.clone();
+    /// assert_matches!(
+    ///     matrix.set(0, 1, BigRational::new_raw(1.into(), 0.into())),
+    ///     Err(LaError::NonFinite {
+    ///         location: NonFiniteLocation::MatrixCell { row: 0, col: 1, .. },
+    ///         origin: NonFiniteOrigin::Input,
+    ///         ..
+    ///     })
+    /// );
+    /// assert_eq!(matrix, before);
+    /// # Ok(())
+    /// # }
+    /// ```
+    ///
     /// # Errors
     /// Returns [`LaError::IndexOutOfBounds`] when `(row, col)` lies outside the
     /// matrix, or [`LaError::NonFinite`] when `value` has a raw zero
@@ -183,6 +213,22 @@ impl<const D: usize> RationalMatrix<D> {
     /// This path clears denominators and reads the sign of the resulting
     /// integer determinant. It does not construct a rational determinant.
     /// For D=0, the empty-product determinant has positive sign.
+    /// Use [`det`](Self::det) when the determinant value is also needed.
+    ///
+    /// # Examples
+    /// ```
+    /// use la_stack::prelude::*;
+    ///
+    /// # fn main() -> Result<(), LaError> {
+    /// let mut matrix = RationalMatrix::<2>::zero();
+    /// assert_eq!(matrix.det_sign(), DeterminantSign::Zero);
+    /// matrix.set(0, 1, BigRational::new(1.into(), 3.into()))?;
+    /// matrix.set(1, 0, BigRational::new(1.into(), 2.into()))?;
+    /// // The exact determinant is -1/6, so its sign is negative.
+    /// assert_eq!(matrix.det_sign(), DeterminantSign::Negative);
+    /// # Ok(())
+    /// # }
+    /// ```
     pub fn det_sign(&self) -> DeterminantSign {
         let (integer_rows, _) = self.integer_rows();
         match det_big_int(integer_rows).sign() {
@@ -197,6 +243,25 @@ impl<const D: usize> RationalMatrix<D> {
     /// Denominators are cleared independently per row. If row `i` uses
     /// positive scale `sᵢ`, the integer determinant is divided by `∏ᵢ sᵢ`.
     /// For D=0, this returns the empty-product determinant `1`.
+    /// Use [`det_sign`](Self::det_sign) when only the sign is needed, or
+    /// [`ExactF64Conversion`](crate::ExactF64Conversion) to convert this result
+    /// under an explicit strict or rounded binary64 contract.
+    ///
+    /// # Examples
+    /// ```
+    /// use la_stack::prelude::*;
+    ///
+    /// # fn main() -> Result<(), LaError> {
+    /// let mut matrix = RationalMatrix::<2>::zero();
+    /// matrix.set(0, 0, BigRational::new(1.into(), 3.into()))?;
+    /// matrix.set(1, 1, BigRational::from_integer(2.into()))?;
+    /// let determinant = matrix.det();
+    /// assert_eq!(determinant, BigRational::new(2.into(), 3.into()));
+    /// // Conversion rounds only after the exact determinant has been computed.
+    /// assert_eq!(determinant.to_rounded_f64()?, 2.0 / 3.0);
+    /// # Ok(())
+    /// # }
+    /// ```
     #[must_use]
     pub fn det(&self) -> BigRational {
         let (integer_rows, row_scales) = self.integer_rows();
@@ -298,6 +363,24 @@ impl<const D: usize> RationalVector<D> {
 
     /// Try to create an exact vector by evaluating a function at every index.
     ///
+    /// The function is evaluated once per entry in increasing index order.
+    /// All entries are generated before validation and canonicalization, as
+    /// in [`try_new`](Self::try_new).
+    ///
+    /// # Examples
+    /// ```
+    /// use la_stack::prelude::*;
+    ///
+    /// # fn main() -> Result<(), LaError> {
+    /// let numerators = [1, 2, 3];
+    /// let rhs = RationalVector::<3>::try_from_fn(|index| {
+    ///     BigRational::new(numerators[index].into(), 2.into())
+    /// })?;
+    /// assert_eq!(rhs.try_to_f64()?.into_array(), [0.5, 1.0, 1.5]);
+    /// # Ok(())
+    /// # }
+    /// ```
+    ///
     /// # Errors
     /// Returns [`LaError::NonFinite`] at the first generated entry whose raw
     /// rational denominator is zero.
@@ -358,6 +441,14 @@ fn integer_at_scale(value: &BigRational, scale: &BigInt) -> BigInt {
 /// forming the larger intermediate `a × b`. The resulting positive scale
 /// clears both denominators without changing determinant sign.
 fn least_common_multiple(lhs: BigInt, rhs: &BigInt) -> BigInt {
+    // Integer entries and repeated denominators leave the current scale unchanged.
+    if rhs.is_one() || lhs == *rhs {
+        return lhs;
+    }
+    // The first non-integer entry establishes the scale without a GCD.
+    if lhs.is_one() {
+        return rhs.clone();
+    }
     let gcd = greatest_common_divisor(lhs.clone(), rhs.clone());
     (lhs / gcd) * rhs
 }
diff --git a/src/rounding.rs b/src/rounding.rs
index f7b1967..1a822a2 100644
--- a/src/rounding.rs
+++ b/src/rounding.rs
@@ -7,18 +7,21 @@
 
 /// Return the exact error in a rounded binary64 sum.
 ///
-/// This is Knuth's `TwoSum` transform. With IEEE-754 round-to-nearest and
-/// gradual underflow, `rounded + error` equals the exact-real sum whenever the
-/// rounded sum is finite.
-/// See `REFERENCES.md` \[8\] for the transform and its error-free arithmetic
-/// analysis. Callers supply finite operands and their rounded sum.
+/// This is `FastTwoSum` with operands ordered by magnitude. With IEEE-754
+/// round-to-nearest and gradual underflow, `rounded + error` equals the
+/// exact-real sum whenever the rounded sum is finite. Ordering prevents an
+/// intermediate overflow even at the finite range boundary; see
+/// `REFERENCES.md` \[17\], Theorem 5.1. Callers supply finite operands and their
+/// rounded sum.
 #[inline]
 pub(crate) const fn two_sum_error(left: f64, right: f64, rounded: f64) -> f64 {
-    let virtual_right = rounded - left;
-    let virtual_left = rounded - virtual_right;
-    let right_error = right - virtual_right;
-    let left_error = left - virtual_left;
-    left_error + right_error
+    let (large, small) = if left.abs() >= right.abs() {
+        (left, right)
+    } else {
+        (right, left)
+    };
+    let virtual_small = rounded - large;
+    small - virtual_small
 }
 
 /// Decompose a nonzero finite binary64 magnitude as `significand × 2^exponent`.
diff --git a/src/scaled_product.rs b/src/scaled_product.rs
index 945516b..c56925d 100644
--- a/src/scaled_product.rs
+++ b/src/scaled_product.rs
@@ -285,6 +285,24 @@ mod tests {
         }
     }
 
+    #[test]
+    fn non_finite_factors_remain_unrepresentable_before_or_after_zero() {
+        for non_finite in [f64::INFINITY, f64::NEG_INFINITY, f64::NAN] {
+            for zero in [0.0, -0.0] {
+                for middle_factors in [[non_finite, zero], [zero, non_finite]] {
+                    let mut product = ScaledProduct::new(false);
+                    product.multiply(1.5);
+                    for factor in middle_factors {
+                        product.multiply(factor);
+                    }
+                    product.multiply(-2.0);
+
+                    assert_eq!(product.finish(), None);
+                }
+            }
+        }
+    }
+
     #[test]
     fn balanced_extreme_factors_do_not_depend_on_storage_order() {
         let mut forward = ScaledProduct::new(false);
diff --git a/src/tolerance.rs b/src/tolerance.rs
index 012a581..fe19e7a 100644
--- a/src/tolerance.rs
+++ b/src/tolerance.rs
@@ -99,6 +99,7 @@ mod tests {
     use approx::assert_abs_diff_eq;
 
     use super::*;
+    use crate::InvalidToleranceReason;
 
     #[test]
     fn default_singular_tol_is_expected() {
@@ -135,7 +136,7 @@ mod tests {
             Tolerance::try_new(-1.0),
             Err(LaError::InvalidTolerance {
                 value: -1.0,
-                reason: crate::InvalidToleranceReason::Negative,
+                reason: InvalidToleranceReason::Negative,
             })
         );
     }
@@ -147,7 +148,7 @@ mod tests {
                 Tolerance::try_new(value),
                 Err(LaError::InvalidTolerance {
                     value: observed,
-                    reason: crate::InvalidToleranceReason::NotFinite,
+                    reason: InvalidToleranceReason::NotFinite,
                 }) if observed.to_bits() == value.to_bits()
             );
         }
diff --git a/src/vector.rs b/src/vector.rs
index ff565d6..f93da7d 100644
--- a/src/vector.rs
+++ b/src/vector.rs
@@ -27,17 +27,18 @@ use crate::{ArithmeticOperation, LaError};
 ///
 /// # Examples
 /// ```
+/// use core::assert_matches;
 /// use la_stack::prelude::*;
 ///
 /// # fn main() -> Result<(), LaError> {
 /// let left = Vector::<2>::try_new([1.0, 2.0])?;
 /// let right = Vector::<2>::try_new([3.0, 4.0])?;
-/// let bounded: ScalarWithErrorBound = left
-///     .dot_with_errbound(&right)?
-///     .expect("small integer products have a binary64 bound");
-/// assert_eq!(bounded.estimate(), 11.0);
-/// assert!(bounded.lower_bound() <= 11.0);
-/// assert!(bounded.upper_bound() >= 11.0);
+/// let certificate = left.dot_with_errbound(&right)?;
+/// assert_eq!(certificate.map(ScalarWithErrorBound::estimate), Some(11.0));
+/// // Preserve None: without a certificate, an exact fallback is needed
+/// // before making a claim about the exact-real result.
+/// let enclosure = certificate.map(|bounded| (bounded.lower_bound(), bounded.upper_bound()));
+/// assert_matches!(enclosure, Some((lower, upper)) if lower <= 11.0 && 11.0 <= upper);
 /// # Ok(())
 /// # }
 /// ```
@@ -82,8 +83,8 @@ impl ScalarWithErrorBound {
 
     /// Construct a validated public result with finite outward endpoints.
     ///
-    /// The `TwoSum` residuals determine whether either rounded endpoint must be
-    /// moved by one binary64 value. Returning `None` instead of publishing an
+    /// Exact addition residuals determine whether either rounded endpoint must
+    /// move by one binary64 value. Returning `None` instead of publishing an
     /// infinite endpoint enforces the public proof-unavailable contract.
     const fn try_new(estimate: f64, absolute_error_bound: f64) -> Option<Self> {
         if !estimate.is_finite() || !absolute_error_bound.is_finite() || absolute_error_bound < 0.0
@@ -502,11 +503,21 @@ impl<const D: usize> Vector<D> {
     /// # fn main() -> Result<(), LaError> {
     /// let left = Vector::<3>::try_new([1.0, 2.0, 3.0])?;
     /// let right = Vector::<3>::try_new([4.0, 5.0, 6.0])?;
-    /// let bounded = left
-    ///     .dot_with_errbound(&right)?
-    ///     .expect("ordinary inputs have a binary64 bound");
-    /// assert_eq!(bounded.estimate(), 32.0);
-    /// assert!(bounded.lower_bound() > 0.0);
+    /// let positive = left.dot_with_errbound(&right)?.and_then(|bounded| {
+    ///     if bounded.lower_bound() > 0.0 {
+    ///         Some(true)
+    ///     } else if bounded.upper_bound() <= 0.0 {
+    ///         Some(false)
+    ///     } else {
+    ///         None // The enclosure cannot establish whether the result is positive.
+    ///     }
+    /// });
+    /// assert_eq!(positive, Some(true));
+    ///
+    /// // Finite inputs can also lack a certificate; use an exact fallback
+    /// // before deciding the sign in this case.
+    /// let tiny = Vector::<2>::try_new([1e-200, 1e-200])?;
+    /// assert_eq!(tiny.dot_with_errbound(&tiny)?, None);
     /// # Ok(())
     /// # }
     /// ```
@@ -574,11 +585,17 @@ impl<const D: usize> Vector<D> {
     /// let axis = Vector::<2>::try_new([2.0, -1.0])?;
     /// let left = Vector::<2>::try_new([4.0, 1.0])?;
     /// let right = Vector::<2>::try_new([1.0, 3.0])?;
-    /// let bounded = axis
-    ///     .dot_difference_with_errbound(&left, &right)?
-    ///     .expect("ordinary inputs have a binary64 bound");
-    /// assert_eq!(bounded.estimate(), 8.0);
-    /// assert!(bounded.lower_bound() > 1.0);
+    /// let separated = axis.dot_difference_with_errbound(&left, &right)?.and_then(|bounded| {
+    ///     if bounded.lower_bound() > 1.0 {
+    ///         Some(true)
+    ///     } else if bounded.upper_bound() <= 1.0 {
+    ///         Some(false)
+    ///     } else {
+    ///         None // An exact fallback is needed to decide this threshold test.
+    ///     }
+    /// });
+    /// // An unavailable certificate also remains None through and_then.
+    /// assert_eq!(separated, Some(true));
     /// # Ok(())
     /// # }
     /// ```
@@ -732,7 +749,10 @@ impl<const D: usize> Vector<D> {
     ///
     /// let large = Vector::<2>::try_new([1.0e200, 1.0e200])?;
     /// assert!(large.norm()?.is_finite());
-    /// assert!(large.norm_squared().is_err());
+    /// assert_eq!(
+    ///     large.norm_squared(),
+    ///     Err(LaError::non_finite_computation_step(ArithmeticOperation::VectorSquaredNorm, 0)),
+    /// );
     /// # Ok(())
     /// # }
     /// ```
@@ -795,6 +815,82 @@ mod tests {
 
     use super::*;
 
+    fn assert_certified_proof_loss_survives_normal_terms<const D: usize>() {
+        let mut left_data = [0.0; D];
+        left_data[0] = f64::MIN_POSITIVE;
+        left_data[D - 1] = 1.0;
+        let mut right_data = [1.0; D];
+        right_data[0] = 0.5;
+        let left = Vector::new(left_data);
+        let right = Vector::new(right_data);
+
+        // The first product is subnormal; the last restores a normal estimate,
+        // but cannot restore the lost certificate for the complete reduction.
+        assert_abs_diff_eq!(left.dot(&right).unwrap(), 1.0, epsilon = 0.0);
+        assert_eq!(left.dot_with_errbound(&right), Ok(None));
+        assert_eq!(
+            left.dot_difference_with_errbound(&right, &Vector::zero()),
+            Ok(None)
+        );
+    }
+
+    fn assert_certified_proof_loss_preserves_later_overflow<const D: usize>() {
+        let mut axis_data = [0.0; D];
+        axis_data[0] = f64::MIN_POSITIVE;
+        axis_data[D - 1] = f64::MAX;
+        let axis = Vector::new(axis_data);
+        let mut left_data = [0.0; D];
+        left_data[0] = 0.5;
+        left_data[D - 1] = 2.0;
+
+        assert_eq!(
+            axis.dot_with_errbound(&Vector::new(left_data)),
+            Err(LaError::non_finite_computation_step(
+                ArithmeticOperation::VectorDotProduct,
+                D - 1,
+            ))
+        );
+        let difference_error = Err(LaError::non_finite_computation_step(
+            ArithmeticOperation::VectorDotDifference,
+            D - 1,
+        ));
+        assert_eq!(
+            axis.dot_difference_with_errbound(&Vector::new(left_data), &Vector::zero()),
+            difference_error
+        );
+
+        // Overflow in the second FMA must also remain observable after the
+        // earlier coordinate has already made certification unavailable.
+        left_data[D - 1] = 0.0;
+        let mut right_data = [0.0; D];
+        right_data[D - 1] = -2.0;
+        assert_eq!(
+            axis.dot_difference_with_errbound(&Vector::new(left_data), &Vector::new(right_data)),
+            difference_error
+        );
+    }
+
+    macro_rules! gen_certified_reduction_sequence_tests {
+        ($d:literal) => {
+            paste! {
+                #[test]
+                fn [<certified_proof_loss_survives_normal_terms_ $d d>]() {
+                    assert_certified_proof_loss_survives_normal_terms::<$d>();
+                }
+
+                #[test]
+                fn [<certified_proof_loss_preserves_later_overflow_ $d d>]() {
+                    assert_certified_proof_loss_preserves_later_overflow::<$d>();
+                }
+            }
+        };
+    }
+
+    gen_certified_reduction_sequence_tests!(2);
+    gen_certified_reduction_sequence_tests!(3);
+    gen_certified_reduction_sequence_tests!(4);
+    gen_certified_reduction_sequence_tests!(5);
+
     macro_rules! gen_vector_tests {
         ($d:literal) => {
             paste! {
@@ -1177,7 +1273,13 @@ mod tests {
         let large = Vector::<2>::new([1.0e200, -1.0e200]);
         let expected = 2.0f64.sqrt() * 1.0e200;
         assert_abs_diff_eq!(large.norm().unwrap(), expected, epsilon = 2.0e184);
-        assert!(large.norm_squared().is_err());
+        assert_eq!(
+            large.norm_squared(),
+            Err(LaError::non_finite_computation_step(
+                ArithmeticOperation::VectorSquaredNorm,
+                0,
+            )),
+        );
 
         let mixed = Vector::<4>::new([1.0e200, 1.0e-200, -f64::from_bits(1), 0.0]);
         assert_eq!(mixed.norm(), Ok(1.0e200));
diff --git a/tests/exact_conversion_boundaries.rs b/tests/exact_conversion_boundaries.rs
index 301dc00..04f131b 100644
--- a/tests/exact_conversion_boundaries.rs
+++ b/tests/exact_conversion_boundaries.rs
@@ -3,11 +3,13 @@
 #![forbid(unsafe_code)]
 #![cfg(feature = "exact")]
 
+use core::array::from_fn;
 use core::cmp::Ordering;
 
-use la_stack::prelude::*;
 use pastey::paste;
 
+use la_stack::prelude::*;
+
 const POSITIVE_ZERO_BITS: u64 = 0;
 const NEGATIVE_ZERO_BITS: u64 = 1_u64 << 63;
 const BELOW_OVERFLOW_MIDPOINT_INCREMENT: f64 = f64::from_bits(1992_u64 << 52); // 2^969
@@ -63,7 +65,7 @@ fn canonical_vector_conversions_preserve_raw_contract<const D: usize>() {
         ),
     ];
     for (value, expected) in success {
-        let raw = std::array::from_fn::<_, D, _>(|_| value.clone());
+        let raw = from_fn::<_, D, _>(|_| value.clone());
         let canonical = RationalVector::try_new(raw.clone()).unwrap();
         for actual in [
             canonical.try_to_f64(),
@@ -97,7 +99,7 @@ fn canonical_vector_conversions_preserve_raw_contract<const D: usize>() {
     ];
     for (value, reason, rounded) in failures {
         for index in 0..D {
-            let raw = std::array::from_fn(|i| match i.cmp(&index) {
+            let raw = from_fn(|i| match i.cmp(&index) {
                 Ordering::Equal => value.clone(),
                 Ordering::Greater => BigRational::from_integer(BigInt::from(1_u8) << 1024_u32),
                 Ordering::Less => raw_rational(3, -6),
diff --git a/tests/prelude_exports.rs b/tests/prelude_exports.rs
index 422e94e..cf8270d 100644
--- a/tests/prelude_exports.rs
+++ b/tests/prelude_exports.rs
@@ -11,6 +11,85 @@ use la_stack::{ERR_COEFF_2, ERR_COEFF_3, ERR_COEFF_4};
 
 const _: [f64; 3] = [ERR_COEFF_2, ERR_COEFF_3, ERR_COEFF_4];
 
+// No Clone, Copy, or Debug: dispatch must accept an ordinary consuming closure.
+struct MoveOnly(usize);
+
+macro_rules! gen_dispatch_capture_tests {
+    ($name:ident, $dispatch:ident, $entry:expr, $max:expr) => {
+        #[test]
+        fn $name() -> Result<(), LaError> {
+            for dimension in 2..=5 {
+                let mut visited = Vec::new();
+                let count = $dispatch!(dimension, |matrix| -> Result<usize, LaError> {
+                    visited.push(matrix.as_rows().len());
+                    Ok(visited.len())
+                })?;
+                assert_eq!(count, 1);
+                assert_eq!(visited, [dimension]);
+
+                let count = $dispatch!(dimension, |mut matrix| -> Result<usize, LaError> {
+                    matrix.set(0, 0, $entry)?;
+                    visited.push(matrix.as_rows().len());
+                    Ok(visited.len())
+                })?;
+                assert_eq!(count, 2);
+                assert_eq!(visited, [dimension, dimension]);
+
+                let captured = MoveOnly(dimension);
+                let moved = $dispatch!(dimension, |matrix| -> Result<MoveOnly, LaError> {
+                    assert_eq!(matrix.as_rows().len(), dimension);
+                    Ok(captured)
+                })?;
+                assert_eq!(moved.0, dimension);
+
+                let captured = MoveOnly(dimension);
+                let moved = $dispatch!(dimension, |mut matrix| -> Result<MoveOnly, LaError> {
+                    matrix.set(0, 0, $entry)?;
+                    assert_eq!(matrix.as_rows().len(), dimension);
+                    Ok(captured)
+                })?;
+                assert_eq!(moved.0, dimension);
+            }
+
+            let mut visited = Vec::new();
+            let rejected = $dispatch!(99usize, |matrix| -> Result<(), LaError> {
+                visited.push(matrix.as_rows().len());
+                Ok(())
+            });
+            assert_eq!(rejected, Err(LaError::unsupported_dimension(99, $max)));
+            assert!(visited.is_empty());
+            let rejected = $dispatch!(99usize, |mut matrix| -> Result<(), LaError> {
+                matrix.set(0, 0, $entry)?;
+                visited.push(matrix.as_rows().len());
+                Ok(())
+            });
+            assert_eq!(rejected, Err(LaError::unsupported_dimension(99, $max)));
+            assert!(visited.is_empty());
+            Ok(())
+        }
+    };
+}
+
+gen_dispatch_capture_tests!(
+    stack_dispatch_accepts_borrowed_and_consuming_captures,
+    try_with_stack_matrix,
+    1.0,
+    MAX_STACK_MATRIX_DISPATCH_DIM
+);
+gen_dispatch_capture_tests!(
+    interval_dispatch_accepts_borrowed_and_consuming_captures,
+    try_with_interval_matrix,
+    Interval::ONE,
+    MAX_INTERVAL_MATRIX_DIM
+);
+#[cfg(feature = "exact")]
+gen_dispatch_capture_tests!(
+    rational_dispatch_accepts_borrowed_and_consuming_captures,
+    try_with_rational_matrix,
+    BigRational::from_integer(1.into()),
+    MAX_RATIONAL_MATRIX_DISPATCH_DIM
+);
+
 #[test]
 fn common_prelude_supports_downstream_composition() -> Result<(), LaError> {
     let matrix = Matrix::<2>::identity();
diff --git a/tests/proptest_interval.rs b/tests/proptest_interval.rs
index 171cb03..5d45962 100644
--- a/tests/proptest_interval.rs
+++ b/tests/proptest_interval.rs
@@ -47,6 +47,16 @@ fn assert_outward_result(
         Ok(interval) => {
             prop_assert!(interval_contains_exact(interval, exact));
             prop_assert!(exact_fits_finite_interval(exact));
+            // A scalar operation needs only the exact point or the two adjacent
+            // floats bracketing it. An unnecessarily wide interval loses proofs.
+            if interval.lower() < interval.upper() {
+                prop_assert_eq!(
+                    exact_f64(interval.lower().next_up()),
+                    exact_f64(interval.upper()),
+                );
+                prop_assert!(exact_f64(interval.lower()) < *exact);
+                prop_assert!(*exact < exact_f64(interval.upper()));
+            }
         }
         Err(LaError::IntervalRangeExhausted { operation, .. }) => {
             prop_assert_eq!(operation, expected_operation);
@@ -81,25 +91,23 @@ fn assert_conclusive_sign_matches(
 fn rational_det<const D: usize>(rows: &[[f64; D]; D]) -> BigRational {
     let mut work: [[BigRational; D]; D] =
         from_fn(|row| from_fn(|column| exact_f64(rows[row][column])));
+    let zero = BigRational::from_integer(0.into());
     let mut negative = false;
 
     for column in 0..D {
-        let mut pivot_row = column;
-        while pivot_row < D && work[pivot_row][column] == BigRational::from_integer(0.into()) {
-            pivot_row += 1;
-        }
-        if pivot_row == D {
-            return BigRational::from_integer(0.into());
-        }
+        let Some(pivot_row) = (column..D).find(|&row| work[row][column] != zero) else {
+            return zero;
+        };
         if pivot_row != column {
             work.swap(pivot_row, column);
             negative = !negative;
         }
 
-        let pivot = work[column][column].clone();
-        let pivot_row = work[column].clone();
-        for row in work.iter_mut().skip(column + 1) {
-            let factor = &row[column] / &pivot;
+        let (pivot_rows, rows_below) = work.split_at_mut(column + 1);
+        let pivot_row = &pivot_rows[column];
+        let pivot = &pivot_row[column];
+        for row in rows_below {
+            let factor = &row[column] / pivot;
             for (entry, pivot_entry) in row.iter_mut().zip(pivot_row.iter()).skip(column) {
                 let reduction = &factor * pivot_entry;
                 *entry -= reduction;
@@ -204,8 +212,8 @@ proptest! {
             square,
             &(exact_f64(left_upper) * exact_f64(left_upper)),
         ));
-        if left.contains(0.0) {
-            prop_assert!(square.contains(0.0));
+        if left_lower <= 0.0 && left_upper >= 0.0 {
+            prop_assert_eq!(square.lower().to_bits(), 0);
         }
     }
 }
@@ -213,6 +221,31 @@ proptest! {
 proptest! {
     #![proptest_config(with_default_cases(512))]
 
+    #[test]
+    fn opposite_sign_extreme_sums_enclose_exact_rational_results(
+        half_ulps in 1_u16..=1024,
+        maximum_offset in prop_oneof![Just(0_u16), 1_u16..=1024],
+        negate in any::<bool>(),
+    ) {
+        let small = -f64::from(half_ulps) * f64::from_bits(1993_u64 << 52);
+        let large = f64::from_bits(f64::MAX.to_bits() - u64::from(maximum_offset));
+        let (small, large) = if negate { (-small, -large) } else { (small, large) };
+        let exact = exact_f64(small) + exact_f64(large);
+
+        for (left, right) in [(small, large), (large, small)] {
+            assert_outward_result(
+                Interval::point(left)?.try_add(&Interval::point(right)?),
+                &exact,
+                ArithmeticOperation::IntervalAddition,
+            )?;
+            assert_outward_result(
+                Interval::try_from_subtraction(left, -right),
+                &exact,
+                ArithmeticOperation::IntervalSubtraction,
+            )?;
+        }
+    }
+
     #[test]
     fn arbitrary_finite_point_operations_are_outward_or_report_true_range_loss(
         left in finite_f64(),
diff --git a/tests/proptest_matrix.rs b/tests/proptest_matrix.rs
index ea73ac2..659ec3b 100644
--- a/tests/proptest_matrix.rs
+++ b/tests/proptest_matrix.rs
@@ -73,14 +73,20 @@ macro_rules! gen_matrix_proptests {
 
                 #[test]
                 fn [<matrix_set_get_in_bounds_ $d d>](
+                    rows in array::[<uniform $d>](array::[<uniform $d>](small_f64())),
                     r in 0usize..$d,
                     c in 0usize..$d,
                     v in small_f64(),
                 ) {
-                    let mut m = Matrix::<$d>::zero();
+                    let mut m = Matrix::<$d>::try_from_rows(rows).unwrap();
+                    let mut expected = rows;
                     prop_assert_eq!(m.set(r, c, v), Ok(()));
+                    expected[r][c] = v;
+                    prop_assert_eq!(m.as_rows(), &expected);
                     assert_abs_diff_eq!(m.get(r, c).unwrap(), v, epsilon = 0.0);
                     prop_assert_eq!(m.set(r, c, -v), Ok(()));
+                    expected[r][c] = -v;
+                    prop_assert_eq!(m.as_rows(), &expected);
                     assert_abs_diff_eq!(m.try_get(r, c).unwrap(), -v, epsilon = 0.0);
                 }
 
diff --git a/tests/proptest_rational.rs b/tests/proptest_rational.rs
index daa61d0..194f579 100644
--- a/tests/proptest_rational.rs
+++ b/tests/proptest_rational.rs
@@ -7,7 +7,7 @@
 use std::array::from_fn;
 
 use pastey::paste;
-use proptest::prelude::*;
+use proptest::{collection, prelude::*};
 
 use la_stack::prelude::*;
 
@@ -52,11 +52,12 @@ fn rational_determinant_gaussian<const D: usize>(mut rows: [[BigRational; D]; D]
             odd_swaps = !odd_swaps;
         }
 
-        let pivot = rows[pivot_col][pivot_col].clone();
-        let pivot_entries = rows[pivot_col].clone();
-        determinant *= &pivot;
-        for row_entries in rows.iter_mut().skip(pivot_col + 1) {
-            let factor = &row_entries[pivot_col] / &pivot;
+        let (pivot_rows, rows_below) = rows.split_at_mut(pivot_col + 1);
+        let pivot_entries = &pivot_rows[pivot_col];
+        let pivot = &pivot_entries[pivot_col];
+        determinant *= pivot;
+        for row_entries in rows_below {
+            let factor = &row_entries[pivot_col] / pivot;
             for (entry, pivot_entry) in row_entries
                 .iter_mut()
                 .zip(pivot_entries.iter())
@@ -92,8 +93,8 @@ macro_rules! gen_rational_properties {
 
                 #[test]
                 fn [<rational_det_sign_value_and_solve_agree_with_reference_ $d d>](
-                    entries in proptest::collection::vec((-5_i16..=5_i16, 1_u8..=9_u8), $d * $d),
-                    solution_entries in proptest::collection::vec((-5_i16..=5_i16, 1_u8..=9_u8), $d),
+                    entries in collection::vec((-5_i16..=5_i16, 1_u8..=9_u8), $d * $d),
+                    solution_entries in collection::vec((-5_i16..=5_i16, 1_u8..=9_u8), $d),
                 ) {
                     let rows: [[BigRational; $d]; $d] = from_fn(|row| {
                         from_fn(|col| {
@@ -121,7 +122,7 @@ macro_rules! gen_rational_properties {
 
                 #[test]
                 fn [<rational_duplicate_row_is_singular_ $d d>](
-                    entries in proptest::collection::vec((-5_i16..=5_i16, 1_u8..=9_u8), $d * $d),
+                    entries in collection::vec((-5_i16..=5_i16, 1_u8..=9_u8), $d * $d),
                 ) {
                     let mut rows: [[BigRational; $d]; $d] = from_fn(|row| {
                         from_fn(|col| {
diff --git a/tests/proptest_vector.rs b/tests/proptest_vector.rs
index a8f8d57..d179fc4 100644
--- a/tests/proptest_vector.rs
+++ b/tests/proptest_vector.rs
@@ -38,6 +38,26 @@ macro_rules! gen_vector_proptests {
                     }
                 }
 
+                #[test]
+                fn [<vector_dot_and_squared_norm_match_integer_oracles_ $d d>](
+                    left in array::[<uniform $d>](-1000i16..=1000),
+                    right in array::[<uniform $d>](-1000i16..=1000),
+                ) {
+                    // At D<=8 these integer products and sums fit in i32 and
+                    // are exact in binary64, independently of the FMA kernel.
+                    let dot: i32 = left.iter().zip(&right)
+                        .map(|(&a, &b)| i32::from(a) * i32::from(b))
+                        .sum();
+                    let squared_norm: i32 = left.iter()
+                        .map(|&value| i32::from(value).pow(2))
+                        .sum();
+                    let a = Vector::<$d>::try_new(left.map(f64::from)).unwrap();
+                    let b = Vector::<$d>::try_new(right.map(f64::from)).unwrap();
+
+                    prop_assert_eq!(a.dot(&b), Ok(f64::from(dot)));
+                    prop_assert_eq!(a.norm_squared(), Ok(f64::from(squared_norm)));
+                }
+
                 #[test]
                 fn [<vector_dot_commutes_and_norm_squared_matches_dot_self_ $d d>](
                     a_arr in array::[<uniform $d>](small_f64()),
diff --git a/tests/regressions.rs b/tests/regressions.rs
index b017126..ec5b624 100644
--- a/tests/regressions.rs
+++ b/tests/regressions.rs
@@ -12,6 +12,31 @@ use proptest::prelude::*;
 #[path = "common/proptest_config.rs"]
 mod proptest_config;
 
+#[test]
+fn interval_extreme_sums_preserve_outward_bounds_in_both_orders() -> Result<(), LaError> {
+    // MAX - 3·2^970 is the exact midpoint between MAX's two predecessors.
+    // Its rounded sum is finite, but unsorted TwoSum can overflow internally.
+    const SMALL: f64 = -3.0 * f64::from_bits(1993_u64 << 52);
+    const DIFFERENCE: Result<Interval, LaError> = Interval::try_from_subtraction(SMALL, -f64::MAX);
+    let upper = f64::MAX.next_down();
+    let expected = Interval::try_new(upper.next_down(), upper)?;
+    assert_eq!(DIFFERENCE, Ok(expected));
+
+    for (small, large, enclosure) in [
+        (SMALL, f64::MAX, expected),
+        (-SMALL, -f64::MAX, expected.negate()),
+    ] {
+        for (left, right) in [(small, large), (large, small)] {
+            assert_eq!(
+                Interval::point(left)?.try_add(&Interval::point(right)?),
+                Ok(enclosure),
+            );
+            assert_eq!(Interval::try_from_subtraction(left, -right), Ok(enclosure));
+        }
+    }
+    Ok(())
+}
+
 /// Pad the independently identified near-overflow pair to a fixed dimension.
 fn norm_boundary_vector<const D: usize>(left: f64, right: f64) -> Vector<D> {
     let mut values = [0.0; D];