num-valid 0.3.3

#![deny(rustdoc::broken_intra_doc_links)]

use crate::functions::Abs;
use getset::Getters;
use num::{Complex, Zero};
use std::ops::{Add, AddAssign, Sub};

//---------------------------------------------------------
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>
        + Abs<Output = RealType>
        + PartialOrd,
{
    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
    };
}

/// 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 {
    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> {
    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::*;

    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 {
        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 {
        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 sum of an iterable object of floating-point (real or scalar) numbers.
///
/// When summing floating-point values in a vector, the standard method of adding values sequentially can lead to precision loss,
/// especially when there are a large number of elements or a mix of very large and very small values.
/// This happens because floating-point arithmetic is not associative due to rounding errors.
/// The most accurate algorithm to sum floating-point numbers typically avoids these precision problems.
///
/// Kahan summation is a popular algorithm that reduces numerical errors when adding a sequence of floating-point numbers.
/// It keeps track of a running compensation for lost low-order bits during summation.
///
/// Neumaier summation is a slight modification of Kahan summation that can handle larger errors better.
/// It uses an extra step to correct for the compensation term if the summation results in a larger round-off error than Kahan’s method can correct.
///
/// # Example
///
/// ```
/// use num_valid::algorithms::neumaier_sum::NeumaierSum;
///
/// let mut neumaier = NeumaierSum::<f64>::zero();
/// neumaier.add(1.0);
/// neumaier.add(1e100);
/// neumaier.add(1.0);
/// neumaier.add(-1e100);
///
/// let sum = neumaier.sum();
/// 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> NeumaierSum<ScalarType>
where
    ScalarType: Clone + Zero + for<'a> Add<&'a ScalarType, Output = ScalarType> + NeumaierAddable,
{
    /// Creates a new instance of the Neumaier compensated sum, initializing the
    /// sum with the provided `value`.
    pub fn new(value: ScalarType) -> Self {
        Self {
            sum_before_compensation: value,
            compensation: ScalarType::zero(),
        }
    }

    /// Creates a new instance of the Neumaier compensated sum, initialized to zero.
    pub fn zero() -> Self {
        Self::new(ScalarType::zero())
    }

    /// Computes and returns the compensated sum.
    pub fn sum(&self) -> ScalarType {
        self.sum_before_compensation.clone() + &self.compensation
    }

    /// Resets the sum to zero.
    pub fn reset(&mut self) {
        self.sum_before_compensation = ScalarType::zero();
        self.compensation = ScalarType::zero();
    }

    /// Adds a `value` to the sum. This method should be called for each value to be summed.
    pub fn add(&mut self, value: ScalarType) {
        <ScalarType as NeumaierAddable>::neumaier_compensated_sum(
            value,
            &mut self.sum_before_compensation,
            &mut self.compensation,
        );
    }

    /// Creates a new instance of the [`NeumaierSum`] object summing the values in the iterable object `values` (sequential_version).
    /// # Example
    ///
    /// ```
    /// use num_valid::algorithms::neumaier_sum::NeumaierSum;
    ///
    /// let values = vec![1.0, 1.0e100, 1.0, -1.0e100];
    ///
    /// let neumaier = NeumaierSum::new_sequential(values);
    /// let sum = neumaier.sum();
    /// println!("Sum: {}", sum);
    /// assert_eq!(sum, 2.0);
    /// ```
    pub fn new_sequential<I>(values: I) -> Self
    where
        I: IntoIterator<Item = ScalarType>,
    {
        let mut neumaier = Self::zero();
        values.into_iter().for_each(|value| {
            neumaier.add(value);
        });
        neumaier
    }
}

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

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

    mod native64 {
        use super::*;

        mod real {
            use super::*;

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

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

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

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

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

                neumaier.reset();
                assert_eq!(neumaier.sum_before_compensation, 0.0);
                assert_eq!(neumaier.compensation, 0.0);
            }

            #[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);
                assert_eq!(neumaier.sum(), 2.0);
                println!("compensated sum = {}", neumaier.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);
                assert_eq!(neumaier.sum(), 1e-100);
                println!("compensated sum = {}", neumaier.sum());
            }
        }

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

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

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

                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>>::zero();

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

                neumaier.add(Complex::new(1.0, 2.0));
                neumaier.add(v);
                neumaier.add(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>>::zero();

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

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

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

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

                neumaier.add(a);
                neumaier.add(v);
                assert_eq!(neumaier.sum_before_compensation, a);
                assert_eq!(neumaier.compensation, v);

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

            #[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);
                assert_eq!(neumaier.sum(), Complex::new(2.0, 4.0));
                println!("compensated sum = {}", neumaier.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);
                assert_eq!(neumaier.sum(), v);
                println!("compensated sum = {}", neumaier.sum());
            }
        }
    }

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

        const PRECISION: u32 = 100;

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

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

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

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

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

                neumaier.add(RealRugStrictFinite::<PRECISION>::try_from_f64(1.0).unwrap());
                neumaier.add(v.clone());
                neumaier.add(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>>::zero();

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

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

                assert_eq!(neumaier.sum_before_compensation, 0.0);
                assert_eq!(&neumaier.compensation, &v);
                assert_eq!(neumaier.sum(), v);
                println!("compensated sum = {}", neumaier.sum());
            }

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

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

                neumaier.add(one.clone());
                neumaier.add(v.clone());

                assert_eq!(&neumaier.sum_before_compensation, &one);
                assert_eq!(&neumaier.compensation, &v);

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

            #[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);
                assert_eq!(
                    neumaier.sum(),
                    RealRugStrictFinite::<PRECISION>::try_from_f64(2.0).unwrap()
                );
                println!("compensated sum = {}", neumaier.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);
                assert_eq!(
                    neumaier.sum(),
                    RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
                        PRECISION,
                        rug::Float::parse("1e-100").unwrap(),
                    ))
                    .expect("valid test value")
                );
                println!("compensated sum = {}", neumaier.sum());
            }
        }

        mod complex {
            use super::*;

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

            #[test]
            fn zero() {
                let neumaier = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::zero();
                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>>::zero();

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

                neumaier.add(
                    ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(1.0, 2.0)).unwrap(),
                );
                neumaier.add(v.clone());
                neumaier.add(
                    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>>::zero();

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

                neumaier.add(
                    ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(1.0, 2.0)).unwrap(),
                );
                neumaier.add(v.clone());
                neumaier.add(
                    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);
                assert_eq!(neumaier.sum(), v);
                println!("compensated sum = {}", neumaier.sum());
            }

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

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

                neumaier.add(one.clone());
                neumaier.add(v.clone());

                assert_eq!(&neumaier.sum_before_compensation, &one);
                assert_eq!(&neumaier.compensation, &v);

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

            #[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);
                assert_eq!(
                    neumaier.sum(),
                    ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(2.0, 4.0)).unwrap()
                );
                println!("compensated sum = {}", neumaier.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);
                assert_eq!(
                    neumaier.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 = {}", neumaier.sum());
            }
        }
    }
}
//------------------------------------------------------------------------------------