boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Fix stack overflow in cohen acceleration from GCC bug

Browse Source
This commit is contained in:
Matt Borland
2023-05-22 13:48:33 +02:00
parent ca31dcfe93
commit faaf4759ac
2 changed files with 23 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
include/boost/math/tools/cohen_acceleration.hpp
View File

@@ -23,17 +23,18 @@ auto cohen_acceleration(G& generator, std::int64_t n = -1)
using std::pow;
using std::ceil;
using std::sqrt;
auto n_ = static_cast<Real>(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_ = static_cast<Real>(ceil(log(4/std::numeric_limits<Real>::epsilon())*0.5672963285532555));
n_ = static_cast<Real>(ceil(log(Real(4)/std::numeric_limits<Real>::epsilon())*Real(0.5672963285532555)));
n = static_cast<std::int64_t>(n_);
}
// d can get huge and overflow if you pick n too large:
auto d = static_cast<Real>(pow(3 + sqrt(Real(8)), n));
d = (d + 1/d)/2;
auto d = static_cast<Real>(pow(Real(3 + sqrt(Real(8))), n_));
d = (d + Real(1)/d)/2;
Real b = -1;
Real c = -d;
Real s = 0;

22
test/cohen_acceleration_test.cpp
View File

@@ -15,6 +15,10 @@ using boost::multiprecision::float128;
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <cmath>
#if __has_include(<stdfloat>)
# include <stdfloat>
#endif
using boost::math::tools::cohen_acceleration;
using boost::multiprecision::cpp_bin_float_100;
using boost::math::constants::pi;
@@ -72,17 +76,29 @@ void test_divergent()
int main()
{
#ifdef __STDCPP_FLOAT32_T__
test_pisq_div12<std::float32_t>();
test_divergent<std::float32_t>();
#else
test_pisq_div12<float>();
test_pisq_div12<double>();
test_pisq_div12<long double>();
test_divergent<float>();
#endif
#ifdef __STDCPP_FLOAT64_T__
test_pisq_div12<std::float64_t>();
test_divergent<std::float64_t>();
#else
test_divergent<double>();
test_pisq_div12<double>();
#endif
test_divergent<long double>();
test_pisq_div12<long double>();
#ifdef BOOST_HAS_FLOAT128
test_pisq_div12<float128>();
test_divergent<float128>();
#endif
return boost::math::test::report_errors();
}