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
// use core::iter::Map;
use itertools::Itertools;
use num::Zero;
// use itertools::TupleWindows;
/// A trait for sequences that can be checked for convergence.
pub trait ConvergentSequence: Iterator<Item = f64> + Sized {
/// Applies Aitken's Δ² process to accelerate the convergence of a sequence.
/// See <https://en.wikipedia.org/wiki/Aitken%27s_delta-squared_process> and
/// <https://en.wikipedia.org/wiki/Shanks_transformation>
///
/// # Returns
///
/// An iterator over the accelerated sequence.
fn aitken(self) -> impl Iterator<Item = f64> {
self.tuple_windows::<(_, _, _)>().filter_map(|(x, y, z)| {
let dx = z - y;
let dx2 = y - x;
let ddx = dx - dx2;
// We can't handle a segment like [2,4,6]
// But e.g. [2, 2, 2] may have already converged
if ddx.is_zero() {
if dx.is_zero() { Some(z) } else { None }
} else {
Some(z - dx.powi(2) / ddx)
}
})
}
/// Finds the limit of the sequence within a given tolerance using Aitken's
/// Δ² process. This should *only* be applied to sequences that are known to
/// converge.
///
/// # Arguments
///
/// * `tol` - The tolerance within which to find the limit.
///
/// # Returns
///
/// The limit of the sequence as a floating-point number.
///
/// # Panics
///
/// Runs forever if the sequence does not converge within the given
/// tolerance.
fn limit(self, tol: f64) -> f64 {
self.aitken()
.aitken()
.aitken()
.aitken()
.tuple_windows::<(_, _)>()
.filter_map(
|(a, b)| {
if (a - b).abs() < tol { Some(b) } else { None }
},
)
.next()
.unwrap()
}
}
impl<T> ConvergentSequence for T where T: Iterator<Item = f64> + Sized {}
#[cfg(test)]
mod tests {
use super::*;
use num::Integer;
#[test]
fn test_aitken_limit() {
let seq = (0..)
.map(|n| {
let sign = if n.is_even() { 1.0 } else { -1.0 };
let val = sign / f64::from(2 * n + 1);
dbg!(val);
val
})
.scan(0.0, |acc, x| {
*acc += x;
Some(*acc)
});
let limit = seq.limit(1e-10);
let pi_over_4 = std::f64::consts::PI / 4.0;
assert!(
(limit - pi_over_4).abs() < 1e-10,
"The limit calculated using Aitken's Δ² process did not converge to π/4 within the tolerance."
);
}
}