num_valid/algorithms/

l2_norm.rs

#![deny(rustdoc::broken_intra_doc_links)]

//! # L2 Norm (Euclidean Norm) Computation
//!
//! This module provides a numerically stable implementation of the L2 norm
//! (Euclidean norm) for slices of [`RealScalar`] values.
//!
//! ## Algorithm: BLAS/LAPACK-Style Incremental Scaling
//!
//! The implementation uses the scaled sum of squares approach, which is the
//! same algorithm used in modern BLAS/LAPACK implementations (e.g., `DNRM2`).
//! This single-pass algorithm avoids overflow and underflow by maintaining
//! the invariant:
//!
//! $$\|x\|_2 = \text{scale} \cdot \sqrt{\text{sumsq}}$$
//!
//! where:
//! - `scale` is the maximum absolute value seen so far
//! - `sumsq` is the sum of squared normalized values: $\sum_i (|x_i| / \text{scale})^2$
//!
//! ### Key Properties
//!
//! | Property | Description |
//! |----------|-------------|
//! | **Single-pass** | O(n) time complexity with one iteration over the data |
//! | **No overflow** | All squared values are ≤ 1 (normalized by max) |
//! | **No underflow** | Scale factor preserves the magnitude of small values |
//! | **Order-independent** | Same result regardless of element ordering |
//! | **Backend-agnostic** | Works with both `f64` and `rug::Float` backends |
//!
//! ### How It Works
//!
//! 1. Track `scale` (max |xᵢ| seen) and `sumsq` (sum of normalized squares)
//! 2. For each element xᵢ:
//!    - If |xᵢ| > scale: rescale sumsq to new scale, then accumulate
//!    - Otherwise: just accumulate (|xᵢ|/scale)²
//! 3. Return scale × √sumsq
//!
//! ## Why Not Naive Implementation?
//!
//! The naive approach `sqrt(sum(x²))` fails for extreme values:
//!
//! ```text
//! Naive overflow:   ||[1e200, 1e200]|| → (1e200)² = Inf → sqrt(Inf) = Inf  ❌
//! Naive underflow:  ||[1e-200, 1e-200]|| → (1e-200)² = 0 → sqrt(0) = 0    ❌
//!
//! Scaled approach:  ||[1e200, 1e200]|| → 1e200 × sqrt(2) ≈ 1.41e200        ✅
//! Scaled approach:  ||[1e-200, 1e-200]|| → 1e-200 × sqrt(2) ≈ 1.41e-200    ✅
//! ```
//!
//! ## Usage Examples
//!
//! ### Basic Usage
//!
//! ```rust
//! use num_valid::{RealNative64StrictFinite, RealScalar, algorithms::l2_norm::l2_norm};
//!
//! let data: Vec<RealNative64StrictFinite> = vec![
//!     RealNative64StrictFinite::from_f64(3.0),
//!     RealNative64StrictFinite::from_f64(4.0),
//! ];
//!
//! let norm = l2_norm(&data);
//! assert_eq!(*norm.as_ref(), 5.0); // 3² + 4² = 25, √25 = 5
//! ```
//!
//! ### Handling Extreme Values
//!
//! ```rust
//! use num_valid::{RealNative64StrictFinite, RealScalar, algorithms::l2_norm::l2_norm};
//!
//! // Values near overflow threshold - naive approach would fail
//! let large: Vec<RealNative64StrictFinite> = vec![
//!     RealNative64StrictFinite::from_f64(1e154),
//!     RealNative64StrictFinite::from_f64(1e154),
//! ];
//! let norm = l2_norm(&large);
//! // Result: √2 × 1e154 ≈ 1.41e154 (no overflow!)
//! assert!((norm.as_ref() / 1e154 - std::f64::consts::SQRT_2).abs() < 1e-10);
//!
//! // Values near underflow threshold
//! let small: Vec<RealNative64StrictFinite> = vec![
//!     RealNative64StrictFinite::from_f64(1e-154),
//!     RealNative64StrictFinite::from_f64(1e-154),
//! ];
//! let norm = l2_norm(&small);
//! // Result: √2 × 1e-154 ≈ 1.41e-154 (no underflow!)
//! assert!((norm.as_ref() / 1e-154 - std::f64::consts::SQRT_2).abs() < 1e-10);
//! ```
//!
//! ### With Arbitrary-Precision Backend
//!
//! ```rust
//! # #[cfg(feature = "rug")]
//! # fn main() -> Result<(), Box<dyn std::error::Error>> {
//! use num_valid::{RealRugStrictFinite, algorithms::l2_norm::l2_norm};
//! use try_create::TryNew;
//!
//! type R = RealRugStrictFinite<200>; // 200-bit precision
//!
//! // Values beyond f64 range (would be infinity in f64)
//! let huge = R::try_new(rug::Float::with_val(200, rug::Float::parse("1e1000")?))?;
//! let data: Vec<R> = vec![huge.clone(), huge.clone()];
//!
//! let norm = l2_norm(&data);
//! // Result: √2 × 1e1000 - computed exactly with arbitrary precision
//! # Ok(())
//! # }
//! # #[cfg(not(feature = "rug"))] fn main() {}
//! ```
//!
//! ## Edge Cases
//!
//! | Input | Result |
//! |-------|--------|
//! | Empty slice `[]` | `0` |
//! | Single element `[x]` | `|x|` |
//! | All zeros `[0, 0, 0]` | `0` |
//! | Contains zeros `[0, 3, 0, 4]` | Same as `[3, 4]` → `5` |
//! | Negative values `[-3, -4]` | Same as `[3, 4]` → `5` |
//!
//! ## Performance Characteristics
//!
//! - **Time complexity**: O(n) - single pass over data
//! - **Space complexity**: O(1) - only two accumulator variables
//! - **Operations per element**: 1 comparison, 1-2 divisions, 1 multiply-add
//!
//! The algorithm performs more divisions than the naive approach, but this
//! tradeoff is necessary for numerical stability. For performance-critical
//! code where values are known to be in a safe range, consider the naive
//! approach with appropriate bounds checking.
//!
//! ## References
//!
//! - Blue, J. L. (1978). "A Portable Fortran Program to Find the Euclidean
//!   Norm of a Vector". ACM Transactions on Mathematical Software, 4(1), 15-23.
//! - Anderson, E. (2017). "Algorithm 978: Safe Scaling in the Level 1 BLAS".
//!   ACM Transactions on Mathematical Software, 44(1), 1-28.
//! - LAPACK Working Note 148: "On Computing LAPACK's XNRM2"
//!
//! [`RealScalar`]: crate::RealScalar

use crate::RealScalar;

/// Computes the L2 norm (Euclidean norm) of a slice of real scalars.
///
/// Uses BLAS/LAPACK-style incremental scaling to prevent overflow and underflow.
/// The algorithm maintains the invariant `||x||₂ = scale × √sumsq` where all
/// accumulated squared values are normalized to the range [0, 1].
///
/// # Arguments
///
/// * `x` - A slice of [`RealScalar`] values
///
/// # Returns
///
/// The L2 norm: $\|x\|_2 = \sqrt{\sum_i x_i^2}$
///
/// # Algorithm Complexity
///
/// - **Time**: O(n) single-pass
/// - **Space**: O(1)
///
/// # Examples
///
/// ```rust
/// use num_valid::{RealNative64StrictFinite, RealScalar, algorithms::l2_norm::l2_norm};
///
/// // Pythagorean triple: 3² + 4² = 5²
/// let v = vec![
///     RealNative64StrictFinite::from_f64(3.0),
///     RealNative64StrictFinite::from_f64(4.0),
/// ];
/// assert_eq!(*l2_norm(&v).as_ref(), 5.0);
///
/// // Empty slice returns zero
/// let empty: Vec<RealNative64StrictFinite> = vec![];
/// assert_eq!(*l2_norm(&empty).as_ref(), 0.0);
///
/// // Works with extreme values (no overflow)
/// let large = vec![
///     RealNative64StrictFinite::from_f64(1e154),
///     RealNative64StrictFinite::from_f64(1e154),
/// ];
/// assert!(l2_norm(&large).as_ref().is_finite());
/// ```
pub fn l2_norm<T: RealScalar>(x: &[T]) -> T {
    let mut scale = T::zero();
    let mut sumsq = T::zero();

    let zero = T::zero();

    for xi in x.iter().cloned() {
        let abs_xi = xi.abs();

        if abs_xi == zero {
            continue;
        }

        if scale < abs_xi {
            // Rescale previous sum to new scale
            if scale > zero {
                let r = scale.clone() / &abs_xi;
                sumsq *= r.pow(2);
            }
            scale = abs_xi.clone();
        }

        // Always accumulate (key insight!)
        let r = abs_xi.clone() / &scale;
        sumsq += r.pow(2);
    }

    if scale == zero {
        zero
    } else {
        scale * sumsq.sqrt()
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::RealNative64StrictFinite;
    use approx::assert_relative_eq;

    /// Helper to create a vector of validated reals from f64 values
    fn vec_f64(vals: &[f64]) -> Vec<RealNative64StrictFinite> {
        vals.iter()
            .map(|&v| RealNative64StrictFinite::from_f64(v))
            .collect()
    }

    mod basic_cases {
        use super::*;

        #[test]
        fn pythagorean_3_4_5() {
            let data = vec_f64(&[3.0, 4.0]);
            let norm = l2_norm(&data);
            assert_relative_eq!(*norm.as_ref(), 5.0, epsilon = 1e-15);
        }

        #[test]
        fn pythagorean_reverse_order() {
            // Same result regardless of order
            let data = vec_f64(&[4.0, 3.0]);
            let norm = l2_norm(&data);
            assert_relative_eq!(*norm.as_ref(), 5.0, epsilon = 1e-15);
        }

        #[test]
        fn unit_vector_x() {
            let data = vec_f64(&[1.0, 0.0, 0.0]);
            let norm = l2_norm(&data);
            assert_relative_eq!(*norm.as_ref(), 1.0, epsilon = 1e-15);
        }

        #[test]
        fn unit_vector_y() {
            let data = vec_f64(&[0.0, 1.0, 0.0]);
            let norm = l2_norm(&data);
            assert_relative_eq!(*norm.as_ref(), 1.0, epsilon = 1e-15);
        }

        #[test]
        fn three_equal_values() {
            // ||[1, 1, 1]|| = sqrt(3)
            let data = vec_f64(&[1.0, 1.0, 1.0]);
            let norm = l2_norm(&data);
            assert_relative_eq!(*norm.as_ref(), 3.0_f64.sqrt(), epsilon = 1e-15);
        }

        #[test]
        fn larger_vector() {
            // ||[1, 2, 3, 4, 5]|| = sqrt(1 + 4 + 9 + 16 + 25) = sqrt(55)
            let data = vec_f64(&[1.0, 2.0, 3.0, 4.0, 5.0]);
            let norm = l2_norm(&data);
            assert_relative_eq!(*norm.as_ref(), 55.0_f64.sqrt(), epsilon = 1e-14);
        }
    }

    mod edge_cases {
        use super::*;

        #[test]
        fn empty_vector() {
            let data: Vec<RealNative64StrictFinite> = vec![];
            let norm = l2_norm(&data);
            assert_eq!(*norm.as_ref(), 0.0);
        }

        #[test]
        fn single_positive_element() {
            let data = vec_f64(&[7.0]);
            let norm = l2_norm(&data);
            assert_relative_eq!(*norm.as_ref(), 7.0, epsilon = 1e-15);
        }

        #[test]
        fn single_negative_element() {
            let data = vec_f64(&[-7.0]);
            let norm = l2_norm(&data);
            assert_relative_eq!(*norm.as_ref(), 7.0, epsilon = 1e-15);
        }

        #[test]
        fn all_zeros() {
            let data = vec_f64(&[0.0, 0.0, 0.0, 0.0]);
            let norm = l2_norm(&data);
            assert_eq!(*norm.as_ref(), 0.0);
        }

        #[test]
        fn zeros_interspersed() {
            let data = vec_f64(&[0.0, 3.0, 0.0, 4.0, 0.0]);
            let norm = l2_norm(&data);
            assert_relative_eq!(*norm.as_ref(), 5.0, epsilon = 1e-15);
        }

        #[test]
        fn negative_values() {
            let data = vec_f64(&[-3.0, -4.0]);
            let norm = l2_norm(&data);
            assert_relative_eq!(*norm.as_ref(), 5.0, epsilon = 1e-15);
        }

        #[test]
        fn mixed_signs() {
            let data = vec_f64(&[-3.0, 4.0]);
            let norm = l2_norm(&data);
            assert_relative_eq!(*norm.as_ref(), 5.0, epsilon = 1e-15);
        }

        #[test]
        fn single_zero() {
            let data = vec_f64(&[0.0]);
            let norm = l2_norm(&data);
            assert_eq!(*norm.as_ref(), 0.0);
        }
    }

    mod numerical_stability {
        use super::*;

        #[test]
        fn large_values_no_overflow() {
            // Values that would overflow with naive x^2 approach
            let scale = 1e154;
            let data = vec_f64(&[3.0 * scale, 4.0 * scale]);
            let norm = l2_norm(&data);
            let expected = 5.0 * scale;
            let rel_err = (*norm.as_ref() - expected).abs() / expected;
            assert!(
                rel_err < 1e-14,
                "Large values: expected {}, got {}, rel_err {}",
                expected,
                *norm.as_ref(),
                rel_err
            );
        }

        #[test]
        fn very_large_values() {
            // Even closer to overflow
            let scale = 1e300;
            let data = vec_f64(&[scale, scale]);
            let norm = l2_norm(&data);
            let expected = scale * 2.0_f64.sqrt();
            let rel_err = (*norm.as_ref() - expected).abs() / expected;
            assert!(
                rel_err < 1e-14,
                "Very large values: expected {}, got {}, rel_err {}",
                expected,
                *norm.as_ref(),
                rel_err
            );
        }

        #[test]
        fn small_values_no_underflow() {
            // Values that would underflow with naive x^2 approach
            let scale = 1e-154;
            let data = vec_f64(&[3.0 * scale, 4.0 * scale]);
            let norm = l2_norm(&data);
            let expected = 5.0 * scale;
            let rel_err = (*norm.as_ref() - expected).abs() / expected;
            assert!(
                rel_err < 1e-14,
                "Small values: expected {}, got {}, rel_err {}",
                expected,
                *norm.as_ref(),
                rel_err
            );
        }

        #[test]
        fn very_small_values() {
            // Even closer to underflow
            let scale = 1e-300;
            let data = vec_f64(&[scale, scale]);
            let norm = l2_norm(&data);
            let expected = scale * 2.0_f64.sqrt();
            let rel_err = (*norm.as_ref() - expected).abs() / expected;
            assert!(
                rel_err < 1e-14,
                "Very small values: expected {}, got {}, rel_err {}",
                expected,
                *norm.as_ref(),
                rel_err
            );
        }

        #[test]
        fn mixed_large_and_small() {
            // Large value dominates
            let data = vec_f64(&[1e150, 1.0, 1e-150]);
            let norm = l2_norm(&data);
            // Norm is approximately 1e150 (small values negligible)
            let rel_err = (*norm.as_ref() - 1e150).abs() / 1e150;
            assert!(
                rel_err < 1e-14,
                "Mixed magnitudes: expected ~1e150, got {}, rel_err {}",
                *norm.as_ref(),
                rel_err
            );
        }

        #[test]
        fn all_same_large_values() {
            // n values of x: ||[x, x, ..., x]|| = |x| * sqrt(n)
            let x = 1e154;
            let n = 100;
            let data: Vec<RealNative64StrictFinite> = (0..n)
                .map(|_| RealNative64StrictFinite::from_f64(x))
                .collect();
            let norm = l2_norm(&data);
            let expected = x * (n as f64).sqrt();
            let rel_err = (*norm.as_ref() - expected).abs() / expected;
            assert!(
                rel_err < 1e-13,
                "100 large values: expected {}, got {}, rel_err {}",
                expected,
                *norm.as_ref(),
                rel_err
            );
        }

        #[test]
        fn all_same_small_values() {
            let x = 1e-154;
            let n = 100;
            let data: Vec<RealNative64StrictFinite> = (0..n)
                .map(|_| RealNative64StrictFinite::from_f64(x))
                .collect();
            let norm = l2_norm(&data);
            let expected = x * (n as f64).sqrt();
            let rel_err = (*norm.as_ref() - expected).abs() / expected;
            assert!(
                rel_err < 1e-13,
                "100 small values: expected {}, got {}, rel_err {}",
                expected,
                *norm.as_ref(),
                rel_err
            );
        }

        #[test]
        fn rescaling_triggered_multiple_times() {
            // Ascending order triggers rescaling at each step
            let data = vec_f64(&[1.0, 2.0, 4.0, 8.0, 16.0]);
            let expected = (1.0 + 4.0 + 16.0 + 64.0 + 256.0_f64).sqrt(); // sqrt(341)
            let norm = l2_norm(&data);
            assert_relative_eq!(*norm.as_ref(), expected, epsilon = 1e-14);
        }

        #[test]
        fn descending_order_no_rescaling() {
            // Descending order: max is found first, no rescaling needed
            let data = vec_f64(&[16.0, 8.0, 4.0, 2.0, 1.0]);
            let expected = (1.0 + 4.0 + 16.0 + 64.0 + 256.0_f64).sqrt();
            let norm = l2_norm(&data);
            assert_relative_eq!(*norm.as_ref(), expected, epsilon = 1e-14);
        }
    }

    mod special_values {
        use super::*;

        #[test]
        fn max_finite_value() {
            // Single f64::MAX should give f64::MAX
            let data = vec_f64(&[f64::MAX]);
            let norm = l2_norm(&data);
            assert_eq!(*norm.as_ref(), f64::MAX);
        }

        #[test]
        fn min_positive_value() {
            // Single f64::MIN_POSITIVE should give f64::MIN_POSITIVE
            let data = vec_f64(&[f64::MIN_POSITIVE]);
            let norm = l2_norm(&data);
            assert_eq!(*norm.as_ref(), f64::MIN_POSITIVE);
        }

        #[test]
        fn epsilon() {
            let data = vec_f64(&[f64::EPSILON]);
            let norm = l2_norm(&data);
            assert_eq!(*norm.as_ref(), f64::EPSILON);
        }
    }

    #[cfg(feature = "rug")]
    mod rug_backend {
        use super::*;
        use crate::functions::{Abs, Sqrt};
        use crate::{Constants, RealRugStrictFinite};
        use num::Zero;
        use try_create::TryNew;

        const PRECISION: u32 = 200;
        type R = RealRugStrictFinite<PRECISION>;

        /// Helper to create a validated rug real from f64
        fn rug_f64(v: f64) -> R {
            R::try_from_f64(v).unwrap()
        }

        /// Helper to create a validated rug real from string (for exact values)
        fn rug_str(s: &str) -> R {
            R::try_new(rug::Float::with_val(
                PRECISION,
                rug::Float::parse(s).unwrap(),
            ))
            .unwrap()
        }

        #[test]
        fn basic_3_4_5() {
            let data: Vec<R> = vec![rug_f64(3.0), rug_f64(4.0)];
            let norm = l2_norm(&data);
            let five = rug_f64(5.0);

            let diff = (norm.clone() - &five).abs();
            assert!(
                diff < R::epsilon(),
                "3-4-5: expected 5, got {}, diff {}",
                norm,
                diff
            );
        }

        #[test]
        fn empty_vector() {
            let data: Vec<R> = vec![];
            let norm = l2_norm(&data);
            assert_eq!(norm, R::zero());
        }

        #[test]
        fn single_element() {
            let data = vec![rug_f64(7.0)];
            let norm = l2_norm(&data);
            let seven = rug_f64(7.0);
            assert_eq!(norm, seven);
        }

        #[test]
        fn single_negative() {
            let data = vec![rug_f64(-7.0)];
            let norm = l2_norm(&data);
            let seven = rug_f64(7.0);
            assert_eq!(norm, seven);
        }

        #[test]
        fn all_zeros() {
            let data: Vec<R> = vec![R::zero(), R::zero(), R::zero()];
            let norm = l2_norm(&data);
            assert_eq!(norm, R::zero());
        }

        #[test]
        fn high_precision_values() {
            // Values that cannot be exactly represented in f64
            let data: Vec<R> = vec![rug_str("1e-100"), rug_str("1e-100")];
            let norm = l2_norm(&data);

            // Expected: sqrt(2) * 1e-100
            let expected = rug_str("1e-100") * rug_f64(2.0).sqrt();
            let diff = (norm.clone() - &expected).abs();

            assert!(
                diff < R::epsilon() * &expected,
                "High precision: expected {}, got {}, diff {}",
                expected,
                norm,
                diff
            );
        }

        #[test]
        fn very_large_exponents() {
            // Rug can handle much larger exponents than f64
            // These values would be infinity in f64
            let large = rug_str("1e1000");
            let data: Vec<R> = vec![large.clone(), large.clone()];
            let norm = l2_norm(&data);

            // Expected: sqrt(2) * 1e1000
            let sqrt2 = rug_f64(2.0).sqrt();
            let expected = large * sqrt2;

            let rel_diff = ((norm.clone() - &expected) / &expected).abs();
            assert!(
                rel_diff < R::epsilon(),
                "Very large exponents: rel_diff = {}",
                rel_diff
            );
        }

        #[test]
        fn very_small_exponents() {
            // Rug can handle much smaller exponents than f64
            // These values would be zero in f64
            let small = rug_str("1e-1000");
            let data: Vec<R> = vec![small.clone(), small.clone()];
            let norm = l2_norm(&data);

            // Expected: sqrt(2) * 1e-1000
            let sqrt2 = rug_f64(2.0).sqrt();
            let expected = small * sqrt2;

            let rel_diff = ((norm.clone() - &expected) / &expected).abs();
            assert!(
                rel_diff < R::epsilon(),
                "Very small exponents: rel_diff = {}",
                rel_diff
            );
        }

        #[test]
        fn mixed_signs_rug() {
            let data: Vec<R> = vec![rug_f64(-3.0), rug_f64(4.0)];
            let norm = l2_norm(&data);
            let five = rug_f64(5.0);

            let diff = (norm.clone() - &five).abs();
            assert!(diff < R::epsilon());
        }

        #[test]
        fn zeros_interspersed_rug() {
            let data: Vec<R> = vec![R::zero(), rug_f64(3.0), R::zero(), rug_f64(4.0), R::zero()];
            let norm = l2_norm(&data);
            let five = rug_f64(5.0);

            let diff = (norm.clone() - &five).abs();
            assert!(diff < R::epsilon());
        }

        #[test]
        fn ascending_order_triggers_rescaling() {
            let data: Vec<R> = vec![
                rug_f64(1.0),
                rug_f64(2.0),
                rug_f64(4.0),
                rug_f64(8.0),
                rug_f64(16.0),
            ];
            // sqrt(1 + 4 + 16 + 64 + 256) = sqrt(341)
            let expected = rug_f64(341.0).sqrt();
            let norm = l2_norm(&data);

            let diff = (norm.clone() - &expected).abs();
            assert!(
                diff < R::epsilon() * &expected,
                "Ascending: expected {}, got {}, diff {}",
                expected,
                norm,
                diff
            );
        }

        #[test]
        fn higher_precision() {
            // Test with even higher precision
            const HIGH_PREC: u32 = 500;
            type HP = RealRugStrictFinite<HIGH_PREC>;

            let hp_f64 = |v: f64| HP::try_from_f64(v).unwrap();

            let data: Vec<HP> = vec![hp_f64(3.0), hp_f64(4.0)];
            let norm = l2_norm(&data);
            let five = hp_f64(5.0);

            let diff = (norm.clone() - &five).abs();
            assert!(diff < HP::epsilon());
        }
    }
}