diff --git a/AGENTS.md b/AGENTS.md index 1177d83..eaaa394 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -94,7 +94,7 @@ invariant over the convenient edit. `Option` instead of relying on `panic!`, `assert!`, `unwrap`, or `expect`. - Borrow by default (`&T`, `&[T]`); return borrowed views when possible. - Type and function names match textbook vocabulary (`Matrix`, `Vector`, - `Lu`, `Ldlt`, `solve`, `det`, `inf_norm`). Avoid Rust-ecosystem + `Lu`, `Ldlt`, `solve`, `det`, `norm_inf`). Avoid Rust-ecosystem abstractions that obscure the math. ### Scientific notation in docs @@ -448,7 +448,7 @@ When creating or updating issues: exact-conversion categories - `src/tolerance.rs`: validated singular-tolerance policy - `src/vector.rs`: `Vector` (`[f64; D]`) - - `src/matrix.rs`: `Matrix` (`[[f64; D]; D]`) + helpers (`get`, `try_get`, `set`, `inf_norm`, `det`, `det_direct`) + - `src/matrix.rs`: `Matrix` (`[[f64; D]; D]`) + helpers (`get`, `try_get`, `set`, `norm_inf`, `det`, `det_direct`) - `src/lu.rs`: `Lu` factorization with partial pivoting (`solve`, `det`) - `src/ldlt.rs`: `Ldlt` factorization without pivoting for exactly symmetric positive-definite matrices (`solve`, `det`) diff --git a/README.md b/README.md index d4ab9a5..0d65a94 100644 --- a/README.md +++ b/README.md @@ -75,6 +75,7 @@ for the algorithms, validity boundaries, and supporting references. - ✅ 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`) +- ✅ Overflow- and underflow-safe Euclidean vector norms (`norm`) - ✅ 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 @@ -565,6 +566,29 @@ 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. +## 📏 Overflow-safe Euclidean norms + +`Vector::norm()` computes the Euclidean norm with a deterministic scaled +sum-of-squares recurrence, so large or subnormal finite coordinates do not fail +merely because their raw squares overflow or underflow. It returns positive zero +for empty and all-zero vectors and reports `LaError::NonFinite` with +`ArithmeticOperation::VectorNorm` only when the exact norm rounds to infinity. +Near the upper range, a fixed-size stack accumulator sums squares exactly and +compares squared rounding midpoints to prevent false or hidden overflow. This +fallback needs no optional dependencies. The general binary64 result remains +approximate and has no certified error bound. + +`Vector::norm_squared()` remains the direct left-to-right FMA sum of squares for +callers that need the squared norm. Its distinct contract deliberately reports +overflow when that square is not finite, even when `norm()` can return a finite +norm. + +**v0.4.6 migration:** `Vector::norm2_sq()` is renamed to `Vector::norm_squared()`, +the unreleased `Vector::norm2()` API is named `Vector::norm()`, and +`Matrix::inf_norm()` is renamed to `Matrix::norm_inf()`. The old method names +are removed; their numerical behavior and error contracts are unchanged by +the renames. `Matrix::norm_inf()` remains the maximum absolute row sum. + ## 🎯 Certified dot products and affine differences `Vector::dot_with_errbound()` evaluates the same left-to-right FMA tree as @@ -698,7 +722,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`, `dot_with_errbound`, `dot_difference_with_errbound`, `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`, `norm`, `norm_squared` | | `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` | diff --git a/REFERENCES.md b/REFERENCES.md index 6b1bbce..49b0237 100644 --- a/REFERENCES.md +++ b/REFERENCES.md @@ -131,6 +131,17 @@ guarantee; worst-case growth and average-case behavior are distinct concerns. See [1-3, 11-12] for stability analysis, finite-precision behavior, and standard algorithmic background. +### Scaled Euclidean vector norm + +`Vector::norm` maintains a scale and a sum of squares relative to that scale, +avoiding raw coordinate squares that would overflow or underflow. This follows +the scaled safe-norm approach described by Blue [15]. The implementation retains +a deterministic coordinate order and documents its binary64 rounding contract; +it does not generally claim a certified error bound or correct rounding. A +fixed-size integer sum of exact binary64 squares handles the upper range, using +the representation and nearest-even rounding model in [9-10] to distinguish +finite results from overflow. + ## References 1. Trefethen, Lloyd N., and Robert S. Schreiber. "Average-case stability of Gaussian elimination." @@ -188,3 +199,6 @@ algorithmic background. *IEEE Std 1788-2015*, 2015: 1–97. [DOI](https://doi.org/10.1109/IEEESTD.2015.7140721) · [IEEE record](https://standards.ieee.org/ieee/1788/4431/) +15. Blue, James L. "A Portable Fortran Program to Find the Euclidean Norm of a Vector." + *ACM Transactions on Mathematical Software* 4.1 (1978): 15–23. + [DOI](https://doi.org/10.1145/355769.355771) diff --git a/benches/common/vs_linalg.rs b/benches/common/vs_linalg.rs index 376490e..0c5cc29 100644 --- a/benches/common/vs_linalg.rs +++ b/benches/common/vs_linalg.rs @@ -4,7 +4,7 @@ use faer::linalg::solvers::{Ldlt as FaerLdlt, PartialPivLu}; use faer::perm::PermRef; -use la_stack::{LaError, Tolerance, Vector}; +use la_stack::{LaError, Matrix, Tolerance, Vector}; use nalgebra::SMatrix; /// Evaluate la-stack's dot product through the ownership contract used by the @@ -33,6 +33,54 @@ pub fn la_stack_dot(left: &Vector, right: &Vector) -> Resu (*left).dot(*right) } +/// Evaluate the squared norm through the current public API. +/// +/// # Errors +/// +/// Returns the library's typed error if squared-norm accumulation overflows. +#[cfg(not(any(la_stack_pre_rational_input_api, la_stack_v0_4_3_api)))] +#[inline] +pub const fn la_stack_norm_squared(vector: &Vector) -> Result { + vector.norm_squared() +} + +/// Evaluate the same squared norm through the pre-v0.4.6 method name. +/// +/// This adapter preserves historical comparisons without adding library aliases. +/// +/// # Errors +/// +/// Returns the selected revision's typed error if squared-norm accumulation overflows. +#[cfg(any(la_stack_pre_rational_input_api, la_stack_v0_4_3_api))] +#[inline] +pub const fn la_stack_norm_squared(vector: &Vector) -> Result { + vector.norm2_sq() +} + +/// Evaluate the matrix infinity norm through the current public API. +/// +/// # Errors +/// +/// Returns the library's typed error if an absolute row sum overflows. +#[cfg(not(any(la_stack_pre_rational_input_api, la_stack_v0_4_3_api)))] +#[inline] +pub const fn la_stack_norm_inf(matrix: &Matrix) -> Result { + matrix.norm_inf() +} + +/// Evaluate the same matrix infinity norm through the pre-v0.4.6 method name. +/// +/// This adapter preserves historical comparisons without adding library aliases. +/// +/// # Errors +/// +/// Returns the selected revision's typed error if an absolute row sum overflows. +#[cfg(any(la_stack_pre_rational_input_api, la_stack_v0_4_3_api))] +#[inline] +pub const fn la_stack_norm_inf(matrix: &Matrix) -> Result { + matrix.inf_norm() +} + /// Parse a tolerance through the constructor exposed by the selected library /// revision. /// @@ -256,6 +304,111 @@ pub fn make_vector_array(offset: f64) -> [f64; D] { data } +/// Build a norm input whose non-zero magnitudes decrease after the first entry. +#[inline] +#[must_use] +pub fn make_norm_descending_array() -> [f64; D] { + std::array::from_fn(|index| { + let magnitude = vector_entry(D - index - 1, 0.0); + if index % 2 == 0 { + magnitude + } else { + -magnitude + } + }) +} + +/// Build a norm input whose entries repeatedly have the same non-zero magnitude. +#[inline] +#[must_use] +pub fn make_norm_repeated_scale_array() -> [f64; D] { + std::array::from_fn(|index| if index % 2 == 0 { 3.0 } else { -3.0 }) +} + +/// Build a norm input with one non-zero entry and otherwise skipped zeros. +#[inline] +#[must_use] +pub fn make_norm_sparse_array() -> [f64; D] { + let mut data = [0.0; D]; + if D > 0 { + data[D / 2] = -3.0; + } + data +} + +/// Build a finite norm input spanning normal, subnormal, and zero magnitudes. +#[inline] +#[must_use] +pub fn make_norm_wide_dynamic_range_array() -> [f64; D] { + const VALUES: [f64; 8] = [ + 1.0e200, + -1.0e-200, + f64::from_bits(1), + 0.0, + -1.0e100, + 1.0e-100, + 1.0, + -1.0e200, + ]; + + std::array::from_fn(|index| VALUES[index % VALUES.len()]) +} + +/// Build the named Euclidean-norm scenario corpus in stable benchmark order. +#[inline] +#[must_use] +pub fn norm_scenarios() -> [(&'static str, [f64; D]); 4] { + [ + ("descending", make_norm_descending_array()), + ("repeated_scale", make_norm_repeated_scale_array()), + ("sparse", make_norm_sparse_array()), + ("wide_dynamic_range", make_norm_wide_dynamic_range_array()), + ] +} + +/// Compute a Euclidean norm by chaining the standard library's binary `hypot`. +/// +/// This is an independent overflow- and underflow-safe reference kernel for the +/// focused `norm2` benchmarks. +#[inline] +#[must_use] +pub fn iterative_hypot(values: &[f64; D]) -> f64 { + values.iter().fold(0.0, |norm, &value| norm.hypot(value)) +} + +/// Reproduce the generalized scaled norm currently owned by Delaunay. +/// +/// D=2 delegates to binary `hypot`; larger dimensions make one pass for the +/// maximum magnitude and a second pass for the scaled sum of squares. +#[inline] +#[must_use] +pub fn delaunay_scaled_norm(values: &[f64; D]) -> f64 { + match D { + 0 => 0.0, + 1 => values[0].abs(), + 2 => values[0].hypot(values[1]), + _ => { + if values.iter().any(|value| value.is_nan()) { + return f64::NAN; + } + if values.iter().any(|value| value.is_infinite()) { + return f64::INFINITY; + } + + let scale = values.iter().map(|value| value.abs()).fold(0.0, f64::max); + if scale == 0.0 { + return 0.0; + } + + let scaled_sum = values.iter().fold(0.0, |sum, &value| { + let ratio = value / scale; + ratio.mul_add(ratio, sum) + }); + scale * scaled_sum.sqrt() + } + } +} + /// Compute nalgebra's matrix infinity norm using la-stack's row-sum convention. #[inline] #[must_use] diff --git a/benches/vs_linalg.rs b/benches/vs_linalg.rs index 1e11b84..13fae4e 100644 --- a/benches/vs_linalg.rs +++ b/benches/vs_linalg.rs @@ -9,6 +9,7 @@ //! - Determinant groups distinguish factorization-inclusive LU, `Matrix::det`, //! precomputed LU determinant queries, and precomputed LDLT/Cholesky queries. //! - Matrix infinity norm is the maximum absolute row sum on all sides. +//! - Euclidean-norm rows include safe reference kernels as well as peer crates. use std::hint::black_box; @@ -27,10 +28,13 @@ mod bench_utils; pub mod vs_linalg_common; use bench_utils::OrAbort; +#[cfg(not(any(la_stack_pre_rational_input_api, la_stack_v0_4_3_api)))] +use vs_linalg_common::norm_scenarios; use vs_linalg_common::{ - PreparedFaerLuDet, faer_det_from_ldlt, la_stack_dot, la_stack_tolerance, - make_balanced_dynamic_range_rows, make_ill_conditioned_matrix_rows, make_matrix_rows, - make_pivoting_matrix_rows, make_vector_array, matrix_entry, nalgebra_inf_norm, vector_entry, + PreparedFaerLuDet, delaunay_scaled_norm, faer_det_from_ldlt, iterative_hypot, la_stack_dot, + la_stack_norm_inf, la_stack_norm_squared, la_stack_tolerance, make_balanced_dynamic_range_rows, + make_ill_conditioned_matrix_rows, make_matrix_rows, make_pivoting_matrix_rows, + make_vector_array, matrix_entry, nalgebra_inf_norm, vector_entry, }; /// Build the deterministic la-stack matrix shared by a benchmark family. @@ -406,7 +410,7 @@ fn register_precomputed_ldlt_determinant_benchmarks( }); } -/// Register vector dot-product and squared-norm benchmarks. +/// Register vector dot-product, squared-norm, and Euclidean-norm benchmarks. fn register_vector_benchmarks(group: &mut BenchmarkGroup<'_, WallTime>) { let v1 = la_vector::(0.0, "la_stack vector construction"); let v2 = la_vector::(1.0, "la_stack vector construction"); @@ -443,7 +447,7 @@ fn register_vector_benchmarks(group: &mut BenchmarkGroup<'_, Wal group.bench_function("la_stack_norm2_sq", |bencher| { bencher.iter(|| { - let result = black_box(&v1).norm2_sq().or_abort("la_stack norm2_sq"); + let result = la_stack_norm_squared(black_box(&v1)).or_abort("la_stack norm_squared"); black_box(result); }); }); @@ -462,6 +466,92 @@ fn register_vector_benchmarks(group: &mut BenchmarkGroup<'_, Wal black_box(result); }); }); + + #[cfg(not(any(la_stack_pre_rational_input_api, la_stack_v0_4_3_api)))] + { + group.bench_function("la_stack_norm2", |bencher| { + bencher.iter(|| { + let result = black_box(&v1).norm().or_abort("la_stack norm"); + black_box(result); + }); + }); + } + + group.bench_function("iterative_f64_hypot", |bencher| { + bencher.iter(|| { + let result = iterative_hypot(black_box(v1.as_array())); + black_box(result); + }); + }); + + group.bench_function("delaunay_scaled_norm", |bencher| { + bencher.iter(|| { + let result = delaunay_scaled_norm(black_box(v1.as_array())); + black_box(result); + }); + }); + + group.bench_function("nalgebra_norm", |bencher| { + bencher.iter(|| { + let result = black_box(&nv1).norm(); + black_box(result); + }); + }); + + group.bench_function("faer_norm_l2", |bencher| { + bencher.iter(|| { + let result = black_box(&fv1).as_mat_ref().norm_l2(); + black_box(result); + }); + }); +} + +/// Validate and register safe Euclidean-norm kernels across branch and range profiles. +#[cfg(not(any(la_stack_pre_rational_input_api, la_stack_v0_4_3_api)))] +fn register_norm_scenario_benchmarks(group: &mut BenchmarkGroup<'_, WallTime>) { + for (scenario, values) in norm_scenarios::() { + let vector = Vector::try_new(values).or_abort("norm scenario vector construction"); + let la_stack = vector.norm().or_abort("la_stack norm scenario validation"); + let hypot = iterative_hypot(&values); + let delaunay = delaunay_scaled_norm(&values); + let scale = la_stack.abs().max(hypot.abs()).max(delaunay.abs()).max(1.0); + assert!( + la_stack.is_finite() + && hypot.is_finite() + && delaunay.is_finite() + && (hypot - la_stack).abs() <= 1.0e-12 * scale + && (delaunay - la_stack).abs() <= 1.0e-12 * scale, + "norm benchmark scenario {scenario} failed setup validation: \ + la_stack={la_stack:?}, hypot={hypot:?}, delaunay={delaunay:?}", + ); + + group.bench_function(format!("la_stack_norm2_scenario_{scenario}"), |bencher| { + bencher.iter(|| { + let result = black_box(&vector).norm().or_abort("la_stack norm scenario"); + black_box(result); + }); + }); + + group.bench_function( + format!("iterative_f64_hypot_norm2_scenario_{scenario}"), + |bencher| { + bencher.iter(|| { + let result = iterative_hypot(black_box(&values)); + black_box(result); + }); + }, + ); + + group.bench_function( + format!("delaunay_scaled_norm_norm2_scenario_{scenario}"), + |bencher| { + bencher.iter(|| { + let result = delaunay_scaled_norm(black_box(&values)); + black_box(result); + }); + }, + ); + } } /// Register matrix infinity-norm benchmarks. @@ -472,7 +562,7 @@ fn register_matrix_norm_benchmarks(group: &mut BenchmarkGroup<'_ group.bench_function("la_stack_inf_norm", |bencher| { bencher.iter(|| { - let result = black_box(&a).inf_norm().or_abort("la_stack inf_norm"); + let result = la_stack_norm_inf(black_box(&a)).or_abort("la_stack norm_inf"); black_box(result); }); }); @@ -585,6 +675,8 @@ macro_rules! define_vs_linalg_benches_for_dim { register_precomputed_lu_determinant_benchmarks::<$d>(&mut group); register_precomputed_ldlt_determinant_benchmarks::<$d>(&mut group); register_vector_benchmarks::<$d>(&mut group); + #[cfg(not(any(la_stack_pre_rational_input_api, la_stack_v0_4_3_api)))] + register_norm_scenario_benchmarks::<$d>(&mut group); register_matrix_norm_benchmarks::<$d>(&mut group); $( $register_stress(&mut group); diff --git a/docs/BENCHMARKING.md b/docs/BENCHMARKING.md index 1bf0603..e80e586 100644 --- a/docs/BENCHMARKING.md +++ b/docs/BENCHMARKING.md @@ -196,7 +196,7 @@ Use local saved baselines when tuning one kernel and comparing several edits against the same starting point. These baselines are local scratch data, not release artifacts. -For example, before optimizing `Matrix::inf_norm`, save a named baseline: +For example, before optimizing `Matrix::norm_inf`, save a named baseline: ```bash just bench-save-baseline inf-norm-before vs_linalg @@ -400,14 +400,33 @@ adapter code computes the agreed mathematical kernel inside the timed closure: | `det_via_lu`, `det_from_lu` | Native `Lu::det` | Native `LU::determinant` | Harness adapter: product of the U diagonal and permutation sign | | `det_from_ldlt` / `det_from_cholesky` | Native `Ldlt::det` | Native `Cholesky::determinant` | Harness adapter: product of the D diagonal | | `dot` | Native `Vector::dot` | Native `dot` | Harness adapter: left-to-right fused multiply-add loop | -| `norm2_sq` | Native `Vector::norm2_sq` | Native `norm_squared` | Native `squared_norm_l2` | -| `inf_norm` | Native `Matrix::inf_norm` | Harness adapter: maximum absolute row sum | Harness adapter: maximum absolute row sum | - -These adapter timings are benchmark-kernel comparisons, not claims about the -speed of an identically named public convenience method in every crate. The -adapter implementation is versioned with the benchmark harness, included in the +| `norm2_sq` | Native `Vector::norm_squared` | Native `norm_squared` | Native `squared_norm_l2` | +| `norm2` | Native `Vector::norm` | Native `norm` | Native `norm_l2` | +| `inf_norm` | Native `Matrix::norm_inf` | Harness adapter: maximum absolute row sum | Harness adapter: maximum absolute row sum | + +The `norm2`, `norm2_sq`, and `inf_norm` benchmark IDs, including scenario suffixes, are +retained for continuity with saved baselines. The current public methods are +`Vector::norm()`, `Vector::norm_squared()`, and `Matrix::norm_inf()`; the historical +benchmark adapters call `norm2_sq()` and `inf_norm()` only when building against +older library releases. + +The `norm2` family also measures iterative `f64::hypot` and Delaunay's existing +dimension-specialized scaled implementation as labeled reference kernels. These +adapter timings are benchmark-kernel comparisons, not claims about the speed of +an identically named public convenience method in every crate. Adapter +implementations are versioned with the benchmark harness, included in the benchmark-contract digest, and covered by the cross-crate input smoke tests. +The ordinary cross-crate norm row uses positive magnitudes in increasing order, +which changes the running scale at every entry in la-stack's one-pass recurrence. +Additional `norm2` rows cover decreasing magnitudes, repeated scales, sparse +vectors, and finite values spanning normal, subnormal, and zero magnitudes. Those +scenario rows compare only la-stack with the two overflow- and underflow-safe +reference kernels; peer methods whose behavior is not contract-equivalent on the +wide-range input are deliberately excluded. Fixture construction and agreement +checks occur before Criterion's measured closures. Run the focused corpus with +`just bench-vs-linalg-quick norm2_scenario`. + All three crates receive equivalent deterministic inputs for a given dimension: - matrix entries come from the same strictly diagonally-dominant generator @@ -426,7 +445,7 @@ All three crates receive equivalent deterministic inputs for a given dimension: closure, applying the same complete-operation protocol to la-stack and nalgebra - borrowed operations receive references through `black_box`; in particular, - `inf_norm` does not copy the matrix inside the measured closure + `Matrix::norm_inf()` does not copy the matrix inside the measured closure Use `iter_batched` only when fixture construction is explicitly outside the scientific quantity being measured. The exclusion must be symmetric across the @@ -436,7 +455,7 @@ the reported kernel time or cross-crate ratio. The integration smoke test `tests/vs_linalg_inputs.rs` reuses the benchmark input helpers and verifies that la-stack, nalgebra, and faer agree on the -determinant, solve, dot, squared-norm, and infinity-norm results for every +determinant, solve, dot, Euclidean-norm, squared-norm, and infinity-norm results for every measured dimension: D=2, 3, 4, 5, 8, 16, 32, and 64. The same focused recipe also tests exact-benchmark range and deterministic-generator configuration: @@ -462,6 +481,7 @@ The main comparable metrics are: - `solve_from_lu` — solve one right-hand side using a precomputed LU factor - `det_from_lu` — compute determinant using a precomputed LU factor - `dot` — vector dot product +- `norm2` — overflow- and underflow-safe Euclidean vector norm - `norm2_sq` — squared Euclidean vector norm - `inf_norm` — matrix infinity norm, implemented as maximum absolute row sum diff --git a/docs/mathematical_basis.md b/docs/mathematical_basis.md index 7d9dadc..54d02c8 100644 --- a/docs/mathematical_basis.md +++ b/docs/mathematical_basis.md @@ -50,9 +50,43 @@ because it enumerates concrete stack types. 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. +error bounds. This includes plain `Vector::dot`, Euclidean and squared norms, +matrix norms, factorizations, and solves. Some kernels use FMA to reduce rounding +steps, but that does not make them exact. + +## Scaled Euclidean vector norm + +`Vector::norm` computes `sqrt(Σᵢ xᵢ²)` with a left-to-right scaled +sum-of-squares recurrence. Its state represents the accumulated squared norm as +`scale² × scaled_sum`. For each nonzero `|xᵢ|`, either `|xᵢ| / scale` is +squared and accumulated, or a larger `|xᵢ|` becomes the new scale and the old +sum is rescaled. Every squared ratio is therefore at most one, which avoids raw +square overflow and scales all-subnormal vectors into a safe range \[15\]. + +Division, FMA, square root, and final rescaling still round in binary64. The +method is deterministic for a fixed coordinate order but does not generally +claim correct rounding or publish a certified error bound. + +Near the upper range, a fixed-size integer accumulator avoids both false and +hidden overflow from the rounded recurrence. Every coordinate square is an +integer multiple of `2^-2148` and is below `2^2048`, so 4196 bits plus +`usize::BITS` carry bits suffice for every representable vector length. The +fallback sums these integer squares exactly, retaining even the smallest +subnormal square. It compares against squared binary64 rounding midpoints to +return the nearest norm, ties to even, without a floating square root. The +overflow midpoint is `f64::MAX + 2^970`; equality rounds to infinity \[9-10\]. + +The fallback is selected from the largest coordinate magnitude, independently +of the rounded norm. With `b = bit_length(D)`, a largest magnitude at most +`2^(1023-b)` gives the conservative L1 upper bound `D × scale < 2^1023`, so +the ordinary recurrence has ample overflow margin. Larger magnitudes use the +exact boundary calculation. Its fixed storage stays on the stack and requires +no `exact` feature. + +A scalar `LaError::NonFinite` tagged with `ArithmeticOperation::VectorNorm` +therefore means the exact norm rounds to infinity. `Vector::norm_squared` +intentionally remains the direct FMA sum `Σᵢ xᵢ²` and may therefore fail even +when `norm` succeeds. ## Outward-rounded interval expressions diff --git a/src/error.rs b/src/error.rs index 6e40661..0c6182c 100644 --- a/src/error.rs +++ b/src/error.rs @@ -51,6 +51,8 @@ pub enum ArithmeticOperation { VectorDotDifference, /// Vector squared-norm calculation. VectorSquaredNorm, + /// Overflow- and underflow-safe vector Euclidean-norm calculation. + VectorNorm, } impl fmt::Display for ArithmeticOperation { @@ -72,6 +74,7 @@ impl fmt::Display for ArithmeticOperation { Self::VectorDotProduct => "vector dot product", Self::VectorDotDifference => "vector dot difference", Self::VectorSquaredNorm => "vector squared norm", + Self::VectorNorm => "vector Euclidean norm", }) } } @@ -915,6 +918,10 @@ mod tests { ArithmeticOperation::VectorSquaredNorm.to_string(), "vector squared norm" ); + assert_eq!( + ArithmeticOperation::VectorNorm.to_string(), + "vector Euclidean norm" + ); assert_eq!( ArithmeticOperation::IntervalAddition.to_string(), "interval addition" diff --git a/src/exact.rs b/src/exact.rs index 7d00b6d..e659dd3 100644 --- a/src/exact.rs +++ b/src/exact.rs @@ -1725,7 +1725,7 @@ impl Matrix { /// [4.0, 5.0, 6.0], /// [7.0, 8.0, 9.0], /// ])?; - /// // This matrix is singular (row 3 = row 1 + row 2 in exact arithmetic). + /// // This matrix is singular (row 3 = 2 × row 2 − row 1 in exact arithmetic). /// assert_eq!(m.det_sign_exact(), DeterminantSign::Zero); /// /// assert_eq!(Matrix::<3>::identity().det_sign_exact(), DeterminantSign::Positive); diff --git a/src/interval.rs b/src/interval.rs index e533d0d..d80449f 100644 --- a/src/interval.rs +++ b/src/interval.rs @@ -2,6 +2,7 @@ //! Outward-rounded intervals and fixed-size interval determinant signs. +use crate::rounding::{compare_product_with_rounded, two_sum_error}; use crate::{ArithmeticOperation, IntervalBound, IntervalOperand, LaError, Matrix}; /// Largest dimension supported by [`IntervalMatrix::det`] and @@ -93,108 +94,6 @@ const fn canonical_zero(value: f64) -> f64 { if value == 0.0 { 0.0 } else { value } } -/// Return the exact error in a rounded binary64 sum. -/// -/// This is Knuth's `TwoSum` transform. With IEEE-754 round-to-nearest and -/// gradual underflow, `rounded + error` equals the exact-real sum whenever the -/// rounded sum is finite. -#[inline] -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 -} - -/// Decompose a nonzero finite binary64 magnitude as `significand × 2^exponent`. -#[inline] -const fn decompose_magnitude(value: f64) -> (u128, i64) { - let magnitude_bits = value.to_bits() & 0x7fff_ffff_ffff_ffff; - let biased_exponent = ((magnitude_bits >> 52) & 0x7ff).cast_signed(); - let fraction = magnitude_bits & 0x000f_ffff_ffff_ffff; - - if biased_exponent == 0 { - (fraction as u128, -1074) - } else { - ( - (fraction | (1_u64 << 52)) as u128, - biased_exponent - 1023 - 52, - ) - } -} - -/// Compare two positive values represented as `significand × 2^exponent`. -#[inline] -const fn compare_binary_magnitudes( - left_significand: u128, - left_exponent: i64, - right_significand: u128, - right_exponent: i64, -) -> i8 { - let left_zeros = left_significand.trailing_zeros() as i64; - let right_zeros = right_significand.trailing_zeros() as i64; - let normalized_left = left_significand >> left_zeros.cast_unsigned(); - let normalized_right = right_significand >> right_zeros.cast_unsigned(); - let normalized_left_exponent = left_exponent + left_zeros; - let normalized_right_exponent = right_exponent + right_zeros; - - let left_top = - normalized_left_exponent + (u128::BITS - normalized_left.leading_zeros() - 1) as i64; - let right_top = - normalized_right_exponent + (u128::BITS - normalized_right.leading_zeros() - 1) as i64; - if left_top < right_top { - return -1; - } - if left_top > right_top { - return 1; - } - - let common_exponent = if normalized_left_exponent < normalized_right_exponent { - normalized_left_exponent - } else { - normalized_right_exponent - }; - let aligned_left = - normalized_left << (normalized_left_exponent - common_exponent).cast_unsigned(); - let aligned_right = - normalized_right << (normalized_right_exponent - common_exponent).cast_unsigned(); - if aligned_left < aligned_right { - -1 - } else if aligned_left > aligned_right { - 1 - } else { - 0 - } -} - -/// Compare the exact-real product `left × right` with its rounded result. -#[inline] -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 }; - } - - let (left_significand, left_exponent) = decompose_magnitude(left); - let (right_significand, right_exponent) = decompose_magnitude(right); - let exact_significand = left_significand * right_significand; - let exact_exponent = left_exponent + right_exponent; - let (rounded_significand, rounded_exponent) = decompose_magnitude(rounded); - let magnitude_relation = compare_binary_magnitudes( - exact_significand, - exact_exponent, - rounded_significand, - rounded_exponent, - ); - - if negative { - -magnitude_relation - } else { - magnitude_relation - } -} - /// Turn a finite rounded sum into the tight adjacent-float enclosure implied by /// its exact `TwoSum` residual. #[inline] diff --git a/src/lib.rs b/src/lib.rs index e881d92..6744f9e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -321,8 +321,10 @@ mod interval; mod ldlt; mod lu; mod matrix; +mod norm; #[cfg(feature = "exact")] mod rational; +mod rounding; mod scaled_product; mod tolerance; mod vector; @@ -1061,7 +1063,7 @@ mod tests { #[test] fn try_with_stack_matrix_converts_unsupported_dimension_error() { let got = try_with_stack_matrix!(9usize, |m| -> Result { - assert_abs_diff_eq!(m.inf_norm()?, 0.0, epsilon = 0.0); + assert_abs_diff_eq!(m.norm_inf()?, 0.0, epsilon = 0.0); Ok(0) }); diff --git a/src/matrix.rs b/src/matrix.rs index de97171..16abf1c 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -486,9 +486,11 @@ impl Matrix { /// Infinity norm (maximum absolute row sum). /// + /// This is the induced matrix L∞ norm: `‖A‖∞ = maxᵢ Σⱼ |Aᵢⱼ|`, not + /// the largest absolute entry. Each row contributes its L₁ norm. + /// /// # Non-finite handling - /// Public constructors and setters reject raw non-finite entries, but - /// `Matrix` values are finite by construction. `inf_norm` returns + /// `Matrix` values are finite by construction. [`Self::norm_inf`] returns /// [`LaError::NonFinite`] with the matrix cell whose addition first makes a /// row sum non-finite. /// @@ -503,7 +505,7 @@ impl Matrix { /// /// # fn main() -> Result<(), LaError> { /// let m = Matrix::<2>::try_from_rows([[1.0, -2.0], [3.0, 4.0]])?; - /// assert!((m.inf_norm()? - 7.0).abs() <= 1e-12); + /// assert!((m.norm_inf()? - 7.0).abs() <= 1e-12); /// /// // Raw NaN entries are rejected with coordinates. /// assert_matches!( @@ -522,7 +524,7 @@ impl Matrix { /// Returns [`LaError::NonFinite`] with matrix coordinates when a row sum /// overflows to NaN or infinity. #[inline] - pub const fn inf_norm(&self) -> Result { + pub const fn norm_inf(&self) -> Result { let mut max_row_sum: f64 = 0.0; let mut r = 0; @@ -536,7 +538,7 @@ impl Matrix { } if !row_sum.is_finite() { cold_path(); - return Err(Self::inf_norm_overflow_error(row, r)); + return Err(Self::norm_inf_overflow_error(row, r)); } if row_sum > max_row_sum { max_row_sum = row_sum; @@ -555,7 +557,7 @@ impl Matrix { /// first column whose addition overflowed; if every earlier prefix is /// finite, the final column is that first failure. #[cold] - const fn inf_norm_overflow_error(row: &[f64; D], row_index: usize) -> LaError { + const fn norm_inf_overflow_error(row: &[f64; D], row_index: usize) -> LaError { let mut row_sum = 0.0; let mut col = 0; let last_col = D.saturating_sub(1); @@ -582,7 +584,7 @@ impl Matrix { /// /// Two entries `self[r][c]` and `self[c][r]` are considered equal (for the /// purposes of symmetry) when - /// `|self[r][c] - self[c][r]| <= rel_tol * max(1.0, inf_norm(self))`. + /// `|self[r][c] - self[c][r]| <= rel_tol * max(1.0, norm_inf(self))`. /// This is a diagnostic predicate for applications that have an /// approximation-specific symmetry threshold. It is not the precondition /// used by [`ldlt`](Self::ldlt), which requires exact mirrored-entry @@ -633,7 +635,7 @@ impl Matrix { /// Iteration order is row-major over the strict upper triangle, so the /// returned indices are the lexicographically smallest such pair. The /// predicate is the same as [`is_symmetric`](Self::is_symmetric): - /// `|self[r][c] - self[c][r]| <= rel_tol * max(1.0, inf_norm(self))`. + /// `|self[r][c] - self[c][r]| <= rel_tol * max(1.0, norm_inf(self))`. /// It is intentionally distinct from the exact equality required by /// [`ldlt`](Self::ldlt). /// @@ -855,7 +857,7 @@ impl Matrix { return Ok(rel_tol); } - if let Ok(norm) = self.inf_norm() { + if let Ok(norm) = self.norm_inf() { let scale = if norm > 1.0 { norm } else { 1.0 }; let eps = rel_tol * scale; if eps.is_finite() { @@ -1804,14 +1806,14 @@ mod tests { #[test] fn []() { let z = Matrix::<$d>::zero(); - assert_abs_diff_eq!(z.inf_norm().unwrap(), 0.0, epsilon = 0.0); + assert_abs_diff_eq!(z.norm_inf().unwrap(), 0.0, epsilon = 0.0); let d = Matrix::<$d>::default(); - assert_abs_diff_eq!(d.inf_norm().unwrap(), 0.0, epsilon = 0.0); + assert_abs_diff_eq!(d.norm_inf().unwrap(), 0.0, epsilon = 0.0); } #[test] - fn []() { + fn []() { let mut rows = [[0.0f64; $d]; $d]; // Row 0 has a smaller absolute row sum. @@ -1825,18 +1827,18 @@ mod tests { } let m = Matrix::<$d>::try_from_rows(rows).unwrap(); - assert_abs_diff_eq!(m.inf_norm().unwrap(), f64::from($d), epsilon = 0.0); + assert_abs_diff_eq!(m.norm_inf().unwrap(), f64::from($d), epsilon = 0.0); } #[test] - fn []() { + fn []() { let mut rows = [[0.0f64; $d]; $d]; rows[$d - 1][0] = f64::MAX; rows[$d - 1][1] = f64::MAX; let m = Matrix::<$d>::try_from_rows(rows).unwrap(); assert_eq!( - m.inf_norm(), + m.norm_inf(), Err(LaError::non_finite_computation_matrix( ArithmeticOperation::MatrixInfinityNorm, $d - 1, @@ -1846,7 +1848,7 @@ mod tests { } #[test] - fn []() { + fn []() { let mut rows = [[0.0f64; $d]; $d]; rows[0][0] = f64::MAX; rows[0][$d - 1] = f64::MAX; @@ -1855,7 +1857,7 @@ mod tests { let m = Matrix::<$d>::try_from_rows(rows).unwrap(); assert_eq!( - m.inf_norm(), + m.norm_inf(), Err(LaError::non_finite_computation_matrix( ArithmeticOperation::MatrixInfinityNorm, 0, @@ -1911,13 +1913,13 @@ mod tests { gen_matrix_tests!(5); #[test] - fn matrix_inf_norm_preserves_left_to_right_row_sum_order() { + fn matrix_norm_inf_preserves_left_to_right_row_sum_order() { let large = 9_007_199_254_740_992.0; let matrix = Matrix::<4>::try_from_rows([[large, 1.0, 1.0, 1.0], [0.0; 4], [0.0; 4], [0.0; 4]]) .unwrap(); - assert_eq!(matrix.inf_norm(), Ok(large)); + assert_eq!(matrix.norm_inf(), Ok(large)); } // === det_direct tests === @@ -2552,27 +2554,27 @@ mod tests { assert_eq!(BOUND, Ok(None)); } - // === inf_norm const-evaluability tests (D = 2..=5) === + // === norm_inf const-evaluability tests (D = 2..=5) === - macro_rules! gen_inf_norm_const_eval_tests { + macro_rules! gen_norm_inf_const_eval_tests { ($d:literal) => { paste! { - /// `Matrix::::inf_norm()` on the identity must const-evaluate + /// `Matrix::::norm_inf()` on the identity must const-evaluate /// to `1.0` for every `D ≥ 1` — each row has a single `1.0` /// entry, so the max absolute row sum is exactly `1.0`. #[test] - fn []() { - const NORM: Result = Matrix::<$d>::identity().inf_norm(); + fn []() { + const NORM: Result = Matrix::<$d>::identity().norm_inf(); assert!((NORM.unwrap() - 1.0).abs() <= 1e-12); } } }; } - gen_inf_norm_const_eval_tests!(2); - gen_inf_norm_const_eval_tests!(3); - gen_inf_norm_const_eval_tests!(4); - gen_inf_norm_const_eval_tests!(5); + gen_norm_inf_const_eval_tests!(2); + gen_norm_inf_const_eval_tests!(3); + gen_norm_inf_const_eval_tests!(4); + gen_norm_inf_const_eval_tests!(5); // === is_symmetric / first_asymmetry (public LDLT preconditions helpers) === @@ -2709,8 +2711,8 @@ mod tests { } #[test] - fn is_symmetric_tolerance_scales_with_inf_norm() { - // Off-diagonal entries differ by 1e-6. With inf_norm ≈ 2e6, the + fn is_symmetric_tolerance_scales_with_norm_inf() { + // Off-diagonal entries differ by 1e-6. With norm_inf ≈ 2e6, the // relative tolerance 1e-12 yields eps ≈ 2e-6, which accepts the gap; // a stricter tol of 1e-15 rejects it. let a = Matrix::<2>::try_from_rows([[1.0e6, 1.0e6 + 1.0e-6], [1.0e6, 1.0e6]]).unwrap(); @@ -2733,7 +2735,7 @@ mod tests { let matrix = Matrix::<5>::try_from_rows(rows).unwrap(); let tolerance = Tolerance::try_new(min_subnormal).unwrap(); - let expected_epsilon = tolerance.get() * matrix.inf_norm().unwrap().max(1.0); + let expected_epsilon = tolerance.get() * matrix.norm_inf().unwrap().max(1.0); assert_eq!(expected_epsilon.to_bits(), 2); assert_eq!(matrix.first_asymmetry(tolerance), Ok(None)); @@ -2746,7 +2748,7 @@ mod tests { Matrix::<2>::try_from_rows([[f64::MAX, f64::MAX], [f64::MAX / 2.0, f64::MAX]]).unwrap(); assert_eq!( - matrix.inf_norm(), + matrix.norm_inf(), Err(LaError::non_finite_computation_matrix( ArithmeticOperation::MatrixInfinityNorm, 0, @@ -2785,7 +2787,7 @@ mod tests { .unwrap(); assert_eq!( - a.inf_norm(), + a.norm_inf(), Err(LaError::non_finite_computation_matrix( ArithmeticOperation::MatrixInfinityNorm, 1, diff --git a/src/norm.rs b/src/norm.rs new file mode 100644 index 0000000..d756738 --- /dev/null +++ b/src/norm.rs @@ -0,0 +1,138 @@ +#![forbid(unsafe_code)] + +//! Exact range-boundary fallback for the otherwise approximate Euclidean norm. + +use core::cmp::Ordering; + +use crate::rounding::decompose_magnitude; +use crate::{ArithmeticOperation, LaError}; + +/// Squared binary64 values are integer multiples of `2^-2148`. +const SQUARE_BASE_EXPONENT: i64 = -2148; +/// Each square is below `2^2048`, requiring 4196 bits at the chosen scale. +/// Summing at most `usize::MAX` squares needs at most `usize::BITS` extra bits. +const SUM_WORDS: usize = (4196 + usize::BITS as usize).div_ceil(u64::BITS as usize); + +/// Non-negative exact sum of squares, in little-endian base-2^64 storage. +/// +/// The fixed capacity covers every finite binary64 vector length representable +/// by `usize`. No term is discarded, including subnormal squares that can break +/// an otherwise exact rounding tie. Representation follows IEEE 754 \[9-10\]. +struct SquareSum { + words: [u64; SUM_WORDS], +} + +impl SquareSum { + const ZERO: Self = Self { + words: [0; SUM_WORDS], + }; + + /// Add one word and propagate its carry within the proven sum capacity. + fn add_word(&mut self, mut index: usize, mut word: u64) { + while word != 0 { + let (sum, carry) = self.words[index].overflowing_add(word); + self.words[index] = sum; + word = u64::from(carry); + index += 1; + } + } + + /// Add `(significand × 2^exponent)²` without rounding. + /// + /// Inputs are binary64 magnitudes or upper-range rounding midpoints, so + /// the significand has at most 54 bits and its square fits `u128`. + /// Their exponents are at least -1074, so all shifts are non-negative; + /// all shifted squares and carries fit the fixed storage above. + #[expect( + clippy::cast_possible_truncation, + clippy::cast_sign_loss, + reason = "extracting low/high words intentionally truncates; finite binary64 square exponents give shifts in 0..=4090" + )] + fn add_square(&mut self, significand: u128, exponent: i64) { + let square = significand * significand; + let shift = (2 * exponent - SQUARE_BASE_EXPONENT) as usize; + let index = shift / u64::BITS as usize; + let offset = shift % u64::BITS as usize; + let low = square as u64; + let high = (square >> u64::BITS) as u64; + + self.add_word(index, low << offset); + if offset == 0 { + self.add_word(index + 1, high); + } else { + self.add_word(index + 1, (high << offset) | (low >> (64 - offset))); + self.add_word(index + 2, high >> (64 - offset)); + } + } + + /// Compare with one squared positive normal binary64 value. + fn compare_square(&self, value: f64) -> Ordering { + let (significand, exponent) = decompose_magnitude(value); + self.compare_scaled_square(significand, exponent) + } + + /// Compare with the squared midpoint above an upper-range normal value. + /// + /// If `lower = m × 2^e`, its next rounding boundary is + /// `(2m + 1) × 2^(e - 1)`, including the boundary above `f64::MAX`. + fn compare_midpoint(&self, lower: f64) -> Ordering { + let (significand, exponent) = decompose_magnitude(lower); + self.compare_scaled_square(2 * significand + 1, exponent - 1) + } + + fn compare_scaled_square(&self, significand: u128, exponent: i64) -> Ordering { + let mut other = Self::ZERO; + other.add_square(significand, exponent); + self.words.iter().rev().cmp(other.words.iter().rev()) + } +} + +/// Round a potentially overflowing norm by comparing exact squared values. +/// +/// `scale` is the largest input magnitude and is normal in this cold path. +/// The overflow midpoint is `f64::MAX + 2^970`. Its tie rounds to infinity; +/// smaller norms round to a finite value. Binary search and an exact squared +/// midpoint comparison then select the nearest finite value, ties to even. +/// This does not require a square-root approximation or optional dependencies. +#[cold] +pub(crate) fn norm_near_overflow( + values: &[f64; D], + scale: f64, +) -> Result { + let mut sum = SquareSum::ZERO; + for &value in values { + if value != 0.0 { + let (significand, exponent) = decompose_magnitude(value); + sum.add_square(significand, exponent); + } + } + + if sum.compare_midpoint(f64::MAX) != Ordering::Less { + return Err(LaError::non_finite_computation_scalar( + ArithmeticOperation::VectorNorm, + )); + } + if sum.compare_square(f64::MAX) != Ordering::Less { + return Ok(f64::MAX); + } + + // The exact norm is at least its largest coordinate. Positive binary64 bit + // patterns are monotonically ordered, and every candidate here is normal. + let mut lower = scale.to_bits(); + let mut upper = f64::MAX.to_bits(); + while upper - lower > 1 { + let middle = lower + (upper - lower) / 2; + match sum.compare_square(f64::from_bits(middle)) { + Ordering::Less => upper = middle, + Ordering::Greater => lower = middle, + Ordering::Equal => return Ok(f64::from_bits(middle)), + } + } + + let rounded = match sum.compare_midpoint(f64::from_bits(lower)) { + Ordering::Less => lower, + Ordering::Greater => upper, + Ordering::Equal => lower + (lower & 1), + }; + Ok(f64::from_bits(rounded)) +} diff --git a/src/rounding.rs b/src/rounding.rs new file mode 100644 index 0000000..a7f47fa --- /dev/null +++ b/src/rounding.rs @@ -0,0 +1,116 @@ +#![forbid(unsafe_code)] + +//! Shared binary64 rounding primitives for certified arithmetic. + +/// Return the exact error in a rounded binary64 sum. +/// +/// This is Knuth's `TwoSum` transform. With IEEE-754 round-to-nearest and +/// gradual underflow, `rounded + error` equals the exact-real sum whenever the +/// rounded sum is finite. +#[inline] +pub(crate) const fn two_sum_error(left: f64, right: f64, rounded: f64) -> f64 { + let virtual_right = rounded - left; + let virtual_left = rounded - virtual_right; + let right_error = right - virtual_right; + let left_error = left - virtual_left; + left_error + right_error +} + +/// Decompose a nonzero finite binary64 magnitude as `significand × 2^exponent`. +#[inline] +pub(crate) const fn decompose_magnitude(value: f64) -> (u128, i64) { + let magnitude_bits = value.to_bits() & 0x7fff_ffff_ffff_ffff; + let biased_exponent = ((magnitude_bits >> 52) & 0x7ff).cast_signed(); + let fraction = magnitude_bits & 0x000f_ffff_ffff_ffff; + + if biased_exponent == 0 { + (fraction as u128, -1074) + } else { + ( + (fraction | (1_u64 << 52)) as u128, + biased_exponent - 1023 - 52, + ) + } +} + +/// Compare two positive values represented as `significand × 2^exponent`. +#[inline] +const fn compare_binary_magnitudes( + left_significand: u128, + left_exponent: i64, + right_significand: u128, + right_exponent: i64, +) -> i8 { + let left_zeros = left_significand.trailing_zeros() as i64; + let right_zeros = right_significand.trailing_zeros() as i64; + let normalized_left = left_significand >> left_zeros.cast_unsigned(); + let normalized_right = right_significand >> right_zeros.cast_unsigned(); + let normalized_left_exponent = left_exponent + left_zeros; + let normalized_right_exponent = right_exponent + right_zeros; + + let left_top = + normalized_left_exponent + (u128::BITS - normalized_left.leading_zeros() - 1) as i64; + let right_top = + normalized_right_exponent + (u128::BITS - normalized_right.leading_zeros() - 1) as i64; + if left_top < right_top { + return -1; + } + if left_top > right_top { + return 1; + } + + let common_exponent = if normalized_left_exponent < normalized_right_exponent { + normalized_left_exponent + } else { + normalized_right_exponent + }; + let aligned_left = + normalized_left << (normalized_left_exponent - common_exponent).cast_unsigned(); + let aligned_right = + normalized_right << (normalized_right_exponent - common_exponent).cast_unsigned(); + if aligned_left < aligned_right { + -1 + } else if aligned_left > aligned_right { + 1 + } else { + 0 + } +} + +/// Compare the exact-real product `left × right` with its rounded result. +#[inline] +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 }; + } + + let (left_significand, left_exponent) = decompose_magnitude(left); + let (right_significand, right_exponent) = decompose_magnitude(right); + let exact_significand = left_significand * right_significand; + let exact_exponent = left_exponent + right_exponent; + let (rounded_significand, rounded_exponent) = decompose_magnitude(rounded); + let magnitude_relation = compare_binary_magnitudes( + exact_significand, + exact_exponent, + rounded_significand, + rounded_exponent, + ); + + if negative { + -magnitude_relation + } else { + magnitude_relation + } +} + +#[cfg(test)] +mod tests { + use super::compare_binary_magnitudes; + + #[test] + fn binary_magnitude_comparison_orders_distinct_top_exponents() { + assert_eq!(compare_binary_magnitudes(1, 1, 1, 0), 1); + assert_eq!(compare_binary_magnitudes(1, 0, 1, 1), -1); + } +} diff --git a/src/vector.rs b/src/vector.rs index ce6d18f..ff565d6 100644 --- a/src/vector.rs +++ b/src/vector.rs @@ -4,7 +4,8 @@ use core::hint::cold_path; -use crate::interval::compare_product_with_rounded; +use crate::norm::norm_near_overflow; +use crate::rounding::{compare_product_with_rounded, two_sum_error}; use crate::{ArithmeticOperation, LaError}; /// A scalar estimate paired with a certified absolute error bound. @@ -134,18 +135,6 @@ 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 -/// 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; - 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. /// /// The rounded estimate continues accumulating after `proof_available` becomes @@ -684,7 +673,7 @@ impl Vector { /// Squared Euclidean norm. /// - /// This is computed as `dot(self, self)`, so `norm2_sq` has the same + /// This is computed as `dot(self, self)`, so `norm_squared` has the same /// `f64` [`mul_add`](f64::mul_add) accumulation behavior as [`dot`](Self::dot). /// Intermediate rounding occurs, and this method does not provide a /// certified absolute rounding bound for the returned squared norm. @@ -697,7 +686,7 @@ impl Vector { /// /// # fn main() -> Result<(), LaError> { /// let v = Vector::<3>::try_new([1.0, 2.0, 3.0])?; - /// assert!((v.norm2_sq()? - 14.0).abs() <= 1e-12); + /// assert!((v.norm_squared()? - 14.0).abs() <= 1e-12); /// # Ok(()) /// # } /// ``` @@ -706,9 +695,88 @@ impl Vector { /// Returns [`LaError::NonFinite`] when the accumulated norm overflows to NaN /// or infinity. #[inline] - pub const fn norm2_sq(&self) -> Result { + pub const fn norm_squared(&self) -> Result { self.dot_with_operation(self, ArithmeticOperation::VectorSquaredNorm) } + + /// Overflow- and underflow-safe Euclidean norm. + /// + /// This computes `sqrt(Σᵢ self[i]²)` with a deterministic left-to-right + /// scaled sum-of-squares recurrence. Each non-zero magnitude is divided by + /// the largest magnitude seen so far before it is squared, so intermediate + /// squares cannot overflow and an all-subnormal vector is scaled into the + /// normal range. See `REFERENCES.md` \[15\]. + /// + /// The divisions, fused multiply-adds, square root, and final rescaling are + /// rounded in binary64. This method does not claim correct rounding or + /// provide a certified absolute error bound. Because [`Vector`] entries are + /// finite by construction, it returns a finite non-negative result unless + /// the exact Euclidean norm rounds outside the finite binary64 range. + /// Near that boundary, a fixed-size stack accumulator sums the coordinate + /// squares exactly and compares squared rounding midpoints. This fallback + /// prevents accumulated roundoff from causing or hiding overflow and does + /// not require the `exact` feature. + /// + /// Unlike [`norm_squared`](Self::norm_squared), this method does not require the + /// squared norm to be representable. For example, the norm of + /// `[1.0e200, 1.0e200]` is finite even though its squared norm is not. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// # fn main() -> Result<(), LaError> { + /// let ordinary = Vector::<2>::try_new([3.0, 4.0])?; + /// assert_eq!(ordinary.norm()?, 5.0); + /// assert_eq!(ordinary.norm_squared()?, 25.0); + /// + /// let large = Vector::<2>::try_new([1.0e200, 1.0e200])?; + /// assert!(large.norm()?.is_finite()); + /// assert!(large.norm_squared().is_err()); + /// # Ok(()) + /// # } + /// ``` + /// + /// # Errors + /// Returns [`LaError::NonFinite`] with + /// [`ArithmeticOperation::VectorNorm`] when the exact Euclidean norm rounds + /// to infinity under round-to-nearest, ties-to-even. + #[inline] + pub fn norm(&self) -> Result { + let mut entries = self.as_array().iter(); + // The first coordinate establishes the scale without a division or FMA. + // A zero (or absent) first coordinate preserves the empty-prefix state. + let mut scale = entries.next().copied().unwrap_or(0.0).abs(); + let mut scaled_sum = 1.0; + + for &entry in entries { + let magnitude = entry.abs(); + if magnitude == 0.0 { + continue; + } + + if scale < magnitude { + let ratio = scale / magnitude; + scaled_sum = (scaled_sum * ratio).mul_add(ratio, 1.0); + scale = magnitude; + } else { + let ratio = magnitude / scale; + scaled_sum = ratio.mul_add(ratio, scaled_sum); + } + } + + // With b = bit_length(D), D < 2^b. If scale <= 2^(1023-b), + // even the L1 upper bound D*scale is below 2^1023. Both the exact + // norm and the rounded recurrence therefore have ample range margin. + // Checking scale, rather than only the computed norm, also catches + // true overflow that rounding in the recurrence could hide. + let dimension_bits = usize::BITS - D.leading_zeros(); + let safe_scale = f64::from_bits(u64::from(2046 - dimension_bits) << 52); + if scale > safe_scale { + return norm_near_overflow(self.as_array(), scale); + } + Ok(scale * scaled_sum.sqrt()) + } } impl Default for Vector { @@ -770,7 +838,7 @@ mod tests { } #[test] - fn []() { + fn []() { // Use black_box to avoid constant-folding/inlining eliminating the actual dot loop, // which can make coverage tools report the mul_add line as uncovered. @@ -800,7 +868,7 @@ mod tests { } acc }; - let expected_norm2_sq = { + let expected_norm_squared = { let mut acc = 0.0; let mut i = 0; while i < $d { @@ -817,8 +885,8 @@ mod tests { // attribution for the loop body. let dot_fn: fn(&Vector<$d>, &Vector<$d>) -> Result = black_box(Vector::<$d>::dot); - let norm2_sq_fn: fn(&Vector<$d>) -> Result = - black_box(Vector::<$d>::norm2_sq); + let norm_squared_fn: fn(&Vector<$d>) -> Result = + black_box(Vector::<$d>::norm_squared); assert_abs_diff_eq!( dot_fn(black_box(&a), black_box(&b)).unwrap(), @@ -826,8 +894,8 @@ mod tests { epsilon = 1e-14 ); assert_abs_diff_eq!( - norm2_sq_fn(black_box(&a)).unwrap(), - expected_norm2_sq, + norm_squared_fn(black_box(&a)).unwrap(), + expected_norm_squared, epsilon = 1e-14 ); } @@ -902,7 +970,7 @@ mod tests { } #[test] - fn []() { + fn []() { let mut a_arr = [1.0f64; $d]; a_arr[0] = f64::MAX; let a = Vector::<$d>::new(a_arr); @@ -933,7 +1001,7 @@ mod tests { )) ); assert_eq!( - a.norm2_sq(), + a.norm_squared(), Err(LaError::non_finite_computation_step( ArithmeticOperation::VectorSquaredNorm, 0, @@ -951,12 +1019,50 @@ mod tests { gen_vector_tests!(3); gen_vector_tests!(4); gen_vector_tests!(5); + gen_vector_tests!(6); + gen_vector_tests!(7); + gen_vector_tests!(8); + + fn known_norm_input() -> ([f64; D], f64) { + let mut data = [0.0; D]; + if D == 1 { + data[0] = -5.0; + } else if D >= 2 { + data[0] = -3.0; + data[1] = 4.0; + } + (data, if D == 0 { 0.0 } else { 5.0 }) + } + + macro_rules! gen_vector_norm_known_answer_tests { + ($d:literal) => { + paste! { + #[test] + fn []() { + let (data, expected) = known_norm_input::<$d>(); + let vector = Vector::<$d>::new(data); + + assert_eq!(vector.norm(), Ok(expected)); + } + } + }; + } + + gen_vector_norm_known_answer_tests!(0); + gen_vector_norm_known_answer_tests!(1); + gen_vector_norm_known_answer_tests!(2); + gen_vector_norm_known_answer_tests!(3); + gen_vector_norm_known_answer_tests!(4); + gen_vector_norm_known_answer_tests!(5); + gen_vector_norm_known_answer_tests!(6); + gen_vector_norm_known_answer_tests!(7); + gen_vector_norm_known_answer_tests!(8); macro_rules! gen_vector_replay_tests { ($d:literal) => { paste! { #[test] - fn []() { + fn []() { let mut dot_lhs = [1.0f64; $d]; dot_lhs[$d - 1] = f64::MAX; let mut dot_rhs = [1.0f64; $d]; @@ -977,7 +1083,7 @@ mod tests { let vector = Vector::<$d>::new(norm_data); assert_eq!( - vector.norm2_sq(), + vector.norm_squared(), Err(LaError::non_finite_computation_step( ArithmeticOperation::VectorSquaredNorm, $d - 1, @@ -994,17 +1100,17 @@ mod tests { gen_vector_replay_tests!(5); macro_rules! gen_vector_const_eval_tests { - ($d:literal, $dot:literal, $norm2_sq:literal) => { + ($d:literal, $dot:literal, $norm_squared:literal) => { paste! { #[test] - fn []() { + fn []() { const DOT: Result = Vector::<$d>::new([1.0; $d]) .dot(&Vector::<$d>::new([2.0; $d])); - const NORM2_SQ: Result = - Vector::<$d>::new([1.0; $d]).norm2_sq(); + const NORM_SQUARED: Result = + Vector::<$d>::new([1.0; $d]).norm_squared(); assert_eq!(DOT, Ok($dot)); - assert_eq!(NORM2_SQ, Ok($norm2_sq)); + assert_eq!(NORM_SQUARED, Ok($norm_squared)); } } }; @@ -1016,10 +1122,10 @@ mod tests { gen_vector_const_eval_tests!(5, 10.0, 5.0); #[test] - fn vector_dot_and_norm2_sq_overflow_const_eval() { + fn vector_dot_and_norm_squared_overflow_const_eval() { const DOT: Result = Vector::<2>::new([f64::MAX; 2]).dot(&Vector::<2>::new([1.0; 2])); - const NORM2_SQ: Result = Vector::<2>::new([f64::MAX; 2]).norm2_sq(); + const NORM_SQUARED: Result = Vector::<2>::new([f64::MAX; 2]).norm_squared(); assert_eq!( DOT, @@ -1029,7 +1135,7 @@ mod tests { )) ); assert_eq!( - NORM2_SQ, + NORM_SQUARED, Err(LaError::non_finite_computation_step( ArithmeticOperation::VectorSquaredNorm, 0, @@ -1038,7 +1144,7 @@ mod tests { } #[test] - fn vector_dot_and_norm2_sq_preserve_fma_and_left_to_right_order() { + fn vector_dot_and_norm_squared_preserve_fma_and_left_to_right_order() { let dot_large = 9_007_199_254_740_992.0; let dot_lhs = Vector::<4>::new([dot_large, 1.0, 1.0, 1.0]); let dot_rhs = Vector::<4>::new([1.0; 4]); @@ -1050,7 +1156,46 @@ mod tests { let norm_large = 134_217_728.0; let vector = Vector::<4>::new([norm_large, 1.0, 1.0, 1.0]); - assert_eq!(vector.norm2_sq(), Ok(norm_large * norm_large)); + assert_eq!(vector.norm_squared(), Ok(norm_large * norm_large)); + } + + #[test] + fn vector_norm_preserves_zero_sign_and_subnormal_magnitudes() { + let signed_zero = Vector::<4>::new([-0.0, 0.0, -0.0, 0.0]); + assert_eq!(signed_zero.norm().unwrap().to_bits(), 0.0f64.to_bits()); + + let least_subnormal = f64::from_bits(1); + let subnormal = Vector::<2>::new([3.0 * least_subnormal, -4.0 * least_subnormal]); + assert_eq!( + subnormal.norm().unwrap().to_bits(), + (5.0 * least_subnormal).to_bits() + ); + } + + #[test] + fn vector_norm_handles_mixed_and_overflowing_magnitudes() { + let large = Vector::<2>::new([1.0e200, -1.0e200]); + let expected = 2.0f64.sqrt() * 1.0e200; + assert_abs_diff_eq!(large.norm().unwrap(), expected, epsilon = 2.0e184); + assert!(large.norm_squared().is_err()); + + let mixed = Vector::<4>::new([1.0e200, 1.0e-200, -f64::from_bits(1), 0.0]); + assert_eq!(mixed.norm(), Ok(1.0e200)); + + let unrepresentable = Vector::<2>::new([f64::MAX, f64::MAX]); + assert_eq!( + unrepresentable.norm(), + Err(LaError::non_finite_computation_scalar( + ArithmeticOperation::VectorNorm, + )) + ); + } + + #[test] + fn vector_norm_accepts_largest_finite_norm() { + let maximum = Vector::<2>::new([f64::MAX, 0.0]); + + assert_eq!(maximum.norm(), Ok(f64::MAX)); } #[test] @@ -1212,7 +1357,7 @@ mod tests { } #[test] - fn vector_dot_and_norm2_sq_report_first_middle_overflowing_step() { + fn vector_dot_and_norm_squared_report_first_middle_overflowing_step() { let dot_lhs = Vector::<3>::new([f64::MAX, f64::MAX, 1.0]); let dot_rhs = Vector::<3>::new([1.0; 3]); assert_eq!( @@ -1226,7 +1371,7 @@ mod tests { let norm_large = 1.0e154; let vector = Vector::<3>::new([norm_large, norm_large, 1.0]); assert_eq!( - vector.norm2_sq(), + vector.norm_squared(), Err(LaError::non_finite_computation_step( ArithmeticOperation::VectorSquaredNorm, 1, @@ -1248,7 +1393,7 @@ mod tests { .unwrap() .unwrap(); assert_abs_diff_eq!(difference_bound.absolute_error_bound(), 0.0, epsilon = 0.0); - assert_eq!(vector.norm2_sq(), Ok(0.0)); + assert_eq!(vector.norm_squared(), Ok(0.0)); } #[test] diff --git a/tests/proptest_exact.rs b/tests/proptest_exact.rs index a3fc1ae..fb026c4 100644 --- a/tests/proptest_exact.rs +++ b/tests/proptest_exact.rs @@ -157,6 +157,39 @@ fn big_rational_dot(left: &[f64; D], right: &[f64; D]) -> BigRat sum } +/// Check an approximate norm against the exact rational squared norm. +/// +/// Squaring nearby binary64 values avoids using a floating-point square root in +/// the oracle. Two steps of latitude account for the exact norm lying on either +/// side of its nearest binary64 value while still tightly constraining the +/// production result in the subnormal and normal regimes. +fn norm_nearby_values_bracket_exact_square(values: &[f64; D]) -> bool { + let Ok(norm) = Vector::::try_new(*values).and_then(|vector| vector.norm()) else { + return false; + }; + let exact_square = big_rational_dot(values, values); + let lower = if norm == 0.0 { + 0.0 + } else { + norm.next_down().next_down() + }; + let upper = norm.next_up().next_up(); + if !upper.is_finite() { + return false; + } + let lower = BigRational::from_f64(lower).expect("finite lower norm neighbor"); + let upper = BigRational::from_f64(upper).expect("finite upper norm neighbor"); + + &lower * &lower <= exact_square && exact_square <= &upper * &upper +} + +#[test] +fn norm_exact_square_oracle_preserves_mixed_scale_rounding_regression() { + let values = [1.0, 0.0, 1.0, 1.0, 0.0, 3.5, 1.0, 1.0]; + + assert!(norm_nearby_values_bracket_exact_square(&values)); +} + /// Evaluate `axis · (left - right)` without rounding coordinate differences. fn big_rational_dot_difference( axis: &[f64; D], @@ -730,6 +763,45 @@ gen_dot_errbound_oracle_proptests!(3); gen_dot_errbound_oracle_proptests!(4); gen_dot_errbound_oracle_proptests!(5); +/// The scaled Euclidean norm must tightly track an independently assembled +/// exact-rational sum of squares. The generator spans signed zero, subnormals, +/// ordinary values, and large finite magnitudes; its `f64::MAX / 4` ceiling +/// keeps the norm finite through D=8. +macro_rules! gen_norm_exact_square_oracle_proptests { + ($d:literal) => { + paste! { + proptest! { + #![proptest_config(with_default_cases(32))] + + #[test] + fn []( + values in array::[](mixed_scale_finite_f64()), + ) { + prop_assert!( + norm_nearby_values_bracket_exact_square(&values), + "D={} norm did not bracket the exact rational squared norm for {values:?}", + $d, + ); + } + } + } + }; +} + +gen_norm_exact_square_oracle_proptests!(1); +gen_norm_exact_square_oracle_proptests!(2); +gen_norm_exact_square_oracle_proptests!(3); +gen_norm_exact_square_oracle_proptests!(4); +gen_norm_exact_square_oracle_proptests!(5); +gen_norm_exact_square_oracle_proptests!(6); +gen_norm_exact_square_oracle_proptests!(7); +gen_norm_exact_square_oracle_proptests!(8); + +#[test] +fn norm_exact_square_oracle_covers_zero_dimension() { + assert!(norm_nearby_values_bracket_exact_square(&[])); +} + /// 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 { diff --git a/tests/proptest_matrix.rs b/tests/proptest_matrix.rs index 7f4aa03..ea73ac2 100644 --- a/tests/proptest_matrix.rs +++ b/tests/proptest_matrix.rs @@ -121,7 +121,7 @@ macro_rules! gen_matrix_proptests { } #[test] - fn []( + fn []( rows in array::[]( array::[](small_f64()), ), @@ -133,7 +133,7 @@ macro_rules! gen_matrix_proptests { .map(|row| row.iter().map(|&x| x.abs()).sum::()) .fold(0.0f64, f64::max); - let actual = m.inf_norm().unwrap(); + let actual = m.norm_inf().unwrap(); assert_abs_diff_eq!(actual, expected, epsilon = 0.0); prop_assert!(actual >= 0.0); } @@ -210,7 +210,7 @@ fn zero_dimension_matrix_obeys_empty_product_and_bounds_contracts() { .. }) )); - assert_eq!(matrix.inf_norm(), Ok(0.0)); + assert_eq!(matrix.norm_inf(), Ok(0.0)); assert_eq!(matrix.det(), Ok(1.0)); assert!( matrix diff --git a/tests/proptest_vector.rs b/tests/proptest_vector.rs index 2a10689..a8f8d57 100644 --- a/tests/proptest_vector.rs +++ b/tests/proptest_vector.rs @@ -39,7 +39,7 @@ macro_rules! gen_vector_proptests { } #[test] - fn []( + fn []( a_arr in array::[](small_f64()), b_arr in array::[](small_f64()), ) { @@ -51,11 +51,23 @@ macro_rules! gen_vector_proptests { assert_abs_diff_eq!(dot_ab, dot_reversed, epsilon = 1e-14); let dot_aa = a.dot(&a).unwrap(); - let norm2_sq = a.norm2_sq().unwrap(); - assert_abs_diff_eq!(norm2_sq, dot_aa, epsilon = 0.0); + let norm_squared = a.norm_squared().unwrap(); + assert_abs_diff_eq!(norm_squared, dot_aa, epsilon = 0.0); // Squared norm is always non-negative for finite inputs. - prop_assert!(norm2_sq >= 0.0); + prop_assert!(norm_squared >= 0.0); + + // The safe norm agrees with an independently ordered chain + // of binary64 hypot operations on this moderate domain. + let hypot_norm = a_arr + .iter() + .fold(0.0f64, |accumulator, &value| accumulator.hypot(value)); + let norm = a.norm().unwrap(); + assert_abs_diff_eq!(norm, hypot_norm, epsilon = 1e-12); + prop_assert!(norm >= 0.0 && norm.is_finite()); + + let negated = Vector::<$d>::try_new(a_arr.map(|value| -value)).unwrap(); + prop_assert_eq!(norm.to_bits(), negated.norm().unwrap().to_bits()); // Dot with zero vector is zero. let z = Vector::<$d>::zero(); @@ -72,3 +84,12 @@ gen_vector_proptests!(2); gen_vector_proptests!(3); gen_vector_proptests!(4); gen_vector_proptests!(5); +gen_vector_proptests!(6); +gen_vector_proptests!(7); +gen_vector_proptests!(8); + +#[test] +fn vector_norm_zero_dimension_is_always_positive_zero() { + let norm = Vector::<0>::zero().norm().unwrap(); + assert_eq!(norm.to_bits(), 0.0f64.to_bits()); +} diff --git a/tests/regressions.rs b/tests/regressions.rs index f4ea74c..b017126 100644 --- a/tests/regressions.rs +++ b/tests/regressions.rs @@ -5,6 +5,246 @@ use la_stack::ERR_COEFF_3; use la_stack::prelude::*; +#[cfg(feature = "exact")] +use proptest::prelude::*; + +#[cfg(feature = "exact")] +#[path = "common/proptest_config.rs"] +mod proptest_config; + +/// Pad the independently identified near-overflow pair to a fixed dimension. +fn norm_boundary_vector(left: f64, right: f64) -> Vector { + let mut values = [0.0; D]; + values[0] = left; + values[1] = right; + Vector::try_new(values).expect("the boundary fixture is finite") +} + +/// Preserve finite answers that the scaled recurrence previously overflowed. +fn assert_norm_boundary_regressions() { + // Expected bits are independently checked with exact squared midpoints in + // the exact-feature tests below, rather than another floating norm kernel. + for (left_bits, right_bits, expected_bits) in [ + ( + 0x7f55_4c98_5f06_f693, + 0x7fef_fffe_3a58_0905, + 0x7fef_ffff_ffff_ffff, + ), + ( + 0x7f61_3404_ea4a_8c14, + 0x7fef_fffb_6032_c601, + 0x7fef_ffff_ffff_fffe, + ), + ( + 0x7f64_7ae1_47ae_147a, + 0x7fef_fff9_7246_996b, + 0x7fef_ffff_ffff_ffff, + ), + ( + 0x7f69_652b_d3c3_6112, + 0x7fef_fff5_ec54_3df0, + 0x7fef_ffff_ffff_ffff, + ), + ] { + let left = f64::from_bits(left_bits); + let right = f64::from_bits(right_bits); + for (left, right) in [(left, right), (right, left), (-left, right), (right, -left)] { + let vector = norm_boundary_vector::(left, right); + #[cfg(feature = "exact")] + assert_high_range_norm_rounding(&vector); + assert_eq!( + vector.norm().unwrap().to_bits(), + expected_bits, + "boundary pair ({left_bits:#x}, {right_bits:#x})" + ); + } + } + + let two_997 = f64::from_bits(2020_u64 << 52); + let two_998 = f64::from_bits(2021_u64 << 52); + let finite = norm_boundary_vector::(f64::MAX, two_997); + assert_eq!(finite.norm(), Ok(f64::MAX)); + #[cfg(feature = "exact")] + assert_high_range_norm_rounding(&finite); + + // The old rounded sqrt could equal 1, hiding true overflow in this case. + for (left, right) in [(f64::MAX, two_998), (f64::MAX, f64::MAX)] { + let vector = norm_boundary_vector::(left, right); + assert_eq!( + vector.norm(), + Err(LaError::non_finite_computation_scalar( + ArithmeticOperation::VectorNorm + )), + ); + #[cfg(feature = "exact")] + assert_high_range_norm_rounding(&vector); + } +} + +fn assert_norm_guard_transition() { + let dimension_bits = usize::BITS - D.leading_zeros(); + let boundary = f64::from_bits(u64::from(2046 - dimension_bits) << 52); + + // Single-coordinate norms are exact on either side of the conservative + // dispatch boundary, regardless of coordinate position or sign. + for magnitude in [boundary.next_down(), boundary, boundary.next_up()] { + for index in 0..D { + for sign in [-1.0, 1.0] { + let mut values = [0.0; D]; + values[index] = sign * magnitude; + let vector = Vector::try_new(values).unwrap(); + assert_eq!(vector.norm(), Ok(magnitude)); + #[cfg(feature = "exact")] + assert_high_range_norm_rounding(&vector); + } + } + } + + // A 3-4-5 triple has its largest coordinate at the boundary but its norm + // above it. Both the recurrence and the exact path must preserve it. + let unit = boundary / 4.0; + for (left, right) in [(3.0, 4.0), (4.0, 3.0), (-3.0, 4.0), (4.0, -3.0)] { + let vector = norm_boundary_vector::(left * unit, right * unit); + assert_eq!(vector.norm(), Ok(5.0 * unit)); + #[cfg(feature = "exact")] + assert_high_range_norm_rounding(&vector); + } +} + +#[test] +fn norm_returns_exact_high_range_pythagorean_result() { + let unit = f64::from_bits(2043_u64 << 52); // 2^1020 + let vector = Vector::<2>::try_new([3.0 * unit, 4.0 * unit]).unwrap(); + + assert_eq!(vector.norm(), Ok(5.0 * unit)); +} + +macro_rules! gen_norm_boundary_regressions { + ($d:literal) => { + pastey::paste! { + #[test] + fn []() { + assert_norm_boundary_regressions::<$d>(); + } + + #[test] + fn []() { + assert_norm_guard_transition::<$d>(); + } + + #[cfg(feature = "exact")] + proptest! { + #![proptest_config(proptest_config::with_default_cases(64))] + + #[test] + fn []( + bits in any::<[u64; $d]>(), + ) { + // Mix upper-range and arbitrary finite coordinates. The + // first coordinate always selects the exact fallback; + // other coordinates exercise widely separated squares, + // dense carry propagation, signs, and true overflow. + let values = std::array::from_fn(|index| { + let raw = bits[index]; + let exponent = if index == 0 { + 2046 + } else if raw & 1 == 0 { + 2043 + ((raw >> 52) % 4) + } else { + (raw >> 52) % 2047 + }; + f64::from_bits((raw & 0x800f_ffff_ffff_ffff) | (exponent << 52)) + }); + let vector = Vector::<$d>::try_new(values).unwrap(); + assert_high_range_norm_rounding(&vector); + } + } + } + }; +} + +gen_norm_boundary_regressions!(2); +gen_norm_boundary_regressions!(3); +gen_norm_boundary_regressions!(4); +gen_norm_boundary_regressions!(5); +gen_norm_boundary_regressions!(6); +gen_norm_boundary_regressions!(7); +gen_norm_boundary_regressions!(8); + +/// Validate a high-range result by squaring exact rational rounding midpoints. +#[cfg(feature = "exact")] +fn assert_high_range_norm_rounding(vector: &Vector) { + let exact = |value| BigRational::from_f64(value).expect("finite oracle input"); + let square: BigRational = vector + .as_array() + .iter() + .map(|&value| exact(value).pow(2)) + .sum(); + let overflow_midpoint = exact(f64::MAX) + BigRational::from_integer(BigInt::from(1) << 970); + match vector.norm() { + Ok(norm) => { + assert!(square < overflow_midpoint.pow(2)); + let lower_midpoint = (exact(norm.next_down()) + exact(norm)) / exact(2.0); + let upper_midpoint = if norm.to_bits() == f64::MAX.to_bits() { + overflow_midpoint + } else { + (exact(norm) + exact(norm.next_up())) / exact(2.0) + }; + let lower_square = lower_midpoint.pow(2); + let upper_square = upper_midpoint.pow(2); + assert!(lower_square <= square && square <= upper_square); + if square == lower_square || square == upper_square { + assert_eq!( + norm.to_bits() & 1, + 0, + "ties must select the even significand" + ); + } + } + Err(error) => { + assert_eq!( + error, + LaError::non_finite_computation_scalar(ArithmeticOperation::VectorNorm) + ); + assert!( + square >= overflow_midpoint.pow(2), + "finite norm misclassified as overflow" + ); + } + } +} + +#[test] +fn norm_upper_range_ties_and_subnormal_tail() { + let unit = f64::from_bits(1993_u64 << 52); // 2^970 + for offset in [1.0, 3.0] { + let k = 2.0_f64.powi(51) + offset; + // A scaled 3-4-5 triple has exact norm 5*k*2^970. Odd 54-bit 5*k + // lies exactly halfway between two binary64 values; the two offsets + // exercise both directions of ties-to-even. + let vector = norm_boundary_vector::<3>((3.0 * k) * unit, (4.0 * k) * unit); + let expected = (5.0 * k) * unit; + assert_eq!(vector.norm().unwrap().to_bits(), expected.to_bits()); + #[cfg(feature = "exact")] + assert_high_range_norm_rounding(&vector); + + let mut values = vector.into_array(); + values[2] = f64::from_bits(1); + let with_tail = Vector::try_new(values).unwrap(); + let expected_with_tail = if offset < 2.0 { + expected.next_up() + } else { + expected + }; + assert_eq!( + with_tail.norm().unwrap().to_bits(), + expected_with_tail.to_bits() + ); + #[cfg(feature = "exact")] + assert_high_range_norm_rounding(&with_tail); + } +} + #[test] #[cfg(feature = "exact")] fn det_exact_f64_preserves_min_positive_subnormal() -> Result<(), LaError> { diff --git a/tests/vs_linalg_inputs.rs b/tests/vs_linalg_inputs.rs index 3848771..5e5aca9 100644 --- a/tests/vs_linalg_inputs.rs +++ b/tests/vs_linalg_inputs.rs @@ -20,9 +20,15 @@ pub mod vs_linalg_common; #[cfg(not(la_stack_v0_4_3_api))] use vs_linalg_common::make_balanced_dynamic_range_rows; use vs_linalg_common::{ - PreparedFaerLuDet, faer_det_from_ldlt, faer_perm_sign, la_stack_dot, la_stack_tolerance, - make_ill_conditioned_matrix_rows, make_matrix_rows, make_pivoting_matrix_rows, - make_vector_array, matrix_entry, nalgebra_inf_norm, vector_entry, + PreparedFaerLuDet, faer_det_from_ldlt, faer_perm_sign, la_stack_dot, la_stack_norm_inf, + la_stack_norm_squared, la_stack_tolerance, make_ill_conditioned_matrix_rows, make_matrix_rows, + make_pivoting_matrix_rows, make_vector_array, matrix_entry, nalgebra_inf_norm, vector_entry, +}; +#[cfg(not(any(la_stack_pre_rational_input_api, la_stack_v0_4_3_api)))] +use vs_linalg_common::{ + delaunay_scaled_norm, iterative_hypot, make_norm_descending_array, + make_norm_repeated_scale_array, make_norm_sparse_array, make_norm_wide_dynamic_range_array, + norm_scenarios, }; /// Assert scalar agreement with a tolerance that scales for larger magnitudes. @@ -169,7 +175,7 @@ fn assert_ldlt_agreement() { ); } -/// Check vector dot-product and squared-norm agreement for one benchmark dimension. +/// Check vector dot-product and norm agreement for one benchmark dimension. fn assert_vector_operation_agreement() { let v1 = Vector::::try_new(make_vector_array::(0.0)) .unwrap_or_else(|err| panic!("la_stack vector construction failed: {err}")); @@ -188,12 +194,67 @@ fn assert_vector_operation_agreement() { } assert_close("faer_dot", fa_dot, la_dot); - let la_norm2_sq = v1 - .norm2_sq() - .unwrap_or_else(|err| panic!("la_stack norm2_sq failed: {err}")); - assert_close("nalgebra_norm_squared", nv1.norm_squared(), la_norm2_sq); - let fa_norm2_sq = fv1.as_mat_ref().squared_norm_l2(); - assert_close("faer_norm2_sq", fa_norm2_sq, la_norm2_sq); + let la_norm_squared = la_stack_norm_squared(&v1) + .unwrap_or_else(|err| panic!("la_stack norm_squared failed: {err}")); + assert_close("nalgebra_norm_squared", nv1.norm_squared(), la_norm_squared); + let fa_norm_squared = fv1.as_mat_ref().squared_norm_l2(); + assert_close("faer_norm2_sq", fa_norm_squared, la_norm_squared); + + #[cfg(not(any(la_stack_pre_rational_input_api, la_stack_v0_4_3_api)))] + { + let la_norm = v1 + .norm() + .unwrap_or_else(|err| panic!("la_stack norm failed: {err}")); + assert_close( + "iterative_f64_hypot", + iterative_hypot(v1.as_array()), + la_norm, + ); + assert_close( + "delaunay_scaled_norm", + delaunay_scaled_norm(v1.as_array()), + la_norm, + ); + assert_close("nalgebra_norm", nv1.norm(), la_norm); + assert_close("faer_norm_l2", fv1.as_mat_ref().norm_l2(), la_norm); + + for (scenario, values) in norm_scenarios::() { + let vector = Vector::::try_new(values).unwrap_or_else(|err| { + panic!("la_stack {scenario} norm scenario construction failed: {err}") + }); + let expected = vector + .norm() + .unwrap_or_else(|err| panic!("la_stack {scenario} norm failed: {err}")); + assert_close(scenario, iterative_hypot(&values), expected); + assert_close(scenario, delaunay_scaled_norm(&values), expected); + } + } +} + +#[cfg(not(any(la_stack_pre_rational_input_api, la_stack_v0_4_3_api)))] +#[test] +fn norm_scenario_inputs_cover_distinct_branch_and_range_profiles() { + let descending = make_norm_descending_array::<8>(); + assert!( + descending + .windows(2) + .all(|pair| pair[0].abs() > pair[1].abs()) + ); + + let repeated = make_norm_repeated_scale_array::<8>(); + assert!( + repeated + .iter() + .all(|entry| entry.abs().to_bits() == 3.0f64.to_bits()) + ); + + let sparse = make_norm_sparse_array::<8>(); + assert_eq!(sparse.iter().filter(|entry| **entry != 0.0).count(), 1); + + let wide = make_norm_wide_dynamic_range_array::<8>(); + assert!(wide.iter().any(|entry| entry.abs() >= 1.0e200)); + assert!(wide.iter().any(|entry| entry.is_subnormal())); + assert!(wide.contains(&0.0)); } #[test] @@ -340,15 +401,14 @@ fn stress_inputs_exercise_pivoting_conditioning_and_scaled_products() { } /// Check matrix infinity-norm agreement for one benchmark dimension. -fn assert_matrix_inf_norm_agreement() { +fn assert_matrix_norm_inf_agreement() { let a = Matrix::::try_from_rows(make_matrix_rows::()) .unwrap_or_else(|err| panic!("la_stack matrix construction failed: {err}")); let na = SMatrix::::from_fn(matrix_entry::); let fa = Mat::::from_fn(D, D, matrix_entry::); - let la_norm = a - .inf_norm() - .unwrap_or_else(|err| panic!("la_stack inf_norm failed: {err}")); + let la_norm = + la_stack_norm_inf(&a).unwrap_or_else(|err| panic!("la_stack norm_inf failed: {err}")); assert_close("nalgebra_inf_norm", nalgebra_inf_norm(&na), la_norm); let mut fa_norm = 0.0; for r in 0..D { @@ -383,7 +443,7 @@ macro_rules! gen_smoke_tests { #[test] fn $norm() { - assert_matrix_inf_norm_agreement::<$d>(); + assert_matrix_norm_inf_agreement::<$d>(); } }; } @@ -393,54 +453,54 @@ gen_smoke_tests!( vs_linalg_lu_agrees_2d, vs_linalg_ldlt_agrees_2d, vs_linalg_vector_operations_agree_2d, - vs_linalg_matrix_inf_norm_agrees_2d + vs_linalg_matrix_norm_inf_agrees_2d ); gen_smoke_tests!( 3, vs_linalg_lu_agrees_3d, vs_linalg_ldlt_agrees_3d, vs_linalg_vector_operations_agree_3d, - vs_linalg_matrix_inf_norm_agrees_3d + vs_linalg_matrix_norm_inf_agrees_3d ); gen_smoke_tests!( 4, vs_linalg_lu_agrees_4d, vs_linalg_ldlt_agrees_4d, vs_linalg_vector_operations_agree_4d, - vs_linalg_matrix_inf_norm_agrees_4d + vs_linalg_matrix_norm_inf_agrees_4d ); gen_smoke_tests!( 5, vs_linalg_lu_agrees_5d, vs_linalg_ldlt_agrees_5d, vs_linalg_vector_operations_agree_5d, - vs_linalg_matrix_inf_norm_agrees_5d + vs_linalg_matrix_norm_inf_agrees_5d ); gen_smoke_tests!( 8, vs_linalg_lu_agrees_8d, vs_linalg_ldlt_agrees_8d, vs_linalg_vector_operations_agree_8d, - vs_linalg_matrix_inf_norm_agrees_8d + vs_linalg_matrix_norm_inf_agrees_8d ); gen_smoke_tests!( 16, vs_linalg_lu_agrees_16d, vs_linalg_ldlt_agrees_16d, vs_linalg_vector_operations_agree_16d, - vs_linalg_matrix_inf_norm_agrees_16d + vs_linalg_matrix_norm_inf_agrees_16d ); gen_smoke_tests!( 32, vs_linalg_lu_agrees_32d, vs_linalg_ldlt_agrees_32d, vs_linalg_vector_operations_agree_32d, - vs_linalg_matrix_inf_norm_agrees_32d + vs_linalg_matrix_norm_inf_agrees_32d ); gen_smoke_tests!( 64, vs_linalg_lu_agrees_64d, vs_linalg_ldlt_agrees_64d, vs_linalg_vector_operations_agree_64d, - vs_linalg_matrix_inf_norm_agrees_64d + vs_linalg_matrix_norm_inf_agrees_64d );