math/tools/cohen_acceleration.hpp

0001 //  (C) Copyright Nick Thompson 2020.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_TOOLS_COHEN_ACCELERATION_HPP
0007 #define BOOST_MATH_TOOLS_COHEN_ACCELERATION_HPP
0008 #include <limits>
0009 #include <cmath>
0010 #include <cstdint>
0011
0012 namespace boost::math::tools {
0013
0014 // Algorithm 1 of https://people.mpim-bonn.mpg.de/zagier/files/exp-math-9/fulltext.pdf
0015 // Convergence Acceleration of Alternating Series: Henri Cohen, Fernando Rodriguez Villegas, and Don Zagier
0016 template<class G>
0017 auto cohen_acceleration(G& generator, std::int64_t n = -1)
0018 {
0019     using Real = decltype(generator());
0020     // This test doesn't pass for float128, sad!
0021     //static_assert(std::is_floating_point_v<Real>, "Real must be a floating point type.");
0022     using std::log;
0023     using std::pow;
0024     using std::ceil;
0025     using std::sqrt;
0026
0027     auto n_ = static_cast<Real>(n);
0028     if (n < 0)
0029     {
0030         // 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).
0031         // 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.)
0032         n_ = static_cast<Real>(ceil(log(Real(4)/std::numeric_limits<Real>::epsilon())*Real(0.5672963285532555)));
0033         n = static_cast<std::int64_t>(n_);
0034     }
0035     // d can get huge and overflow if you pick n too large:
0036     auto d = static_cast<Real>(pow(Real(3 + sqrt(Real(8))), n_));
0037     d = (d + Real(1)/d)/2;
0038     Real b = -1;
0039     Real c = -d;
0040     Real s = 0;
0041     for (Real k = 0; k < n_; ++k) {
0042         c = b - c;
0043         s += c*generator();
0044         b = (k+n_)*(k-n_)*b/((k+Real(1)/Real(2))*(k+1));
0045     }
0046
0047     return s/d;
0048 }
0049
0050 }
0051 #endif