52 lines · plain
1// (C) Copyright Nick Thompson 2020.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_TOOLS_COHEN_ACCELERATION_HPP7#define BOOST_MATH_TOOLS_COHEN_ACCELERATION_HPP8#include <limits>9#include <cmath>10#include <cstdint>11 12namespace boost::math::tools {13 14// Algorithm 1 of https://people.mpim-bonn.mpg.de/zagier/files/exp-math-9/fulltext.pdf15// Convergence Acceleration of Alternating Series: Henri Cohen, Fernando Rodriguez Villegas, and Don Zagier16template<class G>17auto cohen_acceleration(G& generator, std::int64_t n = -1)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 27 auto n_ = static_cast<Real>(n);28 if (n < 0)29 {30 // 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).31 // 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.)32 n_ = static_cast<Real>(ceil(log(Real(4)/std::numeric_limits<Real>::epsilon())*Real(0.5672963285532555)));33 n = static_cast<std::int64_t>(n_);34 }35 // d can get huge and overflow if you pick n too large:36 auto d = static_cast<Real>(pow(Real(3 + sqrt(Real(8))), n_));37 d = (d + Real(1)/d)/2;38 Real b = -1;39 Real c = -d;40 Real s = 0;41 for (Real k = 0; k < n_; ++k) {42 c = b - c;43 s += c*generator();44 b = (k+n_)*(k-n_)*b/((k+Real(1)/Real(2))*(k+1));45 }46 47 return s/d;48}49 50}51#endif52