rustics 1.0.3

simple statistic library for performance analysis
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130

//
//  Copyright 2024 Jonathan L Bertoni
//
//  This code is available under the Berkeley 2-Clause, Berkeley 3-clause,
//  and MIT licenses.
//

// Sort a set by absolute value to try to improve the accuracy of summation.

fn sort(input: &mut [f64]) {
    input.sort_by(|a, b| a.abs().partial_cmp(&b.abs()).unwrap())
}

/// The kbk_sum_sort() function performs a Kahan-Babushka-Klein summation.
/// It has a sort on the inputs in addition to the normal algorithm.
///
/// See the Wikipedia page on Kahan summation for the details of this
/// algorithm.

pub fn kbk_sum_sort(input: &mut [f64]) -> f64 {
    sort(input);
    kbk_sum(input)
}

/// The kbk_sum() function performs a Kahan-Babushka-Klein summation.  It does
/// not sort its inputs.  Currently, it is used by code that has presorted
/// vectors.

pub fn kbk_sum(input: &[f64]) -> f64 {
    let mut sum = 0.0;
    let mut cs  = 0.0;
    let mut ccs = 0.0;

    for addend in input {
        let     addend = *addend;
        let mut t      = sum + addend;

        let c =
            if sum.abs() >=  addend.abs() {
                (sum - t) +  addend
            } else {
                (addend - t) + sum
            };

        sum = t;
        t   = cs + c;

        let cc =
            if cs.abs() >= c.abs() {
                (cs - t) + c
            } else {
                (c - t) + cs
            };

        cs   = t;
        ccs += cc;
    }

    sum + cs + ccs
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn run_tests() {
        let     limit  = 16;
        let mut inputs = Vec::<f64>::new();

        for i in 0..=limit {
            inputs.push(i as f64);
        }

        let result = kbk_sum(&inputs);

        // Compute the expected value of the sum

        let expected = ((limit + 1) * limit) / 2;

        // Now check what we got.

        assert!(result == expected as f64);

        // The version with the sort should match...

        let result = kbk_sum_sort(&mut inputs);
        assert!(result == expected as f64);

        println!("vector sum 1:  {}", result);

        // Now add some negative values to the vector.

        for i in 0..=limit {
            inputs.push(-i as f64);
        }

        let result = kbk_sum(&inputs);

        println!("vector sum 2:  {}", result);

        assert!(result == 0.0);

        // Use the sort...

        let result = kbk_sum_sort(&mut inputs);

        assert!(result == 0.0);

        // Run the test example from Wikipedia:
        //   [ 1, large, -1, large ]

        let large = (10.0 as f64).powi(100);

        inputs.clear();
        inputs.push(1.0);
        inputs.push(large);
        inputs.push(1.0);
        inputs.push(-large);

        let result = kbk_sum(&inputs);

        assert!(result == 2.0);

        let result = kbk_sum_sort(&mut inputs);

        println!("vector sum 3:  {}", result);
        assert!(result == 2.0);
    }
}