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
|
[/
Copyright Nick Thompson 2020
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or copy at
http://www.boost.org/LICENSE_1_0.txt).
]
[section:cohen_acceleration Cohen Acceleration]
[h4 Synopsis]
#include <boost/math/tools/cohen_acceleration.hpp>
namespace boost::math::tools {
template<class G>
auto cohen_acceleration(G& generator, int64_t n = -1);
} // namespaces
The function `cohen_acceleration` rapidly computes the limiting value of an alternating series via a technique developed by [@https://www.johndcook.com/blog/2020/08/06/cohen-acceleration/ Henri Cohen et al].
To compute
[$../equations/alternating_series.svg]
we first define a callable that produces /a/[sub /k/] on the kth call.
For example, suppose we wish to compute
[$../equations/zeta_related_alternating.svg]
First, we need to define a callable which returns the requisite terms:
template<typename Real>
class G {
public:
G(){
k_ = 0;
}
Real operator()() {
k_ += 1;
return 1/(k_*k_);
}
private:
Real k_;
};
Then we pass this into the `cohen_acceleration` function:
auto gen = G<double>();
double computed = cohen_acceleration(gen);
See `cohen_acceleration.cpp` in the `examples` directory for more.
The number of terms consumed is computed from the error model
[$../equations/cohen_acceleration_error_model.svg]
and must be computed /a priori/.
If we read the reference carefully, we notice that this error model is derived under the assumption that the terms /a/[sub /k/] are given as the moments of a positive measure on [0,1].
If this assumption does not hold, then the number of terms chosen by the method is incorrect.
Hence we permit the user to provide a second argument to specify the number of terms:
double computed = cohen_acceleration(gen, 5);
/Nota bene/: When experimenting with this option, we found that adding more terms was no guarantee of additional accuracy, and could not find an example where a user-provided number of terms outperformed the default.
In addition, it is easy to generate intermediates which overflow if we let /n/ grow too large.
Hence we recommend only playing with this parameter to /decrease/ the default number of terms to increase speed.
[heading Performance]
To see that Cohen acceleration is in fact faster than naive summation for the same level of relative accuracy, we can run the `reporting/performance/cohen_acceleration_performance.cpp` file.
This benchmark computes the alternating Basel series discussed above:
```
Running ./reporting/performance/cohen_acceleration_performance.x
Run on (16 X 2300 MHz CPU s)
CPU Caches:
L1 Data 32 KiB (x8)
L1 Instruction 32 KiB (x8)
L2 Unified 256 KiB (x8)
L3 Unified 16384 KiB (x1)
Load Average: 4.13, 3.71, 3.30
-----------------------------------------------------------------
Benchmark Time
-----------------------------------------------------------------
CohenAcceleration<float> 20.7 ns
CohenAcceleration<double> 64.6 ns
CohenAcceleration<long double> 115 ns
NaiveSum<float> 4994 ns
NaiveSum<double> 112803698 ns
NaiveSum<long double> 5009564877 ns
```
In fact not only does the naive sum take orders of magnitude longer to compute, it is less accurate as well.
[heading References]
* Cohen, Henri, Fernando Rodriguez Villegas, and Don Zagier. ['Convergence acceleration of alternating series.] Experimental mathematics 9.1 (2000): 3-12.
[endsect]
|
