]>
Commit | Line | Data |
---|---|---|
20effc67 TL |
1 | // (C) Copyright Nick Thompson 2020. |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. (See accompanying file | |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
5 | ||
6 | #ifndef BOOST_MATH_TOOLS_COHEN_ACCELERATION_HPP | |
7 | #define BOOST_MATH_TOOLS_COHEN_ACCELERATION_HPP | |
8 | #include <limits> | |
9 | #include <cmath> | |
1e59de90 | 10 | #include <cstdint> |
20effc67 TL |
11 | |
12 | namespace boost::math::tools { | |
13 | ||
14 | // Algorithm 1 of https://people.mpim-bonn.mpg.de/zagier/files/exp-math-9/fulltext.pdf | |
15 | // Convergence Acceleration of Alternating Series: Henri Cohen, Fernando Rodriguez Villegas, and Don Zagier | |
16 | template<class G> | |
1e59de90 | 17 | auto cohen_acceleration(G& generator, std::int64_t n = -1) |
20effc67 TL |
18 | { |
19 | using Real = decltype(generator()); | |
20 | // This test doesn't pass for float128, sad! | |
21 | //static_assert(std::is_floating_point_v<Real>, "Real must be a floating point type."); | |
22 | using std::log; | |
23 | using std::pow; | |
24 | using std::ceil; | |
25 | using std::sqrt; | |
26 | Real n_ = n; | |
27 | if (n < 0) | |
28 | { | |
29 | // relative error grows as 2*5.828^-n; take 5.828^-n < eps/4 => -nln(5.828) < ln(eps/4) => n > ln(4/eps)/ln(5.828). | |
30 | // Is there a way to do it rapidly with std::log2? (Yes, of course; but for primitive types it's computed at compile-time anyway.) | |
31 | n_ = ceil(log(4/std::numeric_limits<Real>::epsilon())*0.5672963285532555); | |
1e59de90 | 32 | n = static_cast<std::int64_t>(n_); |
20effc67 TL |
33 | } |
34 | // d can get huge and overflow if you pick n too large: | |
35 | Real d = pow(3 + sqrt(Real(8)), n); | |
36 | d = (d + 1/d)/2; | |
37 | Real b = -1; | |
38 | Real c = -d; | |
39 | Real s = 0; | |
40 | for (Real k = 0; k < n_; ++k) { | |
41 | c = b - c; | |
42 | s += c*generator(); | |
43 | b = (k+n_)*(k-n_)*b/((k+Real(1)/Real(2))*(k+1)); | |
44 | } | |
45 | ||
46 | return s/d; | |
47 | } | |
48 | ||
49 | } | |
50 | #endif |