ftl-numkernel 0.1.0

#![deny(rustdoc::broken_intra_doc_links)]

use crate::{functions::Abs, FunctionsComplexType, FunctionsGeneralType, NumKernel};
use num::Zero;

fn neumaier_sum_and_compensation_real<T: NumKernel>(
    sum_before_compensation: T::RealType,
    value: T::RealType,
) -> (T::RealType, T::RealType) {
    let t = sum_before_compensation.clone() + &value;

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

    (t, compensation_increment)
}

pub trait NeumaierSum<T: NumKernel>: Sized {
    type ScalarType: FunctionsGeneralType<T>;

    /// Creates a new instance of the Neumaier compensated sum with user provided initial `value`.
    /// This is the default constructor.
    fn new(value: Self::ScalarType) -> Self;

    /// Creates a new instance of the Neumaier compensated sum with user initial value set to zero.
    /// This is the default constructor.
    fn zero() -> Self {
        Self::new(Self::ScalarType::zero())
    }

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

    /// Adds a `value` to the sum. This method should be called for each value to be summed.
    fn add(&mut self, value: Self::ScalarType);

    /// Computes and returns the compensated sum.
    fn sum(&self) -> Self::ScalarType;

    /// Resets the sum to zero.
    fn reset(&mut self);
}

//-----------------------------------------------------------------------------------
/// Neumaier compensated sum of an iterable object of floating-point 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 ftl_numkernel::{
///     Native64,
///     neumaier_compensated_sum::{NeumaierSum,NeumaierSumReal}
/// };
///
/// let mut neumaier = NeumaierSumReal::<Native64>::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)]
pub struct NeumaierSumReal<T: NumKernel> {
    /// The sum before the compensation term is added.
    sum_before_compensation: T::RealType,

    /// The compensation term.
    /// This is the correction term that is added to the `sum_before_compensation` to correct for the loss of precision.
    compensation: T::RealType,
}

impl<T: NumKernel> NeumaierSum<T> for NeumaierSumReal<T> {
    type ScalarType = T::RealType;

    /// Creates a new instance of the Neumaier compensated sum with user provided initial `value`.
    /// This is the default constructor.
    fn new(value: T::RealType) -> Self {
        NeumaierSumReal {
            sum_before_compensation: value,
            compensation: T::RealType::zero(),
        }
    }

    /// Adds a `value` to the sum. This method should be called for each value to be summed.
    fn add(&mut self, value: T::RealType) {
        let (t, compensation_increment) =
            neumaier_sum_and_compensation_real::<T>(self.sum_before_compensation.clone(), value);

        self.sum_before_compensation = t;
        self.compensation += compensation_increment;
    }

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

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

//------------------------------------------------------------------------------------
#[derive(Debug, Clone)]
pub struct NeumaierSumComplex<T: NumKernel> {
    /// The sum before the compensation term is added.
    sum_before_compensation: T::ComplexType,

    /// The compensation term.
    /// This is the correction term that is added to the `sum_before_compensation` to correct for the loss of precision.
    compensation: T::ComplexType,
}

impl<T: NumKernel> NeumaierSum<T> for NeumaierSumComplex<T> {
    type ScalarType = T::ComplexType;

    /// Creates a new instance of the Neumaier compensated sum with user provided initial `value`.
    /// This is the default constructor.
    fn new(value: T::ComplexType) -> Self {
        NeumaierSumComplex {
            sum_before_compensation: value,
            compensation: T::ComplexType::zero(),
        }
    }

    /// Adds a `value` to the sum. This method should be called for each value to be summed.
    fn add(&mut self, value: T::ComplexType) {
        {
            // Processing the real part
            let (t, compensation_increment) = neumaier_sum_and_compensation_real::<T>(
                self.sum_before_compensation.real_().clone(),
                value.real_(),
            );

            self.sum_before_compensation.set_real_(t);
            self.compensation.add_real_(&compensation_increment);
        }

        {
            // Processing the imaginary part
            let (t, compensation_increment) = neumaier_sum_and_compensation_real::<T>(
                self.sum_before_compensation.imag_().clone(),
                value.imag_(),
            );

            self.sum_before_compensation.set_imag_(t);
            self.compensation.add_imag_(&compensation_increment);
        }
    }

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

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

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

    mod native64 {
        use super::*;
        use crate::Native64;

        mod real {
            use super::*;

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

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

            #[test]
            fn add() {
                let mut neumaier = NeumaierSumReal::<Native64>::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 = NeumaierSumReal::<Native64>::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 = NeumaierSumReal::<Native64>::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 = NeumaierSumReal::<Native64>::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 = NeumaierSumReal::<Native64>::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 = NeumaierSumComplex::<Native64>::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 = NeumaierSumComplex::<Native64>::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 = NeumaierSumComplex::<Native64>::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.clone());
                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 = NeumaierSumComplex::<Native64>::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 = NeumaierSumComplex::<Native64>::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 = NeumaierSumComplex::<Native64>::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 = NeumaierSumComplex::<Native64>::new_sequential(values);
                assert_eq!(neumaier.sum(), v);
                println!("compensated sum = {}", neumaier.sum());
            }
        }
    }

    #[cfg(feature = "rug")]
    mod rug100 {
        use super::*;
        use crate::Rug;
        use num::One;

        const PRECISION: u32 = 100;

        type Rug100 = Rug<PRECISION>;

        mod real {
            use super::*;
            use crate::{FunctionsRealType, RealRug};

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

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

            #[test]
            fn add() {
                let mut neumaier = NeumaierSumReal::<Rug100>::zero();

                let v = RealRug::<PRECISION>::new(rug::Float::with_val(
                    PRECISION,
                    rug::Float::parse("1e-100").unwrap(),
                ));

                neumaier.add(RealRug::<PRECISION>::try_from_f64_(1.0).unwrap());
                neumaier.add(v.clone());
                neumaier.add(RealRug::<PRECISION>::try_from_f64_(-1.0).unwrap());

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

            #[test]
            fn sum() {
                let mut neumaier = NeumaierSumReal::<Rug100>::zero();

                let v = RealRug::<PRECISION>::new(rug::Float::with_val(
                    PRECISION,
                    rug::Float::parse("1e-100").unwrap(),
                ));

                neumaier.add(RealRug::<PRECISION>::try_from_f64_(1.0).unwrap());
                neumaier.add(v.clone());
                neumaier.add(RealRug::<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 = NeumaierSumReal::<Rug100>::zero();

                let zero = RealRug::<PRECISION>::zero();
                let one = RealRug::<PRECISION>::one();
                let v = RealRug::<PRECISION>::new(rug::Float::with_val(
                    PRECISION,
                    rug::Float::parse("1e-100").unwrap(),
                ));

                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| {
                        RealRug::<PRECISION>::new(rug::Float::with_val(
                            PRECISION,
                            rug::Float::parse(v).unwrap(),
                        ))
                    })
                    .collect::<Vec<_>>();

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

                let neumaier = NeumaierSumReal::<Rug100>::new_sequential(values);
                assert_eq!(
                    neumaier.sum(),
                    RealRug::<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| {
                        RealRug::<PRECISION>::new(rug::Float::with_val(
                            PRECISION,
                            rug::Float::parse(v).unwrap(),
                        ))
                    })
                    .collect::<Vec<_>>();

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

                let neumaier = NeumaierSumReal::<Rug100>::new_sequential(values);
                assert_eq!(
                    neumaier.sum(),
                    RealRug::<PRECISION>::new(rug::Float::with_val(
                        PRECISION,
                        rug::Float::parse("1e-100").unwrap(),
                    ))
                );
                println!("compensated sum = {}", neumaier.sum());
            }
        }

        mod complex {
            use super::*;
            use crate::{ComplexRug, FunctionsComplexType};

            #[test]
            fn new() {
                let one = ComplexRug::<PRECISION>::try_from_f64_(1.0, 0.0).unwrap();
                let neumaier = NeumaierSumComplex::<Rug100>::new(one.clone());
                assert_eq!(neumaier.sum_before_compensation, one);
                assert_eq!(neumaier.compensation, ComplexRug::<PRECISION>::zero());
            }

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

            #[test]
            fn add() {
                let mut neumaier = NeumaierSumComplex::<Rug100>::zero();

                let v = ComplexRug::<PRECISION>::new(rug::Complex::with_val(
                    PRECISION,
                    rug::Complex::parse("(1e-100,2e-100)").unwrap(),
                ));

                neumaier.add(ComplexRug::<PRECISION>::try_from_f64_(1.0, 2.0).unwrap());
                neumaier.add(v.clone());
                neumaier.add(ComplexRug::<PRECISION>::try_from_f64_(-1.0, -2.0).unwrap());

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

            #[test]
            fn sum() {
                let mut neumaier = NeumaierSumComplex::<Rug100>::zero();

                let v = ComplexRug::<PRECISION>::new(rug::Complex::with_val(
                    PRECISION,
                    rug::Complex::parse("(1e-100,2e-100)").unwrap(),
                ));

                neumaier.add(ComplexRug::<PRECISION>::try_from_f64_(1.0, 2.0).unwrap());
                neumaier.add(v.clone());
                neumaier.add(ComplexRug::<PRECISION>::try_from_f64_(-1.0, -2.0).unwrap());

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

            #[test]
            fn reset() {
                let mut neumaier = NeumaierSumComplex::<Rug100>::zero();

                let zero = ComplexRug::<PRECISION>::zero();
                let one = ComplexRug::<PRECISION>::try_from_f64_(1.0, 2.0).unwrap();
                let v = ComplexRug::<PRECISION>::new(rug::Complex::with_val(
                    PRECISION,
                    rug::Complex::parse("(1e-100,2e-100)").unwrap(),
                ));

                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| {
                        ComplexRug::<PRECISION>::new(rug::Complex::with_val(
                            PRECISION,
                            rug::Complex::parse(v).unwrap(),
                        ))
                    })
                    .collect::<Vec<_>>();

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

                let neumaier = NeumaierSumComplex::<Rug100>::new_sequential(values);
                assert_eq!(
                    neumaier.sum(),
                    ComplexRug::<PRECISION>::try_from_f64_(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| {
                        ComplexRug::<PRECISION>::new(rug::Complex::with_val(
                            PRECISION,
                            rug::Complex::parse(v).unwrap(),
                        ))
                    })
                    .collect::<Vec<_>>();

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

                let neumaier = NeumaierSumComplex::<Rug100>::new_sequential(values);
                assert_eq!(
                    neumaier.sum(),
                    ComplexRug::<PRECISION>::new(rug::Complex::with_val(
                        PRECISION,
                        rug::Complex::parse("(1e-100,2e-100)").unwrap(),
                    ))
                );
                println!("compensated sum = {}", neumaier.sum());
            }
        }
    }
}
//------------------------------------------------------------------------------------