num_valid/algorithms/

accumulators.rs

#![deny(rustdoc::broken_intra_doc_links)]

//! Pluggable accumulation strategies for numerical algorithms.
//!
//! This module defines the [`Accumulator`] trait and a set of concrete
//! implementations used internally by
//! [`crate::algorithms::vector_norms`] and available for use in custom
//! algorithms.
//!
//! # Trait hierarchy
//!
//! ```text
//! Accumulator
//! ├── SumAccumulator    (adds rescale_by; used for L1/L2 norms)
//! │   ├── NaiveSum
//! │   └── NeumaierSum
//! ├── MaxAccumulator    (running maximum of real values)
//! ├── MinAccumulator    (running minimum of real values)
//! ├── MaxAbsValueAccumulator  (running maximum of |x|; used for L∞ norm)
//! └── ScaledSumSqAccumulator  (scaled sum-of-squares with dynamic scaling; used for L2 norm)
//! ```
//!
//! # Choosing a summation strategy
//!
//! | Strategy | Speed | Accuracy | Use when |
//! |----------|-------|----------|----------|
//! | [`NaiveSum`] | fastest | O(N·ε) error | data is well-scaled |
//! | [`NeumaierSum`] | ~2× slower | O(ε) error | components span many magnitudes |
//!
//! # Parallel reduction
//!
//! All accumulators support parallel chunked reductions via the default
//! [`Accumulator::new_parallel`] method, which uses rayon `fold` + `reduce`
//! with [`Accumulator::combine`] as the merge operation. `Accumulator::new()`
//! must be the identity element for `combine` — this holds for all types in
//! this module.

use crate::{
    FpScalar, RealScalar,
    core::traits::validation::sealed,
    functions::{Abs, Pow},
    scalars::NonNegativeRealScalar,
};
use getset::Getters;
use num::{Complex, One, Zero};
use rayon::prelude::*;
use std::ops::{Add, AddAssign, Mul, MulAssign};

//---------------------------------------------------------
/// Base trait for incremental accumulation over a sequence of values.
///
/// An accumulator collects items one at a time via [`push`](Self::push) and
/// produces a final result via [`result`](Self::result). Partial accumulators
/// can be merged with [`combine`](Self::combine), enabling chunked and parallel
/// reductions.
///
/// # Correctness requirement for `combine`
///
/// [`new()`](Self::new) must be the identity element for `combine`: for any
/// accumulator `x`, both `combine(new(), x)` and `combine(x, new())` must
/// produce an accumulator equivalent to `x`. This invariant is required by
/// [`new_parallel`](Self::new_parallel).
pub trait Accumulator: Sized {
    /// The type of values pushed into the accumulator.
    type Input;
    /// The type of the final result produced by [`result`](Self::result).
    type Output;

    /// Creates a new, empty accumulator (the identity for [`combine`](Self::combine)).
    fn new() -> Self;

    /// Adds a single value to the accumulator.
    fn push(&mut self, value: Self::Input);

    /// Merges `other` into `self`.
    ///
    /// After this call `self` represents the accumulation of all values pushed
    /// into either `self` or `other`. Used by [`new_parallel`](Self::new_parallel)
    /// to merge per-chunk partial results.
    fn combine(&mut self, other: Self);

    /// Consumes the accumulator and returns the final result.
    fn result(self) -> Self::Output;

    /// Creates a new accumulator by sequentially adding all values from `values`.
    ///
    /// This is the simplest way to build an accumulator from owned input data.
    /// For chunked or parallel reductions, prefer [`Accumulator::combine`]
    /// on partial accumulators.
    ///
    /// # Example
    ///
    /// ```
    /// use num_valid::algorithms::accumulators::*;
    ///
    /// let values = vec![1.0, 1.0e100, 1.0, -1.0e100];
    ///
    /// let neumaier = NeumaierSum::new_sequential(values);
    /// let sum = neumaier.result();
    /// println!("Sum: {}", sum);
    /// assert_eq!(sum, 2.0);
    /// ```
    fn new_sequential<I>(values: I) -> Self
    where
        I: IntoIterator<Item = Self::Input>,
    {
        let mut acc = Self::new();
        values.into_iter().for_each(|v| acc.push(v));
        acc
    }

    /// Builds an accumulator by processing a parallel iterator.
    ///
    /// The iterator is partitioned into chunks by rayon; each chunk is folded
    /// into a partial accumulator with [`push`](Self::push), and the partial
    /// results are merged with [`combine`](Self::combine).
    ///
    /// Requires `Self: Send`.
    ///
    /// Any additional `Send` constraints on iterator items are enforced by the
    /// concrete [`IntoParallelIterator`] implementation.
    ///
    /// # Example
    ///
    /// ```
    /// use num_valid::algorithms::accumulators::*;
    /// use rayon::prelude::*;
    ///
    /// let values = vec![1.0_f64, 2.0, 3.0, 4.0];
    /// let sum = NaiveSum::new_parallel(values.into_par_iter()).result();
    /// assert_eq!(sum, 10.0);
    /// ```
    fn new_parallel<I>(values: I) -> Self
    where
        Self: Send,
        I: IntoParallelIterator<Item = Self::Input>,
    {
        values
            .into_par_iter()
            .fold(
                || Self::new(),
                |mut acc, item| {
                    acc.push(item);
                    acc
                },
            )
            .reduce(
                || Self::new(),
                |mut a, b| {
                    a.combine(b);
                    a
                },
            )
    }
}
//---------------------------------------------------------

//---------------------------------------------------------
/// A generic accumulator for summing a sequence of values with configurable
/// precision/performance trade-offs.
///
/// The trait abstracts over how terms are collected, allowing different
/// accumulation strategies to be swapped in without duplicating the calling
/// code.  It is used internally by the [`crate::algorithms::vector_norms`] vector-norm functions (scaled
/// sum-of-squares loop), but is general enough for any sequential summation
/// task.  Two strategies ship with the crate:
///
/// - **Naive** ([`NaiveSum`]): plain `+=` —
///   fast but accumulates O(N·ε) rounding error.
/// - **Neumaier** ([`NeumaierSum`]): Kahan/Neumaier compensated
///   summation — recovers O(ε) rounding error at a small overhead cost.
///
/// Implement this trait to plug a custom accumulation strategy wherever it is
/// accepted.
pub trait SumAccumulator:
    Accumulator<Input: Zero + Clone, Output = <Self as Accumulator>::Input>
{
    /// Multiply the entire accumulated sum by `r`.
    fn rescale_by(&mut self, r: &Self::Input);
}
//---------------------------------------------------------

//---------------------------------------------------------
/// Naive (plain `+=`) accumulator.
///
/// This strategy keeps a single running sum and updates it with standard
/// floating-point addition. It has minimal overhead and is usually the fastest
/// option, but it may lose small addends when values have very different
/// magnitudes.
///
/// In norm routines this is used by the default APIs such as
/// [`crate::algorithms::vector_norms::vector_norm_l1`],
/// [`crate::algorithms::vector_norms::vector_norm_l2`], and
/// [`crate::algorithms::vector_norms::vector_norm_l2_sq`].
///
/// # Example
///
/// ```
/// use num_valid::algorithms::accumulators::*;
///
/// let mut naive = NaiveSum::<f64>::new();
/// naive.push(1.0);
/// naive.push(2.0);
/// naive.push(3.0);
/// assert_eq!(naive.result(), 6.0);
/// ```
pub struct NaiveSum<T>(T);

impl<T> Accumulator for NaiveSum<T>
where
    T: Zero + Clone + for<'a> Add<&'a T, Output = T> + for<'a> Mul<&'a T, Output = T> + AddAssign,
{
    type Input = T;
    type Output = T;

    #[inline(always)]
    fn new() -> Self {
        NaiveSum(T::zero())
    }

    #[inline(always)]
    fn push(&mut self, r: T) {
        self.0 += r;
    }

    #[inline(always)]
    fn result(self) -> T {
        self.0
    }

    #[inline(always)]
    fn combine(&mut self, other: Self) {
        self.0 += other.0;
    }
}

impl<T> SumAccumulator for NaiveSum<T>
where
    T: Zero
        + Clone
        + for<'a> Add<&'a T, Output = T>
        + for<'a> Mul<&'a T, Output = T>
        + AddAssign
        + for<'a> MulAssign<&'a T>,
{
    #[inline(always)]
    fn rescale_by(&mut self, r: &T) {
        self.0 *= r;
    }
}
//---------------------------------------------------------

//---------------------------------------------------------
#[inline(always)]
fn neumaier_sum_and_compensation_native_real<T>(value: T, sum: &mut T, compensation: &mut T)
where
    T: sealed::NativeReal + AddAssign,
{
    let sum_before_compensation = *sum;
    *sum += value; // here we may lose low-order digits of value if sum is much larger than value or vice versa

    // NOTE: the parenthesis are necessary in order for this algorithm to work correctly.
    *compensation += if sum_before_compensation.abs() >= value.abs() {
        // If sum is bigger, low-order digits of value are lost.
        (sum_before_compensation - *sum) + value
    } else {
        // Else low-order digits of sum are lost.
        (value - *sum) + sum_before_compensation
    };
}

#[inline(always)]
pub(crate) fn neumaier_sum_and_compensation_real<RealType>(
    value: RealType,
    sum: &mut RealType,
    compensation: &mut RealType,
) where
    RealType: RealScalar,
{
    let sum_before_compensation = sum.clone();
    *sum += &value; // here we may lose low-order digits of value if sum is much larger than value or vice versa

    // NOTE: the parenthesis are necessary in order for this algorithm to work correctly.
    *compensation += if sum_before_compensation.clone().abs() >= value.clone().abs() {
        // If sum is bigger, low-order digits of value are lost.
        (sum_before_compensation - &*sum) + value
    } else {
        // Else low-order digits of sum are lost.
        (value - &*sum) + sum_before_compensation
    };
}

/// A trait for types that can participate in Neumaier compensated summation.
///
/// This trait provides the core operation for the Neumaier summation algorithm,
/// which is a numerically stable method for summing floating-point numbers with
/// reduced rounding errors compared to naive summation.
///
/// # Algorithm
///
/// The Neumaier algorithm is an improvement over Kahan summation that handles
/// the case when the new value is larger than the running sum. It maintains
/// a compensation term to capture the low-order bits that are lost during addition.
pub trait NeumaierAddable: Sized {
    /// Performs one step of the Neumaier compensated summation algorithm.
    ///
    /// This method adds `value` to the running `sum` while tracking lost precision
    /// in the `compensation` term. The final result should be computed as
    /// `sum + compensation` after all values have been added.
    ///
    /// # Arguments
    ///
    /// * `value` - The value to add to the sum.
    /// * `sum` - A mutable reference to the running sum.
    /// * `compensation` - A mutable reference to the compensation term that tracks lost precision.
    fn neumaier_compensated_sum(value: Self, sum: &mut Self, compensation: &mut Self);
}

impl<T> NeumaierAddable for T
where
    T: sealed::NativeReal + AddAssign,
{
    #[inline(always)]
    fn neumaier_compensated_sum(value: Self, sum: &mut Self, compensation: &mut Self) {
        neumaier_sum_and_compensation_native_real(value, sum, compensation)
    }
}

impl<T> NeumaierAddable for Complex<T>
where
    T: sealed::NativeReal + AddAssign,
{
    #[inline(always)]
    fn neumaier_compensated_sum(value: Self, sum: &mut Self, compensation: &mut Self) {
        neumaier_sum_and_compensation_native_real(value.re, &mut sum.re, &mut compensation.re);
        neumaier_sum_and_compensation_native_real(value.im, &mut sum.im, &mut compensation.im);
    }
}

#[cfg(feature = "rug")]
mod rug_impls {
    use super::*;

    #[inline(always)]
    fn neumaier_sum_and_compensation_rug_float(
        value: rug::Float,
        sum: &mut rug::Float,
        compensation: &mut rug::Float,
    ) {
        let sum_before_compensation = sum.clone();
        *sum += &value;

        // NOTE: the parenthesis are necessary in order for this algorithm to work correctly.
        *compensation += if sum_before_compensation.clone().abs() >= value.clone().abs() {
            // If sum is bigger, low-order digits of value are lost.
            (sum_before_compensation - &*sum) + value
        } else {
            // Else low-order digits of sum are lost.
            (value - &*sum) + sum_before_compensation
        };
    }

    impl NeumaierAddable for rug::Float {
        #[inline(always)]
        fn neumaier_compensated_sum(value: Self, sum: &mut Self, compensation: &mut Self) {
            neumaier_sum_and_compensation_rug_float(value, sum, compensation)
        }
    }

    impl NeumaierAddable for rug::Complex {
        #[inline(always)]
        fn neumaier_compensated_sum(value: Self, sum: &mut Self, compensation: &mut Self) {
            let (value_real, value_imag) = value.into_real_imag();

            neumaier_sum_and_compensation_rug_float(
                value_real,
                sum.mut_real(),
                compensation.mut_real(),
            );

            neumaier_sum_and_compensation_rug_float(
                value_imag,
                sum.mut_imag(),
                compensation.mut_imag(),
            );
        }
    }
}

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

//------------------------------------------------------------------------------------
/// Neumaier compensated accumulator for sums of floating-point (real or complex) numbers.
///
/// Standard sequential summation loses precision when values of very different
/// magnitudes are added: a small addend can be completely swallowed by a much
/// larger running total.  `NeumaierSum` corrects for this by maintaining a
/// *compensation* term that captures the low-order bits discarded at each step,
/// then folds them back into the result at the end.
///
/// `NeumaierSum<T>` implements [`SumAccumulator`],
/// so it can be used directly as the accumulation strategy in the
/// [`norms`](crate::algorithms::vector_norms) vector-norm functions (e.g.,
/// [`vector_norm_l2_neumaier`](crate::algorithms::vector_norms::vector_norm_l2_neumaier),
/// [`vector_norm_l1_neumaier`](crate::algorithms::vector_norms::vector_norm_l1_neumaier)).
///
/// # Example
///
/// ```
/// use num_valid::algorithms::accumulators::*;
///
/// let mut neumaier = NeumaierSum::<f64>::new();
/// neumaier.push(1.0);
/// neumaier.push(1e100);
/// neumaier.push(1.0);
/// neumaier.push(-1e100);
///
/// let sum = neumaier.result();
/// println!("Sum: {}", sum);
/// assert_eq!(sum, 2.0);
/// ```
///
/// # References
///
/// * [Neumaier Summation](https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements)
#[derive(Debug, Clone, Getters)]
pub struct NeumaierSum<ScalarType> {
    /// The sum before the compensation term is added.
    #[getset(get = "pub")]
    sum_before_compensation: ScalarType,

    /// The compensation term.
    /// This is the correction term that is added to the `sum_before_compensation` to correct for the loss of precision.
    #[getset(get = "pub")]
    compensation: ScalarType,
}

impl<ScalarType> Accumulator for NeumaierSum<ScalarType>
where
    ScalarType: Clone
        + Zero
        + for<'a> Add<&'a ScalarType, Output = ScalarType>
        + for<'a> AddAssign<&'a ScalarType>
        + for<'a> Mul<&'a ScalarType, Output = ScalarType>
        + for<'a> MulAssign<&'a ScalarType>
        + NeumaierAddable,
{
    type Input = ScalarType;
    type Output = ScalarType;

    #[inline(always)]
    fn new() -> Self {
        Self {
            sum_before_compensation: ScalarType::zero(),
            compensation: ScalarType::zero(),
        }
    }

    #[inline(always)]
    fn push(&mut self, value: Self::Input) {
        <ScalarType as NeumaierAddable>::neumaier_compensated_sum(
            value,
            &mut self.sum_before_compensation,
            &mut self.compensation,
        );
    }

    #[inline(always)]
    fn result(self) -> Self::Input {
        self.sum_before_compensation + self.compensation
    }

    #[inline(always)]
    fn combine(&mut self, other: Self) {
        self.push(other.sum_before_compensation);
        self.push(other.compensation);
    }
}

impl<ScalarType> SumAccumulator for NeumaierSum<ScalarType>
where
    ScalarType: Clone
        + Zero
        + for<'a> Add<&'a ScalarType, Output = ScalarType>
        + for<'a> AddAssign<&'a ScalarType>
        + for<'a> Mul<&'a ScalarType, Output = ScalarType>
        + for<'a> MulAssign<&'a ScalarType>
        + NeumaierAddable,
{
    #[inline(always)]
    fn rescale_by(&mut self, r: &Self::Input) {
        // Compute the current true sum before rescaling.
        self.sum_before_compensation += &self.compensation;

        // Rescale the sum and reset the compensation term to zero.
        self.sum_before_compensation *= r;
        self.compensation = ScalarType::zero();
    }
}
//------------------------------------------------------------------------------------

//------------------------------------------------------------------------------------
/// Accumulator that tracks the running maximum of real values.
///
/// Initialises to [`Constants::min_finite`](crate::Constants::min_finite) so that any finite input
/// immediately becomes the maximum. For non-empty inputs this is always the
/// true maximum; for empty inputs [`result`](Accumulator::result) returns
/// `min_finite`.
///
/// # Note on empty input
///
/// For norm computations prefer [`MaxAbsValueAccumulator`], which initialises
/// to zero and is therefore correct for empty vectors (norm = 0).
pub struct MaxAccumulator<RealType: RealScalar> {
    max_value: RealType,
}

impl<RealType: RealScalar> Accumulator for MaxAccumulator<RealType> {
    type Input = RealType;
    type Output = RealType;

    #[inline(always)]
    fn new() -> Self {
        Self {
            max_value: RealType::min_finite(),
        }
    }

    #[inline(always)]
    fn push(&mut self, value: Self::Input) {
        if value > self.max_value {
            self.max_value = value;
        }
    }

    #[inline(always)]
    fn result(self) -> Self::Output {
        self.max_value
    }

    #[inline(always)]
    fn combine(&mut self, other: Self) {
        if other.max_value > self.max_value {
            self.max_value = other.max_value;
        }
    }
}

/// Accumulator that tracks the running minimum of real values.
///
/// Initialises to [`Constants::max_finite`](crate::Constants::max_finite) so that any finite input
/// immediately becomes the minimum. For empty inputs
/// [`result`](Accumulator::result) returns `max_finite`.
pub struct MinAccumulator<RealType: RealScalar> {
    min_value: RealType,
}

impl<RealType: RealScalar> Accumulator for MinAccumulator<RealType> {
    type Input = RealType;
    type Output = RealType;

    #[inline(always)]
    fn new() -> Self {
        Self {
            min_value: RealType::max_finite(),
        }
    }

    #[inline(always)]
    fn push(&mut self, value: Self::Input) {
        if value < self.min_value {
            self.min_value = value;
        }
    }

    #[inline(always)]
    fn result(self) -> Self::Output {
        self.min_value
    }

    #[inline(always)]
    fn combine(&mut self, other: Self) {
        if other.min_value < self.min_value {
            self.min_value = other.min_value;
        }
    }
}
//------------------------------------------------------------------------------------

//------------------------------------------------------------------------------------
/// Accumulator that computes the maximum absolute value over a sequence of
/// scalars: `maxᵢ |xᵢ|`.
///
/// Unlike [`MaxAccumulator`], this type:
/// - accepts any [`FpScalar`] (real **or** complex) as input, applying `abs()`
///   internally, so the modulus is used for complex scalars.
/// - initialises to **zero** rather than `min_finite`, making it correct for
///   empty inputs (result = 0, matching the L∞ norm of an empty vector).
/// - returns a validated [`NonNegativeRealScalar`] from
///   [`result`](Accumulator::result), encoding the non-negativity invariant
///   at the type level.
///
/// This is the accumulator used internally by [`LinfNorm`](crate::algorithms::vector_norms::LinfNorm).
pub struct MaxAbsValueAccumulator<ScalarType: FpScalar> {
    max_abs: NonNegativeRealScalar<ScalarType::RealType>, // internal: RealType for comparison of absolute values
}

impl<ScalarType: FpScalar> Accumulator for MaxAbsValueAccumulator<ScalarType> {
    type Input = ScalarType;
    type Output = NonNegativeRealScalar<ScalarType::RealType>;

    fn new() -> Self {
        Self {
            max_abs: NonNegativeRealScalar::zero(),
        }
    }

    fn push(&mut self, value: ScalarType) {
        let abs_v = value.abs();
        if abs_v > self.max_abs {
            self.max_abs = abs_v;
        }
    }

    fn combine(&mut self, other: Self) {
        if other.max_abs > self.max_abs {
            self.max_abs = other.max_abs;
        }
    }

    fn result(self) -> NonNegativeRealScalar<ScalarType::RealType> {
        self.max_abs
    }
}
//------------------------------------------------------------------------------------

//---------------------------------------------------------
/// Accumulator for computing scaled sum-of-squares: `scale² × Σ(xᵢ/scale)²`.
///
/// This accumulator is used internally for computing vector norms (especially L2),
/// but is general enough for any algorithm that needs to sum squared non-negative values
/// while maintaining numerical stability through dynamic scaling.
///
/// The key invariant is that `scale` tracks the running maximum absolute value,
/// preventing overflow and underflow. When a larger value is encountered, previous
/// accumulated terms are rescaled by `(old_scale / new_scale)²`.
///
/// # Type parameter `SumSqAcc`
///
/// The accumulation strategy for the scaled sum-of-squares term:
/// - [`NaiveSum`]: plain `+=`, fastest
/// - [`NeumaierSum`]: compensated summation for better accuracy
///
/// # Input contract
///
/// The [`Input`](Accumulator::Input) type is [`NonNegativeRealScalar<ScalarType::RealType>`],
/// which explicitly encodes that all pushed values must be non-negative.
/// Callers must validate and wrap values with [`NonNegativeRealScalar`] before pushing.
///
/// # Example
///
/// ```ignore
/// use num_valid::algorithms::accumulators::*;
/// use num_valid::scalars::NonNegativeRealScalar;
///
/// let mut acc = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
/// acc.push(NonNegativeRealScalar::new(3.0));  // validated non-negative
/// acc.push(NonNegativeRealScalar::new(4.0));
/// let result = acc.result();  // 3² + 4² = 25
/// ```
pub struct ScaledSumSqAccumulator<SumSqAcc>
where
    SumSqAcc: SumAccumulator<Input: RealScalar>,
{
    /// - `scale = max |vᵢ|`  (0 if all components are zero)
    scale: NonNegativeRealScalar<SumSqAcc::Input>,
    /// - `sumsq = Σ (|vᵢ|/scale)²`  accumulated via strategy `SumSqAcc`
    sumsq_accumulator: SumSqAcc,
    zero: NonNegativeRealScalar<SumSqAcc::Input>,
}

impl<SumSqAcc> ScaledSumSqAccumulator<SumSqAcc>
where
    SumSqAcc: SumAccumulator<Input: RealScalar>,
{
    fn add_value_and_update_scale(&mut self, value: NonNegativeRealScalar<SumSqAcc::Input>) {
        if value > self.zero {
            let r_sq = if self.scale < value {
                if self.scale > self.zero {
                    // Update scale and rescale the sum-of-squares accumulator
                    //let r = self.scale.as_raw_scalar_ref().clone() / value.as_raw_scalar_ref();
                    //self.scale = value;
                    //
                    // old_scale = self.scale;
                    // r = value / old_scale
                    // self.scale = value;
                    let r = std::mem::replace(&mut self.scale, value).into_inner()
                        / self.scale.as_ref();

                    self.sumsq_accumulator.rescale_by(&r.pow(2));
                } else {
                    self.scale = value;
                }
                SumSqAcc::Input::one()
            } else {
                (value.into_inner() / self.scale.as_ref()).pow(2)
            };
            self.sumsq_accumulator.push(r_sq);
        }
    }

    pub fn result_sqrt(self) -> NonNegativeRealScalar<SumSqAcc::Input> {
        if self.scale == self.zero {
            self.zero
        } else {
            NonNegativeRealScalar::new(self.sumsq_accumulator.result()).sqrt() * self.scale
        }
    }

    /// Build a new accumulator by sequentially pushing non-negative values.
    ///
    /// This is the basic constructor for building from validated non-negative inputs.
    /// For building from arbitrary scalars (which need abs()), use
    /// [`new_sequential_from_scalars`](Self::new_sequential_from_scalars).
    pub fn new_sequential<I>(values: I) -> Self
    where
        I: IntoIterator<Item = NonNegativeRealScalar<SumSqAcc::Input>>,
    {
        let mut acc = Self::new();
        values.into_iter().for_each(|v| acc.push(v));
        acc
    }

    /// Build a new accumulator by sequentially processing scalars, automatically taking abs().
    ///
    /// This is a convenience method for L2 norm computation, where scalar values
    /// (real or complex) are automatically converted to their absolute value
    /// (modulus for complex scalars) before accumulation.
    ///
    /// # Example
    ///
    /// ```ignore
    /// use num_valid::algorithms::accumulators::*;
    ///
    /// let values = [3.0_f64, -4.0, 5.0];
    /// let mut acc = ScaledSumSqAccumulator::<NaiveSum<f64>>::new_sequential_from_scalars(values);
    /// // Accumulates |3| = 3, |-4| = 4, |5| = 5
    /// // Result: √(9 + 16 + 25) = √50
    /// ```
    pub fn new_sequential_from_scalars<I>(values: I) -> Self
    where
        I: IntoIterator<Item: FpScalar<RealType = SumSqAcc::Input>>,
    {
        let mut acc = Self::new();
        values.into_iter().for_each(|v| {
            acc.push(v.abs());
        });
        acc
    }

    /// Build a new accumulator in parallel, automatically taking abs() of each scalar.
    ///
    /// Distributes work across threads using rayon. Each thread accumulates its chunk
    /// independently, then results are merged via [`Accumulator::combine`].
    pub fn new_parallel_from_scalars<I>(values: I) -> Self
    where
        Self: Send,
        I: IntoParallelIterator<Item: FpScalar<RealType = SumSqAcc::Input>>,
    {
        values
            .into_par_iter()
            .fold(Self::new, |mut acc, item| {
                acc.push(item.abs());
                acc
            })
            .reduce(Self::new, |mut a, b| {
                a.combine(b);
                a
            })
    }
}

impl<SumSqAcc> Accumulator for ScaledSumSqAccumulator<SumSqAcc>
where
    SumSqAcc: SumAccumulator<Input: RealScalar>,
{
    /// Input must be already validated as non-negative
    type Input = NonNegativeRealScalar<SumSqAcc::Input>;
    type Output = NonNegativeRealScalar<SumSqAcc::Input>;

    fn new() -> Self {
        Self {
            scale: NonNegativeRealScalar::<SumSqAcc::Input>::zero(),
            sumsq_accumulator: SumSqAcc::new(),
            zero: NonNegativeRealScalar::<SumSqAcc::Input>::zero(),
        }
    }

    /// Adds a non-negative value to the accumulator.
    /// The input must already be validated as non-negative via [`NonNegativeRealScalar`].
    fn push(&mut self, value: Self::Input) {
        self.add_value_and_update_scale(value);
    }

    /// Merges another accumulator, rescaling as needed when scales differ.
    /// Handles all permutations of scale relationships (>, <, ==) and empty accumulators.
    fn combine(&mut self, mut other: Self) {
        if other.scale > self.scale {
            if self.scale > self.zero {
                // old_scale = self.scale
                // self.scale = other.scale
                // r = old_scale / self.scale
                let r = std::mem::replace(&mut self.scale, other.scale).into_inner()
                    / self.scale.as_ref();
                //let r = self.scale.clone().into_inner() / other.scale.as_ref();
                //self.scale = other.scale;
                self.sumsq_accumulator.rescale_by(&r.pow(2));
            } else {
                self.scale = other.scale;
            }
            self.sumsq_accumulator.combine(other.sumsq_accumulator);
        } else if other.scale > self.zero {
            if other.scale < self.scale {
                let r = other.scale.into_inner() / self.scale.as_ref();
                other.sumsq_accumulator.rescale_by(&r.pow(2));
            }
            self.sumsq_accumulator.combine(other.sumsq_accumulator);
        }
        // when both scale==0: zero vector, nothing to do
    }

    /// Returns the squared result: `scale² × sumsq`.
    fn result(self) -> Self::Output {
        if self.scale == self.zero {
            self.zero
        } else {
            NonNegativeRealScalar::new(
                self.sumsq_accumulator.result() * self.scale.into_inner().pow(2),
            )
        }
    }
}
//---------------------------------------------------------

//------------------------------------------------------------------------------------
#[cfg(test)]
mod tests {
    use super::*;

    mod neumaier_sum {
        use super::*;

        mod native {
            use super::*;

            mod real {
                use super::*;

                #[test]
                fn new() {
                    let neumaier = NeumaierSum::<f64>::new();
                    assert_eq!(neumaier.sum_before_compensation, 0.0);
                    assert_eq!(neumaier.compensation, 0.0);
                }

                #[test]
                fn add() {
                    let mut neumaier = NeumaierSum::<f64>::new();
                    neumaier.push(1.0);
                    neumaier.push(1e-16);
                    neumaier.push(-1.0);
                    assert_eq!(neumaier.sum_before_compensation, 0.0);
                    assert_eq!(neumaier.compensation, 1e-16);
                }

                #[test]
                fn sum() {
                    let mut neumaier = NeumaierSum::<f64>::new();
                    neumaier.push(1.0);
                    neumaier.push(1e-16);
                    neumaier.push(-1.0);
                    assert_eq!(neumaier.sum_before_compensation, 0.0);
                    assert_eq!(neumaier.compensation, 1e-16);
                    let sum = neumaier.result();
                    assert_eq!(sum, 1e-16);
                    println!("compensated sum = {}", sum);
                }

                #[test]
                fn sum_big_values() {
                    let values = vec![1.0, 1e100, 1.0, -1e100];
                    let sum = values.iter().sum::<f64>();
                    assert_eq!(sum, 0.0);

                    let neumaier = NeumaierSum::<f64>::new_sequential(values);
                    let sum = neumaier.result();
                    assert_eq!(sum, 2.0);
                    println!("compensated sum = {}", sum);
                }

                #[test]
                fn sum_small_values() {
                    let values = [1.0, 1e-100, -1.0];
                    let sum = values.iter().sum::<f64>();
                    assert_eq!(sum, 0.0);

                    let neumaier = NeumaierSum::<f64>::new_sequential(values);
                    let sum = neumaier.result();
                    assert_eq!(sum, 1e-100);
                    println!("compensated sum = {}", sum);
                }

                #[test]
                fn combine_partial_sums() {
                    let mut left = NaiveSum::<f64>::new();
                    left.push(1.0);
                    left.push(2.0);

                    let mut right = NaiveSum::<f64>::new();
                    right.push(3.0);
                    right.push(4.0);

                    left.combine(right);
                    assert_eq!(left.result(), 10.0);
                }

                #[test]
                fn combine_partial_neumaier_sums() {
                    let mut left = NeumaierSum::<f64>::new();
                    left.push(1.0);
                    left.push(1e-16);

                    let mut right = NeumaierSum::<f64>::new();
                    right.push(2.0);
                    right.push(1e-16);

                    left.combine(right);
                    let sum = left.result();
                    assert!((sum - (3.0 + 2e-16)).abs() < 1e-15);
                }

                #[test]
                fn combine_matches_new_sequential_for_chunks() {
                    let values = [1.0_f64, 1e100, 1.0, -1e100, 1e-16, 2.0];

                    let mut chunk_1 = NaiveSum::<f64>::new();
                    chunk_1.push(values[0]);
                    chunk_1.push(values[1]);

                    let mut chunk_2 = NaiveSum::<f64>::new();
                    chunk_2.push(values[2]);
                    chunk_2.push(values[3]);

                    let mut chunk_3 = NaiveSum::<f64>::new();
                    chunk_3.push(values[4]);
                    chunk_3.push(values[5]);

                    chunk_1.combine(chunk_2);
                    chunk_1.combine(chunk_3);

                    let sequential = NaiveSum::<f64>::new_sequential(values);
                    assert_eq!(chunk_1.result(), sequential.result());
                }

                #[test]
                fn combine_is_order_dependent_for_naive_sum() {
                    let mut left_grouped = NaiveSum::<f64>::new();
                    left_grouped.push(1e16);

                    let mut middle = NaiveSum::<f64>::new();
                    middle.push(1.0);

                    let mut right_last = NaiveSum::<f64>::new();
                    right_last.push(-1e16);

                    left_grouped.combine(middle);
                    left_grouped.combine(right_last);
                    let left_result = left_grouped.result();

                    let mut right_grouped = NaiveSum::<f64>::new();
                    right_grouped.push(1e16);

                    let mut right_last = NaiveSum::<f64>::new();
                    right_last.push(-1e16);

                    let mut middle = NaiveSum::<f64>::new();
                    middle.push(1.0);

                    right_grouped.combine(right_last);
                    right_grouped.combine(middle);
                    let right_result = right_grouped.result();

                    assert_eq!(left_result, 0.0);
                    assert_eq!(right_result, 1.0);
                    assert_ne!(left_result, right_result);
                }

                #[test]
                fn neumaier_combine_matches_new_sequential_for_chunks() {
                    let values = [1.0_f64, 1e100, 1.0, -1e100, 1e-16, 2.0];

                    let mut chunk_1 = NeumaierSum::<f64>::new();
                    chunk_1.push(values[0]);
                    chunk_1.push(values[1]);

                    let mut chunk_2 = NeumaierSum::<f64>::new();
                    chunk_2.push(values[2]);
                    chunk_2.push(values[3]);

                    let mut chunk_3 = NeumaierSum::<f64>::new();
                    chunk_3.push(values[4]);
                    chunk_3.push(values[5]);

                    chunk_1.combine(chunk_2);
                    chunk_1.combine(chunk_3);

                    let sequential = NeumaierSum::<f64>::new_sequential(values);
                    assert_eq!(chunk_1.result(), sequential.result());
                }

                #[test]
                fn neumaier_combine_is_stable_across_orders() {
                    let mut left_order = NeumaierSum::<f64>::new();
                    left_order.push(1e16);

                    let mut middle = NeumaierSum::<f64>::new();
                    middle.push(1.0);

                    let mut right_last = NeumaierSum::<f64>::new();
                    right_last.push(-1e16);

                    left_order.combine(middle);
                    left_order.combine(right_last);
                    let left_result = left_order.result();

                    let mut right_order = NeumaierSum::<f64>::new();
                    right_order.push(1e16);

                    let mut right_last = NeumaierSum::<f64>::new();
                    right_last.push(-1e16);

                    let mut middle = NeumaierSum::<f64>::new();
                    middle.push(1.0);

                    right_order.combine(right_last);
                    right_order.combine(middle);
                    let right_result = right_order.result();

                    assert_eq!(left_result, 1.0);
                    assert_eq!(right_result, 1.0);
                    assert_eq!(left_result, right_result);
                }
            }

            mod complex {
                use super::*;
                use num::Complex;

                #[test]
                fn new() {
                    let neumaier = NeumaierSum::<Complex<f64>>::new();

                    let zero = Complex::new(0.0, 0.0);
                    assert_eq!(&neumaier.sum_before_compensation, &zero);
                    assert_eq!(&neumaier.compensation, &zero);
                }

                #[test]
                fn add() {
                    let mut neumaier = NeumaierSum::<Complex<f64>>::new();

                    let zero = Complex::new(0.0, 0.0);
                    let v = Complex::new(1e-16, 2e-16);

                    neumaier.push(Complex::new(1.0, 2.0));
                    neumaier.push(v);
                    neumaier.push(Complex::new(-1.0, -2.0));

                    assert_eq!(neumaier.sum_before_compensation, zero);
                    assert_eq!(neumaier.compensation, v);
                }

                #[test]
                fn sum() {
                    let mut neumaier = NeumaierSum::<Complex<f64>>::new();

                    let zero = Complex::new(0.0, 0.0);
                    let v = Complex::new(1e-16, 2e-16);

                    neumaier.push(Complex::new(1.0, 2.0));
                    neumaier.push(v);
                    neumaier.push(Complex::new(-1.0, -2.0));
                    assert_eq!(neumaier.sum_before_compensation, zero);
                    assert_eq!(neumaier.compensation, v);
                    let sum = neumaier.result();
                    assert_eq!(sum, v);
                    println!("compensated sum = {}", sum);
                }

                #[test]
                fn sum_big_values() {
                    let values = vec![
                        Complex::new(1.0, 2.0),
                        Complex::new(1e100, 2e100),
                        Complex::new(1.0, 2.0),
                        Complex::new(-1e100, -2e100),
                    ];
                    let sum = values.iter().sum::<Complex<f64>>();
                    assert_eq!(sum, Complex::new(0.0, 0.0));

                    let neumaier = NeumaierSum::<Complex<f64>>::new_sequential(values);
                    let sum = neumaier.result();
                    assert_eq!(sum, Complex::new(2.0, 4.0));
                    println!("compensated sum = {}", sum);
                }

                #[test]
                fn sum_small_values() {
                    let v = Complex::new(1e-100, 2e-100);

                    let values = [Complex::new(1.0, 2.0), v, Complex::new(-1.0, -2.0)];
                    let sum = values.iter().sum::<Complex<f64>>();
                    assert_eq!(sum, Complex::new(0.0, 0.0));

                    let neumaier = NeumaierSum::<Complex<f64>>::new_sequential(values);
                    let sum = neumaier.result();
                    assert_eq!(sum, v);
                    println!("compensated sum = {}", sum);
                }

                #[test]
                fn combine_partial_sums() {
                    let mut left = NaiveSum::<Complex<f64>>::new();
                    left.push(Complex::new(1.0, 2.0));

                    let mut right = NaiveSum::<Complex<f64>>::new();
                    right.push(Complex::new(3.0, 4.0));

                    left.combine(right);
                    assert_eq!(left.result(), Complex::new(4.0, 6.0));
                }

                #[test]
                fn combine_partial_neumaier_sums() {
                    let mut left = NeumaierSum::<Complex<f64>>::new();
                    left.push(Complex::new(1.0, 2.0));
                    left.push(Complex::new(1e-16, 2e-16));

                    let mut right = NeumaierSum::<Complex<f64>>::new();
                    right.push(Complex::new(3.0, 4.0));
                    right.push(Complex::new(1e-16, 2e-16));

                    left.combine(right);
                    assert_eq!(
                        left.result(),
                        Complex::new(4.0, 6.0) + Complex::new(2e-16, 4e-16)
                    );
                }

                #[test]
                fn combine_matches_new_sequential_for_chunks() {
                    let values = [
                        Complex::new(1.0, 2.0),
                        Complex::new(1e100, 2e100),
                        Complex::new(1.0, 2.0),
                        Complex::new(-1e100, -2e100),
                        Complex::new(1e-16, 2e-16),
                        Complex::new(2.0, 3.0),
                    ];

                    let mut chunk_1 = NaiveSum::<Complex<f64>>::new();
                    chunk_1.push(values[0]);
                    chunk_1.push(values[1]);

                    let mut chunk_2 = NaiveSum::<Complex<f64>>::new();
                    chunk_2.push(values[2]);
                    chunk_2.push(values[3]);

                    let mut chunk_3 = NaiveSum::<Complex<f64>>::new();
                    chunk_3.push(values[4]);
                    chunk_3.push(values[5]);

                    chunk_1.combine(chunk_2);
                    chunk_1.combine(chunk_3);

                    let sequential = NaiveSum::<Complex<f64>>::new_sequential(values);
                    assert_eq!(chunk_1.result(), sequential.result());
                }

                #[test]
                fn neumaier_combine_matches_new_sequential_for_chunks() {
                    let values = [
                        Complex::new(1.0, 2.0),
                        Complex::new(1e100, 2e100),
                        Complex::new(1.0, 2.0),
                        Complex::new(-1e100, -2e100),
                        Complex::new(1e-16, 2e-16),
                        Complex::new(2.0, 3.0),
                    ];

                    let mut chunk_1 = NeumaierSum::<Complex<f64>>::new();
                    chunk_1.push(values[0]);
                    chunk_1.push(values[1]);

                    let mut chunk_2 = NeumaierSum::<Complex<f64>>::new();
                    chunk_2.push(values[2]);
                    chunk_2.push(values[3]);

                    let mut chunk_3 = NeumaierSum::<Complex<f64>>::new();
                    chunk_3.push(values[4]);
                    chunk_3.push(values[5]);

                    chunk_1.combine(chunk_2);
                    chunk_1.combine(chunk_3);

                    let sequential = NeumaierSum::<Complex<f64>>::new_sequential(values);
                    assert_eq!(chunk_1.result(), sequential.result());
                }
            }
        }

        #[cfg(feature = "rug")]
        mod rug100 {
            use super::*;
            use crate::{ComplexRugStrictFinite, RealRugStrictFinite};
            use try_create::TryNew;

            const PRECISION: u32 = 100;

            mod real {
                use super::*;
                use crate::RealScalar;

                #[test]
                fn new() {
                    let neumaier = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
                    assert!(neumaier.sum_before_compensation.is_zero());
                    assert!(neumaier.compensation.is_zero());
                }

                #[test]
                fn add() {
                    let mut neumaier = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();

                    let v = RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
                        PRECISION,
                        rug::Float::parse("1e-100").unwrap(),
                    ))
                    .expect("valid test value");

                    neumaier.push(RealRugStrictFinite::<PRECISION>::try_from_f64(1.0).unwrap());
                    neumaier.push(v.clone());
                    neumaier.push(RealRugStrictFinite::<PRECISION>::try_from_f64(-1.0).unwrap());

                    assert_eq!(
                        neumaier.sum_before_compensation,
                        RealRugStrictFinite::<PRECISION>::try_from_f64(0.0).unwrap()
                    );
                    assert_eq!(&neumaier.compensation, &v);
                }

                #[test]
                fn sum() {
                    let mut neumaier = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();

                    let v = RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
                        PRECISION,
                        rug::Float::parse("1e-100").unwrap(),
                    ))
                    .expect("valid test value");

                    neumaier.push(RealRugStrictFinite::<PRECISION>::try_from_f64(1.0).unwrap());
                    neumaier.push(v.clone());
                    neumaier.push(RealRugStrictFinite::<PRECISION>::try_from_f64(-1.0).unwrap());

                    assert!(neumaier.sum_before_compensation.is_zero());
                    assert_eq!(&neumaier.compensation, &v);
                    let sum = neumaier.result();
                    assert_eq!(sum, v);
                    println!("compensated sum = {}", sum);
                }

                #[test]
                fn sum_big_values() {
                    let values = ["1.0", "1e100", "1.0", "-1e100"]
                        .iter()
                        .map(|v| {
                            RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
                                PRECISION,
                                rug::Float::parse(v).unwrap(),
                            ))
                            .expect("valid test value")
                        })
                        .collect::<Vec<_>>();

                    let sum = values
                        .iter()
                        .fold(RealRugStrictFinite::<PRECISION>::zero(), |acc, x| acc + x);
                    assert!(sum.is_zero());

                    let neumaier =
                        NeumaierSum::<RealRugStrictFinite<PRECISION>>::new_sequential(values);
                    let sum = neumaier.result();
                    assert_eq!(
                        sum,
                        RealRugStrictFinite::<PRECISION>::try_from_f64(2.0).unwrap()
                    );
                    println!("compensated sum = {}", sum);
                }

                #[test]
                fn sum_small_values() {
                    let values = ["1.0", "1e-100", "-1.0"]
                        .iter()
                        .map(|v| {
                            RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
                                PRECISION,
                                rug::Float::parse(v).unwrap(),
                            ))
                            .expect("valid test value")
                        })
                        .collect::<Vec<_>>();

                    let sum = values
                        .iter()
                        .fold(RealRugStrictFinite::<PRECISION>::zero(), |acc, x| acc + x);
                    assert_eq!(sum, RealRugStrictFinite::<PRECISION>::zero());

                    let neumaier =
                        NeumaierSum::<RealRugStrictFinite<PRECISION>>::new_sequential(values);
                    let sum = neumaier.result();
                    assert_eq!(
                        sum,
                        RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
                            PRECISION,
                            rug::Float::parse("1e-100").unwrap(),
                        ))
                        .expect("valid test value")
                    );
                    println!("compensated sum = {}", sum);
                }

                #[test]
                fn combine_matches_new_sequential_for_chunks() {
                    let values = ["1.0", "1e100", "1.0", "-1e100", "1e-100", "2.0"]
                        .iter()
                        .map(|v| {
                            RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
                                PRECISION,
                                rug::Float::parse(v).unwrap(),
                            ))
                            .expect("valid test value")
                        })
                        .collect::<Vec<_>>();

                    let mut chunk_1 = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
                    chunk_1.push(values[0].clone());
                    chunk_1.push(values[1].clone());

                    let mut chunk_2 = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
                    chunk_2.push(values[2].clone());
                    chunk_2.push(values[3].clone());

                    let mut chunk_3 = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
                    chunk_3.push(values[4].clone());
                    chunk_3.push(values[5].clone());

                    chunk_1.combine(chunk_2);
                    chunk_1.combine(chunk_3);

                    let sequential =
                        NeumaierSum::<RealRugStrictFinite<PRECISION>>::new_sequential(values);
                    assert_eq!(chunk_1.result(), sequential.result());
                }

                #[test]
                fn neumaier_combine_is_stable_across_orders() {
                    let mut left_order = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
                    left_order.push(RealRugStrictFinite::<PRECISION>::try_from_f64(1e16).unwrap());

                    let mut middle = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
                    middle.push(RealRugStrictFinite::<PRECISION>::try_from_f64(1.0).unwrap());

                    let mut right_last = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
                    right_last.push(RealRugStrictFinite::<PRECISION>::try_from_f64(-1e16).unwrap());

                    left_order.combine(middle);
                    left_order.combine(right_last);
                    let left_result = left_order.result();

                    let mut right_order = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
                    right_order.push(RealRugStrictFinite::<PRECISION>::try_from_f64(1e16).unwrap());

                    let mut right_last = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
                    right_last.push(RealRugStrictFinite::<PRECISION>::try_from_f64(-1e16).unwrap());

                    let mut middle = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
                    middle.push(RealRugStrictFinite::<PRECISION>::try_from_f64(1.0).unwrap());

                    right_order.combine(right_last);
                    right_order.combine(middle);
                    let right_result = right_order.result();

                    assert_eq!(
                        left_result,
                        RealRugStrictFinite::<PRECISION>::try_from_f64(1.0).unwrap()
                    );
                    assert_eq!(
                        right_result,
                        RealRugStrictFinite::<PRECISION>::try_from_f64(1.0).unwrap()
                    );
                    assert_eq!(left_result, right_result);
                }
            }

            mod complex {
                use super::*;

                #[test]
                fn new() {
                    let neumaier = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new();
                    assert_eq!(
                        neumaier.sum_before_compensation,
                        ComplexRugStrictFinite::<PRECISION>::zero()
                    );
                    assert_eq!(
                        neumaier.compensation,
                        ComplexRugStrictFinite::<PRECISION>::zero()
                    );
                }

                #[test]
                fn add() {
                    let mut neumaier = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new();

                    let v = ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
                        PRECISION,
                        rug::Complex::parse("(1e-100,2e-100)").unwrap(),
                    ))
                    .expect("valid test value");

                    neumaier.push(
                        ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(1.0, 2.0))
                            .unwrap(),
                    );
                    neumaier.push(v.clone());
                    neumaier.push(
                        ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(-1.0, -2.0))
                            .unwrap(),
                    );

                    assert_eq!(
                        neumaier.sum_before_compensation,
                        ComplexRugStrictFinite::<PRECISION>::zero()
                    );
                    assert_eq!(&neumaier.compensation, &v);
                }

                #[test]
                fn sum() {
                    let mut neumaier = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new();

                    let v = ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
                        PRECISION,
                        rug::Complex::parse("(1e-100,2e-100)").unwrap(),
                    ))
                    .expect("valid test value");

                    neumaier.push(
                        ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(1.0, 2.0))
                            .unwrap(),
                    );
                    neumaier.push(v.clone());
                    neumaier.push(
                        ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(-1.0, -2.0))
                            .unwrap(),
                    );

                    assert_eq!(
                        neumaier.sum_before_compensation,
                        ComplexRugStrictFinite::<PRECISION>::zero()
                    );
                    assert_eq!(&neumaier.compensation, &v);
                    let sum = neumaier.result();
                    assert_eq!(sum, v);
                    println!("compensated sum = {}", sum);
                }

                #[test]
                fn sum_big_values() {
                    let values = ["(1.0,2.0)", "(1e100,2e100)", "(1.0,2.0)", "(-1e100,-2e100)"]
                        .iter()
                        .map(|v| {
                            ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
                                PRECISION,
                                rug::Complex::parse(v).unwrap(),
                            ))
                            .expect("valid test value")
                        })
                        .collect::<Vec<_>>();

                    let zero = ComplexRugStrictFinite::<PRECISION>::zero();
                    let sum = values.iter().fold(zero.clone(), |acc, x| acc + x);
                    assert_eq!(sum, zero);

                    let neumaier =
                        NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new_sequential(values);
                    let sum = neumaier.result();
                    assert_eq!(
                        sum,
                        ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(2.0, 4.0))
                            .unwrap()
                    );
                    println!("compensated sum = {}", sum);
                }

                #[test]
                fn sum_small_values() {
                    let values = ["(1.0,2.0)", "(1e-100,2e-100)", "(-1.0,-2.0)"]
                        .iter()
                        .map(|v| {
                            ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
                                PRECISION,
                                rug::Complex::parse(v).unwrap(),
                            ))
                            .expect("valid test value")
                        })
                        .collect::<Vec<_>>();

                    let zero = ComplexRugStrictFinite::<PRECISION>::zero();
                    let sum = values.iter().fold(zero.clone(), |acc, x| acc + x);
                    assert_eq!(sum, zero);

                    let neumaier =
                        NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new_sequential(values);
                    let sum = neumaier.result();
                    assert_eq!(
                        sum,
                        ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
                            PRECISION,
                            rug::Complex::parse("(1e-100,2e-100)").unwrap(),
                        ))
                        .expect("valid test value")
                    );
                    println!("compensated sum = {}", sum);
                }

                #[test]
                fn combine_matches_new_sequential_for_chunks() {
                    let values = [
                        "(1.0,2.0)",
                        "(1e100,2e100)",
                        "(1.0,2.0)",
                        "(-1e100,-2e100)",
                        "(1e-100,2e-100)",
                        "(2.0,3.0)",
                    ]
                    .iter()
                    .map(|v| {
                        ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
                            PRECISION,
                            rug::Complex::parse(v).unwrap(),
                        ))
                        .expect("valid test value")
                    })
                    .collect::<Vec<_>>();

                    let mut chunk_1 = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new();
                    chunk_1.push(values[0].clone());
                    chunk_1.push(values[1].clone());

                    let mut chunk_2 = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new();
                    chunk_2.push(values[2].clone());
                    chunk_2.push(values[3].clone());

                    let mut chunk_3 = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new();
                    chunk_3.push(values[4].clone());
                    chunk_3.push(values[5].clone());

                    chunk_1.combine(chunk_2);
                    chunk_1.combine(chunk_3);

                    let sequential =
                        NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new_sequential(values);
                    assert_eq!(chunk_1.result(), sequential.result());
                }
            }
        }
    }

    mod scaled_sum_sq_accumulator_tests {
        use super::*;
        use crate::ScalarCore;

        const EPS: f64 = f64::EPSILON;

        /// Test single value: √(1²) = 1
        #[test]
        fn single_value() {
            let mut acc = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            acc.push(NonNegativeRealScalar::new(3.0));
            let result = acc.result_sqrt();
            assert!((*result.as_raw_scalar_ref() - 3.0).abs() < EPS);
        }

        /// Test 3-4-5 right triangle: √(3² + 4²) = 5
        #[test]
        fn three_four_five_sqrt() {
            let mut acc = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            acc.push(NonNegativeRealScalar::new(3.0));
            acc.push(NonNegativeRealScalar::new(4.0));
            let result = acc.result_sqrt();
            assert!((*result.as_raw_scalar_ref() - 5.0).abs() < EPS);
        }

        /// Test 3-4-5 squared: 3² + 4² = 25
        #[test]
        fn three_four_five_squared() {
            let mut acc = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            acc.push(NonNegativeRealScalar::new(3.0));
            acc.push(NonNegativeRealScalar::new(4.0));
            let result = acc.result();
            assert!((*result.as_raw_scalar_ref() - 25.0).abs() < EPS);
        }

        /// Test zero vector
        #[test]
        fn zero_vector_sqrt() {
            let acc = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            let result = acc.result_sqrt();
            assert_eq!(*result.as_raw_scalar_ref(), 0.0);
        }

        /// Test zero vector squared
        #[test]
        fn zero_vector_squared() {
            let acc = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            let result = acc.result();
            assert_eq!(*result.as_raw_scalar_ref(), 0.0);
        }

        /// Test negative values (abs should handle): √((-3)² + 4²) = 5
        #[test]
        fn negative_values() {
            let mut acc = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            acc.push(NonNegativeRealScalar::new(3.0));
            acc.push(NonNegativeRealScalar::new(4.0));
            let result = acc.result_sqrt();
            assert!((*result.as_raw_scalar_ref() - 5.0).abs() < EPS);
        }

        /// Test combine when other.scale > self.scale (rescale self)
        #[test]
        fn combine_other_scale_larger() {
            // acc1: [1, 1] → scale=1, sumsq=2
            let mut acc1 = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            acc1.push(NonNegativeRealScalar::new(1.0));
            acc1.push(NonNegativeRealScalar::new(1.0));

            // acc2: [100, 100] → scale=100, sumsq=2
            let mut acc2 = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            acc2.push(NonNegativeRealScalar::new(100.0));
            acc2.push(NonNegativeRealScalar::new(100.0));

            // Combined: scale=100, sumsq = (1/100)² * 2 + 2 = 0.0002 + 2 = 2.0002
            // result_sqrt: 100 * √(2.0002) ≈ 141.42
            acc1.combine(acc2);
            let result = acc1.result_sqrt();
            let expected = (1.0_f64 + 1.0 + 100.0 * 100.0 + 100.0 * 100.0).sqrt();
            assert!((*result.as_raw_scalar_ref() - expected).abs() < EPS * expected);
        }

        /// Test combine when self.scale > other.scale (rescale other)
        #[test]
        fn combine_self_scale_larger() {
            // acc1: [100, 100] → scale=100, sumsq=2
            let mut acc1 = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            acc1.push(NonNegativeRealScalar::new(100.0));
            acc1.push(NonNegativeRealScalar::new(100.0));

            // acc2: [1, 1] → scale=1, sumsq=2
            let mut acc2 = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            acc2.push(NonNegativeRealScalar::new(1.0));
            acc2.push(NonNegativeRealScalar::new(1.0));

            // Combined: scale=100, sumsq = 2 + (1/100)² * 2 = 2 + 0.0002
            acc1.combine(acc2);
            let result = acc1.result_sqrt();
            let expected = (100.0_f64 * 100.0 + 100.0 * 100.0 + 1.0 + 1.0).sqrt();
            assert!((*result.as_raw_scalar_ref() - expected).abs() < EPS * expected);
        }

        /// Test combine when both are zero (edge case)
        #[test]
        fn combine_both_zero() {
            let mut acc1 = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            let acc2 = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            acc1.combine(acc2);
            let result = acc1.result_sqrt();
            assert_eq!(*result.as_raw_scalar_ref(), 0.0);
        }

        /// Test combine: one zero, one non-zero
        #[test]
        fn combine_one_zero() {
            // acc1: zero
            let mut acc1 = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();

            // acc2: [3, 4]
            let mut acc2 = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            acc2.push(NonNegativeRealScalar::new(3.0));
            acc2.push(NonNegativeRealScalar::new(4.0));

            acc1.combine(acc2);
            let result = acc1.result_sqrt();
            assert!((*result.as_raw_scalar_ref() - 5.0).abs() < EPS);
        }

        /// Test combine with Neumaier (higher accuracy)
        #[test]
        fn combine_neumaier() {
            // acc1: small values
            let mut acc1 = ScaledSumSqAccumulator::<NeumaierSum<f64>>::new();
            for _ in 0..100 {
                acc1.push(NonNegativeRealScalar::new(0.01));
            }

            // acc2: single large value
            let mut acc2 = ScaledSumSqAccumulator::<NeumaierSum<f64>>::new();
            acc2.push(NonNegativeRealScalar::new(10.0));

            acc1.combine(acc2);
            let result = acc1.result_sqrt();

            // Expected: √(100*0.01² + 10²) = √(0.01 + 100) = √100.01
            let expected = (0.01_f64 + 100.0).sqrt();
            assert!((*result.as_raw_scalar_ref() - expected).abs() < EPS * expected);
        }

        /// Test result() vs result_sqrt(): result_sqrt()² should equal result()
        /// Because result() = scale² × sumsq and result_sqrt() = scale × √sumsq
        /// So (result_sqrt())² = scale² × sumsq = result()
        #[test]
        fn result_vs_result_sqrt() {
            let mut acc1 = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            acc1.push(NonNegativeRealScalar::new(3.0));
            acc1.push(NonNegativeRealScalar::new(4.0));

            let mut acc2 = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            acc2.push(NonNegativeRealScalar::new(3.0));
            acc2.push(NonNegativeRealScalar::new(4.0));

            let result = acc1.result();
            let sqrt_result = acc2.result_sqrt();

            // result_sqrt()² should equal result()
            let sqrt_squared = (*sqrt_result.as_raw_scalar_ref()).powi(2);

            assert!((*result.as_raw_scalar_ref() - sqrt_squared).abs() < EPS);
        }

        /// Test with drastically different scales (numerically sensitive)
        #[test]
        fn extreme_scale_difference() {
            let mut acc1 = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            acc1.push(NonNegativeRealScalar::new(1e-8)); // very small

            let mut acc2 = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            acc2.push(NonNegativeRealScalar::new(1e8)); // very large

            acc1.combine(acc2);
            let result = acc1.result_sqrt();

            // Expected: √((1e-8)² + (1e8)²) ≈ 1e8
            let expected = ((1e-8_f64).powi(2) + (1e8_f64).powi(2)).sqrt();
            let rel_err = (*result.as_raw_scalar_ref() - expected).abs() / expected;
            assert!(rel_err < 1e-10);
        }

        /// Test sequential vs combined: process all data at once vs in chunks
        #[test]
        fn sequential_vs_combined() {
            let data = [1.0_f64, 2.0, 3.0, 4.0, 5.0, 6.0];

            // Sequential: process all at once
            let mut acc_seq = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            for &val in &data {
                acc_seq.push(NonNegativeRealScalar::new(val));
            }
            let result_seq = acc_seq.result_sqrt();

            // Chunked: process in two chunks and combine
            let mut acc_chunk1 = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            for &val in &data[0..3] {
                acc_chunk1.push(NonNegativeRealScalar::new(val));
            }

            let mut acc_chunk2 = ScaledSumSqAccumulator::<NaiveSum<f64>>::new();
            for &val in &data[3..] {
                acc_chunk2.push(NonNegativeRealScalar::new(val));
            }

            acc_chunk1.combine(acc_chunk2);
            let result_combined = acc_chunk1.result_sqrt();

            // Results should match
            assert!(
                (*result_seq.as_raw_scalar_ref() - *result_combined.as_raw_scalar_ref()).abs()
                    < EPS
            );
        }
    }
}