rustics/
sum.rs

1//
2//  Copyright 2024 Jonathan L Bertoni
3//
4//  This code is available under the Berkeley 2-Clause, Berkeley 3-clause,
5//  and MIT licenses.
6//
7
8// Sort a set by absolute value to try to improve the accuracy of summation.
9
10fn sort(input: &mut [f64]) {
11    input.sort_by(|a, b| a.abs().partial_cmp(&b.abs()).unwrap())
12}
13
14/// The kbk_sum_sort() function performs a Kahan-Babushka-Klein summation.
15/// It has a sort on the inputs in addition to the normal algorithm.
16///
17/// See the Wikipedia page on Kahan summation for the details of this
18/// algorithm.
19
20pub fn kbk_sum_sort(input: &mut [f64]) -> f64 {
21    sort(input);
22    kbk_sum(input)
23}
24
25/// The kbk_sum() function performs a Kahan-Babushka-Klein summation.  It does
26/// not sort its inputs.  Currently, it is used by code that has presorted
27/// vectors.
28
29pub fn kbk_sum(input: &[f64]) -> f64 {
30    let mut sum = 0.0;
31    let mut cs  = 0.0;
32    let mut ccs = 0.0;
33
34    for addend in input {
35        let     addend = *addend;
36        let mut t      = sum + addend;
37
38        let c =
39            if sum.abs() >=  addend.abs() {
40                (sum - t) +  addend
41            } else {
42                (addend - t) + sum
43            };
44
45        sum = t;
46        t   = cs + c;
47
48        let cc =
49            if cs.abs() >= c.abs() {
50                (cs - t) + c
51            } else {
52                (c - t) + cs
53            };
54
55        cs   = t;
56        ccs += cc;
57    }
58
59    sum + cs + ccs
60}
61
62#[cfg(test)]
63mod tests {
64    use super::*;
65
66    #[test]
67    fn run_tests() {
68        let     limit  = 16;
69        let mut inputs = Vec::<f64>::new();
70
71        for i in 0..=limit {
72            inputs.push(i as f64);
73        }
74
75        let result = kbk_sum(&inputs);
76
77        // Compute the expected value of the sum
78
79        let expected = ((limit + 1) * limit) / 2;
80
81        // Now check what we got.
82
83        assert!(result == expected as f64);
84
85        // The version with the sort should match...
86
87        let result = kbk_sum_sort(&mut inputs);
88        assert!(result == expected as f64);
89
90        println!("vector sum 1:  {}", result);
91
92        // Now add some negative values to the vector.
93
94        for i in 0..=limit {
95            inputs.push(-i as f64);
96        }
97
98        let result = kbk_sum(&inputs);
99
100        println!("vector sum 2:  {}", result);
101
102        assert!(result == 0.0);
103
104        // Use the sort...
105
106        let result = kbk_sum_sort(&mut inputs);
107
108        assert!(result == 0.0);
109
110        // Run the test example from Wikipedia:
111        //   [ 1, large, -1, large ]
112
113        let large = (10.0 as f64).powi(100);
114
115        inputs.clear();
116        inputs.push(1.0);
117        inputs.push(large);
118        inputs.push(1.0);
119        inputs.push(-large);
120
121        let result = kbk_sum(&inputs);
122
123        assert!(result == 2.0);
124
125        let result = kbk_sum_sort(&mut inputs);
126
127        println!("vector sum 3:  {}", result);
128        assert!(result == 2.0);
129    }
130}
rustics/sum.rs

rustics/
sum.rs
