feat(vector): add certified dot and affine error bounds

Cargo.toml

@@ -86,6 +86,11 @@ name = "interval"
 harness = false
 required-features = [ "bench" ]
 
+[[bench]]
+name = "linear_form"
+harness = false
+required-features = [ "bench" ]
+
 [profile.release]
 lto = "fat"
 codegen-units = 1

README.md

@@ -25,6 +25,8 @@ while keeping the API intentionally small and explicit.
 - `Matrix<const D: usize>` for fixed-size square `f64` matrices backed by `[[f64; D]; D]`
 - `Interval` and `IntervalMatrix<const D: usize>` for outward-rounded,
   proof-bearing determinant filters through D=7
+- `ScalarWithErrorBound` for proof-bearing fixed-vector dot products and
+  affine differences over finite `f64` inputs
 - `RationalVector<const D: usize>` and `RationalMatrix<const D: usize>` for
   exact rational inputs behind the optional `"exact"` feature
 - `Lu<const D: usize>` for LU factorization with partial pivoting (solve + det)
@@ -40,7 +42,9 @@ determinants through D=4. These results remain subject to conditioning and
 binary64 rounding;
 factorization tolerances are rejection thresholds, not accuracy guarantees. For
 D≤4, direct determinants can be paired with a conservative absolute roundoff
-bound when its range preconditions hold.
+bound when its range preconditions hold. Fixed-vector dot products and direct
+affine differences can likewise return a paired estimate and certified absolute
+roundoff bound without enabling arbitrary-precision dependencies.
 
 Derived binary64 expressions can instead be assembled with `Interval`
 subtraction, addition, multiplication, negation, and square. The resulting
@@ -68,8 +72,9 @@ for the algorithms, validity boundaries, and supporting references.
   non-finite classification, exact fallbacks, and reproducibility contract;
   deliberate `f64::mul_add` remains allowed for its defined single-rounding
   semantics
-- ✅ Error-bounded f64 determinant filtering plus optional exact signs
-  (`det_errbound`, `det_sign_exact`)
+- ✅ Error-bounded f64 dot, affine-difference, and determinant filtering plus
+  optional exact signs (`dot_with_errbound`, `dot_difference_with_errbound`,
+  `det_errbound`, `det_sign_exact`)
 - ✅ Outward-rounded interval expressions and division-free determinant signs
   through D=7, with explicit inconclusive evidence
 - ✅ Exact determinant values and linear solves via optional arbitrary-precision
@@ -107,6 +112,8 @@ for current release planning.
   for fixed-size systems
 - You need a cheap, sound interval filter for determinant expressions assembled
   from rounded binary64 operations
+- You need a certified sign or threshold comparison for a fixed-vector dot
+  product or `axis · (left - right)` expression
 - Robust predicates matter for geometry-style workloads near degeneracy
 - You prefer a default build with no runtime dependencies
 
@@ -558,6 +565,62 @@ filter and uses fraction-free Bareiss elimination in `BigInt`.
 Because `Matrix` stores only finite entries, arithmetic range failures in the
 filter are inconclusive rather than errors and the exact fallback is total.
 
+## 🎯 Certified dot products and affine differences
+
+`Vector::dot_with_errbound()` evaluates the same left-to-right FMA tree as
+`Vector::dot()` and returns its estimate together with a certified absolute
+roundoff bound. `Vector::dot_difference_with_errbound()` directly evaluates
+
+```text
+Σᵢ axis[i] × (left[i] - right[i])
+```
+
+as two FMAs per coordinate. It does not first round `left - right` into a new
+`Vector`, so the certificate covers the intended expression over the original
+stored binary64 coordinates.
+
+The opaque `ScalarWithErrorBound` exposes the estimate, absolute error bound,
+and finite outward-rounded lower and upper bounds. Those endpoints support
+positive, negative, and caller-selected threshold proofs:
+
+```rust
+use la_stack::prelude::*;
+
+fn is_separated<const D: usize>(
+    axis: &Vector<D>,
+    left: &Vector<D>,
+    right: &Vector<D>,
+    threshold: f64,
+) -> Result<Option<bool>, LaError> {
+    let Some(value) = axis.dot_difference_with_errbound(left, right)? else {
+        return Ok(None);
+    };
+    if value.lower_bound() > threshold {
+        Ok(Some(true))
+    } else if value.upper_bound() <= threshold {
+        Ok(Some(false))
+    } else {
+        Ok(None)
+    }
+}
+
+# fn main() -> Result<(), LaError> {
+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])?;
+assert_eq!(is_separated(&axis, &left, &right, 1.0)?, Some(true));
+# Ok(())
+# }
+```
+
+An interval that overlaps the threshold is inconclusive, not equal. Likewise,
+`Ok(None)` means gradual underflow or proof-only range exhaustion prevented a
+certificate. A filtered-exact caller should rebuild the same dot or affine
+expression in `BigRational` (available through the `exact` feature) or another
+exact backend. A `LaError::NonFinite` instead reports that the specified FMA
+estimate itself overflowed. These certified bounds describe roundoff in a fixed
+arithmetic tree; they are distinct from user-selected numerical tolerances.
+
 ## 🛡️ Adaptive determinant filtering (D ≤ 4)
 
 `det_direct_with_errbound()` returns a closed-form determinant together with
@@ -635,14 +698,15 @@ out of the common prelude.
 
 | Type | Storage | Purpose | Key methods |
 |---|---|---|---|
-| `Vector<D>` | `[f64; D]` | Finite fixed-length vector for input and computation | `try_new`, `as_array`, `into_array`, `dot`, `norm2_sq` |
+| `Vector<D>` | `[f64; D]` | Finite fixed-length vector for input and computation | `try_new`, `as_array`, `into_array`, `dot`, `dot_with_errbound`, `dot_difference_with_errbound`, `norm2_sq` |
 | `Matrix<D>` | `[[f64; D]; D]` | Finite square matrix for input and computation | See below |
 | `Interval` | Two finite ordered `f64` bounds | Outward-rounded exact-real enclosure | `try_new`, `point`, `try_from_subtraction`, `try_add`, `try_mul`, `negate`, `try_square` |
 | `IntervalMatrix<D>` | `[[Interval; D]; D]` | Division-free determinant enclosure and sign proof through D=7 | `from_rows`, `try_from_point_rows`, `from_matrix`, `det`, `det_sign` |
 | `IntervalDeterminantSign` | enum | Positive, negative, zero, or inconclusive determinant evidence | — |
 | `RationalVector<D>`¹ | `[BigRational; D]` | Exact rational right-hand side and solution | `try_new`, `try_from_fn`, `as_array`, `into_array`, `get` |
 | `RationalMatrix<D>`¹ | `[[BigRational; D]; D]` | Exact rational matrix for determinant and solve operations | `try_from_rows`, `try_from_fn`, `as_rows`, `det_sign`, `det`, `solve` |
 | `DeterminantWithErrorBound` | Opaque validated pair | Paired direct determinant and certified absolute bound | `determinant`, `absolute_error_bound` |
+| `ScalarWithErrorBound` | Opaque validated certificate | Paired scalar estimate, absolute bound, and outward endpoints | `estimate`, `absolute_error_bound`, `lower_bound`, `upper_bound` |
 | `Lu<D>` | Inline factors + permutation | Factorization for solves/det | `solve`, `det` |
 | `Ldlt<D>` | Inline factors | No-pivot SPD factorization for solves/det | `solve`, `det` |
 | `Tolerance` | finite non-negative `f64` | Validated numerical threshold | `try_new`, `get` |
@@ -753,6 +817,11 @@ relative-coordinate lifted determinant signs at D=4 and the maximum supported
 D=7 workload. Run it with `just bench-interval`; fixture validation stays
 outside the timed closures.
 
+The focused `linear_form` Criterion suite compares plain and certified dot
+products and covers both well-separated and inconclusive dot/affine-difference
+filters at D=4. Run it with `just bench-linear-form`; exact small-integer fixture
+expectations are validated outside the timed closures.
+
 <!-- BENCH_TABLE:lu_solve:median:new:BEGIN -->
 
 | D | la-stack median (ns) | nalgebra median (ns) | faer median (ns) | reduction vs nalgebra (point est.) | reduction vs faer (point est.) |

REFERENCES.md

@@ -23,6 +23,19 @@ No generated content was used without human oversight.
 
 ## Linear algebra algorithms
 
+### Certified fixed-vector reductions
+
+`Vector::dot_with_errbound()` and `Vector::dot_difference_with_errbound()` use
+deterministic left-to-right binary64 FMA reductions. When their rounded
+intermediates stay normal or are exact zeros, the standard
+`gamma_n = nu / (1 - nu)` model bounds the absolute forward error by
+`gamma_n Σ |a_i b_i|` (references 9–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
+`axis_i × left_i` and `-axis_i × right_i` FMAs, so its certificate covers the
+original coordinates rather than an already-rounded difference vector.
+
 ### Outward-rounded interval determinant sign
 
 `Interval` uses IEEE-754 round-to-nearest binary64 operations plus adjacent

benches/linear_form.rs

@@ -0,0 +1,110 @@
+#![forbid(unsafe_code)]
+
+//! Criterion coverage for certified dot products and affine differences.
+
+use std::hint::black_box;
+
+use criterion::Criterion;
+
+use la_stack::Vector;
+
+#[path = "common/bench_utils.rs"]
+mod bench_utils;
+use bench_utils::OrAbort;
+
+fn main() {
+    let axis =
+        Vector::<4>::try_new([2.0, -1.0, 3.0, 4.0]).or_abort("well-separated axis construction");
+    let left = Vector::<4>::try_new([4.0, 1.0, 2.0, 3.0])
+        .or_abort("well-separated left vector construction");
+    let right = Vector::<4>::try_new([1.0, 3.0, 0.0, 2.0])
+        .or_abort("well-separated right vector construction");
+    let cancellation_left =
+        Vector::<4>::try_new([1.0, 1.0, 0.0, 0.0]).or_abort("cancellation left construction");
+    let cancellation_right =
+        Vector::<4>::try_new([1.0, -1.0, 7.0, -9.0]).or_abort("cancellation right construction");
+
+    let plain_dot = axis
+        .dot(&left)
+        .or_abort("well-separated plain dot validation");
+    assert_eq!(plain_dot.to_bits(), 25.0_f64.to_bits());
+    let bounded_dot = axis
+        .dot_with_errbound(&left)
+        .or_abort("well-separated bounded dot validation")
+        .or_abort("well-separated bounded dot certificate");
+    assert_eq!(bounded_dot.estimate().to_bits(), plain_dot.to_bits());
+    assert!(bounded_dot.lower_bound() > 0.0);
+
+    let bounded_difference = axis
+        .dot_difference_with_errbound(&left, &right)
+        .or_abort("well-separated bounded difference validation")
+        .or_abort("well-separated bounded difference certificate");
+    assert_eq!(bounded_difference.estimate().to_bits(), 18.0_f64.to_bits());
+    assert!(bounded_difference.lower_bound() > 1.0);
+
+    let inconclusive_dot = cancellation_left
+        .dot_with_errbound(&cancellation_right)
+        .or_abort("inconclusive bounded dot validation")
+        .or_abort("inconclusive bounded dot certificate");
+    assert_eq!(inconclusive_dot.estimate().to_bits(), 0.0_f64.to_bits());
+    assert!(inconclusive_dot.lower_bound() < 0.0 && inconclusive_dot.upper_bound() > 0.0);
+
+    let inconclusive_difference = axis
+        .dot_difference_with_errbound(&left, &left)
+        .or_abort("inconclusive bounded difference validation")
+        .or_abort("inconclusive bounded difference certificate");
+    assert_eq!(
+        inconclusive_difference.estimate().to_bits(),
+        0.0_f64.to_bits()
+    );
+    assert!(
+        inconclusive_difference.lower_bound() < 0.0 && inconclusive_difference.upper_bound() > 0.0
+    );
+
+    let mut criterion = Criterion::default().configure_from_args();
+    {
+        let mut group = criterion.benchmark_group("linear_form_d4");
+        group.bench_function("dot_plain_well_separated", |bencher| {
+            bencher.iter(|| {
+                let result = black_box(&axis)
+                    .dot(black_box(&left))
+                    .or_abort("well-separated plain dot");
+                let _ = black_box(result);
+            });
+        });
+        group.bench_function("dot_bounded_well_separated", |bencher| {
+            bencher.iter(|| {
+                let result = black_box(&axis)
+                    .dot_with_errbound(black_box(&left))
+                    .or_abort("well-separated bounded dot");
+                let _ = black_box(result);
+            });
+        });
+        group.bench_function("dot_bounded_inconclusive", |bencher| {
+            bencher.iter(|| {
+                let result = black_box(&cancellation_left)
+                    .dot_with_errbound(black_box(&cancellation_right))
+                    .or_abort("inconclusive bounded dot");
+                let _ = black_box(result);
+            });
+        });
+        group.bench_function("dot_difference_bounded_well_separated", |bencher| {
+            bencher.iter(|| {
+                let result = black_box(&axis)
+                    .dot_difference_with_errbound(black_box(&left), black_box(&right))
+                    .or_abort("well-separated bounded difference");
+                let _ = black_box(result);
+            });
+        });
+        group.bench_function("dot_difference_bounded_inconclusive", |bencher| {
+            bencher.iter(|| {
+                let result = black_box(&axis)
+                    .dot_difference_with_errbound(black_box(&left), black_box(&left))
+                    .or_abort("inconclusive bounded difference");
+                let _ = black_box(result);
+            });
+        });
+        group.finish();
+    }
+    criterion.final_summary();
+}

docs/BENCHMARKING.md

@@ -29,6 +29,7 @@ the commands measure and where their outputs go.
 | Fast saved-baseline loop | `just bench-save-baseline <name> <suite>` then `just bench-compare <name> <suite> all-benches` |
 | Full crate comparison | `just bench-vs-linalg` |
 | Interval determinant filter | `just bench-interval` |
+| Certified dot/linear-form filter | `just bench-linear-form` |
 | README table and plot | `just performance-release` then `just performance-readme` |
 | Release report | `just performance-release v0.4.5 v0.4.4` |
 | Build docs from retained release inputs | `just performance-doc` |
@@ -64,7 +65,7 @@ promotion in one command.
 
 ## Benchmark Suites
 
-`la-stack` has three Criterion benchmark suites.
+`la-stack` has four 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
@@ -97,6 +98,14 @@ limit. Fixture construction and expected-sign validation occur outside the
 timed closures. This suite is a focused kernel signal; it is not part of the
 release-to-release `vs_linalg` or `exact` report schema.
 
+**`linear_form`** (`benches/linear_form.rs`) measures the default-feature
+certified dot-product and affine-difference filters at D=4. It compares the
+well-separated bounded dot product with the same plain `Vector::dot` input and
+also covers dot and affine-difference cases whose certified intervals overlap
+zero. Fixture construction and exact small-integer expectations are validated
+outside the timed closures. This focused kernel signal is not part of the
+release-to-release report schema.
+
 ## Common Workflows
 
 ### Compare Current Code With The Latest Release

docs/mathematical_basis.md

@@ -48,10 +48,11 @@ fixed dimension” is the intended performance scope. `try_with_stack_matrix!` i
 separately limited to `D = 0..=MAX_STACK_MATRIX_DISPATCH_DIM` (currently 7)
 because it enumerates concrete stack types.
 
-Except for the determinant filter described below, the floating-point APIs do
-not provide certified forward, backward, or absolute error bounds. This includes
-dot products, squared norms, matrix norms, factorizations, and solves. Some
-kernels use FMA to reduce rounding steps, but that does not make them exact.
+Except for the fixed-vector reduction and determinant filters described below,
+the floating-point APIs do not provide certified forward, backward, or absolute
+error bounds. This includes plain `Vector::dot`, squared norms, matrix norms,
+factorizations, and solves. Some kernels use FMA to reduce rounding steps, but
+that does not make them exact.
 
 ## Outward-rounded interval expressions
 
@@ -97,6 +98,51 @@ interval operations, use `IntervalMatrix::det_sign()` as a fast proof, and
 rebuild the expression in `RationalMatrix` or another exact representation when
 the filter is inconclusive or loses range.
 
+## Certified fixed-vector reductions
+
+`Vector::dot_with_errbound()` binds the ordinary left-to-right FMA estimate to
+a certified absolute error bound for the exact-real dot product of the stored
+binary64 inputs. Starting from `s₀ = 0`, its arithmetic tree is
+
+```text
+sᵢ₊₁ = fma(leftᵢ, rightᵢ, sᵢ).
+```
+
+`Vector::dot_difference_with_errbound()` targets the exact-real expression
+
+```text
+Σᵢ axisᵢ(leftᵢ - rightᵢ)
+```
+
+without rounding coordinate differences first. For each coordinate it applies
+`fma(axisᵢ, leftᵢ, s)` followed by `fma(-axisᵢ, rightᵢ, s)`, giving a specified
+`2D`-event tree.
+
+Let `u = 2^-53` be binary64 unit roundoff and
+`γₙ = nu / (1 - nu)`. When every estimate FMA result is normal or an exact
+zero, standard floating-point reduction analysis gives [9-11]
+
+```text
+|estimate - exact value| ≤ γₙ Σⱼ |aⱼbⱼ|,
+```
+
+with `n = D` for a dot product and `n = 2D` for the affine difference. The
+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`.
+
+The relative-error argument is not used across gradual underflow. A nonzero
+product or estimate FMA in the subnormal range, an invalid `γₙ`, or finite-range
+exhaustion in proof-only arithmetic returns `Ok(None)`. This is inconclusive
+evidence, not equality. A non-finite estimate FMA returns `LaError::NonFinite`
+with the failing reduction index and operation. A returned
+`ScalarWithErrorBound` can certify a sign or threshold comparison from its
+lower/upper endpoints; overlap requires an exact fallback that reconstructs the
+same expression over the original inputs. The certificate is a roundoff bound
+for this arithmetic tree, not a numerical tolerance chosen by the caller.
+
 ## Floating-point factorizations
 
 ### LU with partial pivoting