num-valid 0.4.0

#![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)
//! ```
//!
//! # 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,
    functions::{Abs, Sign},
    scalars::NonNegativeRealScalar,
};
use getset::Getters;
use num::{Complex, Zero};
use rayon::prelude::*;
use std::ops::{Add, AddAssign, Mul, MulAssign, Sub, SubAssign};
use try_create::TryNew;

//---------------------------------------------------------
/// 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: FpScalar> Accumulator for NaiveSum<T> {
    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: FpScalar> SumAccumulator for NaiveSum<T> {
    #[inline(always)]
    fn rescale_by(&mut self, r: &T) {
        self.0 *= r;
    }
}
//---------------------------------------------------------

//---------------------------------------------------------
#[inline(always)]
fn neumaier_sum_and_compensation_real<RealType>(
    value: RealType,
    sum: &mut RealType,
    compensation: &mut RealType,
) where
    RealType: Clone
        + Add<RealType, Output = RealType>
        + for<'a> Sub<&'a RealType, Output = RealType>
        + AddAssign
        + for<'a> AddAssign<&'a RealType>
        + for<'a> SubAssign<&'a RealType>
        + Abs<Output = RealType>
        + PartialOrd
        + Sign,
{
    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 NeumaierAddable for f64 {
    #[inline(always)]
    fn neumaier_compensated_sum(value: Self, sum: &mut Self, compensation: &mut Self) {
        neumaier_sum_and_compensation_real(value, sum, compensation)
    }
}

impl NeumaierAddable for Complex<f64> {
    #[inline(always)]
    fn neumaier_compensated_sum(value: Self, sum: &mut Self, compensation: &mut Self) {
        neumaier_sum_and_compensation_real(value.re, &mut sum.re, &mut compensation.re);
        neumaier_sum_and_compensation_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: 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: ScalarType::RealType::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> {
        NonNegativeRealScalar::try_new(self.max_abs)
            .expect("MaxAbsValueAccumulator: max of absolute values is negative (bug)")
    }
}
//------------------------------------------------------------------------------------

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

    mod native64 {
        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_eq!(neumaier.sum_before_compensation, 0.0);
                assert_eq!(neumaier.compensation, 0.0);
            }

            #[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_eq!(neumaier.sum_before_compensation, 0.0);
                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_eq!(sum, 0.0);

                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());
            }
        }
    }
}
//------------------------------------------------------------------------------------