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
//! Shared, auditable floating-point reduction policies.
/// Explicit addition policy for Tensor `sum` and `cumsum`.
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub enum SumMode {
/// Sequential.
Naive,
/// Balanced tree.
Pairwise,
/// Neumaier compensated.
Neumaier,
}
/// Reduces a slice with the named policy.
pub fn sum_f64(v: &[f64], m: SumMode) -> f64 {
match m {
SumMode::Naive => v.iter().sum(),
SumMode::Pairwise => {
if v.len() < 2 {
v.first().copied().unwrap_or(0.)
} else {
let n = v.len() / 2;
sum_f64(&v[..n], m) + sum_f64(&v[n..], m)
}
}
SumMode::Neumaier => {
let (mut s, mut c) = (0., 0.);
for &x in v {
let t = s + x;
c += if s.abs() >= x.abs() {
(s - t) + x
} else {
(x - t) + s
};
s = t
}
s + c
}
}
}
/// Produces prefix sums with the named policy.
pub fn cumsum_f64(v: &[f64], m: SumMode) -> Vec<f64> {
match m {
SumMode::Naive => {
let mut s = 0.;
v.iter()
.map(|&x| {
s += x;
s
})
.collect()
}
SumMode::Pairwise => (1..=v.len()).map(|n| sum_f64(&v[..n], m)).collect(),
SumMode::Neumaier => {
let (mut s, mut c) = (0., 0.);
v.iter()
.map(|&x| {
let t = s + x;
c += if s.abs() >= x.abs() {
(s - t) + x
} else {
(x - t) + s
};
s = t;
s + c
})
.collect()
}
}
}