diff --git a/Cargo.lock b/Cargo.lock
index 9128e81..573b211 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -705,9 +705,9 @@ dependencies = [
 
 [[package]]
 name = "itoa"
-version = "1.0.18"
+version = "1.0.19"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "8f42a60cbdf9a97f5d2305f08a87dc4e09308d1276d28c869c684d7777685682"
+checksum = "3bd1574e617a9a6c1f8f8181bd7e841e1f20e4277607ef28945765188d6284d4"
 
 [[package]]
 name = "js-sys"
diff --git a/README.md b/README.md
index 576ceb7..1942877 100644
--- a/README.md
+++ b/README.md
@@ -238,9 +238,11 @@ for the full contracts.
 
 ¹ Requires `features = ["exact"]`.
 
-For tested usage examples, see [exact dot products and column selection][guide-exact-reductions],
-[interval division and solution enclosures][guide-intervals], and the runnable
-[polynomial determinant example][example-polynomial].
+For tested usage, see the guides to [exact dot products and column selection][guide-exact-reductions]
+and [interval division and solution enclosures][guide-intervals]. Runnable workflows cover
+[column selection][example-column-basis], [interval solves][example-interval-solve],
+[exact rational solves and objectives][example-rational-input], and
+[polynomial determinants][example-polynomial].
 
 [api-reference]: https://docs.rs/la-stack/latest/la_stack/
 [api-guide]: https://docs.rs/la-stack/latest/la_stack/guide/index.html
@@ -457,22 +459,26 @@ for current release planning.
 
 The `examples/` directory contains small, runnable programs:
 
+- **[`column_basis`][example-column-basis]** — select columns of a rectangular matrix and handle partial numerical rank
 - **[`const_det_4x4`][example-const-det]** — compile-time 4×4 determinant via `det_direct()`
 - **[`det_5x5`][example-det]** — determinant of a 5×5 matrix via LU
 - **[`exact_det_3x3`][example-exact-det]** — exact determinant value of a near-singular 3×3 matrix (requires `exact` feature)
 - **[`exact_sign_3x3`][example-exact-sign]** — exact determinant sign of a near-singular 3×3 matrix (requires `exact` feature)
 - **[`exact_solve_3x3`][example-exact-solve]** — exact solve of a near-singular 3×3 system vs f64 LU (requires `exact` feature)
+- **[`interval_solve`][example-interval-solve]** — divide by a zero-free interval, enclose all admitted solutions, and identify when fallback is required
 - **[`ldlt_solve_3x3`][example-ldlt-solve]** — solve a 3×3 symmetric positive definite system via LDLT
 - **[`polynomial_determinant`][example-polynomial]** — exact cancellation, leading coefficients, and explicit work limits (requires `exact` feature)
-- **[`rational_input_5x5`][example-rational-input]** — exact rational solve of a 5×5 system that becomes singular as f64 (requires `exact` feature)
+- **[`rational_input_5x5`][example-rational-input]** — exact solve and rational dot products for a system that becomes singular as f64 (requires `exact`)
 - **[`solve_5x5`][example-solve]** — solve a 5×5 system via LU with partial pivoting
 - **[`vector_angles`][example-vector-angles]** — unsigned angles for fixed vectors and borrowed slices, including extreme magnitudes and subnormal separations
 
+[example-column-basis]: https://github.com/acgetchell/la-stack/blob/main/examples/column_basis.rs
 [example-const-det]: https://github.com/acgetchell/la-stack/blob/main/examples/const_det_4x4.rs
 [example-det]: https://github.com/acgetchell/la-stack/blob/main/examples/det_5x5.rs
 [example-exact-det]: https://github.com/acgetchell/la-stack/blob/main/examples/exact_det_3x3.rs
 [example-exact-sign]: https://github.com/acgetchell/la-stack/blob/main/examples/exact_sign_3x3.rs
 [example-exact-solve]: https://github.com/acgetchell/la-stack/blob/main/examples/exact_solve_3x3.rs
+[example-interval-solve]: https://github.com/acgetchell/la-stack/blob/main/examples/interval_solve.rs
 [example-ldlt-solve]: https://github.com/acgetchell/la-stack/blob/main/examples/ldlt_solve_3x3.rs
 [example-polynomial]: https://github.com/acgetchell/la-stack/blob/main/examples/polynomial_determinant.rs
 [example-rational-input]: https://github.com/acgetchell/la-stack/blob/main/examples/rational_input_5x5.rs
@@ -482,12 +488,15 @@ The `examples/` directory contains small, runnable programs:
 ```bash
 just examples
 # or individually:
+cargo run --locked --example column_basis
 cargo run --example const_det_4x4
 cargo run --example det_5x5
 cargo run --features exact --example exact_det_3x3
 cargo run --features exact --example exact_sign_3x3
 cargo run --features exact --example exact_solve_3x3
+cargo run --locked --example interval_solve
 cargo run --example ldlt_solve_3x3
+cargo run --locked --features exact --example polynomial_determinant
 cargo run --features exact --example rational_input_5x5
 cargo run --example solve_5x5
 cargo run --locked --example vector_angles
diff --git a/docs/code_organization.md b/docs/code_organization.md
index fe28e64..efb29d1 100644
--- a/docs/code_organization.md
+++ b/docs/code_organization.md
@@ -116,10 +116,16 @@ The [Cargo manifest](../Cargo.toml) owns feature and dependency declarations.
 - `tooling/examples.toml` owns the Cargo example package, feature overrides and
   time budgets. The shared runner discovers every binary example from metadata;
   executable numerical assertions stay in `examples/`.
-- `examples/vector_angles.rs` demonstrates fixed-vector and borrowed-slice angles,
-  independent magnitude scaling, and subnormal separations.
+- `examples/column_basis.rs` demonstrates full and partial rectangular column
+  selections with explicit absolute tolerances and a rank check before extracting a minor.
+- `examples/interval_solve.rs` demonstrates zero-free division and reciprocals,
+  an analytical solution enclosure, and the two typed solve fallback cases.
 - `examples/polynomial_determinant.rs` demonstrates exact cancellation, borrowed
   leading coefficients and typed work-budget exhaustion.
+- `examples/rational_input_5x5.rs` demonstrates an exact rational solve and exact
+  objective evaluation with prepared vectors and checked borrowed slices.
+- `examples/vector_angles.rs` demonstrates fixed-vector and borrowed-slice angles,
+  independent magnitude scaling, and subnormal separations.
 
 Use [Testing guidance](dev/testing.md) for dimension coverage and test design,
 and [Documentation guidance](dev/docs.md) for executable-example ownership.
diff --git a/examples/column_basis.rs b/examples/column_basis.rs
new file mode 100644
index 0000000..65d2b24
--- /dev/null
+++ b/examples/column_basis.rs
@@ -0,0 +1,36 @@
+#![forbid(unsafe_code)]
+
+//! Select a provisional column basis from a rectangular matrix.
+//!
+//! See the [exact certification workflow](https://docs.rs/la-stack/latest/la_stack/guide/exact/index.html#exact-reductions-and-provisional-column-selection)
+//! for applications that need a proof about the selected minor.
+
+use std::array::from_fn;
+
+use la_stack::prelude::*;
+
+fn main() -> Result<(), LaError> {
+    // Column 0 is zero, column 2 is twice column 1, and column 3 adds a
+    // nonzero second coordinate. The expected basis uses columns 1 and 3.
+    let rows = [[0.0, 4.0, 8.0, 1.0], [0.0, 0.0, 0.0, 2.0]];
+    let basis = column_basis(&rows, Tolerance::try_new(1e-12)?)?;
+
+    // Check for full numerical row rank before extracting a square minor.
+    // Even at zero tolerance the rank is provisional because elimination rounds.
+    assert_eq!(basis.rank(), rows.len());
+    let columns = basis.columns();
+    assert_eq!(columns, &[1, 3]);
+    let minor = Matrix::<2>::try_from_rows(rows.map(|row| from_fn(|c| row[columns[c]])))?;
+    assert_eq!(minor.det()?.to_bits(), 8.0_f64.to_bits());
+    println!("Selected columns: {columns:?}");
+
+    // The absolute tolerance 2 rejects the second pivot (equal to 2).
+    // The partial selection describes this threshold's decision. The same
+    // input still has a nonzero 2×2 minor with determinant 8.
+    let partial = column_basis(&rows, Tolerance::try_new(2.0)?)?;
+    assert_eq!(partial.rank(), 1);
+    assert_eq!(partial.columns(), &[1]);
+    println!("Selected columns with tolerance 2: {:?}", partial.columns());
+
+    Ok(())
+}
diff --git a/examples/interval_solve.rs b/examples/interval_solve.rs
new file mode 100644
index 0000000..e8abef2
--- /dev/null
+++ b/examples/interval_solve.rs
@@ -0,0 +1,59 @@
+#![forbid(unsafe_code)]
+
+//! Divide by a zero-free interval and enclose the solutions of an interval system.
+
+use la_stack::prelude::*;
+
+fn main() -> Result<(), LaError> {
+    let divisor = Interval::try_new(2.0, 4.0)?;
+    assert_eq!(divisor.try_recip()?, Interval::try_new(0.25, 0.5)?);
+    let quotient = Interval::try_new(4.0, 8.0)?.try_div(&divisor)?;
+    assert_eq!(quotient, Interval::try_new(1.0, 4.0)?);
+
+    // For every a ∈ [2, 4] and b ∈ [6, 10], the equations a*x + y = b
+    // and 2*y = 4 have a unique solution with y = 2 and x = (b - 2)/a.
+    let two = Interval::point(2.0)?;
+    let matrix = IntervalMatrix::from_rows([[divisor, Interval::ONE], [Interval::ZERO, two]]);
+    let rhs = [Interval::try_new(6.0, 10.0)?, Interval::point(4.0)?];
+    let solution = matrix.solve(&rhs)?;
+    assert_eq!(solution, [Interval::try_new(1.0, 4.0)?, two]);
+    println!("Solution enclosure for every admitted system: {solution:?}");
+
+    // Direct division rejects a supplied divisor containing zero. A solve
+    // instead reports when elimination cannot certify a pivot.
+    let crosses_zero = Interval::try_new(-1.0, 1.0)?;
+    assert!(matches!(
+        Interval::ONE.try_div(&crosses_zero),
+        Err(LaError::IntervalZeroDivisor {
+            lower: -1.0,
+            upper: 1.0,
+            ..
+        })
+    ));
+
+    // The exact determinant is 3 × 2^-52. Rounding in elimination can still
+    // leave the final pivot enclosing zero.
+    let regular = IntervalMatrix::try_from_point_rows([[3.0, 3.0], [1.0, 1.0_f64.next_up()]])?;
+    match regular.solve(&[Interval::ONE; 2]) {
+        Err(LaError::IntervalSolveInconclusive { pivot_col, .. }) => {
+            assert_eq!(pivot_col, 1);
+            println!("Pivot {pivot_col} is inconclusive. Use a certified fallback.");
+        }
+        Err(error) => return Err(error),
+        Ok(bounds) => panic!("expected an inconclusive pivot, got {bounds:?}"),
+    }
+
+    // Here the unique real solution is 2 × f64::MAX, outside the finite
+    // endpoint range. A certified fallback needs a wider representation.
+    let small_pivot = IntervalMatrix::try_from_point_rows([[0.5]])?;
+    match small_pivot.solve(&[Interval::point(f64::MAX)?]) {
+        Err(LaError::IntervalRangeExhausted { operation, .. }) => {
+            assert_eq!(operation, ArithmeticOperation::IntervalSolve);
+            println!("Finite endpoint range exhausted. Use a certified fallback.");
+        }
+        Err(error) => return Err(error),
+        Ok(bounds) => panic!("expected finite endpoint range exhaustion, got {bounds:?}"),
+    }
+
+    Ok(())
+}
diff --git a/examples/rational_input_5x5.rs b/examples/rational_input_5x5.rs
index a878fd3..76cb21e 100644
--- a/examples/rational_input_5x5.rs
+++ b/examples/rational_input_5x5.rs
@@ -1,6 +1,8 @@
 #![forbid(unsafe_code)]
 
-//! Solve a rational 5×5 system that becomes singular after conversion to binary64.
+//! Solve a rational 5×5 system and evaluate an exact linear objective.
+//!
+//! The reconstructed binary64 system is singular.
 
 use la_stack::prelude::*;
 
@@ -34,6 +36,32 @@ fn main() -> Result<(), LaError> {
     println!("exact determinant sign: {:?}", matrix.det_sign());
     println!("exact solution: {:?}", exact_solution.as_array());
 
+    let objective = RationalVector::try_new([
+        BigRational::new(1.into(), 3.into()),
+        BigRational::new(1.into(), 2.into()),
+        BigRational::from_integer(0.into()),
+        BigRational::from_integer(0.into()),
+        BigRational::new(1.into(), 4.into()),
+    ])?;
+    // For the known solution [1, -1, 2, 3, 4], the objective is 1/3 - 1/2 + 1.
+    let expected = BigRational::new(5.into(), 6.into());
+    let value = objective.dot(&exact_solution);
+    assert_eq!(value, expected);
+
+    // Borrowed slices support runtime lengths and check shape and denominators.
+    let weights = objective.as_array().as_slice();
+    let values = exact_solution.as_array().as_slice();
+    assert_eq!(rational_dot(weights, values)?, expected);
+    assert!(matches!(
+        rational_dot(weights, &values[..4]),
+        Err(LaError::DimensionMismatch {
+            left: 5,
+            right: 4,
+            ..
+        })
+    ));
+    println!("exact objective: {value}");
+
     let epsilon_f64 = epsilon.try_to_f64()?;
     assert_eq!((1.0 + epsilon_f64).to_bits(), 1.0_f64.to_bits());
     let f64_matrix = Matrix::<5>::try_from_rows([
diff --git a/pyproject.toml b/pyproject.toml
index 27c14ab..dd52575 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -185,7 +185,7 @@ exclude = [ "docs/archive/**", "docs/performance/**", "docs/performance-runs/**"
 path = "README.md"
 pattern = 'https://(?:github\.com/acgetchell/la-stack/(?:blob|raw|tree)/|raw\.githubusercontent\.com/acgetchell/la-stack/)(?P<value>main|v[0-9]+\.[0-9]+\.[0-9]+|[0-9a-f]{7,40})/[^\s)\]>"?]+'
 value = "main"
-count = 42
+count = 44
 exclude = '/docs/(?:assets/bench|performance-runs)/'
 
 [[tool.research-repo-tools.release.rules]]
diff --git a/src/interval.rs b/src/interval.rs
index d389c31..4a794c2 100644
--- a/src/interval.rs
+++ b/src/interval.rs
@@ -2,13 +2,10 @@
 
 //! Outward-rounded intervals, fixed-size determinants, and solution enclosures.
 //!
-//! 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
-//! IEEE 1788 conformance. The
+//! The mathematical basis explains the
 //! [interval construction](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#outward-rounded-interval-expressions)
-//! explains outward endpoints and determinant enclosures.
+//! and the
+//! [solution enclosure proof](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#interval-linear-solves).
 
 use crate::rounding::{compare_product_with_rounded, compare_quotient_with_rounded, two_sum_error};
 use crate::{ArithmeticOperation, IntervalBound, IntervalOperand, LaError, Matrix};
@@ -29,8 +26,12 @@ pub const MAX_INTERVAL_MATRIX_DIM: usize = 7;
 /// zeros are accepted and canonicalized to `+0.0`; subnormal bounds are
 /// retained.
 ///
-/// This is a deliberately small proof-bearing surface, not a general-purpose
-/// interval arithmetic package. Division requires a zero-free divisor.
+/// Division requires a zero-free divisor. For background on interval inclusion,
+/// see \[[20](https://github.com/acgetchell/la-stack/blob/main/REFERENCES.md#ref-20)\]
+/// and the broader IEEE interval standard
+/// \[[14](https://github.com/acgetchell/la-stack/blob/main/REFERENCES.md#ref-14)\].
+/// The [interval construction](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#outward-rounded-interval-expressions)
+/// explains how this crate encloses results with finite binary64 endpoints.
 ///
 /// # Examples
 /// ```
@@ -79,6 +80,12 @@ pub enum IntervalDeterminantSign {
 /// and explicitly reports inconclusive pivots. It has no determinant dimension
 /// limit, uses O(D²) stack storage, and is intended for small systems.
 ///
+/// The Leibniz identity \[[12](https://github.com/acgetchell/la-stack/blob/main/REFERENCES.md#ref-12)\]
+/// and interval matrix arithmetic \[[20](https://github.com/acgetchell/la-stack/blob/main/REFERENCES.md#ref-20)\]
+/// provide the background for the crate's
+/// [determinant enclosure](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#determinant-enclosures-and-sign-evidence)
+/// and [fixed-pivot solve proof](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#interval-linear-solves).
+///
 /// # Examples
 /// ```
 /// use la_stack::prelude::*;
@@ -312,6 +319,9 @@ impl Interval {
     ///
     /// Unlike subtracting first and then calling [`point`](Self::point), this
     /// method preserves the rounding uncertainty introduced by the subtraction.
+    /// The outward endpoint uses a magnitude-ordered `FastTwoSum` residual
+    /// \[[17](https://github.com/acgetchell/la-stack/blob/main/REFERENCES.md#ref-17)\].
+    /// See the [scalar rounding argument](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#scalar-arithmetic-and-finite-bounds).
     ///
     /// # Examples
     /// ```
@@ -382,6 +392,10 @@ impl Interval {
 
     /// Add two intervals with outward rounding.
     ///
+    /// The endpoint sums use `FastTwoSum`
+    /// \[[17](https://github.com/acgetchell/la-stack/blob/main/REFERENCES.md#ref-17)\]
+    /// as described in the [scalar rounding argument](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#scalar-arithmetic-and-finite-bounds).
+    ///
     /// # Examples
     /// ```
     /// use la_stack::prelude::*;
@@ -407,6 +421,10 @@ impl Interval {
     /// For a square of the same represented value, prefer
     /// [`try_square`](Self::try_square), which can give a tighter enclosure.
     ///
+    /// See `TwoProductFMA` \[[18](https://github.com/acgetchell/la-stack/blob/main/REFERENCES.md#ref-18)\]
+    /// and the [product rounding argument](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#scalar-arithmetic-and-finite-bounds)
+    /// for the residual threshold and exact integer comparison near underflow.
+    ///
     /// # Examples
     /// ```
     /// use la_stack::prelude::*;
@@ -429,6 +447,10 @@ impl Interval {
 
     /// Enclose `1 / x` for every real `x` in this interval.
     ///
+    /// This uses the same outward quotient bounds as [`try_div`](Self::try_div).
+    /// See interval inclusion \[[20](https://github.com/acgetchell/la-stack/blob/main/REFERENCES.md#ref-20)\]
+    /// and the [reciprocal construction](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#division-and-reciprocals).
+    ///
     /// # Examples
     /// ```
     /// use la_stack::prelude::*;
@@ -454,6 +476,9 @@ impl Interval {
     /// retain their enclosure, including when ordinary division rounds to zero.
     /// Repeated operands are treated independently: `x / x` need not be `[1, 1]`.
     ///
+    /// See interval inclusion \[[20](https://github.com/acgetchell/la-stack/blob/main/REFERENCES.md#ref-20)\]
+    /// and the [exact quotient comparison](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#division-and-reciprocals).
+    ///
     /// # Examples
     /// ```
     /// use la_stack::prelude::*;
@@ -497,6 +522,11 @@ 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.
     ///
+    /// The endpoint squares share the product rounding kernel, using `TwoProductFMA`
+    /// \[[18](https://github.com/acgetchell/la-stack/blob/main/REFERENCES.md#ref-18)\]
+    /// within the range justified by the
+    /// [scalar rounding argument](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#scalar-arithmetic-and-finite-bounds).
+    ///
     /// # Examples
     /// ```
     /// use la_stack::prelude::*;
@@ -1001,9 +1031,11 @@ impl<const D: usize> IntervalMatrix<D> {
     ///
     /// Uses outward-rounded Gaussian elimination and back substitution, with
     /// row pivots maximizing their minimum distance from zero (first row wins
-    /// ties). No tolerance admits a zero-containing pivot. See `REFERENCES.md`
-    /// \[20\] and the
-    /// [enclosure proof](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#interval-linear-solves).
+    /// ties). Each accepted pivot excludes zero. For interval matrix arithmetic,
+    /// see \[[20](https://github.com/acgetchell/la-stack/blob/main/REFERENCES.md#ref-20)\].
+    /// The [fixed-pivot enclosure proof](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#interval-linear-solves)
+    /// explains why the same row permutation and interval updates cover every
+    /// admitted real system.
     ///
     /// Both shapes are fixed by `D`. All const dimensions are supported, with
     /// O(D³) work and O(D²) inline scratch space intended for small systems;
@@ -1106,12 +1138,13 @@ impl<const D: usize> IntervalMatrix<D> {
 
     /// Enclose the determinant with division-free subset dynamic programming.
     ///
-    /// For each column subset, the DP stores the determinant interval of the
-    /// leading rows and those columns. This evaluates the Leibniz expansion in
-    /// `D × 2^(D-1)` products and additions without choosing or dividing by a
-    /// pivot. The returned interval therefore encloses every exact-real
-    /// determinant represented by the input intervals, subject only to an
-    /// explicit range failure.
+    /// The crate's [subset recurrence](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#determinant-enclosures-and-sign-evidence)
+    /// stores the determinant interval of the leading rows and each column
+    /// subset. It evaluates the Leibniz expansion
+    /// \[[12](https://github.com/acgetchell/la-stack/blob/main/REFERENCES.md#ref-12)\]
+    /// in `D × 2^(D-1)` products and additions without choosing or dividing by
+    /// a pivot. On success the interval encloses every exact-real determinant
+    /// represented by the input intervals.
     ///
     /// The D=0 determinant follows the empty-product convention and is `[1, 1]`.
     ///
@@ -1205,6 +1238,11 @@ impl<const D: usize> IntervalMatrix<D> {
     /// singleton interval `[0, 0]` proves `Zero`; every other overlap with zero
     /// is [`IntervalDeterminantSign::Inconclusive`].
     ///
+    /// This classifies the outward enclosure of the Leibniz determinant
+    /// \[[12](https://github.com/acgetchell/la-stack/blob/main/REFERENCES.md#ref-12)\]
+    /// returned by [`det`](Self::det). See the
+    /// [sign evidence argument](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#determinant-enclosures-and-sign-evidence).
+    ///
     /// # Examples
     /// ```
     /// use la_stack::prelude::*;
diff --git a/src/rounding.rs b/src/rounding.rs
index 05df0be..8ed540d 100644
--- a/src/rounding.rs
+++ b/src/rounding.rs
@@ -2,17 +2,19 @@
 
 //! Shared binary64 rounding primitives for certified arithmetic.
 //!
-//! Representation and rounding use the IEEE 754 model in `REFERENCES.md`
-//! \[9-10\]. The interval and reduction certificates share these primitives.
+//! The interval and reduction certificates share these primitives. They use
+//! IEEE 754 binary64 \[[9](https://github.com/acgetchell/la-stack/blob/main/REFERENCES.md#ref-9)\]
+//! with round-to-nearest, ties-to-even, and gradual underflow.
 
 /// Return the exact error in a rounded binary64 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.
+/// This uses magnitude-ordered `FastTwoSum`
+/// \[[17](https://github.com/acgetchell/la-stack/blob/main/REFERENCES.md#ref-17)\].
+/// With binary64 round-to-nearest, ties-to-even, and gradual underflow,
+/// `rounded + error` in exact arithmetic equals the real sum. Callers supply
+/// finite operands and their finite rounded sum. Magnitude ordering prevents
+/// intermediate overflow under this precondition (Theorem 5.1).
+/// See the [scalar rounding argument](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#scalar-arithmetic-and-finite-bounds).
 #[inline]
 pub(crate) const fn two_sum_error(left: f64, right: f64, rounded: f64) -> f64 {
     let (large, small) = if left.abs() >= right.abs() {
@@ -25,6 +27,9 @@ pub(crate) const fn two_sum_error(left: f64, right: f64, rounded: f64) -> f64 {
 }
 
 /// Decompose a nonzero finite binary64 magnitude as `significand × 2^exponent`.
+///
+/// The implicit leading bit is set only for normal numbers. See the binary64
+/// encoding in \[[9](https://github.com/acgetchell/la-stack/blob/main/REFERENCES.md#ref-9)\], §3.4.
 #[inline]
 pub(crate) const fn decompose_magnitude(value: f64) -> (u128, i64) {
     let magnitude_bits = value.to_bits() & 0x7fff_ffff_ffff_ffff;
@@ -42,6 +47,10 @@ pub(crate) const fn decompose_magnitude(value: f64) -> (u128, i64) {
 }
 
 /// Compare two positive values represented as `significand × 2^exponent`.
+///
+/// This integer comparison uses the binary64 decomposition above. Alignment
+/// is needed only when both highest set bits have the same exponent, so the
+/// aligned significands still fit in `u128`.
 #[inline]
 const fn compare_binary_magnitudes(
     left_significand: u128,
@@ -86,9 +95,12 @@ const fn compare_binary_magnitudes(
 /// Compare the exact-real product `left × right` with its rounded result.
 ///
 /// Nonzero finite operands contribute at most 53 significand bits each, so the
-/// exact product fits in `u128`. Above the residual-underflow threshold,
-/// `TwoProductFMA` gives its exact rounding error; see `REFERENCES.md` \[18\].
-/// Integer comparison handles smaller products, including rounding to zero.
+/// exact significand product fits in `u128`. For `|rounded| ≥ 2^-968`,
+/// `TwoProductFMA` \[[18](https://github.com/acgetchell/la-stack/blob/main/REFERENCES.md#ref-18)\]
+/// gives the exact rounding error. The crate's
+/// [residual threshold argument](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#scalar-arithmetic-and-finite-bounds)
+/// explains the bound and the integer comparison for smaller products,
+/// including rounding to zero.
 /// Callers supply nonzero finite operands and their finite rounded product.
 #[inline]
 pub(crate) const fn compare_product_with_rounded(left: f64, right: f64, rounded: f64) -> i8 {
@@ -146,10 +158,14 @@ const fn compare_product_with_rounded_integer(left: f64, right: f64, rounded: f6
     }
 }
 
-/// Compare the exact quotient with a finite rounded quotient, without relying
-/// on a floating residual that could underflow. The denominator is nonzero;
-/// all three arguments are finite. A 106-bit product suffices to compare
-/// `|rounded| * |denominator|` with `|numerator|` exactly.
+/// Compare the exact quotient with a finite rounded quotient using integers.
+///
+/// The denominator is nonzero and all three arguments are finite. A 106-bit
+/// significand product suffices to compare `|rounded| * |denominator|` with
+/// `|numerator|` exactly. The crate's
+/// [quotient comparison argument](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#division-and-reciprocals)
+/// uses binary64 decomposition and works even when a floating residual would
+/// underflow.
 #[inline]
 pub(crate) const fn compare_quotient_with_rounded(
     numerator: f64,