num_valid/algorithms/
accumulators.rs

1#![deny(rustdoc::broken_intra_doc_links)]
2
3//! Pluggable accumulation strategies for numerical algorithms.
4//!
5//! This module defines the [`Accumulator`] trait and a set of concrete
6//! implementations used internally by
7//! [`crate::algorithms::vector_norms`] and available for use in custom
8//! algorithms.
9//!
10//! # Trait hierarchy
11//!
12//! ```text
13//! Accumulator
14//! ├── SumAccumulator    (adds rescale_by; used for L1/L2 norms)
15//! │   ├── NaiveSum
16//! │   └── NeumaierSum
17//! ├── MaxAccumulator    (running maximum of real values)
18//! ├── MinAccumulator    (running minimum of real values)
19//! └── MaxAbsValueAccumulator  (running maximum of |x|; used for L∞ norm)
20//! ```
21//!
22//! # Choosing a summation strategy
23//!
24//! | Strategy | Speed | Accuracy | Use when |
25//! |----------|-------|----------|----------|
26//! | [`NaiveSum`] | fastest | O(N·ε) error | data is well-scaled |
27//! | [`NeumaierSum`] | ~2× slower | O(ε) error | components span many magnitudes |
28//!
29//! # Parallel reduction
30//!
31//! All accumulators support parallel chunked reductions via the default
32//! [`Accumulator::new_parallel`] method, which uses rayon `fold` + `reduce`
33//! with [`Accumulator::combine`] as the merge operation. `Accumulator::new()`
34//! must be the identity element for `combine` — this holds for all types in
35//! this module.
36
37use crate::{
38    FpScalar, RealScalar,
39    functions::{Abs, Sign},
40    scalars::NonNegativeRealScalar,
41};
42use getset::Getters;
43use num::{Complex, Zero};
44use rayon::prelude::*;
45use std::ops::{Add, AddAssign, Mul, MulAssign, Sub, SubAssign};
46use try_create::TryNew;
47
48//---------------------------------------------------------
49/// Base trait for incremental accumulation over a sequence of values.
50///
51/// An accumulator collects items one at a time via [`push`](Self::push) and
52/// produces a final result via [`result`](Self::result). Partial accumulators
53/// can be merged with [`combine`](Self::combine), enabling chunked and parallel
54/// reductions.
55///
56/// # Correctness requirement for `combine`
57///
58/// [`new()`](Self::new) must be the identity element for `combine`: for any
59/// accumulator `x`, both `combine(new(), x)` and `combine(x, new())` must
60/// produce an accumulator equivalent to `x`. This invariant is required by
61/// [`new_parallel`](Self::new_parallel).
62pub trait Accumulator: Sized {
63    /// The type of values pushed into the accumulator.
64    type Input;
65    /// The type of the final result produced by [`result`](Self::result).
66    type Output;
67
68    /// Creates a new, empty accumulator (the identity for [`combine`](Self::combine)).
69    fn new() -> Self;
70
71    /// Adds a single value to the accumulator.
72    fn push(&mut self, value: Self::Input);
73
74    /// Merges `other` into `self`.
75    ///
76    /// After this call `self` represents the accumulation of all values pushed
77    /// into either `self` or `other`. Used by [`new_parallel`](Self::new_parallel)
78    /// to merge per-chunk partial results.
79    fn combine(&mut self, other: Self);
80
81    /// Consumes the accumulator and returns the final result.
82    fn result(self) -> Self::Output;
83
84    /// Creates a new accumulator by sequentially adding all values from `values`.
85    ///
86    /// This is the simplest way to build an accumulator from owned input data.
87    /// For chunked or parallel reductions, prefer [`Accumulator::combine`]
88    /// on partial accumulators.
89    ///
90    /// # Example
91    ///
92    /// ```
93    /// use num_valid::algorithms::accumulators::*;
94    ///
95    /// let values = vec![1.0, 1.0e100, 1.0, -1.0e100];
96    ///
97    /// let neumaier = NeumaierSum::new_sequential(values);
98    /// let sum = neumaier.result();
99    /// println!("Sum: {}", sum);
100    /// assert_eq!(sum, 2.0);
101    /// ```
102    fn new_sequential<I>(values: I) -> Self
103    where
104        I: IntoIterator<Item = Self::Input>,
105    {
106        let mut acc = Self::new();
107        values.into_iter().for_each(|v| acc.push(v));
108        acc
109    }
110
111    /// Builds an accumulator by processing a parallel iterator.
112    ///
113    /// The iterator is partitioned into chunks by rayon; each chunk is folded
114    /// into a partial accumulator with [`push`](Self::push), and the partial
115    /// results are merged with [`combine`](Self::combine).
116    ///
117    /// Requires `Self: Send`.
118    ///
119    /// Any additional `Send` constraints on iterator items are enforced by the
120    /// concrete [`IntoParallelIterator`] implementation.
121    ///
122    /// # Example
123    ///
124    /// ```
125    /// use num_valid::algorithms::accumulators::*;
126    /// use rayon::prelude::*;
127    ///
128    /// let values = vec![1.0_f64, 2.0, 3.0, 4.0];
129    /// let sum = NaiveSum::new_parallel(values.into_par_iter()).result();
130    /// assert_eq!(sum, 10.0);
131    /// ```
132    fn new_parallel<I>(values: I) -> Self
133    where
134        Self: Send,
135        I: IntoParallelIterator<Item = Self::Input>,
136    {
137        values
138            .into_par_iter()
139            .fold(
140                || Self::new(),
141                |mut acc, item| {
142                    acc.push(item);
143                    acc
144                },
145            )
146            .reduce(
147                || Self::new(),
148                |mut a, b| {
149                    a.combine(b);
150                    a
151                },
152            )
153    }
154}
155//---------------------------------------------------------
156
157//---------------------------------------------------------
158/// A generic accumulator for summing a sequence of values with configurable
159/// precision/performance trade-offs.
160///
161/// The trait abstracts over how terms are collected, allowing different
162/// accumulation strategies to be swapped in without duplicating the calling
163/// code.  It is used internally by the [`crate::algorithms::vector_norms`] vector-norm functions (scaled
164/// sum-of-squares loop), but is general enough for any sequential summation
165/// task.  Two strategies ship with the crate:
166///
167/// - **Naive** ([`NaiveSum`]): plain `+=` —
168///   fast but accumulates O(N·ε) rounding error.
169/// - **Neumaier** ([`NeumaierSum`]): Kahan/Neumaier compensated
170///   summation — recovers O(ε) rounding error at a small overhead cost.
171///
172/// Implement this trait to plug a custom accumulation strategy wherever it is
173/// accepted.
174pub trait SumAccumulator:
175    Accumulator<Input: Zero + Clone, Output = <Self as Accumulator>::Input>
176{
177    /// Multiply the entire accumulated sum by `r`.
178    fn rescale_by(&mut self, r: &Self::Input);
179}
180//---------------------------------------------------------
181
182//---------------------------------------------------------
183/// Naive (plain `+=`) accumulator.
184///
185/// This strategy keeps a single running sum and updates it with standard
186/// floating-point addition. It has minimal overhead and is usually the fastest
187/// option, but it may lose small addends when values have very different
188/// magnitudes.
189///
190/// In norm routines this is used by the default APIs such as
191/// [`crate::algorithms::vector_norms::vector_norm_l1`],
192/// [`crate::algorithms::vector_norms::vector_norm_l2`], and
193/// [`crate::algorithms::vector_norms::vector_norm_l2_sq`].
194///
195/// # Example
196///
197/// ```
198/// use num_valid::algorithms::accumulators::*;
199///
200/// let mut naive = NaiveSum::<f64>::new();
201/// naive.push(1.0);
202/// naive.push(2.0);
203/// naive.push(3.0);
204/// assert_eq!(naive.result(), 6.0);
205/// ```
206pub struct NaiveSum<T>(T);
207
208impl<T: FpScalar> Accumulator for NaiveSum<T> {
209    type Input = T;
210    type Output = T;
211
212    #[inline(always)]
213    fn new() -> Self {
214        NaiveSum(T::zero())
215    }
216
217    #[inline(always)]
218    fn push(&mut self, r: T) {
219        self.0 += r;
220    }
221
222    #[inline(always)]
223    fn result(self) -> T {
224        self.0
225    }
226
227    #[inline(always)]
228    fn combine(&mut self, other: Self) {
229        self.0 += other.0;
230    }
231}
232
233impl<T: FpScalar> SumAccumulator for NaiveSum<T> {
234    #[inline(always)]
235    fn rescale_by(&mut self, r: &T) {
236        self.0 *= r;
237    }
238}
239//---------------------------------------------------------
240
241//---------------------------------------------------------
242#[inline(always)]
243fn neumaier_sum_and_compensation_real<RealType>(
244    value: RealType,
245    sum: &mut RealType,
246    compensation: &mut RealType,
247) where
248    RealType: Clone
249        + Add<RealType, Output = RealType>
250        + for<'a> Sub<&'a RealType, Output = RealType>
251        + AddAssign
252        + for<'a> AddAssign<&'a RealType>
253        + for<'a> SubAssign<&'a RealType>
254        + Abs<Output = RealType>
255        + PartialOrd
256        + Sign,
257{
258    let sum_before_compensation = sum.clone();
259    *sum += &value; // here we may lose low-order digits of value if sum is much larger than value or vice versa
260
261    // NOTE: the parenthesis are necessary in order for this algorithm to work correctly.
262    *compensation += if sum_before_compensation.clone().abs() >= value.clone().abs() {
263        // If sum is bigger, low-order digits of value are lost.
264        (sum_before_compensation - &*sum) + value
265    } else {
266        // Else low-order digits of sum are lost.
267        (value - &*sum) + sum_before_compensation
268    };
269}
270
271/// A trait for types that can participate in Neumaier compensated summation.
272///
273/// This trait provides the core operation for the Neumaier summation algorithm,
274/// which is a numerically stable method for summing floating-point numbers with
275/// reduced rounding errors compared to naive summation.
276///
277/// # Algorithm
278///
279/// The Neumaier algorithm is an improvement over Kahan summation that handles
280/// the case when the new value is larger than the running sum. It maintains
281/// a compensation term to capture the low-order bits that are lost during addition.
282pub trait NeumaierAddable: Sized {
283    /// Performs one step of the Neumaier compensated summation algorithm.
284    ///
285    /// This method adds `value` to the running `sum` while tracking lost precision
286    /// in the `compensation` term. The final result should be computed as
287    /// `sum + compensation` after all values have been added.
288    ///
289    /// # Arguments
290    ///
291    /// * `value` - The value to add to the sum.
292    /// * `sum` - A mutable reference to the running sum.
293    /// * `compensation` - A mutable reference to the compensation term that tracks lost precision.
294    fn neumaier_compensated_sum(value: Self, sum: &mut Self, compensation: &mut Self);
295}
296
297impl NeumaierAddable for f64 {
298    #[inline(always)]
299    fn neumaier_compensated_sum(value: Self, sum: &mut Self, compensation: &mut Self) {
300        neumaier_sum_and_compensation_real(value, sum, compensation)
301    }
302}
303
304impl NeumaierAddable for Complex<f64> {
305    #[inline(always)]
306    fn neumaier_compensated_sum(value: Self, sum: &mut Self, compensation: &mut Self) {
307        neumaier_sum_and_compensation_real(value.re, &mut sum.re, &mut compensation.re);
308        neumaier_sum_and_compensation_real(value.im, &mut sum.im, &mut compensation.im);
309    }
310}
311
312#[cfg(feature = "rug")]
313mod rug_impls {
314    use super::*;
315
316    #[inline(always)]
317    fn neumaier_sum_and_compensation_rug_float(
318        value: rug::Float,
319        sum: &mut rug::Float,
320        compensation: &mut rug::Float,
321    ) {
322        let sum_before_compensation = sum.clone();
323        *sum += &value;
324
325        // NOTE: the parenthesis are necessary in order for this algorithm to work correctly.
326        *compensation += if sum_before_compensation.clone().abs() >= value.clone().abs() {
327            // If sum is bigger, low-order digits of value are lost.
328            (sum_before_compensation - &*sum) + value
329        } else {
330            // Else low-order digits of sum are lost.
331            (value - &*sum) + sum_before_compensation
332        };
333    }
334
335    impl NeumaierAddable for rug::Float {
336        #[inline(always)]
337        fn neumaier_compensated_sum(value: Self, sum: &mut Self, compensation: &mut Self) {
338            neumaier_sum_and_compensation_rug_float(value, sum, compensation)
339        }
340    }
341
342    impl NeumaierAddable for rug::Complex {
343        #[inline(always)]
344        fn neumaier_compensated_sum(value: Self, sum: &mut Self, compensation: &mut Self) {
345            let (value_real, value_imag) = value.into_real_imag();
346
347            neumaier_sum_and_compensation_rug_float(
348                value_real,
349                sum.mut_real(),
350                compensation.mut_real(),
351            );
352
353            neumaier_sum_and_compensation_rug_float(
354                value_imag,
355                sum.mut_imag(),
356                compensation.mut_imag(),
357            );
358        }
359    }
360}
361
362//-----------------------------------------------------------------------------------
363
364//------------------------------------------------------------------------------------
365/// Neumaier compensated accumulator for sums of floating-point (real or complex) numbers.
366///
367/// Standard sequential summation loses precision when values of very different
368/// magnitudes are added: a small addend can be completely swallowed by a much
369/// larger running total.  `NeumaierSum` corrects for this by maintaining a
370/// *compensation* term that captures the low-order bits discarded at each step,
371/// then folds them back into the result at the end.
372///
373/// `NeumaierSum<T>` implements [`SumAccumulator`],
374/// so it can be used directly as the accumulation strategy in the
375/// [`norms`](crate::algorithms::vector_norms) vector-norm functions (e.g.,
376/// [`vector_norm_l2_neumaier`](crate::algorithms::vector_norms::vector_norm_l2_neumaier),
377/// [`vector_norm_l1_neumaier`](crate::algorithms::vector_norms::vector_norm_l1_neumaier)).
378///
379/// # Example
380///
381/// ```
382/// use num_valid::algorithms::accumulators::*;
383///
384/// let mut neumaier = NeumaierSum::<f64>::new();
385/// neumaier.push(1.0);
386/// neumaier.push(1e100);
387/// neumaier.push(1.0);
388/// neumaier.push(-1e100);
389///
390/// let sum = neumaier.result();
391/// println!("Sum: {}", sum);
392/// assert_eq!(sum, 2.0);
393/// ```
394///
395/// # References
396///
397/// * [Neumaier Summation](https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements)
398#[derive(Debug, Clone, Getters)]
399pub struct NeumaierSum<ScalarType> {
400    /// The sum before the compensation term is added.
401    #[getset(get = "pub")]
402    sum_before_compensation: ScalarType,
403
404    /// The compensation term.
405    /// This is the correction term that is added to the `sum_before_compensation` to correct for the loss of precision.
406    #[getset(get = "pub")]
407    compensation: ScalarType,
408}
409
410impl<ScalarType> Accumulator for NeumaierSum<ScalarType>
411where
412    ScalarType: Clone
413        + Zero
414        + for<'a> Add<&'a ScalarType, Output = ScalarType>
415        + for<'a> AddAssign<&'a ScalarType>
416        + for<'a> Mul<&'a ScalarType, Output = ScalarType>
417        + for<'a> MulAssign<&'a ScalarType>
418        + NeumaierAddable,
419{
420    type Input = ScalarType;
421    type Output = ScalarType;
422
423    #[inline(always)]
424    fn new() -> Self {
425        Self {
426            sum_before_compensation: ScalarType::zero(),
427            compensation: ScalarType::zero(),
428        }
429    }
430
431    #[inline(always)]
432    fn push(&mut self, value: Self::Input) {
433        <ScalarType as NeumaierAddable>::neumaier_compensated_sum(
434            value,
435            &mut self.sum_before_compensation,
436            &mut self.compensation,
437        );
438    }
439
440    #[inline(always)]
441    fn result(self) -> Self::Input {
442        self.sum_before_compensation + self.compensation
443    }
444
445    #[inline(always)]
446    fn combine(&mut self, other: Self) {
447        self.push(other.sum_before_compensation);
448        self.push(other.compensation);
449    }
450}
451
452impl<ScalarType> SumAccumulator for NeumaierSum<ScalarType>
453where
454    ScalarType: Clone
455        + Zero
456        + for<'a> Add<&'a ScalarType, Output = ScalarType>
457        + for<'a> AddAssign<&'a ScalarType>
458        + for<'a> Mul<&'a ScalarType, Output = ScalarType>
459        + for<'a> MulAssign<&'a ScalarType>
460        + NeumaierAddable,
461{
462    #[inline(always)]
463    fn rescale_by(&mut self, r: &Self::Input) {
464        // Compute the current true sum before rescaling.
465        self.sum_before_compensation += &self.compensation;
466
467        // Rescale the sum and reset the compensation term to zero.
468        self.sum_before_compensation *= r;
469        self.compensation = ScalarType::zero();
470    }
471}
472//------------------------------------------------------------------------------------
473
474//------------------------------------------------------------------------------------
475/// Accumulator that tracks the running maximum of real values.
476///
477/// Initialises to [`Constants::min_finite`](crate::Constants::min_finite) so that any finite input
478/// immediately becomes the maximum. For non-empty inputs this is always the
479/// true maximum; for empty inputs [`result`](Accumulator::result) returns
480/// `min_finite`.
481///
482/// # Note on empty input
483///
484/// For norm computations prefer [`MaxAbsValueAccumulator`], which initialises
485/// to zero and is therefore correct for empty vectors (norm = 0).
486pub struct MaxAccumulator<RealType: RealScalar> {
487    max_value: RealType,
488}
489
490impl<RealType: RealScalar> Accumulator for MaxAccumulator<RealType> {
491    type Input = RealType;
492    type Output = RealType;
493
494    #[inline(always)]
495    fn new() -> Self {
496        Self {
497            max_value: RealType::min_finite(),
498        }
499    }
500
501    #[inline(always)]
502    fn push(&mut self, value: Self::Input) {
503        if value > self.max_value {
504            self.max_value = value;
505        }
506    }
507
508    #[inline(always)]
509    fn result(self) -> Self::Output {
510        self.max_value
511    }
512
513    #[inline(always)]
514    fn combine(&mut self, other: Self) {
515        if other.max_value > self.max_value {
516            self.max_value = other.max_value;
517        }
518    }
519}
520
521/// Accumulator that tracks the running minimum of real values.
522///
523/// Initialises to [`Constants::max_finite`](crate::Constants::max_finite) so that any finite input
524/// immediately becomes the minimum. For empty inputs
525/// [`result`](Accumulator::result) returns `max_finite`.
526pub struct MinAccumulator<RealType: RealScalar> {
527    min_value: RealType,
528}
529
530impl<RealType: RealScalar> Accumulator for MinAccumulator<RealType> {
531    type Input = RealType;
532    type Output = RealType;
533
534    #[inline(always)]
535    fn new() -> Self {
536        Self {
537            min_value: RealType::max_finite(),
538        }
539    }
540
541    #[inline(always)]
542    fn push(&mut self, value: Self::Input) {
543        if value < self.min_value {
544            self.min_value = value;
545        }
546    }
547
548    #[inline(always)]
549    fn result(self) -> Self::Output {
550        self.min_value
551    }
552
553    #[inline(always)]
554    fn combine(&mut self, other: Self) {
555        if other.min_value < self.min_value {
556            self.min_value = other.min_value;
557        }
558    }
559}
560//------------------------------------------------------------------------------------
561
562//------------------------------------------------------------------------------------
563/// Accumulator that computes the maximum absolute value over a sequence of
564/// scalars: `maxᵢ |xᵢ|`.
565///
566/// Unlike [`MaxAccumulator`], this type:
567/// - accepts any [`FpScalar`] (real **or** complex) as input, applying `abs()`
568///   internally, so the modulus is used for complex scalars.
569/// - initialises to **zero** rather than `min_finite`, making it correct for
570///   empty inputs (result = 0, matching the L∞ norm of an empty vector).
571/// - returns a validated [`NonNegativeRealScalar`] from
572///   [`result`](Accumulator::result), encoding the non-negativity invariant
573///   at the type level.
574///
575/// This is the accumulator used internally by [`LinfNorm`](crate::algorithms::vector_norms::LinfNorm).
576pub struct MaxAbsValueAccumulator<ScalarType: FpScalar> {
577    max_abs: ScalarType::RealType, // internal: RealType for comparison of absolute values
578}
579
580impl<ScalarType: FpScalar> Accumulator for MaxAbsValueAccumulator<ScalarType> {
581    type Input = ScalarType;
582    type Output = NonNegativeRealScalar<ScalarType::RealType>;
583
584    fn new() -> Self {
585        Self {
586            max_abs: ScalarType::RealType::zero(),
587        }
588    }
589
590    fn push(&mut self, value: ScalarType) {
591        let abs_v = value.abs();
592        if abs_v > self.max_abs {
593            self.max_abs = abs_v;
594        }
595    }
596
597    fn combine(&mut self, other: Self) {
598        if other.max_abs > self.max_abs {
599            self.max_abs = other.max_abs;
600        }
601    }
602
603    fn result(self) -> NonNegativeRealScalar<ScalarType::RealType> {
604        NonNegativeRealScalar::try_new(self.max_abs)
605            .expect("MaxAbsValueAccumulator: max of absolute values is negative (bug)")
606    }
607}
608//------------------------------------------------------------------------------------
609
610//------------------------------------------------------------------------------------
611#[cfg(test)]
612mod tests_neumaier_sum {
613    use super::*;
614
615    mod native64 {
616        use super::*;
617
618        mod real {
619            use super::*;
620
621            #[test]
622            fn new() {
623                let neumaier = NeumaierSum::<f64>::new();
624                assert_eq!(neumaier.sum_before_compensation, 0.0);
625                assert_eq!(neumaier.compensation, 0.0);
626            }
627
628            #[test]
629            fn add() {
630                let mut neumaier = NeumaierSum::<f64>::new();
631                neumaier.push(1.0);
632                neumaier.push(1e-16);
633                neumaier.push(-1.0);
634                assert_eq!(neumaier.sum_before_compensation, 0.0);
635                assert_eq!(neumaier.compensation, 1e-16);
636            }
637
638            #[test]
639            fn sum() {
640                let mut neumaier = NeumaierSum::<f64>::new();
641                neumaier.push(1.0);
642                neumaier.push(1e-16);
643                neumaier.push(-1.0);
644                assert_eq!(neumaier.sum_before_compensation, 0.0);
645                assert_eq!(neumaier.compensation, 1e-16);
646                let sum = neumaier.result();
647                assert_eq!(sum, 1e-16);
648                println!("compensated sum = {}", sum);
649            }
650
651            #[test]
652            fn sum_big_values() {
653                let values = vec![1.0, 1e100, 1.0, -1e100];
654                let sum = values.iter().sum::<f64>();
655                assert_eq!(sum, 0.0);
656
657                let neumaier = NeumaierSum::<f64>::new_sequential(values);
658                let sum = neumaier.result();
659                assert_eq!(sum, 2.0);
660                println!("compensated sum = {}", sum);
661            }
662
663            #[test]
664            fn sum_small_values() {
665                let values = [1.0, 1e-100, -1.0];
666                let sum = values.iter().sum::<f64>();
667                assert_eq!(sum, 0.0);
668
669                let neumaier = NeumaierSum::<f64>::new_sequential(values);
670                let sum = neumaier.result();
671                assert_eq!(sum, 1e-100);
672                println!("compensated sum = {}", sum);
673            }
674
675            #[test]
676            fn combine_partial_sums() {
677                let mut left = NaiveSum::<f64>::new();
678                left.push(1.0);
679                left.push(2.0);
680
681                let mut right = NaiveSum::<f64>::new();
682                right.push(3.0);
683                right.push(4.0);
684
685                left.combine(right);
686                assert_eq!(left.result(), 10.0);
687            }
688
689            #[test]
690            fn combine_partial_neumaier_sums() {
691                let mut left = NeumaierSum::<f64>::new();
692                left.push(1.0);
693                left.push(1e-16);
694
695                let mut right = NeumaierSum::<f64>::new();
696                right.push(2.0);
697                right.push(1e-16);
698
699                left.combine(right);
700                let sum = left.result();
701                assert!((sum - (3.0 + 2e-16)).abs() < 1e-15);
702            }
703
704            #[test]
705            fn combine_matches_new_sequential_for_chunks() {
706                let values = [1.0_f64, 1e100, 1.0, -1e100, 1e-16, 2.0];
707
708                let mut chunk_1 = NaiveSum::<f64>::new();
709                chunk_1.push(values[0]);
710                chunk_1.push(values[1]);
711
712                let mut chunk_2 = NaiveSum::<f64>::new();
713                chunk_2.push(values[2]);
714                chunk_2.push(values[3]);
715
716                let mut chunk_3 = NaiveSum::<f64>::new();
717                chunk_3.push(values[4]);
718                chunk_3.push(values[5]);
719
720                chunk_1.combine(chunk_2);
721                chunk_1.combine(chunk_3);
722
723                let sequential = NaiveSum::<f64>::new_sequential(values);
724                assert_eq!(chunk_1.result(), sequential.result());
725            }
726
727            #[test]
728            fn combine_is_order_dependent_for_naive_sum() {
729                let mut left_grouped = NaiveSum::<f64>::new();
730                left_grouped.push(1e16);
731
732                let mut middle = NaiveSum::<f64>::new();
733                middle.push(1.0);
734
735                let mut right_last = NaiveSum::<f64>::new();
736                right_last.push(-1e16);
737
738                left_grouped.combine(middle);
739                left_grouped.combine(right_last);
740                let left_result = left_grouped.result();
741
742                let mut right_grouped = NaiveSum::<f64>::new();
743                right_grouped.push(1e16);
744
745                let mut right_last = NaiveSum::<f64>::new();
746                right_last.push(-1e16);
747
748                let mut middle = NaiveSum::<f64>::new();
749                middle.push(1.0);
750
751                right_grouped.combine(right_last);
752                right_grouped.combine(middle);
753                let right_result = right_grouped.result();
754
755                assert_eq!(left_result, 0.0);
756                assert_eq!(right_result, 1.0);
757                assert_ne!(left_result, right_result);
758            }
759
760            #[test]
761            fn neumaier_combine_matches_new_sequential_for_chunks() {
762                let values = [1.0_f64, 1e100, 1.0, -1e100, 1e-16, 2.0];
763
764                let mut chunk_1 = NeumaierSum::<f64>::new();
765                chunk_1.push(values[0]);
766                chunk_1.push(values[1]);
767
768                let mut chunk_2 = NeumaierSum::<f64>::new();
769                chunk_2.push(values[2]);
770                chunk_2.push(values[3]);
771
772                let mut chunk_3 = NeumaierSum::<f64>::new();
773                chunk_3.push(values[4]);
774                chunk_3.push(values[5]);
775
776                chunk_1.combine(chunk_2);
777                chunk_1.combine(chunk_3);
778
779                let sequential = NeumaierSum::<f64>::new_sequential(values);
780                assert_eq!(chunk_1.result(), sequential.result());
781            }
782
783            #[test]
784            fn neumaier_combine_is_stable_across_orders() {
785                let mut left_order = NeumaierSum::<f64>::new();
786                left_order.push(1e16);
787
788                let mut middle = NeumaierSum::<f64>::new();
789                middle.push(1.0);
790
791                let mut right_last = NeumaierSum::<f64>::new();
792                right_last.push(-1e16);
793
794                left_order.combine(middle);
795                left_order.combine(right_last);
796                let left_result = left_order.result();
797
798                let mut right_order = NeumaierSum::<f64>::new();
799                right_order.push(1e16);
800
801                let mut right_last = NeumaierSum::<f64>::new();
802                right_last.push(-1e16);
803
804                let mut middle = NeumaierSum::<f64>::new();
805                middle.push(1.0);
806
807                right_order.combine(right_last);
808                right_order.combine(middle);
809                let right_result = right_order.result();
810
811                assert_eq!(left_result, 1.0);
812                assert_eq!(right_result, 1.0);
813                assert_eq!(left_result, right_result);
814            }
815        }
816
817        mod complex {
818            use super::*;
819            use num::Complex;
820
821            #[test]
822            fn new() {
823                let neumaier = NeumaierSum::<Complex<f64>>::new();
824
825                let zero = Complex::new(0.0, 0.0);
826                assert_eq!(&neumaier.sum_before_compensation, &zero);
827                assert_eq!(&neumaier.compensation, &zero);
828            }
829
830            #[test]
831            fn add() {
832                let mut neumaier = NeumaierSum::<Complex<f64>>::new();
833
834                let zero = Complex::new(0.0, 0.0);
835                let v = Complex::new(1e-16, 2e-16);
836
837                neumaier.push(Complex::new(1.0, 2.0));
838                neumaier.push(v);
839                neumaier.push(Complex::new(-1.0, -2.0));
840
841                assert_eq!(neumaier.sum_before_compensation, zero);
842                assert_eq!(neumaier.compensation, v);
843            }
844
845            #[test]
846            fn sum() {
847                let mut neumaier = NeumaierSum::<Complex<f64>>::new();
848
849                let zero = Complex::new(0.0, 0.0);
850                let v = Complex::new(1e-16, 2e-16);
851
852                neumaier.push(Complex::new(1.0, 2.0));
853                neumaier.push(v);
854                neumaier.push(Complex::new(-1.0, -2.0));
855                assert_eq!(neumaier.sum_before_compensation, zero);
856                assert_eq!(neumaier.compensation, v);
857                let sum = neumaier.result();
858                assert_eq!(sum, v);
859                println!("compensated sum = {}", sum);
860            }
861
862            #[test]
863            fn sum_big_values() {
864                let values = vec![
865                    Complex::new(1.0, 2.0),
866                    Complex::new(1e100, 2e100),
867                    Complex::new(1.0, 2.0),
868                    Complex::new(-1e100, -2e100),
869                ];
870                let sum = values.iter().sum::<Complex<f64>>();
871                assert_eq!(sum, Complex::new(0.0, 0.0));
872
873                let neumaier = NeumaierSum::<Complex<f64>>::new_sequential(values);
874                let sum = neumaier.result();
875                assert_eq!(sum, Complex::new(2.0, 4.0));
876                println!("compensated sum = {}", sum);
877            }
878
879            #[test]
880            fn sum_small_values() {
881                let v = Complex::new(1e-100, 2e-100);
882
883                let values = [Complex::new(1.0, 2.0), v, Complex::new(-1.0, -2.0)];
884                let sum = values.iter().sum::<Complex<f64>>();
885                assert_eq!(sum, Complex::new(0.0, 0.0));
886
887                let neumaier = NeumaierSum::<Complex<f64>>::new_sequential(values);
888                let sum = neumaier.result();
889                assert_eq!(sum, v);
890                println!("compensated sum = {}", sum);
891            }
892
893            #[test]
894            fn combine_partial_sums() {
895                let mut left = NaiveSum::<Complex<f64>>::new();
896                left.push(Complex::new(1.0, 2.0));
897
898                let mut right = NaiveSum::<Complex<f64>>::new();
899                right.push(Complex::new(3.0, 4.0));
900
901                left.combine(right);
902                assert_eq!(left.result(), Complex::new(4.0, 6.0));
903            }
904
905            #[test]
906            fn combine_partial_neumaier_sums() {
907                let mut left = NeumaierSum::<Complex<f64>>::new();
908                left.push(Complex::new(1.0, 2.0));
909                left.push(Complex::new(1e-16, 2e-16));
910
911                let mut right = NeumaierSum::<Complex<f64>>::new();
912                right.push(Complex::new(3.0, 4.0));
913                right.push(Complex::new(1e-16, 2e-16));
914
915                left.combine(right);
916                assert_eq!(
917                    left.result(),
918                    Complex::new(4.0, 6.0) + Complex::new(2e-16, 4e-16)
919                );
920            }
921
922            #[test]
923            fn combine_matches_new_sequential_for_chunks() {
924                let values = [
925                    Complex::new(1.0, 2.0),
926                    Complex::new(1e100, 2e100),
927                    Complex::new(1.0, 2.0),
928                    Complex::new(-1e100, -2e100),
929                    Complex::new(1e-16, 2e-16),
930                    Complex::new(2.0, 3.0),
931                ];
932
933                let mut chunk_1 = NaiveSum::<Complex<f64>>::new();
934                chunk_1.push(values[0]);
935                chunk_1.push(values[1]);
936
937                let mut chunk_2 = NaiveSum::<Complex<f64>>::new();
938                chunk_2.push(values[2]);
939                chunk_2.push(values[3]);
940
941                let mut chunk_3 = NaiveSum::<Complex<f64>>::new();
942                chunk_3.push(values[4]);
943                chunk_3.push(values[5]);
944
945                chunk_1.combine(chunk_2);
946                chunk_1.combine(chunk_3);
947
948                let sequential = NaiveSum::<Complex<f64>>::new_sequential(values);
949                assert_eq!(chunk_1.result(), sequential.result());
950            }
951
952            #[test]
953            fn neumaier_combine_matches_new_sequential_for_chunks() {
954                let values = [
955                    Complex::new(1.0, 2.0),
956                    Complex::new(1e100, 2e100),
957                    Complex::new(1.0, 2.0),
958                    Complex::new(-1e100, -2e100),
959                    Complex::new(1e-16, 2e-16),
960                    Complex::new(2.0, 3.0),
961                ];
962
963                let mut chunk_1 = NeumaierSum::<Complex<f64>>::new();
964                chunk_1.push(values[0]);
965                chunk_1.push(values[1]);
966
967                let mut chunk_2 = NeumaierSum::<Complex<f64>>::new();
968                chunk_2.push(values[2]);
969                chunk_2.push(values[3]);
970
971                let mut chunk_3 = NeumaierSum::<Complex<f64>>::new();
972                chunk_3.push(values[4]);
973                chunk_3.push(values[5]);
974
975                chunk_1.combine(chunk_2);
976                chunk_1.combine(chunk_3);
977
978                let sequential = NeumaierSum::<Complex<f64>>::new_sequential(values);
979                assert_eq!(chunk_1.result(), sequential.result());
980            }
981        }
982    }
983
984    #[cfg(feature = "rug")]
985    mod rug100 {
986        use super::*;
987        use crate::{ComplexRugStrictFinite, RealRugStrictFinite};
988        use try_create::TryNew;
989
990        const PRECISION: u32 = 100;
991
992        mod real {
993            use super::*;
994            use crate::RealScalar;
995
996            #[test]
997            fn new() {
998                let neumaier = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
999                assert_eq!(neumaier.sum_before_compensation, 0.0);
1000                assert_eq!(neumaier.compensation, 0.0);
1001            }
1002
1003            #[test]
1004            fn add() {
1005                let mut neumaier = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
1006
1007                let v = RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
1008                    PRECISION,
1009                    rug::Float::parse("1e-100").unwrap(),
1010                ))
1011                .expect("valid test value");
1012
1013                neumaier.push(RealRugStrictFinite::<PRECISION>::try_from_f64(1.0).unwrap());
1014                neumaier.push(v.clone());
1015                neumaier.push(RealRugStrictFinite::<PRECISION>::try_from_f64(-1.0).unwrap());
1016
1017                assert_eq!(
1018                    neumaier.sum_before_compensation,
1019                    RealRugStrictFinite::<PRECISION>::try_from_f64(0.0).unwrap()
1020                );
1021                assert_eq!(&neumaier.compensation, &v);
1022            }
1023
1024            #[test]
1025            fn sum() {
1026                let mut neumaier = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
1027
1028                let v = RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
1029                    PRECISION,
1030                    rug::Float::parse("1e-100").unwrap(),
1031                ))
1032                .expect("valid test value");
1033
1034                neumaier.push(RealRugStrictFinite::<PRECISION>::try_from_f64(1.0).unwrap());
1035                neumaier.push(v.clone());
1036                neumaier.push(RealRugStrictFinite::<PRECISION>::try_from_f64(-1.0).unwrap());
1037
1038                assert_eq!(neumaier.sum_before_compensation, 0.0);
1039                assert_eq!(&neumaier.compensation, &v);
1040                let sum = neumaier.result();
1041                assert_eq!(sum, v);
1042                println!("compensated sum = {}", sum);
1043            }
1044
1045            #[test]
1046            fn sum_big_values() {
1047                let values = ["1.0", "1e100", "1.0", "-1e100"]
1048                    .iter()
1049                    .map(|v| {
1050                        RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
1051                            PRECISION,
1052                            rug::Float::parse(v).unwrap(),
1053                        ))
1054                        .expect("valid test value")
1055                    })
1056                    .collect::<Vec<_>>();
1057
1058                let sum = values
1059                    .iter()
1060                    .fold(RealRugStrictFinite::<PRECISION>::zero(), |acc, x| acc + x);
1061                assert_eq!(sum, 0.0);
1062
1063                let neumaier =
1064                    NeumaierSum::<RealRugStrictFinite<PRECISION>>::new_sequential(values);
1065                let sum = neumaier.result();
1066                assert_eq!(
1067                    sum,
1068                    RealRugStrictFinite::<PRECISION>::try_from_f64(2.0).unwrap()
1069                );
1070                println!("compensated sum = {}", sum);
1071            }
1072
1073            #[test]
1074            fn sum_small_values() {
1075                let values = ["1.0", "1e-100", "-1.0"]
1076                    .iter()
1077                    .map(|v| {
1078                        RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
1079                            PRECISION,
1080                            rug::Float::parse(v).unwrap(),
1081                        ))
1082                        .expect("valid test value")
1083                    })
1084                    .collect::<Vec<_>>();
1085
1086                let sum = values
1087                    .iter()
1088                    .fold(RealRugStrictFinite::<PRECISION>::zero(), |acc, x| acc + x);
1089                assert_eq!(sum, RealRugStrictFinite::<PRECISION>::zero());
1090
1091                let neumaier =
1092                    NeumaierSum::<RealRugStrictFinite<PRECISION>>::new_sequential(values);
1093                let sum = neumaier.result();
1094                assert_eq!(
1095                    sum,
1096                    RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
1097                        PRECISION,
1098                        rug::Float::parse("1e-100").unwrap(),
1099                    ))
1100                    .expect("valid test value")
1101                );
1102                println!("compensated sum = {}", sum);
1103            }
1104
1105            #[test]
1106            fn combine_matches_new_sequential_for_chunks() {
1107                let values = ["1.0", "1e100", "1.0", "-1e100", "1e-100", "2.0"]
1108                    .iter()
1109                    .map(|v| {
1110                        RealRugStrictFinite::<PRECISION>::try_new(rug::Float::with_val(
1111                            PRECISION,
1112                            rug::Float::parse(v).unwrap(),
1113                        ))
1114                        .expect("valid test value")
1115                    })
1116                    .collect::<Vec<_>>();
1117
1118                let mut chunk_1 = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
1119                chunk_1.push(values[0].clone());
1120                chunk_1.push(values[1].clone());
1121
1122                let mut chunk_2 = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
1123                chunk_2.push(values[2].clone());
1124                chunk_2.push(values[3].clone());
1125
1126                let mut chunk_3 = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
1127                chunk_3.push(values[4].clone());
1128                chunk_3.push(values[5].clone());
1129
1130                chunk_1.combine(chunk_2);
1131                chunk_1.combine(chunk_3);
1132
1133                let sequential =
1134                    NeumaierSum::<RealRugStrictFinite<PRECISION>>::new_sequential(values);
1135                assert_eq!(chunk_1.result(), sequential.result());
1136            }
1137
1138            #[test]
1139            fn neumaier_combine_is_stable_across_orders() {
1140                let mut left_order = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
1141                left_order.push(RealRugStrictFinite::<PRECISION>::try_from_f64(1e16).unwrap());
1142
1143                let mut middle = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
1144                middle.push(RealRugStrictFinite::<PRECISION>::try_from_f64(1.0).unwrap());
1145
1146                let mut right_last = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
1147                right_last.push(RealRugStrictFinite::<PRECISION>::try_from_f64(-1e16).unwrap());
1148
1149                left_order.combine(middle);
1150                left_order.combine(right_last);
1151                let left_result = left_order.result();
1152
1153                let mut right_order = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
1154                right_order.push(RealRugStrictFinite::<PRECISION>::try_from_f64(1e16).unwrap());
1155
1156                let mut right_last = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
1157                right_last.push(RealRugStrictFinite::<PRECISION>::try_from_f64(-1e16).unwrap());
1158
1159                let mut middle = NeumaierSum::<RealRugStrictFinite<PRECISION>>::new();
1160                middle.push(RealRugStrictFinite::<PRECISION>::try_from_f64(1.0).unwrap());
1161
1162                right_order.combine(right_last);
1163                right_order.combine(middle);
1164                let right_result = right_order.result();
1165
1166                assert_eq!(
1167                    left_result,
1168                    RealRugStrictFinite::<PRECISION>::try_from_f64(1.0).unwrap()
1169                );
1170                assert_eq!(
1171                    right_result,
1172                    RealRugStrictFinite::<PRECISION>::try_from_f64(1.0).unwrap()
1173                );
1174                assert_eq!(left_result, right_result);
1175            }
1176        }
1177
1178        mod complex {
1179            use super::*;
1180
1181            #[test]
1182            fn new() {
1183                let neumaier = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new();
1184                assert_eq!(
1185                    neumaier.sum_before_compensation,
1186                    ComplexRugStrictFinite::<PRECISION>::zero()
1187                );
1188                assert_eq!(
1189                    neumaier.compensation,
1190                    ComplexRugStrictFinite::<PRECISION>::zero()
1191                );
1192            }
1193
1194            #[test]
1195            fn add() {
1196                let mut neumaier = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new();
1197
1198                let v = ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
1199                    PRECISION,
1200                    rug::Complex::parse("(1e-100,2e-100)").unwrap(),
1201                ))
1202                .expect("valid test value");
1203
1204                neumaier.push(
1205                    ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(1.0, 2.0)).unwrap(),
1206                );
1207                neumaier.push(v.clone());
1208                neumaier.push(
1209                    ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(-1.0, -2.0))
1210                        .unwrap(),
1211                );
1212
1213                assert_eq!(
1214                    neumaier.sum_before_compensation,
1215                    ComplexRugStrictFinite::<PRECISION>::zero()
1216                );
1217                assert_eq!(&neumaier.compensation, &v);
1218            }
1219
1220            #[test]
1221            fn sum() {
1222                let mut neumaier = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new();
1223
1224                let v = ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
1225                    PRECISION,
1226                    rug::Complex::parse("(1e-100,2e-100)").unwrap(),
1227                ))
1228                .expect("valid test value");
1229
1230                neumaier.push(
1231                    ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(1.0, 2.0)).unwrap(),
1232                );
1233                neumaier.push(v.clone());
1234                neumaier.push(
1235                    ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(-1.0, -2.0))
1236                        .unwrap(),
1237                );
1238
1239                assert_eq!(
1240                    neumaier.sum_before_compensation,
1241                    ComplexRugStrictFinite::<PRECISION>::zero()
1242                );
1243                assert_eq!(&neumaier.compensation, &v);
1244                let sum = neumaier.result();
1245                assert_eq!(sum, v);
1246                println!("compensated sum = {}", sum);
1247            }
1248
1249            #[test]
1250            fn sum_big_values() {
1251                let values = ["(1.0,2.0)", "(1e100,2e100)", "(1.0,2.0)", "(-1e100,-2e100)"]
1252                    .iter()
1253                    .map(|v| {
1254                        ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
1255                            PRECISION,
1256                            rug::Complex::parse(v).unwrap(),
1257                        ))
1258                        .expect("valid test value")
1259                    })
1260                    .collect::<Vec<_>>();
1261
1262                let zero = ComplexRugStrictFinite::<PRECISION>::zero();
1263                let sum = values.iter().fold(zero.clone(), |acc, x| acc + x);
1264                assert_eq!(sum, zero);
1265
1266                let neumaier =
1267                    NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new_sequential(values);
1268                let sum = neumaier.result();
1269                assert_eq!(
1270                    sum,
1271                    ComplexRugStrictFinite::<PRECISION>::try_from(Complex::new(2.0, 4.0)).unwrap()
1272                );
1273                println!("compensated sum = {}", sum);
1274            }
1275
1276            #[test]
1277            fn sum_small_values() {
1278                let values = ["(1.0,2.0)", "(1e-100,2e-100)", "(-1.0,-2.0)"]
1279                    .iter()
1280                    .map(|v| {
1281                        ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
1282                            PRECISION,
1283                            rug::Complex::parse(v).unwrap(),
1284                        ))
1285                        .expect("valid test value")
1286                    })
1287                    .collect::<Vec<_>>();
1288
1289                let zero = ComplexRugStrictFinite::<PRECISION>::zero();
1290                let sum = values.iter().fold(zero.clone(), |acc, x| acc + x);
1291                assert_eq!(sum, zero);
1292
1293                let neumaier =
1294                    NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new_sequential(values);
1295                let sum = neumaier.result();
1296                assert_eq!(
1297                    sum,
1298                    ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
1299                        PRECISION,
1300                        rug::Complex::parse("(1e-100,2e-100)").unwrap(),
1301                    ))
1302                    .expect("valid test value")
1303                );
1304                println!("compensated sum = {}", sum);
1305            }
1306
1307            #[test]
1308            fn combine_matches_new_sequential_for_chunks() {
1309                let values = [
1310                    "(1.0,2.0)",
1311                    "(1e100,2e100)",
1312                    "(1.0,2.0)",
1313                    "(-1e100,-2e100)",
1314                    "(1e-100,2e-100)",
1315                    "(2.0,3.0)",
1316                ]
1317                .iter()
1318                .map(|v| {
1319                    ComplexRugStrictFinite::<PRECISION>::try_new(rug::Complex::with_val(
1320                        PRECISION,
1321                        rug::Complex::parse(v).unwrap(),
1322                    ))
1323                    .expect("valid test value")
1324                })
1325                .collect::<Vec<_>>();
1326
1327                let mut chunk_1 = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new();
1328                chunk_1.push(values[0].clone());
1329                chunk_1.push(values[1].clone());
1330
1331                let mut chunk_2 = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new();
1332                chunk_2.push(values[2].clone());
1333                chunk_2.push(values[3].clone());
1334
1335                let mut chunk_3 = NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new();
1336                chunk_3.push(values[4].clone());
1337                chunk_3.push(values[5].clone());
1338
1339                chunk_1.combine(chunk_2);
1340                chunk_1.combine(chunk_3);
1341
1342                let sequential =
1343                    NeumaierSum::<ComplexRugStrictFinite<PRECISION>>::new_sequential(values);
1344                assert_eq!(chunk_1.result(), sequential.result());
1345            }
1346        }
1347    }
1348}
1349//------------------------------------------------------------------------------------
num_valid/algorithms/accumulators.rs

num_valid/algorithms/
accumulators.rs
