num_valid/algorithms/
neumaier_sum.rs

1#![deny(rustdoc::broken_intra_doc_links)]
2
3use crate::functions::Abs;
4use getset::Getters;
5use num::{Complex, Zero};
6use std::ops::{Add, AddAssign, Sub};
7
8//---------------------------------------------------------
9fn neumaier_sum_and_compensation_real<RealType>(
10    value: RealType,
11    sum: &mut RealType,
12    compensation: &mut RealType,
13) where
14    RealType: Clone
15        + Add<RealType, Output = RealType>
16        + for<'a> Sub<&'a RealType, Output = RealType>
17        + AddAssign
18        + for<'a> AddAssign<&'a RealType>
19        + Abs<Output = RealType>
20        + PartialOrd,
21{
22    let sum_before_compensation = sum.clone();
23    *sum += &value;
24
25    // NOTE: the parenthesis are necessary in order for this algorithm to work correctly.
26    *compensation += if sum_before_compensation.clone().abs() >= value.clone().abs() {
27        // If sum is bigger, low-order digits of value are lost.
28        (sum_before_compensation - &*sum) + value
29    } else {
30        // Else low-order digits of sum are lost.
31        (value - &*sum) + sum_before_compensation
32    };
33}
34
35/// A trait for types that can participate in Neumaier compensated summation.
36///
37/// This trait provides the core operation for the Neumaier summation algorithm,
38/// which is a numerically stable method for summing floating-point numbers with
39/// reduced rounding errors compared to naive summation.
40///
41/// # Algorithm
42///
43/// The Neumaier algorithm is an improvement over Kahan summation that handles
44/// the case when the new value is larger than the running sum. It maintains
45/// a compensation term to capture the low-order bits that are lost during addition.
46pub trait NeumaierAddable: Sized {
47    /// Performs one step of the Neumaier compensated summation algorithm.
48    ///
49    /// This method adds `value` to the running `sum` while tracking lost precision
50    /// in the `compensation` term. The final result should be computed as
51    /// `sum + compensation` after all values have been added.
52    ///
53    /// # Arguments
54    ///
55    /// * `value` - The value to add to the sum.
56    /// * `sum` - A mutable reference to the running sum.
57    /// * `compensation` - A mutable reference to the compensation term that tracks lost precision.
58    fn neumaier_compensated_sum(value: Self, sum: &mut Self, compensation: &mut Self);
59}
60
61impl NeumaierAddable for f64 {
62    fn neumaier_compensated_sum(value: Self, sum: &mut Self, compensation: &mut Self) {
63        neumaier_sum_and_compensation_real(value, sum, compensation)
64    }
65}
66
67impl NeumaierAddable for Complex<f64> {
68    fn neumaier_compensated_sum(value: Self, sum: &mut Self, compensation: &mut Self) {
69        neumaier_sum_and_compensation_real(value.re, &mut sum.re, &mut compensation.re);
70        neumaier_sum_and_compensation_real(value.im, &mut sum.im, &mut compensation.im);
71    }
72}
73
74#[cfg(feature = "rug")]
75mod rug_impls {
76    use super::*;
77
78    fn neumaier_sum_and_compensation_rug_float(
79        value: rug::Float,
80        sum: &mut rug::Float,
81        compensation: &mut rug::Float,
82    ) {
83        let sum_before_compensation = sum.clone();
84        *sum += &value;
85
86        // NOTE: the parenthesis are necessary in order for this algorithm to work correctly.
87        *compensation += if sum_before_compensation.clone().abs() >= value.clone().abs() {
88            // If sum is bigger, low-order digits of value are lost.
89            (sum_before_compensation - &*sum) + value
90        } else {
91            // Else low-order digits of sum are lost.
92            (value - &*sum) + sum_before_compensation
93        };
94    }
95
96    impl NeumaierAddable for rug::Float {
97        fn neumaier_compensated_sum(value: Self, sum: &mut Self, compensation: &mut Self) {
98            neumaier_sum_and_compensation_rug_float(value, sum, compensation)
99        }
100    }
101
102    impl NeumaierAddable for rug::Complex {
103        fn neumaier_compensated_sum(value: Self, sum: &mut Self, compensation: &mut Self) {
104            let (value_real, value_imag) = value.into_real_imag();
105
106            neumaier_sum_and_compensation_rug_float(
107                value_real,
108                sum.mut_real(),
109                compensation.mut_real(),
110            );
111
112            neumaier_sum_and_compensation_rug_float(
113                value_imag,
114                sum.mut_imag(),
115                compensation.mut_imag(),
116            );
117        }
118    }
119}
120
121//-----------------------------------------------------------------------------------
122
123//------------------------------------------------------------------------------------
124/// Neumaier compensated sum of an iterable object of floating-point (real or scalar) numbers.
125///
126/// When summing floating-point values in a vector, the standard method of adding values sequentially can lead to precision loss,
127/// especially when there are a large number of elements or a mix of very large and very small values.
128/// This happens because floating-point arithmetic is not associative due to rounding errors.
129/// The most accurate algorithm to sum floating-point numbers typically avoids these precision problems.
130///
131/// Kahan summation is a popular algorithm that reduces numerical errors when adding a sequence of floating-point numbers.
132/// It keeps track of a running compensation for lost low-order bits during summation.
133///
134/// Neumaier summation is a slight modification of Kahan summation that can handle larger errors better.
135/// 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.
136///
137/// # Example
138///
139/// ```
140/// use num_valid::algorithms::neumaier_sum::NeumaierSum;
141///
142/// let mut neumaier = NeumaierSum::<f64>::zero();
143/// neumaier.add(1.0);
144/// neumaier.add(1e100);
145/// neumaier.add(1.0);
146/// neumaier.add(-1e100);
147///
148/// let sum = neumaier.sum();
149/// println!("Sum: {}", sum);
150/// assert_eq!(sum, 2.0);
151/// ```
152///
153/// # References
154///
155/// * [Neumaier Summation](https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements)
156#[derive(Debug, Clone, Getters)]
157pub struct NeumaierSum<ScalarType> {
158    /// The sum before the compensation term is added.
159    #[getset(get = "pub")]
160    sum_before_compensation: ScalarType,
161
162    /// The compensation term.
163    /// This is the correction term that is added to the `sum_before_compensation` to correct for the loss of precision.
164    #[getset(get = "pub")]
165    compensation: ScalarType,
166}
167
168impl<ScalarType> NeumaierSum<ScalarType>
169where
170    ScalarType: Clone + Zero + for<'a> Add<&'a ScalarType, Output = ScalarType> + NeumaierAddable,
171{
172    /// Creates a new instance of the Neumaier compensated sum, initializing the
173    /// sum with the provided `value`.
174    pub fn new(value: ScalarType) -> Self {
175        Self {
176            sum_before_compensation: value,
177            compensation: ScalarType::zero(),
178        }
179    }
180
181    /// Creates a new instance of the Neumaier compensated sum, initialized to zero.
182    pub fn zero() -> Self {
183        Self::new(ScalarType::zero())
184    }
185
186    /// Computes and returns the compensated sum.
187    pub fn sum(&self) -> ScalarType {
188        self.sum_before_compensation.clone() + &self.compensation
189    }
190
191    /// Resets the sum to zero.
192    pub fn reset(&mut self) {
193        self.sum_before_compensation = ScalarType::zero();
194        self.compensation = ScalarType::zero();
195    }
196
197    /// Adds a `value` to the sum. This method should be called for each value to be summed.
198    pub fn add(&mut self, value: ScalarType) {
199        <ScalarType as NeumaierAddable>::neumaier_compensated_sum(
200            value,
201            &mut self.sum_before_compensation,
202            &mut self.compensation,
203        );
204    }
205
206    /// Creates a new instance of the [`NeumaierSum`] object summing the values in the iterable object `values` (sequential_version).
207    /// # Example
208    ///
209    /// ```
210    /// use num_valid::algorithms::neumaier_sum::NeumaierSum;
211    ///
212    /// let values = vec![1.0, 1.0e100, 1.0, -1.0e100];
213    ///
214    /// let neumaier = NeumaierSum::new_sequential(values);
215    /// let sum = neumaier.sum();
216    /// println!("Sum: {}", sum);
217    /// assert_eq!(sum, 2.0);
218    /// ```
219    pub fn new_sequential<I>(values: I) -> Self
220    where
221        I: IntoIterator<Item = ScalarType>,
222    {
223        let mut neumaier = Self::zero();
224        values.into_iter().for_each(|value| {
225            neumaier.add(value);
226        });
227        neumaier
228    }
229}
230
231//------------------------------------------------------------------------------------
232
233//------------------------------------------------------------------------------------
234#[cfg(test)]
235mod tests_neumaier_sum {
236    use super::*;
237
238    mod native64 {
239        use super::*;
240
241        mod real {
242            use super::*;
243
244            #[test]
245            fn new() {
246                let neumaier = NeumaierSum::new(1.0);
247                assert_eq!(neumaier.sum_before_compensation, 1.0);
248                assert_eq!(neumaier.compensation, 0.0);
249            }
250
251            #[test]
252            fn zero() {
253                let neumaier = NeumaierSum::<f64>::zero();
254                assert_eq!(neumaier.sum_before_compensation, 0.0);
255                assert_eq!(neumaier.compensation, 0.0);
256            }
257
258            #[test]
259            fn add() {
260                let mut neumaier = NeumaierSum::<f64>::zero();
261                neumaier.add(1.0);
262                neumaier.add(1e-16);
263                neumaier.add(-1.0);
264                assert_eq!(neumaier.sum_before_compensation, 0.0);
265                assert_eq!(neumaier.compensation, 1e-16);
266            }
267
268            #[test]
269            fn sum() {
270                let mut neumaier = NeumaierSum::<f64>::zero();
271                neumaier.add(1.0);
272                neumaier.add(1e-16);
273                neumaier.add(-1.0);
274                assert_eq!(neumaier.sum_before_compensation, 0.0);
275                assert_eq!(neumaier.compensation, 1e-16);
276                assert_eq!(neumaier.sum(), 1e-16);
277                println!("compensated sum = {}", neumaier.sum());
278            }
279
280            #[test]
281            fn reset() {
282                let mut neumaier = NeumaierSum::<f64>::zero();
283                neumaier.add(1.0);
284                neumaier.add(1e-16);
285                assert_eq!(neumaier.sum_before_compensation, 1.0);
286                assert_eq!(neumaier.compensation, 1e-16);
287
288                neumaier.reset();
289                assert_eq!(neumaier.sum_before_compensation, 0.0);
290                assert_eq!(neumaier.compensation, 0.0);
291            }
292
293            #[test]
294            fn sum_big_values() {
295                let values = vec![1.0, 1e100, 1.0, -1e100];
296                let sum = values.iter().sum::<f64>();
297                assert_eq!(sum, 0.0);
298
299                let neumaier = NeumaierSum::<f64>::new_sequential(values);
300                assert_eq!(neumaier.sum(), 2.0);
301                println!("compensated sum = {}", neumaier.sum());
302            }
303
304            #[test]
305            fn sum_small_values() {
306                let values = [1.0, 1e-100, -1.0];
307                let sum = values.iter().sum::<f64>();
308                assert_eq!(sum, 0.0);
309
310                let neumaier = NeumaierSum::<f64>::new_sequential(values);
311                assert_eq!(neumaier.sum(), 1e-100);
312                println!("compensated sum = {}", neumaier.sum());
313            }
314        }
315
316        mod complex {
317            use super::*;
318            use num::Complex;
319
320            #[test]
321            fn new() {
322                let neumaier = NeumaierSum::<Complex<f64>>::new(Complex::new(1.0, 2.0));
323                assert_eq!(neumaier.sum_before_compensation, Complex::new(1.0, 2.0));
324                assert_eq!(neumaier.compensation, Complex::new(0.0, 0.0));
325            }
326
327            #[test]
328            fn zero() {
329                let neumaier = NeumaierSum::<Complex<f64>>::zero();
330
331                let zero = Complex::new(0.0, 0.0);
332                assert_eq!(&neumaier.sum_before_compensation, &zero);
333                assert_eq!(&neumaier.compensation, &zero);
334            }
335
336            #[test]
337            fn add() {
338                let mut neumaier = NeumaierSum::<Complex<f64>>::zero();
339
340                let zero = Complex::new(0.0, 0.0);
341                let v = Complex::new(1e-16, 2e-16);
342
343                neumaier.add(Complex::new(1.0, 2.0));
344                neumaier.add(v);
345                neumaier.add(Complex::new(-1.0, -2.0));
346
347                assert_eq!(neumaier.sum_before_compensation, zero);
348                assert_eq!(neumaier.compensation, v);
349            }
350
351            #[test]
352            fn sum() {
353                let mut neumaier = NeumaierSum::<Complex<f64>>::zero();
354
355                let zero = Complex::new(0.0, 0.0);
356                let v = Complex::new(1e-16, 2e-16);
357
358                neumaier.add(Complex::new(1.0, 2.0));
359                neumaier.add(v);
360                neumaier.add(Complex::new(-1.0, -2.0));
361                assert_eq!(neumaier.sum_before_compensation, zero);
362                assert_eq!(neumaier.compensation, v);
363                assert_eq!(neumaier.sum(), v);
364                println!("compensated sum = {}", neumaier.sum());
365            }
366
367            #[test]
368            fn reset() {
369                let mut neumaier = NeumaierSum::<Complex<f64>>::zero();
370
371                let zero = Complex::new(0.0, 0.0);
372                let a = Complex::new(1.0, 2.0);
373                let v = Complex::new(1e-16, 2e-16);
374
375                neumaier.add(a);
376                neumaier.add(v);
377                assert_eq!(neumaier.sum_before_compensation, a);
378                assert_eq!(neumaier.compensation, v);
379
380                neumaier.reset();
381                assert_eq!(neumaier.sum_before_compensation, zero);
382                assert_eq!(neumaier.compensation, zero);
383            }
384
385            #[test]
386            fn sum_big_values() {
387                let values = vec![
388                    Complex::new(1.0, 2.0),
389                    Complex::new(1e100, 2e100),
390                    Complex::new(1.0, 2.0),
391                    Complex::new(-1e100, -2e100),
392                ];
393                let sum = values.iter().sum::<Complex<f64>>();
394                assert_eq!(sum, Complex::new(0.0, 0.0));
395
396                let neumaier = NeumaierSum::<Complex<f64>>::new_sequential(values);
397                assert_eq!(neumaier.sum(), Complex::new(2.0, 4.0));
398                println!("compensated sum = {}", neumaier.sum());
399            }
400
401            #[test]
402            fn sum_small_values() {
403                let v = Complex::new(1e-100, 2e-100);
404
405                let values = [Complex::new(1.0, 2.0), v, Complex::new(-1.0, -2.0)];
406                let sum = values.iter().sum::<Complex<f64>>();
407                assert_eq!(sum, Complex::new(0.0, 0.0));
408
409                let neumaier = NeumaierSum::<Complex<f64>>::new_sequential(values);
410                assert_eq!(neumaier.sum(), v);
411                println!("compensated sum = {}", neumaier.sum());
412            }
413        }
414    }
415
416    #[cfg(feature = "rug")]
417    mod rug100 {
418        use super::*;
419        use crate::{ComplexRugStrictFinite, RealRugStrictFinite};
420        use num::One;
421        use try_create::TryNew;
422
423        const PRECISION: u32 = 100;
424
425        mod real {
426            use super::*;
427            use crate::RealScalar;
428
429            #[test]
430            fn new() {
431                let neumaier = NeumaierSum::new(RealRugStrictFinite::<PRECISION>::one());
432                assert_eq!(neumaier.sum_before_compensation, 1.0);
433                assert_eq!(neumaier.compensation, 0.0);
434            }
435
436            #[test]
437            fn zero() {
438                let neumaier = NeumaierSum::<RealRugStrictFinite<PRECISION>>::zero();
439                assert_eq!(neumaier.sum_before_compensation, 0.0);
440                assert_eq!(neumaier.compensation, 0.0);
441            }
442
443            #[test]
444            fn add() {
445                let mut neumaier = NeumaierSum::<RealRugStrictFinite<PRECISION>>::zero();
446
447                let v = RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
448                    PRECISION,
449                    rug::Float::parse("1e-100").unwrap(),
450                ))
451                .expect("valid test value");
452
453                neumaier.add(RealRugStrictFinite::<PRECISION>::try_from_f64(1.0).unwrap());
454                neumaier.add(v.clone());
455                neumaier.add(RealRugStrictFinite::<PRECISION>::try_from_f64(-1.0).unwrap());
456
457                assert_eq!(
458                    neumaier.sum_before_compensation,
459                    RealRugStrictFinite::<PRECISION>::try_from_f64(0.0).unwrap()
460                );
461                assert_eq!(&neumaier.compensation, &v);
462            }
463
464            #[test]
465            fn sum() {
466                let mut neumaier = NeumaierSum::<RealRugStrictFinite<PRECISION>>::zero();
467
468                let v = RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
469                    PRECISION,
470                    rug::Float::parse("1e-100").unwrap(),
471                ))
472                .expect("valid test value");
473
474                neumaier.add(RealRugStrictFinite::<PRECISION>::try_from_f64(1.0).unwrap());
475                neumaier.add(v.clone());
476                neumaier.add(RealRugStrictFinite::<PRECISION>::try_from_f64(-1.0).unwrap());
477
478                assert_eq!(neumaier.sum_before_compensation, 0.0);
479                assert_eq!(&neumaier.compensation, &v);
480                assert_eq!(neumaier.sum(), v);
481                println!("compensated sum = {}", neumaier.sum());
482            }
483
484            #[test]
485            fn reset() {
486                let mut neumaier = NeumaierSum::<RealRugStrictFinite<PRECISION>>::zero();
487
488                let zero = RealRugStrictFinite::<PRECISION>::zero();
489                let one = RealRugStrictFinite::<PRECISION>::one();
490                let v = RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
491                    PRECISION,
492                    rug::Float::parse("1e-100").unwrap(),
493                ))
494                .expect("valid test value");
495
496                neumaier.add(one.clone());
497                neumaier.add(v.clone());
498
499                assert_eq!(&neumaier.sum_before_compensation, &one);
500                assert_eq!(&neumaier.compensation, &v);
501
502                neumaier.reset();
503                assert_eq!(&neumaier.sum_before_compensation, &zero);
504                assert_eq!(&neumaier.compensation, &zero);
505            }
506
507            #[test]
508            fn sum_big_values() {
509                let values = ["1.0", "1e100", "1.0", "-1e100"]
510                    .iter()
511                    .map(|v| {
512                        RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
513                            PRECISION,
514                            rug::Float::parse(v).unwrap(),
515                        ))
516                        .expect("valid test value")
517                    })
518                    .collect::<Vec<_>>();
519
520                let sum = values
521                    .iter()
522                    .fold(RealRugStrictFinite::<PRECISION>::zero(), |acc, x| acc + x);
523                assert_eq!(sum, 0.0);
524
525                let neumaier =
526                    NeumaierSum::<RealRugStrictFinite<PRECISION>>::new_sequential(values);
527                assert_eq!(
528                    neumaier.sum(),
529                    RealRugStrictFinite::<PRECISION>::try_from_f64(2.0).unwrap()
530                );
531                println!("compensated sum = {}", neumaier.sum());
532            }
533
534            #[test]
535            fn sum_small_values() {
536                let values = ["1.0", "1e-100", "-1.0"]
537                    .iter()
538                    .map(|v| {
539                        RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
540                            PRECISION,
541                            rug::Float::parse(v).unwrap(),
542                        ))
543                        .expect("valid test value")
544                    })
545                    .collect::<Vec<_>>();
546
547                let sum = values
548                    .iter()
549                    .fold(RealRugStrictFinite::<PRECISION>::zero(), |acc, x| acc + x);
550                assert_eq!(sum, RealRugStrictFinite::<PRECISION>::zero());
551
552                let neumaier =
553                    NeumaierSum::<RealRugStrictFinite<PRECISION>>::new_sequential(values);
554                assert_eq!(
555                    neumaier.sum(),
556                    RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
557                        PRECISION,
558                        rug::Float::parse("1e-100").unwrap(),
559                    ))
560                    .expect("valid test value")
561                );
562                println!("compensated sum = {}", neumaier.sum());
563            }
564        }
565
566        mod complex {
567            use super::*;
568
569            #[test]
570            fn new() {
571                let one = ComplexRugStrictFinite::<PRECISION>::one();
572                let neumaier = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new(one.clone());
573                assert_eq!(neumaier.sum_before_compensation, one);
574                assert_eq!(
575                    neumaier.compensation,
576                    ComplexRugStrictFinite::<PRECISION>::zero()
577                );
578            }
579
580            #[test]
581            fn zero() {
582                let neumaier = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::zero();
583                assert_eq!(
584                    neumaier.sum_before_compensation,
585                    ComplexRugStrictFinite::<PRECISION>::zero()
586                );
587                assert_eq!(
588                    neumaier.compensation,
589                    ComplexRugStrictFinite::<PRECISION>::zero()
590                );
591            }
592
593            #[test]
594            fn add() {
595                let mut neumaier = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::zero();
596
597                let v = ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
598                    PRECISION,
599                    rug::Complex::parse("(1e-100,2e-100)").unwrap(),
600                ))
601                .expect("valid test value");
602
603                neumaier.add(
604                    ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(1.0, 2.0)).unwrap(),
605                );
606                neumaier.add(v.clone());
607                neumaier.add(
608                    ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(-1.0, -2.0))
609                        .unwrap(),
610                );
611
612                assert_eq!(
613                    neumaier.sum_before_compensation,
614                    ComplexRugStrictFinite::<PRECISION>::zero()
615                );
616                assert_eq!(&neumaier.compensation, &v);
617            }
618
619            #[test]
620            fn sum() {
621                let mut neumaier = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::zero();
622
623                let v = ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
624                    PRECISION,
625                    rug::Complex::parse("(1e-100,2e-100)").unwrap(),
626                ))
627                .expect("valid test value");
628
629                neumaier.add(
630                    ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(1.0, 2.0)).unwrap(),
631                );
632                neumaier.add(v.clone());
633                neumaier.add(
634                    ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(-1.0, -2.0))
635                        .unwrap(),
636                );
637
638                assert_eq!(
639                    neumaier.sum_before_compensation,
640                    ComplexRugStrictFinite::<PRECISION>::zero()
641                );
642                assert_eq!(&neumaier.compensation, &v);
643                assert_eq!(neumaier.sum(), v);
644                println!("compensated sum = {}", neumaier.sum());
645            }
646
647            #[test]
648            fn reset() {
649                let mut neumaier = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::zero();
650
651                let zero = ComplexRugStrictFinite::<PRECISION>::zero();
652                let one =
653                    ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(1.0, 2.0)).unwrap();
654                let v = ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
655                    PRECISION,
656                    rug::Complex::parse("(1e-100,2e-100)").unwrap(),
657                ))
658                .expect("valid test value");
659
660                neumaier.add(one.clone());
661                neumaier.add(v.clone());
662
663                assert_eq!(&neumaier.sum_before_compensation, &one);
664                assert_eq!(&neumaier.compensation, &v);
665
666                neumaier.reset();
667                assert_eq!(&neumaier.sum_before_compensation, &zero);
668                assert_eq!(&neumaier.compensation, &zero);
669            }
670
671            #[test]
672            fn sum_big_values() {
673                let values = ["(1.0,2.0)", "(1e100,2e100)", "(1.0,2.0)", "(-1e100,-2e100)"]
674                    .iter()
675                    .map(|v| {
676                        ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
677                            PRECISION,
678                            rug::Complex::parse(v).unwrap(),
679                        ))
680                        .expect("valid test value")
681                    })
682                    .collect::<Vec<_>>();
683
684                let zero = ComplexRugStrictFinite::<PRECISION>::zero();
685                let sum = values.iter().fold(zero.clone(), |acc, x| acc + x);
686                assert_eq!(sum, zero);
687
688                let neumaier =
689                    NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new_sequential(values);
690                assert_eq!(
691                    neumaier.sum(),
692                    ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(2.0, 4.0)).unwrap()
693                );
694                println!("compensated sum = {}", neumaier.sum());
695            }
696
697            #[test]
698            fn sum_small_values() {
699                let values = ["(1.0,2.0)", "(1e-100,2e-100)", "(-1.0,-2.0)"]
700                    .iter()
701                    .map(|v| {
702                        ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
703                            PRECISION,
704                            rug::Complex::parse(v).unwrap(),
705                        ))
706                        .expect("valid test value")
707                    })
708                    .collect::<Vec<_>>();
709
710                let zero = ComplexRugStrictFinite::<PRECISION>::zero();
711                let sum = values.iter().fold(zero.clone(), |acc, x| acc + x);
712                assert_eq!(sum, zero);
713
714                let neumaier =
715                    NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new_sequential(values);
716                assert_eq!(
717                    neumaier.sum(),
718                    ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
719                        PRECISION,
720                        rug::Complex::parse("(1e-100,2e-100)").unwrap(),
721                    ))
722                    .expect("valid test value")
723                );
724                println!("compensated sum = {}", neumaier.sum());
725            }
726        }
727    }
728}
729//------------------------------------------------------------------------------------
num_valid/algorithms/neumaier_sum.rs

num_valid/algorithms/
neumaier_sum.rs
