// (C) Copyright Nick Thompson 2020.
|
// Use, modification and distribution are subject to 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)
|
|
#ifndef BOOST_MATH_TOOLS_COHEN_ACCELERATION_HPP
|
#define BOOST_MATH_TOOLS_COHEN_ACCELERATION_HPP
|
#include <limits>
|
#include <cmath>
|
|
namespace boost::math::tools {
|
|
// Algorithm 1 of https://people.mpim-bonn.mpg.de/zagier/files/exp-math-9/fulltext.pdf
|
// Convergence Acceleration of Alternating Series: Henri Cohen, Fernando Rodriguez Villegas, and Don Zagier
|
template<class G>
|
auto cohen_acceleration(G& generator, int64_t n = -1)
|
{
|
using Real = decltype(generator());
|
// This test doesn't pass for float128, sad!
|
//static_assert(std::is_floating_point_v<Real>, "Real must be a floating point type.");
|
using std::log;
|
using std::pow;
|
using std::ceil;
|
using std::sqrt;
|
Real n_ = n;
|
if (n < 0)
|
{
|
// 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).
|
// 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.)
|
n_ = ceil(log(4/std::numeric_limits<Real>::epsilon())*0.5672963285532555);
|
n = static_cast<int64_t>(n_);
|
}
|
// d can get huge and overflow if you pick n too large:
|
Real d = pow(3 + sqrt(Real(8)), n);
|
d = (d + 1/d)/2;
|
Real b = -1;
|
Real c = -d;
|
Real s = 0;
|
for (Real k = 0; k < n_; ++k) {
|
c = b - c;
|
s += c*generator();
|
b = (k+n_)*(k-n_)*b/((k+Real(1)/Real(2))*(k+1));
|
}
|
|
return s/d;
|
}
|
|
}
|
#endif
|
