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

num_valid/
neumaier_compensated_sum.rs
