num_valid/algorithms/

vector_norms.rs

//! Functions and utilities for computing vector norms.
//!
//! This module provides public utilities for vector-norm computation on validated
//! scalars (`ScalarType: FpScalar`), supporting both real
//! ([`RealScalar`]) and complex
//! ([`ComplexScalar`](crate::ComplexScalar)) scalars.
//!
//! # Norm type structs
//!
//! Three zero-sized structs identify norm kinds and implement the [`VectorNorm`]
//! and [`ParallelVectorNorm`] traits:
//!
//! | Struct | Formula | Accumulator |
//! |--------|---------|-------------|
//! | [`L1Norm<Acc>`] | `‖v‖₁ = Σᵢ \|vᵢ\|` | [`SumAccumulator`] |
//! | [`L2Norm<Acc>`] | `‖v‖₂ = √(Σᵢ \|vᵢ\|²)` | [`SumAccumulator`] |
//! | [`LinfNorm<Acc>`] | `‖v‖∞ = maxᵢ \|vᵢ\|` | [`MaxAbsValueAccumulator`] |
//!
//! The type parameter `Acc` selects the summation strategy:
//! - `NaiveSum<R>` — fast plain `+=` (default for most functions)
//! - `NeumaierSum<R>` — compensated summation for higher accuracy
//!
//! For the two most common combinations, convenience type aliases are provided:
//!
//! | Alias | Expands to |
//! |-------|------------|
//! | [`L1NormNaiveSum<T>`] | `L1Norm<NaiveSum<T::RealType>>` |
//! | [`L1NormNeumaierSum<T>`] | `L1Norm<NeumaierSum<T::RealType>>` |
//! | [`L2NormNaiveSum<T>`] | `L2Norm<NaiveSum<T::RealType>>` |
//! | [`L2NormNeumaierSum<T>`] | `L2Norm<NeumaierSum<T::RealType>>` |
//!
//! ## Direct usage
//!
//! The structs can be used directly via their static methods — useful when you
//! need a specific norm+strategy combination without going through the free
//! functions. The type aliases make the type annotations shorter:
//!
//! ```
//! use num_valid::algorithms::vector_norms::{
//!     L1Norm, L2Norm, LinfNorm, VectorNorm, ParallelVectorNorm,
//!     L1NormNaiveSum, L2NormNaiveSum, L2NormNeumaierSum,
//! };
//! use num_valid::algorithms::accumulators::{NaiveSum, NeumaierSum};
//!
//! let v = [3.0_f64, -4.0, 5.0];
//!
//! // Using a type alias — equivalent to L1Norm::<NaiveSum<f64>>::compute_vector_norm(&v)
//! let l1 = L1NormNaiveSum::<f64>::compute_vector_norm(&v);
//! assert_eq!(*l1.as_ref(), 12.0);
//!
//! // L2 norm with Neumaier compensated summation (higher accuracy)
//! let l2 = L2NormNeumaierSum::<f64>::compute_vector_norm(&v);
//! assert!((*l2.as_ref() - 50.0_f64.sqrt()).abs() < 1e-12);
//!
//! // L∞ norm (default marker type `Acc = ()`)
//! let linf = LinfNorm::<()>::compute_vector_norm(&v);
//! assert_eq!(*linf.as_ref(), 5.0);
//!
//! // L∞ norm in parallel
//! let linf_par = LinfNorm::<()>::compute_vector_norm_par(&v);
//! assert_eq!(linf_par, linf);
//!
//! // L2 squared norm (not a norm, but avoids the sqrt)
//! let l2_sq = L2NormNaiveSum::<f64>::compute_vector_norm_sq(&v);
//! assert_eq!(*l2_sq.as_ref(), 50.0);
//! ```
//!
//! # Trait-based API
//!
//! [`VectorNorm<ScalarType>`] and [`ParallelVectorNorm<T>`] provide a generic
//! interface that lets algorithms be parameterised over the norm type:
//!
//! ```
//! use num_valid::algorithms::vector_norms::{VectorNorm, L2Norm};
//! use num_valid::algorithms::accumulators::NaiveSum;
//! use num_valid::FpScalar;
//! use num_valid::scalars::NonNegativeRealScalar;
//!
//! fn compute<T: FpScalar, N: VectorNorm<T>>(data: &[T]) -> NonNegativeRealScalar<T::RealType> {
//!     N::compute_vector_norm(data)
//! }
//!
//! let v = [3.0_f64, 4.0];
//! let norm = compute::<_, L2Norm<NaiveSum<f64>>>(&v);
//! assert!((*norm.as_ref() - 5.0).abs() < 1e-12);
//! ```
//!
//! # Free functions
//!
//! Ergonomic wrappers for the most common use cases. All variants accept both
//! slice (`&[ScalarType]`) and owned-iterator (`_iter`-suffixed) inputs.
//! All return [`NonNegativeRealScalar<ScalarType::RealType>`].
//!
//! For real scalars `RealType = ScalarType`; for complex scalars `RealType` is
//! the underlying real component type — the norm is always real-valued because
//! each element contributes via its modulus `|·|`.
//!
//! When the input is borrowed (e.g. `slice.iter()`), chain `.cloned()` before
//! passing to an `_iter`-suffixed function.
//!
//! ## L2 (Euclidean) norms — `‖v‖₂ = √(Σᵢ |vᵢ|²)`
//!
//! | Function | Formula | Summation |
//! |----------|---------|-----------|
//! | [`vector_norm_l2`] | `‖v‖₂` | naive |
//! | [`vector_norm_l2_iter`] | `‖v‖₂` | naive |
//! | [`vector_norm_l2_neumaier`] | `‖v‖₂` | Neumaier compensated |
//! | [`vector_norm_l2_neumaier_iter`] | `‖v‖₂` | Neumaier compensated |
//! | [`vector_norm_l2_sq`] | `‖v‖₂²` | naive |
//! | [`vector_norm_l2_sq_iter`] | `‖v‖₂²` | naive |
//! | [`vector_norm_l2_sq_neumaier`] | `‖v‖₂²` | Neumaier compensated |
//! | [`vector_norm_l2_sq_neumaier_iter`] | `‖v‖₂²` | Neumaier compensated |
//!
//! All L2 variants use a **scaled sum-of-squares** algorithm (equivalent to
//! LAPACK `dnrm2`): the running `scale = max |vᵢ|` ensures that no intermediate
//! value overflows or underflows. The `_neumaier` variants additionally apply
//! compensated summation to recover low-order bits lost when components span
//! many orders of magnitude.
//!
//! Note: `‖v‖₂²` is not a norm (it violates the triangle inequality) and is
//! therefore intentionally absent from the [`VectorNorm`] trait. It is exposed
//! as a computational convenience and as a method on [`L2Norm`].
//!
//! ## L1 (Manhattan) norm — `‖v‖₁ = Σᵢ |vᵢ|`
//!
//! | Function | Formula | Summation |
//! |----------|---------|-----------|
//! | [`vector_norm_l1`] | `‖v‖₁` | naive |
//! | [`vector_norm_l1_iter`] | `‖v‖₁` | naive |
//! | [`vector_norm_l1_neumaier`] | `‖v‖₁` | Neumaier compensated |
//! | [`vector_norm_l1_neumaier_iter`] | `‖v‖₁` | Neumaier compensated |
//!
//! ## L∞ (Chebyshev) norm — `‖v‖∞ = maxᵢ |vᵢ|`
//!
//! | Function | Execution |
//! |----------|-----------|
//! | [`vector_norm_linf`] | sequential, slice |
//! | [`vector_norm_linf_iter`] | sequential, iterator |
//! | [`vector_norm_linf_par`] | parallel, slice |
//! | [`vector_norm_linf_par_iter`] | parallel, iterator |
//!
//! ## Advanced API — custom accumulator
//!
//! For full control over the summation strategy use the accumulator-generic
//! functions. These accept any [`SumAccumulator`] as a type parameter:
//! - [`vector_norm_l1_with_accumulator`] / [`vector_norm_l1_iter_with_accumulator`]
//! - [`vector_norm_l2_with_accumulator`] / [`vector_norm_l2_iter_with_accumulator`]
//! - [`vector_norm_l2_sq_with_accumulator`] / [`vector_norm_l2_sq_iter_with_accumulator`]

#![deny(rustdoc::broken_intra_doc_links)]

use crate::{
    FpScalar, RealScalar,
    algorithms::accumulators::{
        Accumulator, MaxAbsValueAccumulator, NaiveSum, NeumaierSum, ScaledSumSqAccumulator,
        SumAccumulator,
    },
    scalars::NonNegativeRealScalar,
};
use rayon::prelude::*;
use std::marker::PhantomData;

//---------------------------------------------------------
/// Identifies the Manhattan (L1) norm, `‖v‖₁ = Σᵢ |vᵢ|`.
///
/// The type parameter `Acc` selects the summation strategy for absolute values.
/// Typical choices:
/// - `NaiveSum<R>` — plain `+=`, fastest option.
/// - `NeumaierSum<R>` — compensated summation, more accurate when components
///   span many orders of magnitude.
///
/// `L1Norm<Acc>` implements [`VectorNorm`] and [`ParallelVectorNorm`] for any
/// [`FpScalar`] whose `RealType` matches `Acc::Input`.
pub struct L1Norm<Acc>(PhantomData<Acc>);

/// Type alias for [`L1Norm`] using [`NaiveSum`] as summation strategy.
///
/// `L1NormNaiveSum<T>` expands to `L1Norm<NaiveSum<NonNegativeRealScalar<T::RealType>>>`. Use this
/// when performance is the priority and the data is well-scaled.
pub type L1NormNaiveSum<ScalarType> =
    L1Norm<NaiveSum<NonNegativeRealScalar<<ScalarType as FpScalar>::RealType>>>;

/// Type alias for [`L1Norm`] using [`NeumaierSum`] as summation strategy.
///
/// `L1NormNeumaierSum<T>` expands to `L1Norm<NeumaierSum<NonNegativeRealScalar<T::RealType>>>`. Use
/// this when components span many orders of magnitude and accuracy matters.
pub type L1NormNeumaierSum<ScalarType> =
    L1Norm<NeumaierSum<NonNegativeRealScalar<<ScalarType as FpScalar>::RealType>>>;

//---------------------------------------------------------

//---------------------------------------------------------
/// Identifies the Euclidean (L2) norm, `‖v‖₂ = √(Σᵢ |vᵢ|²)`.
///
/// The type parameter `Acc` selects the summation strategy for the scaled
/// sum-of-squares loop. Typical choices:
/// - `NaiveSum<R>` — plain `+=`, fastest option.
/// - `NeumaierSum<R>` — compensated summation for higher accuracy when
///   components span many orders of magnitude.
///
/// `L2Norm<Acc>` implements [`VectorNorm`] and [`ParallelVectorNorm`] for any
/// [`FpScalar`] whose `RealType` matches `Acc::Input`.
///
/// In addition to the standard norm, `L2Norm` exposes the squared Euclidean
/// norm `‖v‖₂² = Σᵢ |vᵢ|²` as a separate operation (avoiding the square root)
/// via [`L2Norm::compute_vector_norm_sq`] and its iterator/parallel variants.
/// Note that `‖v‖₂²` is not a norm — it violates the triangle inequality — and
/// is therefore not part of the [`VectorNorm`] trait.
pub struct L2Norm<Acc>(PhantomData<Acc>);

/// Type alias for [`L2Norm`] using [`NaiveSum`] as summation strategy.
///
/// `L2NormNaiveSum<T>` expands to `L2Norm<NaiveSum<T::RealType>>`. Use this
/// when performance is the priority and the data is well-scaled.
pub type L2NormNaiveSum<ScalarType> = L2Norm<NaiveSum<<ScalarType as FpScalar>::RealType>>;

/// Type alias for [`L2Norm`] using [`NeumaierSum`] as summation strategy.
///
/// `L2NormNeumaierSum<T>` expands to `L2Norm<NeumaierSum<T::RealType>>`. Use
/// this when components span many orders of magnitude and accuracy matters.
pub type L2NormNeumaierSum<ScalarType> = L2Norm<NeumaierSum<<ScalarType as FpScalar>::RealType>>;

impl<Acc: SumAccumulator<Input: RealScalar>> L2Norm<Acc> {
    /// Computes the squared Euclidean (L2) norm `‖v‖₂²` of a slice.
    ///
    /// Equivalent to [`compute_vector_norm_sq_iter`](Self::compute_vector_norm_sq_iter)
    /// called with `v.iter().cloned()`.
    ///
    /// This avoids the square root present in [`VectorNorm::compute_vector_norm`]
    /// and is useful when only relative comparisons of norms are needed.
    pub fn compute_vector_norm_sq<T>(v: &[T]) -> NonNegativeRealScalar<Acc::Output>
    where
        T: FpScalar<RealType = Acc::Output>,
    {
        Self::compute_vector_norm_sq_iter(v.iter().cloned())
    }

    /// Computes the squared Euclidean (L2) norm `‖v‖₂²` from an owned iterator.
    ///
    /// Uses the same scaled sum-of-squares algorithm as [`VectorNorm`] but
    /// returns `scale² × sumsq` directly without taking the square root.
    pub fn compute_vector_norm_sq_iter<I>(v: I) -> NonNegativeRealScalar<Acc::Output>
    where
        I: IntoIterator<Item: FpScalar<RealType = Acc::Output>>,
    {
        ScaledSumSqAccumulator::<Acc>::new_sequential_from_scalars(v).result()
    }
}

impl<Acc> L2Norm<Acc>
where
    Acc: SumAccumulator<Input: RealScalar> + Send,
{
    /// Computes the squared Euclidean (L2) norm `‖v‖₂²` of a slice in parallel.
    ///
    /// Equivalent to [`compute_vector_norm_sq_par_iter`](Self::compute_vector_norm_sq_par_iter)
    /// called with `v.par_iter().cloned()`.
    pub fn compute_vector_norm_sq_par<T>(v: &[T]) -> NonNegativeRealScalar<Acc::Output>
    where
        T: FpScalar<RealType = Acc::Output>,
    {
        Self::compute_vector_norm_sq_par_iter(v.par_iter().cloned())
    }

    /// Computes the squared Euclidean (L2) norm `‖v‖₂²` from a parallel iterator.
    ///
    /// Chunks are processed in parallel and merged via
    /// [`Accumulator::combine`],
    /// preserving the numerical properties of the chosen `Acc` strategy across chunks.
    pub fn compute_vector_norm_sq_par_iter<I>(v: I) -> NonNegativeRealScalar<Acc::Output>
    where
        I: IntoParallelIterator<Item: FpScalar<RealType = Acc::Output>>,
    {
        ScaledSumSqAccumulator::<Acc>::new_parallel_from_scalars(v).result()
    }
}

//---------------------------------------------------------

//pub struct L2SqNorm<Strategy>(PhantomData<Strategy>);

//---------------------------------------------------------
/// Identifies the Chebyshev (L∞) norm, `‖v‖∞ = maxᵢ |vᵢ|`.
///
/// Unlike [`L1Norm`] and [`L2Norm`], `LinfNorm` does not use a summation
/// strategy: it is a single linear scan for the maximum absolute value.
///
/// The `Acc` parameter is currently a marker type (default `()`) and is not
/// used by the vector implementation itself. It exists to keep the API
/// extensible and to align type-level dispatch with norm families where a
/// dedicated accumulator may be useful in related contexts (for example,
/// matrix norms).
///
/// `LinfNorm` implements both [`VectorNorm`] and [`ParallelVectorNorm`] for any
/// [`FpScalar`]. The parallel variant uses
/// [`MaxAbsValueAccumulator`]
/// with rayon `fold` + `reduce`, which is numerically exact (no rounding).
pub struct LinfNorm<Acc = ()>(PhantomData<Acc>);
//---------------------------------------------------------

//---------------------------------------------------------
/// Sequential vector norm computation.
///
/// Implemented by [`L1Norm`], [`L2Norm`], and [`LinfNorm`]. Generic code that
/// accepts any norm should bound `N: VectorNorm<T>`.
///
/// All implementations accept both real and complex scalars: for complex types
/// each element contributes via its modulus `|zᵢ|`, so the result is always
/// a real-valued [`NonNegativeRealScalar`].
///
/// # Methods
///
/// - [`compute_vector_norm`](VectorNorm::compute_vector_norm) — slice input
///   (delegates to `compute_vector_norm_iter` by default).
/// - [`compute_vector_norm_iter`](VectorNorm::compute_vector_norm_iter) — owned
///   iterator input; **must** be implemented by each norm type.
pub trait VectorNorm<ScalarType: FpScalar> {
    /// Computes the norm of a slice of scalars.
    ///
    /// The default implementation calls [`compute_vector_norm_iter`](Self::compute_vector_norm_iter)
    /// with `v.iter().cloned()`.
    fn compute_vector_norm(v: &[ScalarType]) -> NonNegativeRealScalar<ScalarType::RealType> {
        Self::compute_vector_norm_iter(v.iter().cloned())
    }

    /// Computes the norm from an owned-item iterator.
    ///
    /// Iterator-based APIs require **owned** items. When the input is borrowed
    /// (e.g. from a slice), chain `.iter().cloned()` before calling this method.
    fn compute_vector_norm_iter<I>(v: I) -> NonNegativeRealScalar<ScalarType::RealType>
    where
        I: IntoIterator<Item = ScalarType>;
}

/// Parallel vector norm computation.
///
/// Extends [`VectorNorm`] with parallel variants that distribute work across
/// threads using rayon. Chunks are accumulated independently and merged via
/// [`Accumulator::combine`].
///
/// Requires `T: Send + Sync`.
///
/// # Methods
///
/// - [`compute_vector_norm_par`](ParallelVectorNorm::compute_vector_norm_par) — parallel slice.
/// - [`compute_vector_norm_par_iter`](ParallelVectorNorm::compute_vector_norm_par_iter) —
///   parallel iterator; **must** be implemented by each norm type.
pub trait ParallelVectorNorm<T: FpScalar + Send + Sync>: VectorNorm<T> {
    /// Computes the norm of a slice in parallel.
    ///
    /// The default implementation calls
    /// [`compute_vector_norm_par_iter`](Self::compute_vector_norm_par_iter)
    /// with `v.par_iter().cloned()`.
    fn compute_vector_norm_par(v: &[T]) -> NonNegativeRealScalar<T::RealType> {
        Self::compute_vector_norm_par_iter(v.par_iter().cloned())
    }

    /// Computes the norm from an owned parallel iterator.
    fn compute_vector_norm_par_iter<I: IntoParallelIterator<Item = T>>(
        iter: I,
    ) -> NonNegativeRealScalar<T::RealType>;
}
//---------------------------------------------------------

//---------------------------------------------------------
impl<ScalarType, SumAcc> VectorNorm<ScalarType> for L2Norm<SumAcc>
where
    ScalarType: FpScalar,
    SumAcc: SumAccumulator<Input = ScalarType::RealType>,
{
    fn compute_vector_norm_iter<I>(v: I) -> NonNegativeRealScalar<ScalarType::RealType>
    where
        I: IntoIterator<Item = ScalarType>,
    {
        ScaledSumSqAccumulator::<SumAcc>::new_sequential_from_scalars(v).result_sqrt()
    }
}

impl<T, SumAcc> ParallelVectorNorm<T> for L2Norm<SumAcc>
where
    T: FpScalar,
    SumAcc: SumAccumulator<Input = T::RealType> + Send,
{
    fn compute_vector_norm_par_iter<I: IntoParallelIterator<Item = T>>(
        iter: I,
    ) -> NonNegativeRealScalar<T::RealType> {
        ScaledSumSqAccumulator::<SumAcc>::new_parallel_from_scalars(iter).result_sqrt()
    }
}
//---------------------------------------------------------

//---------------------------------------------------------
impl<ScalarType, SumAcc> VectorNorm<ScalarType> for L1Norm<SumAcc>
where
    ScalarType: FpScalar,
    SumAcc: SumAccumulator<Input = NonNegativeRealScalar<ScalarType::RealType>>,
{
    fn compute_vector_norm_iter<I>(v: I) -> NonNegativeRealScalar<ScalarType::RealType>
    where
        I: IntoIterator<Item = ScalarType>,
    {
        SumAcc::new_sequential(v.into_iter().map(|x| x.abs())).result()
    }
}

impl<T, SumAcc> ParallelVectorNorm<T> for L1Norm<SumAcc>
where
    T: FpScalar,
    SumAcc: SumAccumulator<Input = NonNegativeRealScalar<T::RealType>> + Send,
{
    fn compute_vector_norm_par_iter<I: IntoParallelIterator<Item = T>>(
        iter: I,
    ) -> NonNegativeRealScalar<T::RealType> {
        SumAcc::new_parallel(iter.into_par_iter().map(|x| x.abs())).result()
    }
}

//---------------------------------------------------------

//---------------------------------------------------------
impl<ScalarType> VectorNorm<ScalarType> for LinfNorm
where
    ScalarType: FpScalar,
{
    fn compute_vector_norm_iter<I>(v: I) -> NonNegativeRealScalar<ScalarType::RealType>
    where
        I: IntoIterator<Item = ScalarType>,
    {
        MaxAbsValueAccumulator::new_sequential(v).result()
    }
}

impl<ScalarType> ParallelVectorNorm<ScalarType> for LinfNorm
where
    ScalarType: FpScalar,
{
    fn compute_vector_norm_par_iter<I>(v: I) -> NonNegativeRealScalar<ScalarType::RealType>
    where
        I: IntoParallelIterator<Item = ScalarType>,
    {
        MaxAbsValueAccumulator::new_parallel(v).result()
    }
}
//---------------------------------------------------------

// ── Public L2 norm functions ─────────────────────────────────────────────────

/// Accumulator-driven implementation of the Euclidean (L2) norm.
///
/// This helper computes `vector_norm_l2` using the provided accumulation strategy
/// `SumSqAcc` for the scaled sum-of-squares term.
///
/// The final value is guaranteed to be non-negative by construction and is
/// wrapped into [`NonNegativeRealScalar`].
///
/// # Panics
///
/// Panics only if the computed value is negative, which indicates an internal
/// invariant violation.
#[inline(always)]
pub fn vector_norm_l2_with_accumulator<ScalarType, SumSqAcc>(
    v: &[ScalarType],
) -> NonNegativeRealScalar<ScalarType::RealType>
where
    ScalarType: FpScalar,
    SumSqAcc: SumAccumulator<Input = ScalarType::RealType>,
{
    L2Norm::<SumSqAcc>::compute_vector_norm(v)
}

/// Iterator-based accumulator-driven implementation of the Euclidean (L2) norm.
///
/// This is the iterator counterpart of [`vector_norm_l2_with_accumulator`] and accepts
/// any `IntoIterator<Item = ScalarType>`.
///
/// Use this when your inputs are already produced as an iterator chain.
#[inline(always)]
pub fn vector_norm_l2_iter_with_accumulator<ScalarType, SumSqAcc, I>(
    v: I,
) -> NonNegativeRealScalar<ScalarType::RealType>
where
    ScalarType: FpScalar,
    SumSqAcc: SumAccumulator<Input = ScalarType::RealType>,
    I: IntoIterator<Item = ScalarType>,
{
    L2Norm::<SumSqAcc>::compute_vector_norm_iter(v)
}

/// Accumulator-driven implementation of the squared Euclidean (L2) norm.
///
/// This function computes `vector_norm_l2_sq` using the provided accumulation
/// strategy `SumSqAcc` for the scaled sum-of-squares term.
///
/// Use this API when you need explicit control over how terms are accumulated,
/// for example to compare [`NaiveSum`] with [`NeumaierSum`] or a custom
/// accumulator implementation.
///
/// The final value is guaranteed to be non-negative by construction and is
/// wrapped into [`NonNegativeRealScalar`].
///
/// # Examples
///
/// ```
/// use num_valid::algorithms::vector_norms::vector_norm_l2_sq_with_accumulator;
/// use num_valid::algorithms::accumulators::{NaiveSum, NeumaierSum};
///
/// let v = [3.0_f64, 4.0];
/// let naive = vector_norm_l2_sq_with_accumulator::<_, NaiveSum<_>>(&v);
/// let compensated = vector_norm_l2_sq_with_accumulator::<_, NeumaierSum<_>>(&v);
/// assert_eq!(*naive.as_ref(), 25.0);
/// assert_eq!(naive, compensated);
/// ```
///
/// # Panics
///
/// Panics only if the computed value is negative, which indicates an internal
/// invariant violation.
#[inline(always)]
pub fn vector_norm_l2_sq_with_accumulator<ScalarType, SumSqAcc>(
    v: &[ScalarType],
) -> NonNegativeRealScalar<ScalarType::RealType>
where
    ScalarType: FpScalar,
    SumSqAcc: SumAccumulator<Input = ScalarType::RealType>,
{
    L2Norm::<SumSqAcc>::compute_vector_norm_sq(v)
}

/// Iterator-based accumulator-driven implementation of the squared Euclidean
/// (L2) norm.
///
/// This is the iterator counterpart of [`vector_norm_l2_sq_with_accumulator`] and
/// accepts any `IntoIterator<Item = ScalarType>`.
#[inline(always)]
pub fn vector_norm_l2_sq_iter_with_accumulator<ScalarType, SumSqAcc, I>(
    v: I,
) -> NonNegativeRealScalar<ScalarType::RealType>
where
    ScalarType: FpScalar,
    SumSqAcc: SumAccumulator<Input = ScalarType::RealType>,
    I: IntoIterator<Item = ScalarType>,
{
    L2Norm::<SumSqAcc>::compute_vector_norm_sq_iter(v)
}

/// Accumulator-driven implementation of the Manhattan (L1) norm.
///
/// This function computes `vector_norm_l1` by summing absolute values with the
/// accumulation strategy `SumAbsAcc`.
///
/// Use this API when you need explicit control over how absolute values are
/// accumulated, for example to compare [`NaiveSum`] with [`NeumaierSum`] or a
/// custom accumulator implementation.
///
/// The final value is guaranteed to be non-negative by construction and is
/// wrapped into [`NonNegativeRealScalar`].
///
/// # Examples
///
/// ```
/// use num_valid::algorithms::vector_norms::vector_norm_l1_with_accumulator;
/// use num_valid::algorithms::accumulators::{NaiveSum, NeumaierSum};
///
/// let v = [3.0_f64, 4.0];
/// let naive = vector_norm_l1_with_accumulator::<_, NaiveSum<_>>(&v);
/// let compensated = vector_norm_l1_with_accumulator::<_, NeumaierSum<_>>(&v);
/// assert_eq!(*naive.as_ref(), 7.0);
/// assert_eq!(naive, compensated);
/// ```
///
/// # Panics
///
/// Panics only if the computed value is negative, which indicates an internal
/// invariant violation.
#[inline(always)]
pub fn vector_norm_l1_with_accumulator<ScalarType, SumAcc>(
    a: &[ScalarType],
) -> NonNegativeRealScalar<ScalarType::RealType>
where
    ScalarType: FpScalar,
    SumAcc: SumAccumulator<Input = NonNegativeRealScalar<ScalarType::RealType>>,
{
    L1Norm::<SumAcc>::compute_vector_norm(a)
}

/// Iterator-based accumulator-driven implementation of the Manhattan (L1) norm.
///
/// This is the iterator counterpart of [`vector_norm_l1_with_accumulator`] and accepts
/// any `IntoIterator<Item = ScalarType>`.
#[inline(always)]
pub fn vector_norm_l1_iter_with_accumulator<ScalarType, SumAcc, I>(
    a: I,
) -> NonNegativeRealScalar<ScalarType::RealType>
where
    ScalarType: FpScalar,
    SumAcc: SumAccumulator<Input = NonNegativeRealScalar<ScalarType::RealType>>,
    I: IntoIterator<Item = ScalarType>,
{
    L1Norm::<SumAcc>::compute_vector_norm_iter(a)
}

/// Computes the Euclidean (L2) norm of a slice of scalars.
///
/// Returns ‖**v**‖₂ = √(Σᵢ |vᵢ|²), where |·| denotes the modulus (absolute
/// value for reals, complex modulus for complex scalars).
///
/// # Numerical Stability
///
/// The result is computed using a **scaled sum-of-squares** algorithm: the
/// running `scale` equals the largest |vᵢ| seen so far and `sumsq`
/// accumulates (|vᵢ| / scale)².  The final result is `scale × √sumsq`.  This
/// avoids overflow and underflow of intermediate values.
///
/// The algorithm is equivalent to the one used in LAPACK's `dnrm2` routine.
/// Note: scaling alone does not prevent the loss of small addends when they are
/// absorbed into a larger running total (absorption error); use [`vector_norm_l2_neumaier`]
/// when that precision is required.
///
/// Equivalent to calling the internal accumulator-driven implementation with
/// [`NaiveSum`].
///
/// # Examples
///
/// ```
/// use num_valid::algorithms::vector_norms::vector_norm_l2;
///
/// // 3-4-5 right triangle: norm of (3, 4) is 5
/// let v = [3.0_f64, 4.0];
/// assert!((*vector_norm_l2(&v).as_ref() - 5.0).abs() < 1e-12);
///
/// // Zero vector
/// let z = [0.0_f64; 3];
/// assert_eq!(*vector_norm_l2(&z).as_ref(), 0.0);
///
/// // Unit axis vector has norm 1
/// let e1 = [1.0_f64, 0.0, 0.0];
/// assert!((*vector_norm_l2(&e1).as_ref() - 1.0).abs() < 1e-15);
/// ```
///
/// Complex scalars are also supported; `|zᵢ|` is the complex modulus:
///
/// ```
/// use num_valid::algorithms::vector_norms::vector_norm_l2;
/// use num_valid::ComplexNative64StrictFinite;
/// use num::Complex;
/// use try_create::TryNew;
///
/// // ‖[3+0i, 0+4i]‖₂ = √(3² + 4²) = 5
/// let v = [
///     ComplexNative64StrictFinite::try_new(Complex::new(3.0_f64, 0.0)).unwrap(),
///     ComplexNative64StrictFinite::try_new(Complex::new(0.0_f64, 4.0)).unwrap(),
/// ];
/// let result = vector_norm_l2(&v).into_inner();
/// assert!((*result.as_ref() - 5.0_f64).abs() < 1e-12);
/// ```
#[inline(always)]
pub fn vector_norm_l2<ScalarType>(v: &[ScalarType]) -> NonNegativeRealScalar<ScalarType::RealType>
where
    ScalarType: FpScalar,
{
    L2NormNaiveSum::<ScalarType>::compute_vector_norm(v)
}

/// Computes the Euclidean (L2) norm from an owned-item iterator.
///
/// This is the iterator counterpart of [`vector_norm_l2`].
///
/// # Examples
///
/// ```
/// use num_valid::algorithms::vector_norms::vector_norm_l2_iter;
///
/// let v = vec![3.0_f64, 4.0];
/// assert!((*vector_norm_l2_iter(v).as_ref() - 5.0).abs() < 1e-12);
/// ```
#[inline(always)]
pub fn vector_norm_l2_iter<ScalarType, I>(v: I) -> NonNegativeRealScalar<ScalarType::RealType>
where
    ScalarType: FpScalar,
    I: IntoIterator<Item = ScalarType>,
{
    L2NormNaiveSum::<ScalarType>::compute_vector_norm_iter(v)
}

/// Computes the squared Euclidean (L2) norm of a slice of scalars.
///
/// Returns ‖**v**‖₂² = Σᵢ |vᵢ|², where |·| denotes the modulus.
///
/// # Numerical Stability
///
/// The result is computed using a **scaled sum-of-squares** algorithm: the
/// running `scale` equals the largest |vᵢ| seen so far and `sumsq`
/// accumulates (|vᵢ| / scale)².  The final result is `scale² × sumsq`.  This
/// avoids overflow and underflow of intermediate values.
///
/// Note: `scale²` must not overflow the maximum representable value of the
/// underlying real type. For vectors where accuracy of small squared terms
/// matters, prefer [`vector_norm_l2_sq_neumaier`].
///
/// Equivalent to calling [`vector_norm_l2_sq_with_accumulator`] with [`NaiveSum`].
///
/// # Examples
///
/// ```
/// use num_valid::algorithms::vector_norms::vector_norm_l2_sq;
///
/// // 3-4-5 right triangle: ‖(3, 4)‖₂² == 9 + 16 == 25
/// let v = [3.0_f64, 4.0];
/// assert_eq!(*vector_norm_l2_sq(&v).as_ref(), 25.0);
///
/// // Zero vector
/// let z = [0.0_f64; 3];
/// assert_eq!(*vector_norm_l2_sq(&z).as_ref(), 0.0);
///
/// // Unit axis vector has squared norm 1
/// let e1 = [1.0_f64, 0.0, 0.0];
/// assert_eq!(*vector_norm_l2_sq(&e1).as_ref(), 1.0);
/// ```
#[inline(always)]
pub fn vector_norm_l2_sq<ScalarType: FpScalar>(
    v: &[ScalarType],
) -> NonNegativeRealScalar<ScalarType::RealType> {
    L2NormNaiveSum::<ScalarType>::compute_vector_norm_sq(v)
}

/// Computes the squared Euclidean (L2) norm from an owned-item iterator.
///
/// This is the iterator counterpart of [`vector_norm_l2_sq`].
#[inline(always)]
pub fn vector_norm_l2_sq_iter<ScalarType, I>(v: I) -> NonNegativeRealScalar<ScalarType::RealType>
where
    ScalarType: FpScalar,
    I: IntoIterator<Item = ScalarType>,
{
    L2NormNaiveSum::<ScalarType>::compute_vector_norm_sq_iter(v)
}

/// Computes the squared Euclidean (L2) norm of a slice of scalars using
/// **scaling + Neumaier compensated summation**.
///
/// Returns ‖**v**‖₂² = Σᵢ |vᵢ|², where |·| denotes the modulus.
///
/// # Numerical Accuracy
///
/// This is an improved variant of [`vector_norm_l2_sq`] that combines two complementary
/// techniques:
///
/// 1. **Scaling** (same as [`vector_norm_l2_sq`]): tracks `scale = max |vᵢ|` and
///    accumulates `Σ (|vᵢ|/scale)²`, so that no individual term overflows or
///    underflows before being squared.
///
/// 2. **Neumaier compensated summation**: replaces the plain `sumsq += r²` with
///    Kahan/Neumaier compensation, capturing low-order bits that are silently
///    discarded in the naive loop. This is particularly beneficial when a large
///    dominant component (near `scale`) causes the unit-in-last-place of `sumsq`
///    to exceed many small squared terms.
///
/// **When rescaling occurs** (a new maximum is found mid-iteration), the
/// accumulated Neumaier sum is multiplied by `r²` and the compensation term is
/// reset to zero. Rescaling is rare in practice (occurs at most log(N) times for
/// random data), so the accuracy benefit of Neumaier is preserved for the bulk
/// of the summation.
///
/// # Special Cases
///
/// | Input | Result |
/// |-------|--------|
/// | all-zero vector | `0.0` |
/// | single dominant component | full Neumaier accuracy for the remaining terms |
///
/// # Examples
///
/// ```
/// use num_valid::algorithms::vector_norms::vector_norm_l2_sq_neumaier;
///
/// // 3-4-5 right triangle: ‖(3, 4)‖₂² == 25
/// let v = [3.0_f64, 4.0];
/// assert_eq!(*vector_norm_l2_sq_neumaier(&v).as_ref(), 25.0);
///
/// // Zero vector
/// let z = [0.0_f64; 3];
/// assert_eq!(*vector_norm_l2_sq_neumaier(&z).as_ref(), 0.0);
///
/// // Recovers small squared terms lost by naive accumulation:
/// // v = [1.0, 1e-9 × 10_000]; true ‖v‖₂² = 1.0 + 1e-14
/// const N: usize = 10_000;
/// let mut v_large: [f64; N + 1] = [1e-9; N + 1];
/// v_large[0] = 1.0;
/// let result = *vector_norm_l2_sq_neumaier(&v_large).as_ref();
/// let expected = 1.0_f64 + N as f64 * 1e-18;
/// assert!((result - expected).abs() < 1e-15 * expected);
/// ```
#[inline(always)]
pub fn vector_norm_l2_sq_neumaier<ScalarType: FpScalar>(
    v: &[ScalarType],
) -> NonNegativeRealScalar<ScalarType::RealType> {
    L2NormNeumaierSum::<ScalarType>::compute_vector_norm_sq(v)
}

/// Computes the squared Euclidean (L2) norm from an owned-item iterator using
/// scaling plus Neumaier compensated summation.
///
/// This is the iterator counterpart of [`vector_norm_l2_sq_neumaier`].
#[inline(always)]
pub fn vector_norm_l2_sq_neumaier_iter<ScalarType, I>(
    v: I,
) -> NonNegativeRealScalar<ScalarType::RealType>
where
    ScalarType: FpScalar,
    I: IntoIterator<Item = ScalarType>,
{
    L2NormNeumaierSum::<ScalarType>::compute_vector_norm_sq_iter(v)
}

/// Computes the Euclidean (L2) norm of a slice of scalars using
/// **scaling + Neumaier compensated summation**.
///
/// Returns ‖**v**‖₂ = √(Σᵢ |vᵢ|²), where |·| denotes the modulus.
///
/// # Numerical Accuracy
///
/// This is an improved variant of [`vector_norm_l2`] that combines two complementary
/// techniques:
///
/// 1. **Scaling** (same as [`vector_norm_l2`]): tracks `scale = max |vᵢ|` and
///    accumulates `Σ (|vᵢ|/scale)²`, so that no individual term overflows or
///    underflows before being squared.
///
/// 2. **Neumaier compensated summation**: replaces the plain `sumsq += r²` with
///    Kahan/Neumaier compensation, capturing low-order bits that are silently
///    discarded in the naive loop. The final result is `scale × √(neumaier.result())`.
///
/// **When rescaling occurs** (a new maximum is found mid-iteration), the
/// accumulated Neumaier sum is multiplied by `r²` and the compensation term is
/// reset to zero. Rescaling is rare in practice, so the Neumaier accuracy benefit
/// is preserved for the bulk of the summation.
///
/// # Special Cases
///
/// | Input | Result |
/// |-------|--------|
/// | all-zero vector | `0.0` |
/// | single dominant component | full Neumaier accuracy for the remaining terms |
///
/// # Examples
///
/// ```
/// use num_valid::algorithms::vector_norms::vector_norm_l2_neumaier;
///
/// // 3-4-5 right triangle: ‖(3, 4)‖₂ == 5
/// let v = [3.0_f64, 4.0];
/// assert!((*vector_norm_l2_neumaier(&v).as_ref() - 5.0).abs() < 1e-12);
///
/// // Zero vector
/// let z = [0.0_f64; 3];
/// assert_eq!(*vector_norm_l2_neumaier(&z).as_ref(), 0.0);
///
/// // Recovers small squared terms lost by naive accumulation:
/// // v = [1.0, 1e-9 × 10_000]; true ‖v‖₂ = sqrt(1.0 + 1e-14) ≈ 1.0 + 5e-15
/// const N: usize = 10_000;
/// let mut v_large: [f64; N + 1] = [1e-9; N + 1];
/// v_large[0] = 1.0;
/// let result = *vector_norm_l2_neumaier(&v_large).as_ref();
/// let expected = (1.0_f64 + N as f64 * 1e-18).sqrt();
/// assert!((result - expected).abs() < 1e-15 * expected);
/// ```
#[inline(always)]
pub fn vector_norm_l2_neumaier<ScalarType: FpScalar>(
    v: &[ScalarType],
) -> NonNegativeRealScalar<ScalarType::RealType> {
    L2NormNeumaierSum::<ScalarType>::compute_vector_norm(v)
}

/// Computes the Euclidean (L2) norm from an owned-item iterator using scaling
/// plus Neumaier compensated summation.
///
/// This is the iterator counterpart of [`vector_norm_l2_neumaier`].
#[inline(always)]
pub fn vector_norm_l2_neumaier_iter<ScalarType, I>(
    v: I,
) -> NonNegativeRealScalar<ScalarType::RealType>
where
    ScalarType: FpScalar,
    I: IntoIterator<Item = ScalarType>,
{
    L2NormNeumaierSum::<ScalarType>::compute_vector_norm_iter(v)
}

/// Computes the Manhattan (L1) norm of a slice of scalars.
///
/// Returns ‖**v**‖₁ = Σᵢ |vᵢ|, where |·| denotes the modulus.
///
/// # Numerical Accuracy
///
/// The absolute values are summed with a **naive `+=` loop**.
/// For typical inputs this is exact up to the rounding of each individual
/// `abs()` call.  When components span many orders of magnitude a large
/// running total may absorb small addends (absorption error), silently
/// discarding their low-order bits.  Use [`vector_norm_l1_neumaier`] when that
/// extra precision matters.
///
/// Equivalent to calling [`vector_norm_l1_with_accumulator`] with [`NaiveSum`].
///
/// # Examples
///
/// ```
/// use num_valid::algorithms::vector_norms::vector_norm_l1;
///
/// // L1 norm of (3, 4) is 3 + 4 = 7
/// let v = [3.0_f64, 4.0];
/// assert_eq!(*vector_norm_l1(&v).as_ref(), 7.0);
///
/// // Zero vector
/// let z = [0.0_f64; 3];
/// assert_eq!(*vector_norm_l1(&z).as_ref(), 0.0);
///
/// // Unit axis vector has L1 norm 1
/// let e1 = [1.0_f64, 0.0, 0.0];
/// assert_eq!(*vector_norm_l1(&e1).as_ref(), 1.0);
/// ```
///
/// For complex scalars the L1 norm sums the complex moduli:
///
/// ```
/// use num_valid::algorithms::vector_norms::vector_norm_l1;
/// use num_valid::ComplexNative64StrictFinite;
/// use num::Complex;
/// use try_create::TryNew;
///
/// // ‖[3+0i, 0+4i]‖₁ = |3+0i| + |0+4i| = 3 + 4 = 7
/// let v = [
///     ComplexNative64StrictFinite::try_new(Complex::new(3.0_f64, 0.0)).unwrap(),
///     ComplexNative64StrictFinite::try_new(Complex::new(0.0_f64, 4.0)).unwrap(),
/// ];
/// let result = vector_norm_l1(&v).into_inner();
/// assert!((*result.as_ref() - 7.0_f64).abs() < 1e-12);
/// ```
#[inline(always)]
pub fn vector_norm_l1<ScalarType: FpScalar>(
    a: &[ScalarType],
) -> NonNegativeRealScalar<ScalarType::RealType> {
    L1NormNaiveSum::<ScalarType>::compute_vector_norm(a)
}

/// Computes the Manhattan (L1) norm from an owned-item iterator.
///
/// This is the iterator counterpart of [`vector_norm_l1`].
#[inline(always)]
pub fn vector_norm_l1_iter<ScalarType, I>(a: I) -> NonNegativeRealScalar<ScalarType::RealType>
where
    ScalarType: FpScalar,
    I: IntoIterator<Item = ScalarType>,
{
    L1NormNaiveSum::<ScalarType>::compute_vector_norm_iter(a)
}

/// Computes the Manhattan (L1) norm of a slice of scalars using
/// **Neumaier compensated summation**.
///
/// Returns ‖**v**‖₁ = Σᵢ |vᵢ|, where |·| denotes the modulus.
///
/// # Numerical Accuracy
///
/// The absolute values are accumulated using the **Neumaier compensated
/// summation** algorithm, which corrects for round-off errors that can
/// accumulate in a naïve loop.  For typical coordinate values the improvement
/// is small, but for vectors whose components span many orders of magnitude
/// the compensation captures the round-off lost when a small addend is absorbed
/// by a much larger running total (absorption error).  Catastrophic cancellation
/// — the loss of significant bits by subtracting nearly-equal numbers — cannot
/// occur in a sum of absolute values, so the benefit here is specifically the
/// recovery of those absorbed low-order bits.
///
/// # Examples
///
/// ```
/// use num_valid::algorithms::vector_norms::vector_norm_l1_neumaier;
///
/// // L1 norm of (3, 4) is 3 + 4 = 7
/// let v = [3.0_f64, 4.0];
/// assert_eq!(*vector_norm_l1_neumaier(&v).as_ref(), 7.0);
///
/// // Zero vector
/// let z = [0.0_f64; 3];
/// assert_eq!(*vector_norm_l1_neumaier(&z).as_ref(), 0.0);
///
/// // Unit axis vector has L1 norm 1
/// let e1 = [1.0_f64, 0.0, 0.0];
/// assert_eq!(*vector_norm_l1_neumaier(&e1).as_ref(), 1.0);
/// ```
#[inline(always)]
pub fn vector_norm_l1_neumaier<ScalarType: FpScalar>(
    a: &[ScalarType],
) -> NonNegativeRealScalar<ScalarType::RealType> {
    L1NormNeumaierSum::<ScalarType>::compute_vector_norm(a)
}

/// Computes the Manhattan (L1) norm from an owned-item iterator using
/// Neumaier compensated summation.
///
/// This is the iterator counterpart of [`vector_norm_l1_neumaier`].
#[inline(always)]
pub fn vector_norm_l1_neumaier_iter<ScalarType, I>(
    a: I,
) -> NonNegativeRealScalar<ScalarType::RealType>
where
    ScalarType: FpScalar,
    I: IntoIterator<Item = ScalarType>,
{
    L1NormNeumaierSum::<ScalarType>::compute_vector_norm_iter(a)
}

/// Computes the Chebyshev (L∞) norm of a slice of scalars.
///
/// Returns ‖**v**‖∞ = maxᵢ |vᵢ|, where |·| denotes the modulus.
///
/// # Algorithm
///
/// A single linear scan keeps the running maximum absolute value.  The result
/// requires no square root and no compensation — it is numerically exact up to
/// a single call to `abs()` on each component.
///
/// # Examples
///
/// ```
/// use num_valid::algorithms::vector_norms::vector_norm_linf;
///
/// // L∞ norm of (3, 4) = max(3, 4) = 4
/// let v = [3.0_f64, 4.0];
/// assert_eq!(*vector_norm_linf(&v).as_ref(), 4.0);
///
/// // Zero vector
/// let z = [0.0_f64; 3];
/// assert_eq!(*vector_norm_linf(&z).as_ref(), 0.0);
///
/// // Unit axis vector has L∞ norm 1
/// let e1 = [1.0_f64, 0.0, 0.0];
/// assert_eq!(*vector_norm_linf(&e1).as_ref(), 1.0);
/// ```
///
/// For complex scalars `|·|` is the complex modulus:
///
/// ```
/// use num_valid::algorithms::vector_norms::vector_norm_linf;
/// use num_valid::ComplexNative64StrictFinite;
/// use num::Complex;
/// use try_create::TryNew;
///
/// // ‖[3+0i, 0+4i]‖∞ = max(|3+0i|, |0+4i|) = max(3, 4) = 4
/// let v = [
///     ComplexNative64StrictFinite::try_new(Complex::new(3.0_f64, 0.0)).unwrap(),
///     ComplexNative64StrictFinite::try_new(Complex::new(0.0_f64, 4.0)).unwrap(),
/// ];
/// let result = vector_norm_linf(&v).into_inner();
/// assert!((*result.as_ref() - 4.0_f64).abs() < 1e-12);
/// ```
#[inline(always)]
pub fn vector_norm_linf<ScalarType: FpScalar>(
    a: &[ScalarType],
) -> NonNegativeRealScalar<ScalarType::RealType> {
    LinfNorm::<()>::compute_vector_norm(a)
}

/// Computes the Chebyshev (L∞) norm of a slice of scalars in parallel.
///
/// Returns ‖**v**‖∞ = maxᵢ |vᵢ|, where |·| denotes the modulus.
///
/// This is the parallel counterpart of [`vector_norm_linf`].
///
/// # Examples
///
/// ```
/// use num_valid::algorithms::vector_norms::vector_norm_linf_par;
///
/// let v = [3.0_f64, 4.0, -2.0];
/// assert_eq!(*vector_norm_linf_par(&v).as_ref(), 4.0);
/// ```
#[inline(always)]
pub fn vector_norm_linf_par<ScalarType: FpScalar>(
    a: &[ScalarType],
) -> NonNegativeRealScalar<ScalarType::RealType> {
    LinfNorm::<()>::compute_vector_norm_par(a)
}

/// Computes the Chebyshev (L∞) norm from an owned-item iterator.
///
/// This is the iterator counterpart of [`vector_norm_linf`].
#[inline(always)]
pub fn vector_norm_linf_iter<ScalarType: FpScalar, I>(
    a: I,
) -> NonNegativeRealScalar<ScalarType::RealType>
where
    I: IntoIterator<Item = ScalarType>,
{
    LinfNorm::<()>::compute_vector_norm_iter(a)
}

/// Computes the infinity norm (maximum absolute value) of a parallel iterator of values.
///
/// This function computes the infinity norm, which is the maximum absolute value of the elements
/// in the given parallel iterator. It uses the `Abs` trait to compute the absolute value of each
/// element and the `Max` trait to find the maximum value.
///
/// This is the parallel iterator counterpart of [`vector_norm_linf_iter`].
///
/// # Examples
///
/// ```
/// use num_valid::algorithms::vector_norms::vector_norm_linf_par_iter;
///
/// let v = vec![3.0_f64, 4.0, -2.0];
/// assert_eq!(*vector_norm_linf_par_iter(v).as_ref(), 4.0);
/// ```
#[inline(always)]
pub fn vector_norm_linf_par_iter<ScalarType: FpScalar, I>(
    a: I,
) -> NonNegativeRealScalar<ScalarType::RealType>
where
    I: IntoParallelIterator<Item = ScalarType>,
{
    LinfNorm::<()>::compute_vector_norm_par_iter(a)
}

// ================================================================================================
// Tests
// ================================================================================================

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

    mod vector_norm_l2_tests {
        use super::*;

        const EPS: f64 = 1e-12;

        /// 3-4-5: norm of (3, 4) is 5.
        #[test]
        fn three_four_five() {
            let v = [3.0_f64, 4.0];
            assert!((*vector_norm_l2(&v).as_ref() - 5.0).abs() < EPS);
        }

        /// Zero vector has norm 0.
        #[test]
        fn zero_vector() {
            let v = [0.0_f64; 4];
            assert_eq!(*vector_norm_l2(&v).as_ref(), 0.0);
        }

        /// 1-D: norm is |v[0]|.
        #[test]
        fn one_dimension_positive() {
            let v = [5.0_f64];
            assert!((*vector_norm_l2(&v).as_ref() - 5.0).abs() < EPS);
        }

        /// 1-D: norm of negative scalar is its absolute value.
        #[test]
        fn one_dimension_negative() {
            let v = [-7.0_f64];
            assert!((*vector_norm_l2(&v).as_ref() - 7.0).abs() < EPS);
        }

        /// Unit axis vectors always have norm 1.
        #[test]
        fn unit_axis_vectors_3d() {
            for axis in 0..3 {
                let mut e = [0.0_f64; 3];
                e[axis] = 1.0;
                assert!((*vector_norm_l2(&e).as_ref() - 1.0).abs() < EPS);
            }
        }

        /// Known 3-D norm: sqrt(1² + 2² + 2²) = 3.
        #[test]
        fn known_3d_norm() {
            let v = [1.0_f64, 2.0, 2.0];
            assert!((*vector_norm_l2(&v).as_ref() - 3.0).abs() < EPS);
        }

        /// Negating a vector preserves its norm.
        #[test]
        fn negation_preserves_norm() {
            let v = [3.0_f64, -4.0, 5.0];
            let neg_v = [-3.0_f64, 4.0, -5.0];
            assert_eq!(vector_norm_l2(&v), vector_norm_l2(&neg_v));
        }

        /// Numerical stability: component near f64::MAX must not overflow.
        #[test]
        fn numerical_stability_large_values() {
            let huge = f64::MAX / 2.0;
            let v = [huge, 0.0_f64];
            let n = *vector_norm_l2(&v).as_ref();
            assert!(n.is_finite(), "norm should be finite, got {n}");
            assert!((n / huge - 1.0).abs() < 1e-10);
        }

        /// Numerical stability: subnormal component must not underflow to zero.
        #[test]
        fn numerical_stability_tiny_values() {
            let tiny = f64::MIN_POSITIVE;
            let v = [tiny, 0.0_f64];
            let n = *vector_norm_l2(&v).as_ref();
            assert!(n.is_finite());
            assert!(n > 0.0);
        }
    }

    mod vector_norm_l2_sq_tests {
        use super::*;
        use crate::RealNative64StrictFinite;

        const EPS: f64 = 1e-16;

        /// 3-4-5: vector_norm_l2_sq([3, 4]) == 3² + 4² == 25.
        #[test]
        fn three_four_five() {
            let v = [3.0_f64, 4.0];
            assert_eq!(*vector_norm_l2_sq(&v).as_ref(), 25.0);
        }

        /// Zero vector: squared norm is 0.
        #[test]
        fn zero_vector() {
            let v = [0.0_f64; 4];
            assert_eq!(*vector_norm_l2_sq(&v).as_ref(), 0.0);
        }

        /// 1-D positive: vector_norm_l2_sq([x]) == x².
        #[test]
        fn one_dimension_positive() {
            let v = [5.0_f64];
            assert_eq!(*vector_norm_l2_sq(&v).as_ref(), 25.0);
        }

        /// 1-D negative: vector_norm_l2_sq([x]) == x² (sign disappears).
        #[test]
        fn one_dimension_negative() {
            let v = [-7.0_f64];
            assert_eq!(*vector_norm_l2_sq(&v).as_ref(), 49.0);
        }

        /// Unit axis vectors always have squared norm 1.
        #[test]
        fn unit_axis_vectors_3d() {
            for axis in 0..3 {
                let mut e = [0.0_f64; 3];
                e[axis] = 1.0;
                assert_eq!(*vector_norm_l2_sq(&e).as_ref(), 1.0);
            }
        }

        /// Known 3-D: 1² + 2² + 2² = 1 + 4 + 4 = 9.
        #[test]
        fn known_3d_vector_norm_sq() {
            let v = [1.0_f64, 2.0, 2.0];
            assert_eq!(*vector_norm_l2_sq(&v).as_ref(), 9.0);
        }

        /// Negating a vector preserves its squared norm.
        #[test]
        fn negation_preserves_vector_norm_sq() {
            let v = [3.0_f64, -4.0, 5.0];
            let neg_v = [-3.0_f64, 4.0, -5.0];
            assert_eq!(vector_norm_l2_sq(&v), vector_norm_l2_sq(&neg_v));
        }

        /// vector_norm_l2_sq == vector_norm_l2² (consistency with vector_norm_l2).
        #[test]
        fn consistent_with_vector_norm_l2() {
            let v = [2.0_f64, 3.0, 6.0];
            let sq = *vector_norm_l2_sq(&v).as_ref();
            let n = *vector_norm_l2(&v).as_ref();
            assert_eq!(sq, n * n);
        }

        /// Numerical stability: component below sqrt(f64::MAX) ≈ 1.34e154 must not overflow.
        /// Note: vector_norm_l2_sq computes scale² internally; values above sqrt(f64::MAX) will overflow.
        #[test]
        fn numerical_stability_large_values() {
            // Safe upper bound: scale must satisfy scale² < f64::MAX
            let large = 1.0e150_f64;
            let v = [large, 0.0_f64];
            let sq = *vector_norm_l2_sq(&v).as_ref();
            assert!(sq.is_finite(), "squared norm should be finite, got {sq}");
            assert_eq!(sq / (large * large), 1.0);
        }

        /// Numerical stability: component above sqrt(f64::MIN_POSITIVE) ≈ 1.05e-154 must not underflow.
        /// Note: vector_norm_l2_sq computes scale² internally; values below sqrt(f64::MIN_POSITIVE) will underflow to 0.
        #[test]
        fn numerical_stability_tiny_values() {
            // Safe lower bound: scale must satisfy scale² > 0 (no underflow)
            let small = 1.0e-100_f64;
            let v = [small, 0.0_f64];
            let sq = *vector_norm_l2_sq(&v).as_ref();
            assert!(sq.is_finite());
            assert!(sq > 0.0);
        }

        /// Overflow is expected when the mathematically correct squared norm is
        /// not representable (e.g. f64::MAX squared).
        ///
        /// We use strict validated scalars so overflow is rejected in both debug
        /// and release builds.
        #[test]
        #[should_panic]
        fn overflow_panics_slice_api() {
            let v = [RealNative64StrictFinite::from_f64(f64::MAX)];
            let _ = vector_norm_l2_sq(&v);
        }

        /// Iterator API must fail in the same way when final scale² overflows.
        #[test]
        #[should_panic]
        fn overflow_panics_iter_api() {
            let v = vec![RealNative64StrictFinite::from_f64(f64::MAX)];
            let _ = vector_norm_l2_sq_iter(v);
        }

        /// Accumulator-driven API must also panic on non-finite final output.
        #[test]
        #[should_panic]
        fn overflow_panics_accumulator_api() {
            let v = [RealNative64StrictFinite::from_f64(f64::MAX)];
            let _ = vector_norm_l2_sq_with_accumulator::<_, NaiveSum<_>>(&v);
        }

        /// vector_norm_l2_sq is always ≥ vector_norm_linf² (since ‖v‖∞ ≤ ‖v‖₂).
        #[test]
        fn geq_vector_norm_linf_squared() {
            let v = [3.0_f64, -4.0, 5.0];
            let sq = *vector_norm_l2_sq(&v).as_ref();
            let linf = *vector_norm_linf(&v).as_ref();
            assert!(sq >= linf * linf - EPS);
        }

        /// Documents the inaccuracy of naive accumulation for small squared terms.
        ///
        /// v = [1.0, δ, δ, …, δ] with δ = 1e-9, N = 10_000.
        /// True ‖v‖₂² = 1.0 + N·δ² = 1.0 + 1e-14.
        ///
        /// The scaled algorithm keeps `sumsq ≈ 1.0` throughout; each δ² = 1e-18
        /// is well below the floating-point precision of 1.0 (ε ≈ 2.2e-16), so
        /// `sumsq += δ²` rounds back to `sumsq` unchanged and every small term is
        /// silently discarded. The final result equals exactly 1.0, an error of
        /// 1e-14 — four orders of magnitude above machine epsilon.
        ///
        /// With Neumaier compensated summation the δ² corrections would be captured
        /// in the running compensation term and the result would satisfy
        /// `(result - expected).abs() < ε·expected ≈ 2.2e-16`.
        #[test]
        fn naive_accumulation_loses_small_squared_terms() {
            const N: usize = 10_000;
            let mut v = [1e-9_f64; N + 1];
            v[0] = 1.0;

            let result = *vector_norm_l2_sq(&v).as_ref();
            let expected = 1.0_f64 + N as f64 * 1e-18; // = 1.0 + 1e-14

            // Current: every δ² = 1e-18 < ε is lost → result collapses to 1.0.
            assert_eq!(result, 1.0_f64);

            // The relative error equals the entire small-terms contribution (~1e-14),
            // far larger than machine epsilon (~2.2e-16).
            let relative_error = (result - expected).abs() / expected;
            assert!(
                relative_error > 1e-15,
                "relative error {relative_error:.2e} should be significant \
                 (would be ~ε ≈ 2.2e-16 with Neumaier-sum)"
            );

            // With Neumaier the following assertion would hold instead:
            // assert!((result - expected).abs() < 1e-15 * expected);
        }
    }

    mod vector_norm_l2_sq_neumaier_tests {
        use super::*;

        const EPS: f64 = f64::EPSILON;

        /// 3-4-5: vector_norm_l2_sq_neumaier([3, 4]) == 3² + 4² == 25.
        #[test]
        fn three_four_five() {
            let v = [3.0_f64, 4.0];
            assert_eq!(*vector_norm_l2_sq_neumaier(&v).as_ref(), 25.0);
        }

        /// Zero vector → 0.
        #[test]
        fn zero_vector() {
            let v = [0.0_f64; 4];
            assert_eq!(*vector_norm_l2_sq_neumaier(&v).as_ref(), 0.0);
        }

        /// 1-D: vector_norm_l2_sq_neumaier([x]) == x².
        #[test]
        fn one_dimension() {
            let v = [7.0_f64];
            assert_eq!(*vector_norm_l2_sq_neumaier(&v).as_ref(), 49.0);
        }

        /// Matches vector_norm_l2_sq on a generic 3-D vector.
        #[test]
        fn matches_naive_for_normal_values() {
            let v = [1.0_f64, 2.0, 3.0];
            let neumaier = *vector_norm_l2_sq_neumaier(&v).as_ref();
            let naive = *vector_norm_l2_sq(&v).as_ref();
            assert!((neumaier - naive).abs() < EPS);
        }

        /// Neumaier version recovers small squared terms that naive accumulation loses.
        ///
        /// v = [1.0, δ × N] where δ = 1e-9. True ‖v‖₂² = 1.0 + N·δ² = 1.0 + 1e-14.
        /// Naive: every δ² = 1e-18 < ε is silently dropped → result = 1.0.
        /// Neumaier: compensation accumulates all N·δ² → result ≈ 1.0 + 1e-14.
        #[test]
        fn recovers_small_squared_terms() {
            const N: usize = 10_000;
            let mut v = [1e-9_f64; N + 1];
            v[0] = 1.0;

            let result = *vector_norm_l2_sq_neumaier(&v).as_ref();
            let expected = 1.0_f64 + N as f64 * 1e-18; // = 1.0 + 1e-14

            // The Neumaier-compensated result must be within machine-epsilon of expected.
            let relative_error = (result - expected).abs() / expected;
            assert!(
                relative_error < EPS,
                "relative error {relative_error:.2e} should be tiny with Neumaier-sum"
            );

            // Demonstrate the contrast with naive accumulation.
            let naive = *vector_norm_l2_sq(&v).as_ref();
            assert_eq!(
                naive, 1.0_f64,
                "naive should collapse to 1.0 (known limitation)"
            );
            assert!(
                result > naive,
                "Neumaier result {result} should exceed naive {naive}"
            );
        }

        /// Numerical stability: large-valued input does not overflow.
        #[test]
        fn numerical_stability_large_values() {
            let v = [1e150_f64, 1e150, 1e150];
            let result = *vector_norm_l2_sq_neumaier(&v).as_ref();
            let expected = 3.0 * 1e300_f64;
            let relative_error = (result - expected).abs() / expected;
            assert!(relative_error < EPS, "large values: {relative_error:.2e}");
        }

        /// Numerical stability: tiny-valued input does not underflow to zero.
        #[test]
        fn numerical_stability_tiny_values() {
            let v = [1e-100_f64, 1e-100, 1e-100];
            let result = *vector_norm_l2_sq_neumaier(&v).as_ref();
            let expected = 3.0e-200_f64;
            let relative_error = (result - expected).abs() / expected;
            assert!(relative_error < EPS, "tiny values: {relative_error:.2e}");
        }

        /// Consistent with vector_norm_l2: vector_norm_l2_sq_neumaier(v) ≈ vector_norm_l2(v)².
        #[test]
        fn consistent_with_vector_norm_l2() {
            let v = [1.5_f64, 2.5, 3.5, 4.5];
            let sq = *vector_norm_l2_sq_neumaier(&v).as_ref();
            let l2 = *vector_norm_l2(&v).as_ref();
            let err = (sq - l2 * l2).abs();
            println!("sq - l2² = {}", err);
            assert_eq!(err, 0.7105427357601002e-14);
        }
    }

    mod vector_norm_l2_neumaier_tests {
        use super::*;

        const EPS: f64 = f64::EPSILON;

        /// 3-4-5 right triangle: ‖(3, 4)‖₂ == 5.
        #[test]
        fn three_four_five() {
            let v = [3.0_f64, 4.0];
            assert_eq!(*vector_norm_l2_neumaier(&v).as_ref(), 5.0);
        }

        /// Zero vector → 0.
        #[test]
        fn zero_vector() {
            let v = [0.0_f64; 4];
            assert_eq!(*vector_norm_l2_neumaier(&v).as_ref(), 0.0);
        }

        /// 1-D: ‖[x]‖₂ == |x|.
        #[test]
        fn one_dimension() {
            let v = [7.0_f64];
            assert_eq!(*vector_norm_l2_neumaier(&v).as_ref(), 7.0);
        }

        /// Matches vector_norm_l2 on ordinary values (no dominant small-term cancellation).
        #[test]
        fn matches_naive_for_normal_values() {
            let v = [1.0_f64, 2.0, 3.0];
            let neumaier = *vector_norm_l2_neumaier(&v).as_ref();
            let naive = *vector_norm_l2(&v).as_ref();
            assert!((neumaier - naive).abs() < EPS);
        }

        /// Neumaier version recovers small squared terms that naive accumulation loses.
        ///
        /// v = [1.0, δ × N] where δ = 1e-9. True ‖v‖₂ = √(1.0 + N·δ²) = √(1.0 + 1e-14).
        /// Naive: sumsq collapses to 1.0 → norm = 1.0 (all small terms lost).
        /// Neumaier: sumsq ≈ 1.0 + 1e-14 → norm ≈ √(1.0 + 1e-14) ≈ 1.0 + 5e-15.
        #[test]
        fn recovers_small_squared_terms() {
            const N: usize = 10_000;
            let mut v = [1e-9_f64; N + 1];
            v[0] = 1.0;

            let result = *vector_norm_l2_neumaier(&v).as_ref();
            let expected = (1.0_f64 + N as f64 * 1e-18).sqrt();

            let relative_error = (result - expected).abs() / expected;
            assert!(
                relative_error < 1e-12,
                "relative error {relative_error:.2e} should be tiny with Neumaier-sum"
            );

            // Contrast with naive: vector_norm_l2 loses the small terms, giving exactly 1.0.
            let naive = *vector_norm_l2(&v).as_ref();
            assert_eq!(
                naive, 1.0_f64,
                "naive should collapse to 1.0 (known limitation)"
            );
            assert!(
                result > naive,
                "Neumaier result {result} should exceed naive {naive}"
            );
        }

        /// Numerical stability: large-valued input does not overflow.
        #[test]
        fn numerical_stability_large_values() {
            let v = [1e150_f64, 1e150, 1e150];
            let result = *vector_norm_l2_neumaier(&v).as_ref();
            let expected = (3.0_f64).sqrt() * 1e150;
            let relative_error = (result - expected).abs() / expected;
            assert!(relative_error < EPS, "large values: {relative_error:.2e}");
        }

        /// Numerical stability: tiny-valued input does not underflow to zero.
        #[test]
        fn numerical_stability_tiny_values() {
            let v = [1e-150_f64, 1e-150, 1e-150];
            let result = *vector_norm_l2_neumaier(&v).as_ref();
            let expected = (3.0_f64).sqrt() * 1e-150;
            let relative_error = (result - expected).abs() / expected;
            assert!(relative_error < EPS, "tiny values: {relative_error:.2e}");
        }

        /// Consistent with vector_norm_l2_sq_neumaier: vector_norm_l2_neumaier(v)² ≈ vector_norm_l2_sq_neumaier(v).
        #[test]
        fn consistent_with_vector_norm_l2_sq_neumaier() {
            let v = [1.5_f64, 2.5, 3.5, 4.5];
            let norm = *vector_norm_l2_neumaier(&v).as_ref();
            let sq = *vector_norm_l2_sq_neumaier(&v).as_ref();
            assert!(
                (norm * norm - sq).abs() < 100. * EPS,
                "norm²={} sq={}",
                norm * norm,
                sq
            );
        }

        /// Unit-axis vectors have norm 1.
        #[test]
        fn unit_axis_vectors() {
            let e1 = [1.0_f64, 0.0, 0.0];
            let e2 = [0.0_f64, 1.0, 0.0];
            let e3 = [0.0_f64, 0.0, 1.0];
            for e in [&e1, &e2, &e3] {
                let n = *vector_norm_l2_neumaier(e).as_ref();
                assert!((n - 1.0).abs() < EPS, "unit axis: {n}");
            }
        }
    }

    mod vector_norm_l1_tests {
        use super::*;

        const EPS: f64 = 1e-12;

        /// L1 norm of (3, 4) is 3 + 4 = 7.
        #[test]
        fn three_plus_four() {
            let v = [3.0_f64, 4.0];
            assert_eq!(*vector_norm_l1(&v).as_ref(), 7.0);
        }

        /// Zero vector has L1 norm 0.
        #[test]
        fn zero_vector() {
            let v = [0.0_f64; 4];
            assert_eq!(*vector_norm_l1(&v).as_ref(), 0.0);
        }

        /// 1-D: norm is |v[0]|.
        #[test]
        fn one_dimension_positive() {
            let v = [5.0_f64];
            assert!((*vector_norm_l1(&v).as_ref() - 5.0).abs() < EPS);
        }

        /// 1-D: norm of negative scalar is its absolute value.
        #[test]
        fn one_dimension_negative() {
            let v = [-7.0_f64];
            assert!((*vector_norm_l1(&v).as_ref() - 7.0).abs() < EPS);
        }

        /// Unit axis vectors always have L1 norm 1.
        #[test]
        fn unit_axis_vectors_3d() {
            for axis in 0..3 {
                let mut e = [0.0_f64; 3];
                e[axis] = 1.0;
                assert_eq!(*vector_norm_l1(&e).as_ref(), 1.0);
            }
        }

        /// Known 3-D: |1| + |2| + |2| = 5.
        #[test]
        fn known_3d_norm() {
            let v = [1.0_f64, 2.0, 2.0];
            assert!((*vector_norm_l1(&v).as_ref() - 5.0).abs() < EPS);
        }

        /// Negating a vector preserves its L1 norm.
        #[test]
        fn negation_preserves_norm() {
            let v = [3.0_f64, -4.0, 5.0];
            let neg_v = [-3.0_f64, 4.0, -5.0];
            assert_eq!(vector_norm_l1(&v), vector_norm_l1(&neg_v));
        }

        /// L1 norm is the sum of absolute values: mixed signs.
        #[test]
        fn mixed_signs() {
            let v = [-1.0_f64, 2.0, -3.0, 4.0];
            assert!((*vector_norm_l1(&v).as_ref() - 10.0).abs() < EPS);
        }

        /// L1 ≥ L2 for any vector.
        #[test]
        fn l1_geq_l2() {
            let v = [3.0_f64, -4.0, 5.0];
            assert!(vector_norm_l1(&v).as_ref() >= vector_norm_l2(&v).as_ref());
        }

        /// Numerical stability: large components must not overflow.
        #[test]
        fn numerical_stability_large_values() {
            let huge = f64::MAX / 4.0;
            // Two equal components: L1 norm = 2 × huge
            let v = [huge, huge];
            let n = *vector_norm_l1(&v).as_ref();
            assert!(n.is_finite(), "norm should be finite, got {n}");
            assert!((n / (2.0 * huge) - 1.0).abs() < 1e-10);
        }

        /// Numerical stability: subnormal component must not underflow to zero.
        #[test]
        fn numerical_stability_tiny_values() {
            let tiny = f64::MIN_POSITIVE;
            let v = [tiny, 0.0_f64];
            let n = *vector_norm_l1(&v).as_ref();
            assert!(n.is_finite());
            assert!(n > 0.0);
        }
    }

    mod vector_norm_l1_neumaier_tests {
        use super::*;

        const EPS: f64 = 1e-12;

        /// L1 norm of (3, 4) is 3 + 4 = 7.
        #[test]
        fn three_plus_four() {
            let v = [3.0_f64, 4.0];
            assert_eq!(*vector_norm_l1_neumaier(&v).as_ref(), 7.0);
        }

        /// Zero vector has L1 norm 0.
        #[test]
        fn zero_vector() {
            let v = [0.0_f64; 4];
            assert_eq!(*vector_norm_l1_neumaier(&v).as_ref(), 0.0);
        }

        /// 1-D: norm is |v[0]|.
        #[test]
        fn one_dimension_positive() {
            let v = [5.0_f64];
            assert!((*vector_norm_l1_neumaier(&v).as_ref() - 5.0).abs() < EPS);
        }

        /// 1-D: norm of negative scalar is its absolute value.
        #[test]
        fn one_dimension_negative() {
            let v = [-7.0_f64];
            assert!((*vector_norm_l1_neumaier(&v).as_ref() - 7.0).abs() < EPS);
        }

        /// Unit axis vectors always have L1 norm 1.
        #[test]
        fn unit_axis_vectors_3d() {
            for axis in 0..3 {
                let mut e = [0.0_f64; 3];
                e[axis] = 1.0;
                assert_eq!(*vector_norm_l1_neumaier(&e).as_ref(), 1.0);
            }
        }

        /// Known 3-D: |1| + |2| + |2| = 5.
        #[test]
        fn known_3d_norm() {
            let v = [1.0_f64, 2.0, 2.0];
            assert!((*vector_norm_l1_neumaier(&v).as_ref() - 5.0).abs() < EPS);
        }

        /// Negating a vector preserves its L1 norm.
        #[test]
        fn negation_preserves_norm() {
            let v = [3.0_f64, -4.0, 5.0];
            let neg_v = [-3.0_f64, 4.0, -5.0];
            assert_eq!(vector_norm_l1_neumaier(&v), vector_norm_l1_neumaier(&neg_v));
        }

        /// L1 norm is the sum of absolute values: mixed signs.
        #[test]
        fn mixed_signs() {
            let v = [-1.0_f64, 2.0, -3.0, 4.0];
            assert!((*vector_norm_l1_neumaier(&v).as_ref() - 10.0).abs() < EPS);
        }

        /// L1 ≥ L2 for any vector.
        #[test]
        fn l1_geq_l2() {
            let v = [3.0_f64, -4.0, 5.0];
            assert!(vector_norm_l1_neumaier(&v).as_ref() >= vector_norm_l2(&v).as_ref());
        }

        /// Numerical stability: large components must not overflow.
        #[test]
        fn numerical_stability_large_values() {
            let huge = f64::MAX / 4.0;
            // Two equal components: L1 norm = 2 × huge
            let v = [huge, huge];
            let n = *vector_norm_l1_neumaier(&v).as_ref();
            assert!(n.is_finite(), "norm should be finite, got {n}");
            assert!((n / (2.0 * huge) - 1.0).abs() < 1e-10);
        }

        /// Numerical stability: subnormal component must not underflow to zero.
        #[test]
        fn numerical_stability_tiny_values() {
            let tiny = f64::MIN_POSITIVE;
            let v = [tiny, 0.0_f64];
            let n = *vector_norm_l1_neumaier(&v).as_ref();
            assert!(n.is_finite());
            assert!(n > 0.0);
        }
    }

    mod vector_norm_linf_tests {
        use super::*;

        const EPS: f64 = 1e-12;

        #[test]
        fn parallel_slice_matches_sequential() {
            let v = [3.0_f64, -4.0, 5.0, -2.5, 1.0];
            assert_eq!(vector_norm_linf_par(&v), vector_norm_linf(&v));
        }

        #[test]
        fn parallel_iter_matches_sequential_iter() {
            let v = vec![3.0_f64, -4.0, 5.0, -2.5, 1.0];
            let par = vector_norm_linf_par_iter(v.clone());
            let seq = vector_norm_linf_iter(v);
            assert_eq!(par, seq);
        }

        /// L∞ norm of (3, 4) = max(3, 4) = 4.
        #[test]
        fn max_of_three_four() {
            let v = [3.0_f64, 4.0];
            assert_eq!(*vector_norm_linf(&v).as_ref(), 4.0);
        }

        /// Zero vector has L∞ norm 0.
        #[test]
        fn zero_vector() {
            let v = [0.0_f64; 4];
            assert_eq!(*vector_norm_linf(&v).as_ref(), 0.0);
        }

        /// 1-D: norm is |v[0]|.
        #[test]
        fn one_dimension_positive() {
            let v = [5.0_f64];
            assert!((*vector_norm_linf(&v).as_ref() - 5.0).abs() < EPS);
        }

        /// 1-D: norm of negative scalar is its absolute value.
        #[test]
        fn one_dimension_negative() {
            let v = [-7.0_f64];
            assert!((*vector_norm_linf(&v).as_ref() - 7.0).abs() < EPS);
        }

        /// Unit axis vectors always have L∞ norm 1.
        #[test]
        fn unit_axis_vectors_3d() {
            for axis in 0..3 {
                let mut e = [0.0_f64; 3];
                e[axis] = 1.0;
                assert_eq!(*vector_norm_linf(&e).as_ref(), 1.0);
            }
        }

        /// Known 3-D: max(|1|, |-2|, |3|) = 3.
        #[test]
        fn known_3d_norm() {
            let v = [1.0_f64, -2.0, 3.0];
            assert!((*vector_norm_linf(&v).as_ref() - 3.0).abs() < EPS);
        }

        /// Negating a vector preserves its L∞ norm.
        #[test]
        fn negation_preserves_norm() {
            let v = [3.0_f64, -4.0, 5.0];
            let neg_v = [-3.0_f64, 4.0, -5.0];
            assert_eq!(vector_norm_linf(&v), vector_norm_linf(&neg_v));
        }

        /// The dominant component determines the result.
        #[test]
        fn dominant_component() {
            let v = [1.0_f64, 0.001, 100.0, 0.5];
            assert!((*vector_norm_linf(&v).as_ref() - 100.0).abs() < EPS);
        }

        /// L∞ ≤ L2 for any vector.
        #[test]
        fn linf_leq_l2() {
            let v = [3.0_f64, -4.0, 5.0];
            assert!(vector_norm_linf(&v).as_ref() <= vector_norm_l2(&v).as_ref());
        }

        /// L∞ ≤ L1 for any vector.
        #[test]
        fn linf_leq_l1() {
            let v = [3.0_f64, -4.0, 5.0];
            assert!(vector_norm_linf(&v).as_ref() <= vector_norm_l1(&v).as_ref());
        }

        /// Numerical stability: large component must remain finite.
        #[test]
        fn numerical_stability_large_values() {
            let huge = f64::MAX / 2.0;
            let v = [1.0_f64, huge];
            let n = *vector_norm_linf(&v).as_ref();
            assert!(n.is_finite(), "norm should be finite, got {n}");
            assert!((n / huge - 1.0).abs() < 1e-10);
        }

        /// Numerical stability: subnormal component must not underflow to zero.
        #[test]
        fn numerical_stability_tiny_values() {
            let tiny = f64::MIN_POSITIVE;
            let v = [tiny, 0.0_f64];
            let n = *vector_norm_linf(&v).as_ref();
            assert!(n.is_finite());
            assert!(n > 0.0);
        }
    }

    mod other_vector_norm_l2_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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                assert_eq!(*norm.as_ref(), 0.0);
            }

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

            #[test]
            fn single_negative_element() {
                let data = vec_f64(&[-7.0]);
                let norm = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                assert_relative_eq!(*norm.as_ref(), 5.0, epsilon = 1e-15);
            }

            #[test]
            fn single_zero() {
                let data = vec_f64(&[0.0]);
                let norm = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                // 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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                assert_eq!(*norm.as_ref(), f64::MIN_POSITIVE);
            }

            #[test]
            fn epsilon() {
                let data = vec_f64(&[f64::EPSILON]);
                let norm = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                let five = rug_f64(5.0);

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

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

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

            #[test]
            fn single_negative() {
                let data = vec![rug_f64(-7.0)];
                let norm = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();
                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 = vector_norm_l2(&data).into_inner();

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

                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 = vector_norm_l2(&data).into_inner();

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

                let rel_diff = ((norm.clone() - &expected) / &expected).abs().into_inner();
                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 = vector_norm_l2(&data).into_inner();

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

                let rel_diff = ((norm.clone() - &expected) / &expected).abs().into_inner();
                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 = vector_norm_l2(&data).into_inner();
                let five = rug_f64(5.0);

                let diff = (norm.clone() - &five).abs().into_inner();
                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 = vector_norm_l2(&data).into_inner();
                let five = rug_f64(5.0);

                let diff = (norm.clone() - &five).abs().into_inner();
                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 = vector_norm_l2(&data).into_inner();

                let diff = (norm.clone() - &expected).abs().into_inner();
                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 = vector_norm_l2(&data).into_inner();
                let five = hp_f64(5.0);

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

    // ── Complex scalar tests ─────────────────────────────────────────────────

    mod complex_types {
        use super::*;
        use crate::{ComplexNative64StrictFinite, RealNative64StrictFinite};
        use num::Complex;
        use try_create::TryNew;

        /// Helper: build a `ComplexNative64StrictFinite` from two f64 components.
        fn c(re: f64, im: f64) -> ComplexNative64StrictFinite {
            ComplexNative64StrictFinite::try_new(Complex::new(re, im)).unwrap()
        }

        /// Helper: extract the inner f64 from a `NonNegativeRealScalar<RealNative64StrictFinite>`.
        fn inner_f64(x: crate::scalars::NonNegativeRealScalar<RealNative64StrictFinite>) -> f64 {
            *x.into_inner().as_ref()
        }

        mod vector_norm_l2_complex {
            use super::*;

            /// ‖[3+0i, 0+4i]‖₂ = √(3² + 4²) = 5
            #[test]
            fn components_on_axes() {
                let v = [c(3.0, 0.0), c(0.0, 4.0)];
                assert!((inner_f64(vector_norm_l2(&v)) - 5.0).abs() < 1e-12);
            }

            /// ‖[3+4i]‖₂ = |3+4i| = 5
            #[test]
            fn single_complex() {
                let v = [c(3.0, 4.0)];
                assert!((inner_f64(vector_norm_l2(&v)) - 5.0).abs() < 1e-12);
            }

            /// ‖[0+0i, 0+0i]‖₂ = 0
            #[test]
            fn zero_vector() {
                let v = [c(0.0, 0.0), c(0.0, 0.0)];
                assert_eq!(inner_f64(vector_norm_l2(&v)), 0.0);
            }

            /// Negating imaginary part does not change the norm.
            #[test]
            fn conjugate_preserves_norm() {
                let v = [c(1.0, 2.0), c(3.0, -4.0)];
                let vc = [c(1.0, -2.0), c(3.0, 4.0)];
                assert!(
                    (inner_f64(vector_norm_l2(&v)) - inner_f64(vector_norm_l2(&vc))).abs() < 1e-14
                );
            }

            /// vector_norm_l2 and vector_norm_l2_neumaier agree on ordinary values.
            #[test]
            fn matches_neumaier() {
                let v = [c(1.0, 2.0), c(3.0, -1.0), c(-1.0, 4.0)];
                let diff =
                    (inner_f64(vector_norm_l2(&v)) - inner_f64(vector_norm_l2_neumaier(&v))).abs();
                assert!(diff < 1e-14);
            }
        }

        mod vector_norm_l2_sq_complex {
            use super::*;

            /// ‖[3+0i, 0+4i]‖₂² = 9 + 16 = 25
            #[test]
            fn components_on_axes() {
                let v = [c(3.0, 0.0), c(0.0, 4.0)];
                assert!((inner_f64(vector_norm_l2_sq(&v)) - 25.0).abs() < 1e-12);
            }

            /// ‖[3+4i]‖₂² = |3+4i|² = 25
            #[test]
            fn single_complex() {
                let v = [c(3.0, 4.0)];
                assert!((inner_f64(vector_norm_l2_sq(&v)) - 25.0).abs() < 1e-12);
            }

            /// vector_norm_l2_sq(v) ≈ vector_norm_l2(v)²
            #[test]
            fn consistent_with_vector_norm_l2() {
                let v = [c(1.0, 2.0), c(3.0, 4.0)];
                let l2 = inner_f64(vector_norm_l2(&v));
                let sq = inner_f64(vector_norm_l2_sq(&v));
                assert!((sq - l2 * l2).abs() < 1e-12);
            }

            /// vector_norm_l2_sq and vector_norm_l2_sq_neumaier agree on ordinary values.
            #[test]
            fn matches_neumaier() {
                let v = [c(3.0, 0.0), c(0.0, 4.0)];
                let diff = (inner_f64(vector_norm_l2_sq(&v))
                    - inner_f64(vector_norm_l2_sq_neumaier(&v)))
                .abs();
                assert!(diff < 1e-14);
            }
        }

        mod vector_norm_l2_sq_neumaier_complex {
            use super::*;

            #[test]
            fn zero_vector() {
                let v = [c(0.0, 0.0), c(0.0, 0.0)];
                assert_eq!(inner_f64(vector_norm_l2_sq_neumaier(&v)), 0.0);
            }

            #[test]
            fn matches_naive() {
                let v = [c(2.0, 1.0), c(0.0, 3.0), c(-1.0, -1.0)];
                let diff = (inner_f64(vector_norm_l2_sq_neumaier(&v))
                    - inner_f64(vector_norm_l2_sq(&v)))
                .abs();
                assert!(diff < 1e-14);
            }
        }

        mod vector_norm_l2_neumaier_complex {
            use super::*;

            #[test]
            fn zero_vector() {
                let v = [c(0.0, 0.0)];
                assert_eq!(inner_f64(vector_norm_l2_neumaier(&v)), 0.0);
            }

            #[test]
            fn matches_naive() {
                let v = [c(2.0, 1.0), c(0.0, 3.0), c(-1.0, -1.0)];
                let diff =
                    (inner_f64(vector_norm_l2_neumaier(&v)) - inner_f64(vector_norm_l2(&v))).abs();
                assert!(diff < 1e-14);
            }
        }

        mod vector_norm_l1_complex {
            use super::*;

            /// ‖[3+0i, 0+4i]‖₁ = |3+0i| + |0+4i| = 3 + 4 = 7
            #[test]
            fn components_on_axes() {
                let v = [c(3.0, 0.0), c(0.0, 4.0)];
                assert!((inner_f64(vector_norm_l1(&v)) - 7.0).abs() < 1e-12);
            }

            /// ‖[3+4i]‖₁ = |3+4i| = 5
            #[test]
            fn single_complex() {
                let v = [c(3.0, 4.0)];
                assert!((inner_f64(vector_norm_l1(&v)) - 5.0).abs() < 1e-12);
            }

            /// ‖v‖₁ ≥ ‖v‖₂ for multi-element vectors (triangle inequality corollary)
            #[test]
            fn l1_geq_l2() {
                let v = [c(3.0, 0.0), c(0.0, 4.0)];
                assert!(inner_f64(vector_norm_l1(&v)) >= inner_f64(vector_norm_l2(&v)) - 1e-14);
            }

            #[test]
            fn zero_vector() {
                let v = [c(0.0, 0.0), c(0.0, 0.0)];
                assert_eq!(inner_f64(vector_norm_l1(&v)), 0.0);
            }

            /// vector_norm_l1 and vector_norm_l1_neumaier agree on ordinary values.
            #[test]
            fn matches_neumaier() {
                let v = [c(1.0, 2.0), c(3.0, 4.0), c(5.0, 6.0)];
                let diff =
                    (inner_f64(vector_norm_l1(&v)) - inner_f64(vector_norm_l1_neumaier(&v))).abs();
                assert!(diff < 1e-14);
            }
        }

        mod vector_norm_l1_neumaier_complex {
            use super::*;

            #[test]
            fn components_on_axes() {
                let v = [c(3.0, 0.0), c(0.0, 4.0)];
                assert!((inner_f64(vector_norm_l1_neumaier(&v)) - 7.0).abs() < 1e-12);
            }

            #[test]
            fn zero_vector() {
                let v = [c(0.0, 0.0)];
                assert_eq!(inner_f64(vector_norm_l1_neumaier(&v)), 0.0);
            }
        }

        mod vector_norm_linf_complex {
            use super::*;

            /// ‖[3+0i, 0+4i]‖∞ = max(3, 4) = 4
            #[test]
            fn maximum_modulus() {
                let v = [c(3.0, 0.0), c(0.0, 4.0)];
                assert!((inner_f64(vector_norm_linf(&v)) - 4.0).abs() < 1e-12);
            }

            /// ‖[3+4i, 0+5i]‖∞ = max(5, 5) = 5 (equal moduli)
            #[test]
            fn equal_moduli() {
                let v = [c(3.0, 4.0), c(0.0, 5.0)];
                assert!((inner_f64(vector_norm_linf(&v)) - 5.0).abs() < 1e-12);
            }

            /// ‖v‖∞ ≤ ‖v‖₁ always holds
            #[test]
            fn linf_leq_l1() {
                let v = [c(1.0, 2.0), c(3.0, 4.0)];
                assert!(inner_f64(vector_norm_linf(&v)) <= inner_f64(vector_norm_l1(&v)) + 1e-14);
            }

            /// ‖v‖∞ ≤ ‖v‖₂ ≤ √n · ‖v‖∞ (standard norm inequality)
            #[test]
            fn linf_leq_l2() {
                let v = [c(1.0, 2.0), c(3.0, 4.0)];
                let linf = inner_f64(vector_norm_linf(&v));
                let l2 = inner_f64(vector_norm_l2(&v));
                let n = v.len() as f64;
                assert!(linf <= l2 + 1e-14);
                assert!(l2 <= n.sqrt() * linf + 1e-14);
            }

            #[test]
            fn zero_vector() {
                let v = [c(0.0, 0.0), c(0.0, 0.0)];
                assert_eq!(inner_f64(vector_norm_linf(&v)), 0.0);
            }

            #[test]
            fn parallel_matches_sequential() {
                let v = [c(1.0, 2.0), c(3.0, 4.0), c(-2.0, -1.0)];
                let seq = inner_f64(vector_norm_linf(&v));
                let par = inner_f64(vector_norm_linf_par(&v));
                assert!((seq - par).abs() < 1e-14);
            }
        }

        mod api_variants {
            use super::*;

            type L1Naive = NaiveSum<NonNegativeRealScalar<f64>>;
            type L2Naive = NaiveSum<f64>;

            #[test]
            fn exercises_iterator_and_parallel_variants() {
                let values = [3.0_f64, -4.0, 5.0];

                assert_eq!(
                    *L2Norm::<L2Naive>::compute_vector_norm_sq_par(&values).as_ref(),
                    50.0
                );
                assert_eq!(
                    *L2Norm::<L2Naive>::compute_vector_norm_sq_par_iter(values.par_iter().cloned())
                        .as_ref(),
                    50.0
                );

                assert_eq!(
                    *L2Norm::<L2Naive>::compute_vector_norm_iter(values.iter().cloned()).as_ref(),
                    50.0_f64.sqrt()
                );
                assert_eq!(
                    *L2Norm::<L2Naive>::compute_vector_norm_par(&values).as_ref(),
                    50.0_f64.sqrt()
                );
                assert_eq!(
                    *L2Norm::<L2Naive>::compute_vector_norm_par_iter(values.par_iter().cloned())
                        .as_ref(),
                    50.0_f64.sqrt()
                );

                assert_eq!(
                    *L1Norm::<L1Naive>::compute_vector_norm_iter(values.iter().cloned()).as_ref(),
                    12.0
                );
                assert_eq!(
                    *L1Norm::<L1Naive>::compute_vector_norm_par(&values).as_ref(),
                    12.0
                );
                assert_eq!(
                    *L1Norm::<L1Naive>::compute_vector_norm_par_iter(values.par_iter().cloned())
                        .as_ref(),
                    12.0
                );

                assert_eq!(
                    *LinfNorm::<()>::compute_vector_norm_iter(values.iter().cloned()).as_ref(),
                    5.0
                );
                assert_eq!(
                    *LinfNorm::<()>::compute_vector_norm_par(&values).as_ref(),
                    5.0
                );
                assert_eq!(
                    *LinfNorm::<()>::compute_vector_norm_par_iter(values.par_iter().cloned())
                        .as_ref(),
                    5.0
                );
            }

            #[test]
            fn exercises_free_iterator_and_accumulator_helpers() {
                let values = [3.0_f64, -4.0, 5.0];

                assert_eq!(
                    *vector_norm_l2_with_accumulator::<_, L2Naive>(&values).as_ref(),
                    50.0_f64.sqrt()
                );
                assert_eq!(
                    *vector_norm_l2_iter_with_accumulator::<_, L2Naive, _>(values.iter().cloned())
                        .as_ref(),
                    50.0_f64.sqrt()
                );
                assert_eq!(
                    *vector_norm_l2_sq_with_accumulator::<_, L2Naive>(&values).as_ref(),
                    50.0
                );
                assert_eq!(
                    *vector_norm_l2_sq_iter_with_accumulator::<_, L2Naive, _>(
                        values.iter().cloned()
                    )
                    .as_ref(),
                    50.0
                );
                assert_eq!(
                    *vector_norm_l1_with_accumulator::<_, L1Naive>(&values).as_ref(),
                    12.0
                );
                assert_eq!(
                    *vector_norm_l1_iter_with_accumulator::<_, L1Naive, _>(values.iter().cloned())
                        .as_ref(),
                    12.0
                );

                assert_eq!(
                    *vector_norm_l2_iter(values.iter().cloned()).as_ref(),
                    50.0_f64.sqrt()
                );
                assert_eq!(
                    *vector_norm_l2_sq_iter(values.iter().cloned()).as_ref(),
                    50.0
                );
                assert_eq!(
                    *vector_norm_l2_sq_neumaier_iter(values.iter().cloned()).as_ref(),
                    50.0
                );
                assert_eq!(
                    *vector_norm_l2_neumaier_iter(values.iter().cloned()).as_ref(),
                    50.0_f64.sqrt()
                );
                assert_eq!(*vector_norm_l1_iter(values.iter().cloned()).as_ref(), 12.0);
                assert_eq!(
                    *vector_norm_l1_neumaier_iter(values.iter().cloned()).as_ref(),
                    12.0
                );
                assert_eq!(*vector_norm_linf_iter(values.iter().cloned()).as_ref(), 5.0);
                assert_eq!(
                    *vector_norm_linf_par_iter(values.par_iter().cloned()).as_ref(),
                    5.0
                );
            }
        }
    }
}