feat(vector): add overflow-safe Euclidean norms

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<const D: usize>` (`[f64; D]`)
-  - `src/matrix.rs`: `Matrix<const D: usize>` (`[[f64; D]; D]`) + helpers (`get`, `try_get`, `set`, `inf_norm`, `det`, `det_direct`)
+  - `src/matrix.rs`: `Matrix<const D: usize>` (`[[f64; D]; D]`) + helpers (`get`, `try_get`, `set`, `norm_inf`, `det`, `det_direct`)
   - `src/lu.rs`: `Lu<const D: usize>` factorization with partial pivoting (`solve`, `det`)
   - `src/ldlt.rs`: `Ldlt<const D: usize>` factorization without pivoting for exactly
     symmetric positive-definite matrices (`solve`, `det`)

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<D>` | `[f64; D]` | Finite fixed-length vector for input and computation | `try_new`, `as_array`, `into_array`, `dot`, `dot_with_errbound`, `dot_difference_with_errbound`, `norm2_sq` |
+| `Vector<D>` | `[f64; D]` | Finite fixed-length vector for input and computation | `try_new`, `as_array`, `into_array`, `dot`, `dot_with_errbound`, `dot_difference_with_errbound`, `norm`, `norm_squared` |
 | `Matrix<D>` | `[[f64; D]; D]` | Finite square matrix for input and computation | See below |
 | `Interval` | Two finite ordered `f64` bounds | Outward-rounded exact-real enclosure | `try_new`, `point`, `try_from_subtraction`, `try_add`, `try_mul`, `negate`, `try_square` |
 | `IntervalMatrix<D>` | `[[Interval; D]; D]` | Division-free determinant enclosure and sign proof through D=7 | `from_rows`, `try_from_point_rows`, `from_matrix`, `det`, `det_sign` |

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)

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<const D: usize>(left: &Vector<D>, right: &Vector<D>) -> 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<const D: usize>(vector: &Vector<D>) -> Result<f64, LaError> {
+    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<const D: usize>(vector: &Vector<D>) -> Result<f64, LaError> {
+    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<const D: usize>(matrix: &Matrix<D>) -> Result<f64, LaError> {
+    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<const D: usize>(matrix: &Matrix<D>) -> Result<f64, LaError> {
+    matrix.inf_norm()
+}
+
 /// Parse a tolerance through the constructor exposed by the selected library
 /// revision.
 ///
@@ -256,6 +304,111 @@ pub fn make_vector_array<const D: usize>(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<const D: usize>() -> [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<const D: usize>() -> [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<const D: usize>() -> [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<const D: usize>() -> [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<const D: usize>() -> [(&'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<const D: usize>(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<const D: usize>(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]