From 973e662bfa1bb85631f6826034f128b181a326c4 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Fri, 4 Sep 2026 09:07:36 -0700 Subject: [PATCH 1/2] feat(vector): add certified dot and affine error bounds MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Add proof-bearing scalar certificates for dot products and unrounded axis · (left - right) reductions. - Expose outward bounds for sign and threshold filtering, with inconclusive results when proof conditions fail. - Document and benchmark the deterministic FMA error model. Closes #220 --- Cargo.lock | 8 +- Cargo.toml | 5 + README.md | 77 ++++- REFERENCES.md | 13 + benches/linear_form.rs | 110 +++++++ docs/BENCHMARKING.md | 11 +- docs/mathematical_basis.md | 54 +++- justfile | 5 + src/error.rs | 7 + src/interval.rs | 2 +- src/lib.rs | 40 ++- src/vector.rs | 571 +++++++++++++++++++++++++++++++++++++ tests/prelude_exports.rs | 6 + tests/proptest_exact.rs | 130 +++++++++ 14 files changed, 1021 insertions(+), 18 deletions(-) create mode 100644 benches/linear_form.rs diff --git a/Cargo.lock b/Cargo.lock index 29c0ba4..98a4ff5 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -108,9 +108,9 @@ checksum = "37b2a672a2cb129a2e41c10b1224bb368f9f37a2b16b612598138befd7b37eb5" [[package]] name = "cc" -version = "1.4.4" +version = "1.4.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0ad534f4357a5264cce5019c989cf66a4f0dc4e0d1b1d15f8aacec0ff7360273" +checksum = "005ec2760ca554fae18df7a11195552ec576cd665632a881bc011d5bb2fd4d80" dependencies = [ "find-msvc-tools", "shlex", @@ -371,9 +371,9 @@ checksum = "da7c62ceae207dd37ea5b845da6a0696c799f85e97da1ab5b7910be3c1c80223" [[package]] name = "find-msvc-tools" -version = "0.1.11" +version = "0.1.12" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d45db016d36b838f563236e9193d0ee6ce38f3f68b6c94e914b4929c96bbb890" +checksum = "3e0f1c7c3a72c66fd80abe965175f7523475c0489a87d3ff9d6e8c87d87a9d2d" [[package]] name = "fnv" diff --git a/Cargo.toml b/Cargo.toml index 1096b86..a07212f 100644 --- a/Cargo.toml +++ b/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 diff --git a/README.md b/README.md index 4592779..d4ab9a5 100644 --- a/README.md +++ b/README.md @@ -25,6 +25,8 @@ while keeping the API intentionally small and explicit. - `Matrix` for fixed-size square `f64` matrices backed by `[[f64; D]; D]` - `Interval` and `IntervalMatrix` 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` and `RationalMatrix` for exact rational inputs behind the optional `"exact"` feature - `Lu` 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( + axis: &Vector, + left: &Vector, + right: &Vector, + threshold: f64, +) -> Result, 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,7 +698,7 @@ out of the common prelude. | Type | Storage | Purpose | Key methods | |---|---|---|---| -| `Vector` | `[f64; D]` | Finite fixed-length vector for input and computation | `try_new`, `as_array`, `into_array`, `dot`, `norm2_sq` | +| `Vector` | `[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` | `[[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` | `[[Interval; D]; D]` | Division-free determinant enclosure and sign proof through D=7 | `from_rows`, `try_from_point_rows`, `from_matrix`, `det`, `det_sign` | @@ -643,6 +706,7 @@ out of the common prelude. | `RationalVector`¹ | `[BigRational; D]` | Exact rational right-hand side and solution | `try_new`, `try_from_fn`, `as_array`, `into_array`, `get` | | `RationalMatrix`¹ | `[[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` | Inline factors + permutation | Factorization for solves/det | `solve`, `det` | | `Ldlt` | 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. + | D | la-stack median (ns) | nalgebra median (ns) | faer median (ns) | reduction vs nalgebra (point est.) | reduction vs faer (point est.) | diff --git a/REFERENCES.md b/REFERENCES.md index 8cffdee..6b1bbce 100644 --- a/REFERENCES.md +++ b/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 diff --git a/benches/linear_form.rs b/benches/linear_form.rs new file mode 100644 index 0000000..46fe751 --- /dev/null +++ b/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(); +} diff --git a/docs/BENCHMARKING.md b/docs/BENCHMARKING.md index b1abaa2..1bf0603 100644 --- a/docs/BENCHMARKING.md +++ b/docs/BENCHMARKING.md @@ -29,6 +29,7 @@ the commands measure and where their outputs go. | Fast saved-baseline loop | `just bench-save-baseline ` then `just bench-compare 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 diff --git a/docs/mathematical_basis.md b/docs/mathematical_basis.md index c20259a..7d9dadc 100644 --- a/docs/mathematical_basis.md +++ b/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 diff --git a/justfile b/justfile index 60bc919..679e5e1 100644 --- a/justfile +++ b/justfile @@ -274,6 +274,10 @@ bench-exact: bench-interval: cargo bench --locked --features bench --bench interval +# Run the certified dot-product and affine-difference benchmark suite. +bench-linear-form: + cargo bench --locked --features bench --bench linear_form + # Run the cheaper latest measurements used for latest-vs-last reports. bench-latest: bench-vs-linalg-la-stack bench-exact @@ -506,6 +510,7 @@ help-workflows: @echo " just bench-latest-vs-last # Run latest and compare against last" @echo " just bench-exact # Run exact-arithmetic benchmarks" @echo " just bench-interval # Run interval determinant benchmarks" + @echo " just bench-linear-form # Run certified linear-form benchmarks" @echo " just bench-save-last # Save full baseline as 'last'" @echo " just bench-vs-linalg # Run vs_linalg bench (optional filter)" @echo " just bench-vs-linalg-la-stack # Run la-stack rows from vs_linalg" diff --git a/src/error.rs b/src/error.rs index 803c030..c299c6e 100644 --- a/src/error.rs +++ b/src/error.rs @@ -47,6 +47,8 @@ pub enum ArithmeticOperation { IntervalDeterminant, /// Vector dot-product calculation. VectorDotProduct, + /// Dot product with an unrounded vector difference. + VectorDotDifference, /// Vector squared-norm calculation. VectorSquaredNorm, } @@ -68,6 +70,7 @@ impl fmt::Display for ArithmeticOperation { Self::IntervalSquare => "interval square", Self::IntervalDeterminant => "interval determinant", Self::VectorDotProduct => "vector dot product", + Self::VectorDotDifference => "vector dot difference", Self::VectorSquaredNorm => "vector squared norm", }) } @@ -904,6 +907,10 @@ mod tests { ArithmeticOperation::VectorDotProduct.to_string(), "vector dot product" ); + assert_eq!( + ArithmeticOperation::VectorDotDifference.to_string(), + "vector dot difference" + ); assert_eq!( ArithmeticOperation::VectorSquaredNorm.to_string(), "vector squared norm" diff --git a/src/interval.rs b/src/interval.rs index 302e2cb..e533d0d 100644 --- a/src/interval.rs +++ b/src/interval.rs @@ -170,7 +170,7 @@ const fn compare_binary_magnitudes( /// Compare the exact-real product `left × right` with its rounded result. #[inline] -const fn compare_product_with_rounded(left: f64, right: f64, rounded: f64) -> i8 { +pub(crate) const fn compare_product_with_rounded(left: f64, right: f64, rounded: f64) -> i8 { let negative = left.is_sign_negative() != right.is_sign_negative(); if rounded == 0.0 { return if negative { -1 } else { 1 }; diff --git a/src/lib.rs b/src/lib.rs index 7bd3398..e881d92 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -111,6 +111,37 @@ mod readme_doctests { /// ``` fn interval_determinant_example() {} + /// ```rust + /// use la_stack::prelude::*; + /// + /// fn is_separated( + /// axis: &Vector, + /// left: &Vector, + /// right: &Vector, + /// threshold: f64, + /// ) -> Result, 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(()) + /// # } + /// ``` + fn certified_linear_form_example() {} + #[cfg(feature = "exact")] /// ```rust /// use la_stack::prelude::*; @@ -479,7 +510,7 @@ pub use ldlt::Ldlt; pub use lu::Lu; pub use matrix::{DeterminantWithErrorBound, Matrix}; pub use tolerance::{DEFAULT_SINGULAR_TOL, Tolerance}; -pub use vector::Vector; +pub use vector::{ScalarWithErrorBound, Vector}; /// Fallibly dispatch a runtime dimension to a concrete stack-allocated matrix. /// @@ -735,7 +766,8 @@ macro_rules! try_with_rational_matrix { /// /// This prelude re-exports the primary types and common constants: [`Matrix`], /// [`DeterminantWithErrorBound`], [`Interval`], [`IntervalMatrix`], -/// [`IntervalDeterminantSign`], [`Vector`], [`Lu`], [`Ldlt`], [`Tolerance`], +/// [`IntervalDeterminantSign`], [`ScalarWithErrorBound`], [`Vector`], [`Lu`], +/// [`Ldlt`], [`Tolerance`], /// and [`LaError`]. Its typed /// error categories include [`ArithmeticOperation`], [`FactorizationKind`], /// [`IntervalBound`], [`IntervalOperand`], [`InvalidToleranceReason`], @@ -799,8 +831,8 @@ pub mod prelude { Interval, IntervalBound, IntervalDeterminantSign, IntervalMatrix, IntervalOperand, InvalidToleranceReason, LaError, Ldlt, Lu, MAX_INTERVAL_MATRIX_DIM, MAX_STACK_MATRIX_DISPATCH_DIM, Matrix, NonFiniteLocation, NonFiniteOrigin, - PositiveSemidefiniteViolation, SingularityReason, Tolerance, UnrepresentableReason, Vector, - try_with_interval_matrix, try_with_stack_matrix, + PositiveSemidefiniteViolation, ScalarWithErrorBound, SingularityReason, Tolerance, + UnrepresentableReason, Vector, try_with_interval_matrix, try_with_stack_matrix, }; #[cfg(feature = "exact")] diff --git a/src/vector.rs b/src/vector.rs index 5db4bdd..3fe1c4f 100644 --- a/src/vector.rs +++ b/src/vector.rs @@ -4,8 +4,267 @@ use core::hint::cold_path; +use crate::interval::compare_product_with_rounded; use crate::{ArithmeticOperation, LaError}; +/// A scalar estimate paired with a certified absolute error bound. +/// +/// Values of this type are produced by [`Vector::dot_with_errbound`] and +/// [`Vector::dot_difference_with_errbound`]. The exact-real value of the +/// corresponding expression over the stored binary64 inputs lies between +/// [`lower_bound`](Self::lower_bound) and [`upper_bound`](Self::upper_bound), +/// and differs from [`estimate`](Self::estimate) by at most +/// [`absolute_error_bound`](Self::absolute_error_bound). +/// +/// The bound certifies floating-point roundoff in one specified arithmetic +/// tree. It is not a caller-selected numerical tolerance and does not classify +/// an interval containing zero as equality. +#[must_use] +#[non_exhaustive] +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct ScalarWithErrorBound { + estimate: f64, + absolute_error_bound: f64, + lower_bound: f64, + upper_bound: f64, +} + +impl ScalarWithErrorBound { + /// Return the rounded scalar estimate. + #[inline] + #[must_use] + pub const fn estimate(self) -> f64 { + self.estimate + } + + /// Return the certified absolute error bound. + #[inline] + #[must_use] + pub const fn absolute_error_bound(self) -> f64 { + self.absolute_error_bound + } + + /// Return a finite outward-rounded lower bound on the exact-real value. + #[inline] + #[must_use] + pub const fn lower_bound(self) -> f64 { + self.lower_bound + } + + /// Return a finite outward-rounded upper bound on the exact-real value. + #[inline] + #[must_use] + pub const fn upper_bound(self) -> f64 { + self.upper_bound + } + + /// Construct a certificate only when both outward endpoints remain finite. + const fn try_new(estimate: f64, absolute_error_bound: f64) -> Option { + if !estimate.is_finite() || !absolute_error_bound.is_finite() || absolute_error_bound < 0.0 + { + return None; + } + + if absolute_error_bound == 0.0 { + return Some(Self { + estimate, + absolute_error_bound: 0.0, + lower_bound: estimate, + upper_bound: estimate, + }); + } + + let lower_rounded = estimate - absolute_error_bound; + let upper_rounded = estimate + absolute_error_bound; + if !lower_rounded.is_finite() || !upper_rounded.is_finite() { + return None; + } + + let lower_error = two_sum_error(estimate, -absolute_error_bound, lower_rounded); + let upper_error = two_sum_error(estimate, absolute_error_bound, upper_rounded); + if !lower_error.is_finite() || !upper_error.is_finite() { + return None; + } + + let lower_bound = if lower_error < 0.0 { + lower_rounded.next_down() + } else { + lower_rounded + }; + let upper_bound = if upper_error > 0.0 { + upper_rounded.next_up() + } else { + upper_rounded + }; + if !lower_bound.is_finite() || !upper_bound.is_finite() { + return None; + } + + Some(Self { + estimate, + absolute_error_bound, + lower_bound, + upper_bound, + }) + } +} + +/// Return the exact residual of a finite rounded binary64 addition. +/// +/// This is Knuth's `TwoSum` transform. It is used only to round the published +/// certificate endpoints outward; it is independent of the reduction bound. +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 +} + +/// State for one certified left-to-right FMA reduction. +#[derive(Clone, Copy, Debug, PartialEq)] +struct CertifiedReduction { + estimate: f64, + magnitude_upper: f64, + proof_available: bool, +} + +impl CertifiedReduction { + const ZERO: Self = Self { + estimate: 0.0, + magnitude_upper: 0.0, + proof_available: true, + }; + + /// Add one exact-real product through one rounded FMA. + const fn add_product( + mut self, + left: f64, + right: f64, + operation: ArithmeticOperation, + index: usize, + ) -> Result { + let prior = self.estimate; + let estimate = left.mul_add(right, prior); + if !estimate.is_finite() { + cold_path(); + return Err(LaError::non_finite_computation_step(operation, index)); + } + + if self.proof_available { + self.proof_available = estimate.is_normal() + || (estimate == 0.0 && Self::fma_result_is_exact_zero(left, right, prior)); + } + if self.proof_available { + match Self::add_product_magnitude_upper(self.magnitude_upper, left, right) { + Some(magnitude_upper) => self.magnitude_upper = magnitude_upper, + None => self.proof_available = false, + } + } + self.estimate = estimate; + Ok(self) + } + + /// Return whether `left × right + addend` is exactly zero. + const fn fma_result_is_exact_zero(left: f64, right: f64, addend: f64) -> bool { + if left == 0.0 || right == 0.0 { + return addend == 0.0; + } + + let rounded_product = left * right; + let rounded_bits = rounded_product.to_bits(); + let negated_addend_bits = (-addend).to_bits(); + let same_rounded_value = rounded_bits == negated_addend_bits + || (rounded_bits << 1 == 0 && negated_addend_bits << 1 == 0); + rounded_product.is_finite() + && same_rounded_value + && compare_product_with_rounded(left, right, rounded_product) == 0 + } + + /// Add an upward-rounded bound on `|left × right|` to the magnitude sum. + const fn add_product_magnitude_upper( + magnitude_upper: f64, + left: f64, + right: f64, + ) -> Option { + if left == 0.0 || right == 0.0 { + return Some(magnitude_upper); + } + + let left_magnitude = left.abs(); + let right_magnitude = right.abs(); + let rounded_product = left_magnitude * right_magnitude; + if !rounded_product.is_normal() { + return None; + } + + let product_upper = + if compare_product_with_rounded(left_magnitude, right_magnitude, rounded_product) > 0 { + rounded_product.next_up() + } else { + rounded_product + }; + if !product_upper.is_finite() { + return None; + } + + if magnitude_upper == 0.0 { + return Some(product_upper); + } + let rounded_sum = magnitude_upper + product_upper; + if !rounded_sum.is_finite() { + return None; + } + let sum_upper = rounded_sum.next_up(); + if sum_upper.is_finite() { + Some(sum_upper) + } else { + None + } + } + + /// Finish the reduction with an upward-rounded `gamma_n` error bound. + #[expect( + clippy::cast_precision_loss, + reason = "a usable gamma requires a term count below 2^53, where the cast is exact" + )] + const fn finish(self, term_count: Option) -> Option { + if !self.proof_available { + return None; + } + if self.magnitude_upper == 0.0 { + return ScalarWithErrorBound::try_new(self.estimate, 0.0); + } + + let Some(term_count) = term_count else { + return None; + }; + let scaled_roundoff = (term_count as f64) * (f64::EPSILON / 2.0); + if !scaled_roundoff.is_finite() || scaled_roundoff >= 1.0 { + return None; + } + + // The count conversion, multiplication by 2^-53, and subtraction from + // one are exact throughout the usable range. Round the division and + // final multiplication upward to retain a certified upper bound. + let gamma = scaled_roundoff / (1.0 - scaled_roundoff); + let gamma_upper = gamma.next_up(); + if !gamma_upper.is_finite() { + return None; + } + let rounded_bound = gamma_upper * self.magnitude_upper; + if !rounded_bound.is_finite() { + return None; + } + let absolute_error_bound = if rounded_bound == 0.0 { + 0.0 + } else { + rounded_bound.next_up() + }; + ScalarWithErrorBound::try_new(self.estimate, absolute_error_bound) + } +} + /// Finite fixed-size vector of length `D`, stored inline. /// /// Public construction rejects NaN and infinity through [`try_new`](Self::try_new), @@ -180,6 +439,145 @@ impl Vector { self.dot_with_operation(other, ArithmeticOperation::VectorDotProduct) } + /// Dot product with a certified absolute roundoff bound. + /// + /// The estimate uses the deterministic left-to-right recurrence + /// `s[0] = 0` and `s[i + 1] = self[i].mul_add(other[i], s[i])`. When the + /// relative-error model is valid, the returned certificate bounds the + /// difference between `s[D]` and the exact-real expression + /// `Σᵢ self[i] × other[i]` over the stored binary64 inputs. + /// + /// The bound is `gamma_D × Σᵢ |self[i] × other[i]|`, where + /// `gamma_D = D u / (1 - D u)` and `u = 2^-53`. The magnitude sum and the + /// published bound are rounded upward. See `REFERENCES.md` \[9-11\]. + /// + /// `Ok(None)` means no certificate is available because a nonzero product + /// or FMA result entered the subnormal range, the reduction dimension made + /// `gamma_D` invalid, or a proof-only magnitude/bound calculation exhausted + /// the finite binary64 range. It does not mean the exact dot product is + /// zero. Unlike a user-selected tolerance, a returned error bound describes + /// rounding in this specific arithmetic tree. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// # 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 certificate = left + /// .dot_with_errbound(&right)? + /// .expect("ordinary inputs have a binary64 certificate"); + /// assert_eq!(certificate.estimate(), 32.0); + /// assert!(certificate.lower_bound() > 0.0); + /// # Ok(()) + /// # } + /// ``` + /// + /// # Errors + /// Returns [`LaError::NonFinite`] with the first failing reduction index + /// when the FMA estimate overflows to NaN or infinity. + #[inline] + pub const fn dot_with_errbound( + &self, + other: &Self, + ) -> Result, LaError> { + let left = self.as_array(); + let right = other.as_array(); + let mut reduction = CertifiedReduction::ZERO; + let mut i = 0; + while i < D { + reduction = match reduction.add_product( + left[i], + right[i], + ArithmeticOperation::VectorDotProduct, + i, + ) { + Ok(reduction) => reduction, + Err(error) => return Err(error), + }; + i += 1; + } + Ok(reduction.finish(Some(D))) + } + + /// Certified dot product with an unrounded vector difference. + /// + /// This evaluates the exact-real expression + /// `Σᵢ self[i] × (left[i] - right[i])` without first rounding + /// `left - right` into a [`Vector`]. Its deterministic arithmetic tree is + /// + /// ```text + /// s[0] = 0 + /// s[2i + 1] = self[i].mul_add(left[i], s[2i]) + /// s[2i + 2] = (-self[i]).mul_add(right[i], s[2i + 1]). + /// ``` + /// + /// A returned certificate therefore includes all `2D` FMA rounding events + /// in that tree and bounds the intended expression over the original + /// binary64 coordinates. Its lower and upper endpoints can certify a sign + /// or separation from a caller's threshold. An overlapping endpoint range + /// remains inconclusive and should trigger the caller's exact fallback. + /// + /// `Ok(None)` has the same proof-unavailable meaning as in + /// [`dot_with_errbound`](Self::dot_with_errbound), including gradual + /// underflow and proof-only range exhaustion. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// # 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])?; + /// let certificate = axis + /// .dot_difference_with_errbound(&left, &right)? + /// .expect("ordinary inputs have a binary64 certificate"); + /// assert_eq!(certificate.estimate(), 8.0); + /// assert!(certificate.lower_bound() > 1.0); + /// # Ok(()) + /// # } + /// ``` + /// + /// # Errors + /// Returns [`LaError::NonFinite`] with the first failing coordinate index + /// when an FMA estimate overflows to NaN or infinity. + #[inline] + pub const fn dot_difference_with_errbound( + &self, + left: &Self, + right: &Self, + ) -> Result, LaError> { + let axis = self.as_array(); + let left = left.as_array(); + let right = right.as_array(); + let mut reduction = CertifiedReduction::ZERO; + let mut i = 0; + while i < D { + reduction = match reduction.add_product( + axis[i], + left[i], + ArithmeticOperation::VectorDotDifference, + i, + ) { + Ok(reduction) => reduction, + Err(error) => return Err(error), + }; + reduction = match reduction.add_product( + -axis[i], + right[i], + ArithmeticOperation::VectorDotDifference, + i, + ) { + Ok(reduction) => reduction, + Err(error) => return Err(error), + }; + i += 1; + } + Ok(reduction.finish(D.checked_mul(2))) + } + /// Accumulate a dot product while retaining the public operation that owns it. const fn dot_with_operation( &self, @@ -379,6 +777,38 @@ mod tests { ); } + #[test] + fn []() { + let mut left_data = [0.0; $d]; + let mut right_data = [0.0; $d]; + let left_values = [1.0, 2.0, 3.0, 4.0, 5.0]; + let right_values = [2.0, 3.0, 4.0, 5.0, 6.0]; + for (destination, source) in left_data.iter_mut().zip(left_values) { + *destination = source; + } + for (destination, source) in right_data.iter_mut().zip(right_values) { + *destination = source; + } + let left = Vector::<$d>::new(left_data); + let right = Vector::<$d>::new(right_data); + let zero = Vector::<$d>::zero(); + + let dot = left.dot(&right).unwrap(); + let dot_certificate = left.dot_with_errbound(&right).unwrap().unwrap(); + assert_abs_diff_eq!(dot_certificate.estimate(), dot, epsilon = 0.0); + assert!(dot_certificate.absolute_error_bound() >= 0.0); + assert!(dot_certificate.lower_bound() <= dot); + assert!(dot <= dot_certificate.upper_bound()); + + let difference_certificate = left + .dot_difference_with_errbound(&right, &zero) + .unwrap() + .unwrap(); + assert_abs_diff_eq!(difference_certificate.estimate(), dot, epsilon = 0.0); + assert!(difference_certificate.lower_bound() <= dot); + assert!(dot <= difference_certificate.upper_bound()); + } + #[test] fn []() { for value in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] { @@ -433,6 +863,20 @@ mod tests { 0, )) ); + assert_eq!( + a.dot_with_errbound(&b), + Err(LaError::non_finite_computation_step( + ArithmeticOperation::VectorDotProduct, + 0, + )) + ); + assert_eq!( + a.dot_difference_with_errbound(&b, &Vector::zero()), + Err(LaError::non_finite_computation_step( + ArithmeticOperation::VectorDotDifference, + 0, + )) + ); assert_eq!( a.norm2_sq(), Err(LaError::non_finite_computation_step( @@ -554,6 +998,122 @@ mod tests { assert_eq!(vector.norm2_sq(), Ok(norm_large * norm_large)); } + #[test] + fn certified_dot_preserves_fma_estimate_and_withholds_range_exhausted_bound() { + let left = Vector::<2>::new([f64::MAX, f64::MAX]); + let right = Vector::<2>::new([-1.0, 2.0]); + + assert_eq!(left.dot(&right), Ok(f64::MAX)); + assert_eq!(left.dot_with_errbound(&right), Ok(None)); + } + + #[test] + fn certified_bounds_distinguish_conclusive_and_inconclusive_results() { + let conclusive = Vector::<2>::new([1.0, 2.0]) + .dot_with_errbound(&Vector::new([3.0, 4.0])) + .unwrap() + .unwrap(); + assert_abs_diff_eq!(conclusive.estimate(), 11.0, epsilon = 0.0); + assert!(conclusive.lower_bound() > 1.0); + + let inconclusive = Vector::<2>::new([1.0, 1.0]) + .dot_with_errbound(&Vector::new([1.0, -1.0])) + .unwrap() + .unwrap(); + assert_abs_diff_eq!(inconclusive.estimate(), 0.0, epsilon = 0.0); + assert!(inconclusive.absolute_error_bound() > 0.0); + assert!(inconclusive.lower_bound() < 0.0); + assert!(inconclusive.upper_bound() > 0.0); + } + + #[test] + fn certified_zero_and_signed_zero_have_an_exact_zero_bound() { + let left = Vector::<3>::new([-0.0, 0.0, -0.0]); + let right = Vector::<3>::new([f64::MAX, -1.0, f64::MIN_POSITIVE]); + let certificate = left.dot_with_errbound(&right).unwrap().unwrap(); + + assert_abs_diff_eq!(certificate.estimate(), 0.0, epsilon = 0.0); + assert_abs_diff_eq!(certificate.absolute_error_bound(), 0.0, epsilon = 0.0); + assert_abs_diff_eq!(certificate.lower_bound(), 0.0, epsilon = 0.0); + assert_abs_diff_eq!(certificate.upper_bound(), 0.0, epsilon = 0.0); + } + + #[test] + fn certified_reductions_withhold_bounds_for_subnormal_products() { + let tiny = Vector::<1>::new([f64::MIN_POSITIVE]); + let half = Vector::<1>::new([0.5]); + assert_eq!(tiny.dot_with_errbound(&half), Ok(None)); + + let min_subnormal = Vector::<1>::new([f64::from_bits(1)]); + assert_eq!( + min_subnormal.dot_with_errbound(&Vector::new([1.0])), + Ok(None) + ); + assert_eq!( + tiny.dot_difference_with_errbound(&half, &Vector::zero()), + Ok(None) + ); + } + + #[test] + fn certified_reduction_detects_fma_cancellation_below_subnormal_range() { + let near_sqrt_min = f64::from_bits((512_u64 << 52) | 1); + let rounded_product = near_sqrt_min * near_sqrt_min; + assert!(rounded_product.is_normal()); + assert_abs_diff_eq!( + near_sqrt_min.mul_add(near_sqrt_min, -rounded_product), + 0.0, + epsilon = 0.0 + ); + + let left = Vector::<2>::new([-rounded_product, near_sqrt_min]); + let right = Vector::<2>::new([1.0, near_sqrt_min]); + assert_eq!(left.dot(&right), Ok(0.0)); + assert_eq!(left.dot_with_errbound(&right), Ok(None)); + } + + #[test] + fn certified_dot_handles_mixed_normal_magnitudes() { + let left = Vector::<2>::new([1.0e100, 1.0e-100]); + let right = Vector::<2>::new([1.0e-100, 1.0e100]); + let certificate = left.dot_with_errbound(&right).unwrap().unwrap(); + + assert_abs_diff_eq!(certificate.estimate(), 2.0, epsilon = 0.0); + assert!(certificate.lower_bound() <= 2.0); + assert!(certificate.upper_bound() >= 2.0); + } + + #[test] + fn certified_dot_difference_does_not_round_coordinates_first() { + let scale = 18_014_398_509_481_984.0; + let axis = Vector::<1>::new([scale]); + let left = Vector::<1>::new([1.0]); + let right = Vector::<1>::new([1.0 / scale]); + assert_abs_diff_eq!(left.as_array()[0] - right.as_array()[0], 1.0, epsilon = 0.0); + + let certificate = axis + .dot_difference_with_errbound(&left, &right) + .unwrap() + .unwrap(); + assert_abs_diff_eq!(certificate.estimate(), scale, epsilon = 0.0); + assert!(certificate.lower_bound() < scale); + assert!(certificate.upper_bound() >= scale); + } + + #[test] + fn certified_reductions_are_const_evaluable() { + const DOT: Result, LaError> = + Vector::<2>::new([1.0, 2.0]).dot_with_errbound(&Vector::<2>::new([3.0, 4.0])); + const DIFFERENCE: Result, LaError> = + Vector::<2>::new([2.0, -1.0]).dot_difference_with_errbound( + &Vector::<2>::new([4.0, 1.0]), + &Vector::<2>::new([1.0, 3.0]), + ); + + assert_abs_diff_eq!(DOT.unwrap().unwrap().estimate(), 11.0, epsilon = 0.0); + assert_abs_diff_eq!(DIFFERENCE.unwrap().unwrap().estimate(), 8.0, epsilon = 0.0); + } + #[test] fn vector_dot_and_norm2_sq_report_first_middle_overflowing_step() { let dot_lhs = Vector::<3>::new([f64::MAX, f64::MAX, 1.0]); @@ -584,6 +1144,17 @@ mod tests { assert!(vector.as_array().is_empty()); assert!(vector.into_array().is_empty()); assert_eq!(vector.dot(&Vector::zero()), Ok(0.0)); + let dot_certificate = vector.dot_with_errbound(&Vector::zero()).unwrap().unwrap(); + assert_abs_diff_eq!(dot_certificate.absolute_error_bound(), 0.0, epsilon = 0.0); + let difference_certificate = vector + .dot_difference_with_errbound(&Vector::zero(), &Vector::zero()) + .unwrap() + .unwrap(); + assert_abs_diff_eq!( + difference_certificate.absolute_error_bound(), + 0.0, + epsilon = 0.0 + ); assert_eq!(vector.norm2_sq(), Ok(0.0)); } } diff --git a/tests/prelude_exports.rs b/tests/prelude_exports.rs index 35b4b8c..422e94e 100644 --- a/tests/prelude_exports.rs +++ b/tests/prelude_exports.rs @@ -22,6 +22,12 @@ fn common_prelude_supports_downstream_composition() -> Result<(), LaError> { }; assert_abs_diff_eq!(estimate.determinant(), 1.0, epsilon = 0.0); assert!(estimate.absolute_error_bound() >= 0.0); + let dot_estimate: ScalarWithErrorBound = vector + .dot_with_errbound(&vector)? + .expect("ordinary vector inputs must have a certified dot bound"); + assert_abs_diff_eq!(dot_estimate.estimate(), 5.0, epsilon = 0.0); + assert!(dot_estimate.lower_bound() <= 5.0); + assert!(dot_estimate.upper_bound() >= 5.0); let lu: Lu<2> = matrix.lu(tolerance)?; let ldlt: Ldlt<2> = matrix.ldlt(tolerance)?; diff --git a/tests/proptest_exact.rs b/tests/proptest_exact.rs index b4ee7b0..a3fc1ae 100644 --- a/tests/proptest_exact.rs +++ b/tests/proptest_exact.rs @@ -146,6 +146,50 @@ fn big_rational_matvec( }) } +/// Evaluate a dot product over the exact rational values of binary64 inputs. +fn big_rational_dot(left: &[f64; D], right: &[f64; D]) -> BigRational { + let mut sum = BigRational::from_integer(BigInt::from(0)); + for (&left, &right) in left.iter().zip(right) { + let left = BigRational::from_f64(left).expect("finite f64 converts exactly"); + let right = BigRational::from_f64(right).expect("finite f64 converts exactly"); + sum += left * right; + } + sum +} + +/// Evaluate `axis · (left - right)` without rounding coordinate differences. +fn big_rational_dot_difference( + axis: &[f64; D], + left: &[f64; D], + right: &[f64; D], +) -> BigRational { + let mut sum = BigRational::from_integer(BigInt::from(0)); + for ((&axis, &left), &right) in axis.iter().zip(left).zip(right) { + let axis = BigRational::from_f64(axis).expect("finite f64 converts exactly"); + let left = BigRational::from_f64(left).expect("finite f64 converts exactly"); + let right = BigRational::from_f64(right).expect("finite f64 converts exactly"); + sum += axis * (left - right); + } + sum +} + +/// Check both published forms of a scalar certificate against an exact oracle. +fn scalar_certificate_contains_exact( + certificate: ScalarWithErrorBound, + exact: &BigRational, +) -> bool { + let estimate = BigRational::from_f64(certificate.estimate()) + .expect("a scalar certificate has a finite estimate"); + let error_bound = BigRational::from_f64(certificate.absolute_error_bound()) + .expect("a scalar certificate has a finite error bound"); + let lower = BigRational::from_f64(certificate.lower_bound()) + .expect("a scalar certificate has a finite lower bound"); + let upper = BigRational::from_f64(certificate.upper_bound()) + .expect("a scalar certificate has a finite upper bound"); + + (estimate - exact).abs() <= error_bound && lower <= *exact && *exact <= upper +} + /// Compute an exact determinant via the Leibniz permutation expansion. /// /// This is intentionally independent from the production Bareiss core. It is @@ -645,6 +689,92 @@ gen_det_errbound_leibniz_oracle_proptests!(2); gen_det_errbound_leibniz_oracle_proptests!(3); gen_det_errbound_leibniz_oracle_proptests!(4); +/// Certified dot products must enclose an independent exact-rational sum of +/// products over the stored binary64 values. +macro_rules! gen_dot_errbound_oracle_proptests { + ($d:literal) => { + paste! { + proptest! { + #![proptest_config(with_default_cases(64))] + + #[test] + fn []( + left in array::[]( + (-50i16..=50i16).prop_map(|value| f64::from(value) / 10.0) + ), + right in array::[]( + (-50i16..=50i16).prop_map(|value| f64::from(value) / 10.0) + ), + ) { + let left_vector = Vector::<$d>::try_new(left).unwrap(); + let right_vector = Vector::<$d>::try_new(right).unwrap(); + let certificate = left_vector + .dot_with_errbound(&right_vector) + .unwrap() + .expect("moderate inputs stay in the certified range"); + let exact = big_rational_dot(&left, &right); + + prop_assert!( + scalar_certificate_contains_exact(certificate, &exact), + "D={} dot certificate {certificate:?} did not contain {exact}", + $d, + ); + } + } + } + }; +} + +gen_dot_errbound_oracle_proptests!(2); +gen_dot_errbound_oracle_proptests!(3); +gen_dot_errbound_oracle_proptests!(4); +gen_dot_errbound_oracle_proptests!(5); + +/// The affine-difference certificate is checked against the exact expression +/// over the original coordinates, never an already-rounded `left - right`. +macro_rules! gen_dot_difference_errbound_oracle_proptests { + ($d:literal) => { + paste! { + proptest! { + #![proptest_config(with_default_cases(64))] + + #[test] + fn []( + axis in array::[]( + (-20i16..=20i16).prop_map(|value| f64::from(value) / 10.0) + ), + left in array::[]( + (-50i16..=50i16).prop_map(|value| f64::from(value) / 10.0) + ), + right in array::[]( + (-50i16..=50i16).prop_map(|value| f64::from(value) / 10.0) + ), + ) { + let axis_vector = Vector::<$d>::try_new(axis).unwrap(); + let left_vector = Vector::<$d>::try_new(left).unwrap(); + let right_vector = Vector::<$d>::try_new(right).unwrap(); + let certificate = axis_vector + .dot_difference_with_errbound(&left_vector, &right_vector) + .unwrap() + .expect("moderate inputs stay in the certified range"); + let exact = big_rational_dot_difference(&axis, &left, &right); + + prop_assert!( + scalar_certificate_contains_exact(certificate, &exact), + "D={} affine certificate {certificate:?} did not contain {exact}", + $d, + ); + } + } + } + }; +} + +gen_dot_difference_errbound_oracle_proptests!(2); +gen_dot_difference_errbound_oracle_proptests!(3); +gen_dot_difference_errbound_oracle_proptests!(4); +gen_dot_difference_errbound_oracle_proptests!(5); + /// Exercise the determinant certificate with independently mixed per-entry /// exponents spanning zero, subnormal, tiny normal, ordinary, and large finite /// regimes. `det_sign_exact` must always match the independent Leibniz oracle; From d8d21dbf76d7a4fa3ba15f37a6086179e7394976 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Fri, 4 Sep 2026 10:36:50 -0700 Subject: [PATCH 2/2] fix(codeql): avoid false certificate logging alerts - Rename numerical-bound locals so CodeQL does not mistake them for sensitive certificate data. - Document finite bound invariants, the affine error formula, and typed failure contexts. - Make proof-range and second-FMA overflow expectations explicit. --- src/error.rs | 2 +- src/vector.rs | 195 ++++++++++++++++++++++++++++++++++++++------------ 2 files changed, 151 insertions(+), 46 deletions(-) diff --git a/src/error.rs b/src/error.rs index c299c6e..6e40661 100644 --- a/src/error.rs +++ b/src/error.rs @@ -47,7 +47,7 @@ pub enum ArithmeticOperation { IntervalDeterminant, /// Vector dot-product calculation. VectorDotProduct, - /// Dot product with an unrounded vector difference. + /// `axis · (left - right)` without first rounding the vector difference. VectorDotDifference, /// Vector squared-norm calculation. VectorSquaredNorm, diff --git a/src/vector.rs b/src/vector.rs index 3fe1c4f..ce6d18f 100644 --- a/src/vector.rs +++ b/src/vector.rs @@ -19,6 +19,27 @@ use crate::{ArithmeticOperation, LaError}; /// The bound certifies floating-point roundoff in one specified arithmetic /// tree. It is not a caller-selected numerical tolerance and does not classify /// an interval containing zero as equality. +/// +/// Callers cannot construct this type directly. Every value has a finite +/// estimate, a finite non-negative absolute error bound, and finite ordered +/// endpoints. +/// +/// # Examples +/// ``` +/// 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); +/// # Ok(()) +/// # } +/// ``` #[must_use] #[non_exhaustive] #[derive(Clone, Copy, Debug, PartialEq)] @@ -58,7 +79,11 @@ impl ScalarWithErrorBound { self.upper_bound } - /// Construct a certificate only when both outward endpoints remain finite. + /// 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 + /// infinite endpoint enforces the public proof-unavailable contract. const fn try_new(estimate: f64, absolute_error_bound: f64) -> Option { if !estimate.is_finite() || !absolute_error_bound.is_finite() || absolute_error_bound < 0.0 { @@ -112,7 +137,7 @@ impl ScalarWithErrorBound { /// Return the exact residual of a finite rounded binary64 addition. /// /// This is Knuth's `TwoSum` transform. It is used only to round the published -/// certificate endpoints outward; it is independent of the reduction bound. +/// bound endpoints outward; it is independent of the reduction bound. const fn two_sum_error(left: f64, right: f64, rounded: f64) -> f64 { let virtual_right = rounded - left; let virtual_left = rounded - virtual_right; @@ -122,6 +147,12 @@ const fn two_sum_error(left: f64, right: f64, rounded: f64) -> f64 { } /// State for one certified left-to-right FMA reduction. +/// +/// The rounded estimate continues accumulating after `proof_available` becomes +/// false. This lets the public methods distinguish a non-finite estimate +/// (`Err`) from a finite estimate whose proof arithmetic is unavailable +/// (`Ok(None)`). `magnitude_upper` encloses the sum of exact product +/// magnitudes while that proof remains available. #[derive(Clone, Copy, Debug, PartialEq)] struct CertifiedReduction { estimate: f64, @@ -137,6 +168,10 @@ impl CertifiedReduction { }; /// Add one exact-real product through one rounded FMA. + /// + /// A non-finite estimate becomes a typed public error. Underflow or range + /// loss confined to proof construction instead clears `proof_available`, + /// preserving the finite estimate for the eventual `Ok(None)` result. const fn add_product( mut self, left: f64, @@ -166,6 +201,10 @@ impl CertifiedReduction { } /// Return whether `left × right + addend` is exactly zero. + /// + /// This distinguishes exact cancellation from a nonzero value rounded to + /// zero. Only exact zero is admissible under the relative-error model used + /// by the public certified reductions. const fn fma_result_is_exact_zero(left: f64, right: f64, addend: f64) -> bool { if left == 0.0 || right == 0.0 { return addend == 0.0; @@ -182,6 +221,10 @@ impl CertifiedReduction { } /// Add an upward-rounded bound on `|left × right|` to the magnitude sum. + /// + /// `None` means a nonzero product was not normal or the product/sum could + /// not be enclosed by a finite binary64 value, so the public result must be + /// proof-unavailable. const fn add_product_magnitude_upper( magnitude_upper: f64, left: f64, @@ -224,6 +267,11 @@ impl CertifiedReduction { } /// Finish the reduction with an upward-rounded `gamma_n` error bound. + /// + /// Returns `None` when an earlier proof step failed, the term count cannot + /// support the relative-error model, or the final bound/endpoints cannot + /// remain finite. Otherwise the result satisfies every public + /// [`ScalarWithErrorBound`] invariant. #[expect( clippy::cast_precision_loss, reason = "a usable gamma requires a term count below 2^53, where the cast is exact" @@ -465,18 +513,19 @@ impl Vector { /// # 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 certificate = left + /// let bounded = left /// .dot_with_errbound(&right)? - /// .expect("ordinary inputs have a binary64 certificate"); - /// assert_eq!(certificate.estimate(), 32.0); - /// assert!(certificate.lower_bound() > 0.0); + /// .expect("ordinary inputs have a binary64 bound"); + /// assert_eq!(bounded.estimate(), 32.0); + /// assert!(bounded.lower_bound() > 0.0); /// # Ok(()) /// # } /// ``` /// /// # Errors - /// Returns [`LaError::NonFinite`] with the first failing reduction index - /// when the FMA estimate overflows to NaN or infinity. + /// Returns [`LaError::NonFinite`] with the first failing reduction index and + /// [`ArithmeticOperation::VectorDotProduct`] when an FMA estimate becomes + /// non-finite. #[inline] pub const fn dot_with_errbound( &self, @@ -515,8 +564,13 @@ impl Vector { /// /// A returned certificate therefore includes all `2D` FMA rounding events /// in that tree and bounds the intended expression over the original - /// binary64 coordinates. Its lower and upper endpoints can certify a sign - /// or separation from a caller's threshold. An overlapping endpoint range + /// binary64 coordinates. When available, its absolute bound is + /// `gamma_2D × Σᵢ (|self[i] × left[i]| + |self[i] × right[i]|)`, where + /// `gamma_2D = 2D u / (1 - 2D u)` and `u = 2^-53`; every magnitude and the + /// final bound are rounded upward. Its + /// [`lower_bound`](ScalarWithErrorBound::lower_bound) and + /// [`upper_bound`](ScalarWithErrorBound::upper_bound) can certify a sign or + /// separation from a caller's threshold. An overlapping endpoint range /// remains inconclusive and should trigger the caller's exact fallback. /// /// `Ok(None)` has the same proof-unavailable meaning as in @@ -531,18 +585,19 @@ impl Vector { /// 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 certificate = axis + /// let bounded = axis /// .dot_difference_with_errbound(&left, &right)? - /// .expect("ordinary inputs have a binary64 certificate"); - /// assert_eq!(certificate.estimate(), 8.0); - /// assert!(certificate.lower_bound() > 1.0); + /// .expect("ordinary inputs have a binary64 bound"); + /// assert_eq!(bounded.estimate(), 8.0); + /// assert!(bounded.lower_bound() > 1.0); /// # Ok(()) /// # } /// ``` /// /// # Errors /// Returns [`LaError::NonFinite`] with the first failing coordinate index - /// when an FMA estimate overflows to NaN or infinity. + /// and [`ArithmeticOperation::VectorDotDifference`] when either FMA for + /// that coordinate produces a non-finite estimate. #[inline] pub const fn dot_difference_with_errbound( &self, @@ -794,19 +849,19 @@ mod tests { let zero = Vector::<$d>::zero(); let dot = left.dot(&right).unwrap(); - let dot_certificate = left.dot_with_errbound(&right).unwrap().unwrap(); - assert_abs_diff_eq!(dot_certificate.estimate(), dot, epsilon = 0.0); - assert!(dot_certificate.absolute_error_bound() >= 0.0); - assert!(dot_certificate.lower_bound() <= dot); - assert!(dot <= dot_certificate.upper_bound()); + let dot_bound = left.dot_with_errbound(&right).unwrap().unwrap(); + assert_abs_diff_eq!(dot_bound.estimate(), dot, epsilon = 0.0); + assert!(dot_bound.absolute_error_bound() >= 0.0); + assert!(dot_bound.lower_bound() <= dot); + assert!(dot <= dot_bound.upper_bound()); - let difference_certificate = left + let difference_bound = left .dot_difference_with_errbound(&right, &zero) .unwrap() .unwrap(); - assert_abs_diff_eq!(difference_certificate.estimate(), dot, epsilon = 0.0); - assert!(difference_certificate.lower_bound() <= dot); - assert!(dot <= difference_certificate.upper_bound()); + assert_abs_diff_eq!(difference_bound.estimate(), dot, epsilon = 0.0); + assert!(difference_bound.lower_bound() <= dot); + assert!(dot <= difference_bound.upper_bound()); } #[test] @@ -1007,6 +1062,29 @@ mod tests { assert_eq!(left.dot_with_errbound(&right), Ok(None)); } + #[test] + fn certified_dot_withholds_bound_when_finite_endpoints_cannot_be_published() { + let maximum = Vector::<1>::new([f64::MAX]); + let one = Vector::<1>::new([1.0]); + + assert_eq!(maximum.dot(&one), Ok(f64::MAX)); + assert_eq!(maximum.dot_with_errbound(&one), Ok(None)); + } + + #[test] + fn certified_dot_withholds_bound_when_magnitude_sum_exhausts_range() { + let maximum = Vector::<2>::new([f64::MAX, f64::MAX]); + + for factor in [0.5, 0.75] { + let cancelling = Vector::<2>::new([factor, -factor]); + let estimate = maximum + .dot(&cancelling) + .expect("the cancelling FMA estimate must remain finite"); + assert!(estimate.is_finite()); + assert_eq!(maximum.dot_with_errbound(&cancelling), Ok(None)); + } + } + #[test] fn certified_bounds_distinguish_conclusive_and_inconclusive_results() { let conclusive = Vector::<2>::new([1.0, 2.0]) @@ -1030,12 +1108,12 @@ mod tests { fn certified_zero_and_signed_zero_have_an_exact_zero_bound() { let left = Vector::<3>::new([-0.0, 0.0, -0.0]); let right = Vector::<3>::new([f64::MAX, -1.0, f64::MIN_POSITIVE]); - let certificate = left.dot_with_errbound(&right).unwrap().unwrap(); + let bounded = left.dot_with_errbound(&right).unwrap().unwrap(); - assert_abs_diff_eq!(certificate.estimate(), 0.0, epsilon = 0.0); - assert_abs_diff_eq!(certificate.absolute_error_bound(), 0.0, epsilon = 0.0); - assert_abs_diff_eq!(certificate.lower_bound(), 0.0, epsilon = 0.0); - assert_abs_diff_eq!(certificate.upper_bound(), 0.0, epsilon = 0.0); + assert_abs_diff_eq!(bounded.estimate(), 0.0, epsilon = 0.0); + assert_abs_diff_eq!(bounded.absolute_error_bound(), 0.0, epsilon = 0.0); + assert_abs_diff_eq!(bounded.lower_bound(), 0.0, epsilon = 0.0); + assert_abs_diff_eq!(bounded.upper_bound(), 0.0, epsilon = 0.0); } #[test] @@ -1049,6 +1127,10 @@ mod tests { min_subnormal.dot_with_errbound(&Vector::new([1.0])), Ok(None) ); + assert_eq!( + min_subnormal.dot_with_errbound(&Vector::new([0.5])), + Ok(None) + ); assert_eq!( tiny.dot_difference_with_errbound(&half, &Vector::zero()), Ok(None) @@ -1076,11 +1158,11 @@ mod tests { fn certified_dot_handles_mixed_normal_magnitudes() { let left = Vector::<2>::new([1.0e100, 1.0e-100]); let right = Vector::<2>::new([1.0e-100, 1.0e100]); - let certificate = left.dot_with_errbound(&right).unwrap().unwrap(); + let bounded = left.dot_with_errbound(&right).unwrap().unwrap(); - assert_abs_diff_eq!(certificate.estimate(), 2.0, epsilon = 0.0); - assert!(certificate.lower_bound() <= 2.0); - assert!(certificate.upper_bound() >= 2.0); + assert_abs_diff_eq!(bounded.estimate(), 2.0, epsilon = 0.0); + assert!(bounded.lower_bound() <= 2.0); + assert!(bounded.upper_bound() >= 2.0); } #[test] @@ -1091,13 +1173,13 @@ mod tests { let right = Vector::<1>::new([1.0 / scale]); assert_abs_diff_eq!(left.as_array()[0] - right.as_array()[0], 1.0, epsilon = 0.0); - let certificate = axis + let bounded = axis .dot_difference_with_errbound(&left, &right) .unwrap() .unwrap(); - assert_abs_diff_eq!(certificate.estimate(), scale, epsilon = 0.0); - assert!(certificate.lower_bound() < scale); - assert!(certificate.upper_bound() >= scale); + assert_abs_diff_eq!(bounded.estimate(), scale, epsilon = 0.0); + assert!(bounded.lower_bound() < scale); + assert!(bounded.upper_bound() >= scale); } #[test] @@ -1114,6 +1196,21 @@ mod tests { assert_abs_diff_eq!(DIFFERENCE.unwrap().unwrap().estimate(), 8.0, epsilon = 0.0); } + #[test] + fn certified_dot_difference_reports_second_fma_overflow() { + let axis = Vector::<2>::new([1.0, f64::MAX]); + let left = Vector::<2>::new([0.0, 1.0]); + let right = Vector::<2>::new([0.0, -1.0]); + + assert_eq!( + axis.dot_difference_with_errbound(&left, &right), + Err(LaError::non_finite_computation_step( + ArithmeticOperation::VectorDotDifference, + 1, + )) + ); + } + #[test] fn vector_dot_and_norm2_sq_report_first_middle_overflowing_step() { let dot_lhs = Vector::<3>::new([f64::MAX, f64::MAX, 1.0]); @@ -1144,17 +1241,25 @@ mod tests { assert!(vector.as_array().is_empty()); assert!(vector.into_array().is_empty()); assert_eq!(vector.dot(&Vector::zero()), Ok(0.0)); - let dot_certificate = vector.dot_with_errbound(&Vector::zero()).unwrap().unwrap(); - assert_abs_diff_eq!(dot_certificate.absolute_error_bound(), 0.0, epsilon = 0.0); - let difference_certificate = vector + let dot_bound = vector.dot_with_errbound(&Vector::zero()).unwrap().unwrap(); + assert_abs_diff_eq!(dot_bound.absolute_error_bound(), 0.0, epsilon = 0.0); + let difference_bound = vector .dot_difference_with_errbound(&Vector::zero(), &Vector::zero()) .unwrap() .unwrap(); - assert_abs_diff_eq!( - difference_certificate.absolute_error_bound(), - 0.0, - epsilon = 0.0 - ); + assert_abs_diff_eq!(difference_bound.absolute_error_bound(), 0.0, epsilon = 0.0); assert_eq!(vector.norm2_sq(), Ok(0.0)); } + + #[test] + fn certified_dot_withholds_bound_when_exact_product_exceeds_finite_range() { + let left_value = f64::from_bits(0x7fe3_0319_b612_3729); + let right_value = f64::from_bits(0x3ffa_ee21_bf46_bc00); + let left = Vector::<1>::new([left_value]); + let right = Vector::<1>::new([right_value]); + + assert!(left_value.mul_add(right_value, -f64::MAX) > 0.0); + assert_eq!(left.dot(&right), Ok(f64::MAX)); + assert_eq!(left.dot_with_errbound(&right), Ok(None)); + } }