ftl_numkernel/
neumaier_compensated_sum.rs

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

ftl_numkernel/
neumaier_compensated_sum.rs
