From faaf4759aca2db615269caf725006695775ff20c Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Mon, 22 May 2023 13:48:33 +0200
Subject: [PATCH] Fix stack overflow in cohen acceleration from GCC bug

---
 .../boost/math/tools/cohen_acceleration.hpp   |  7 +++---
 test/cohen_acceleration_test.cpp              | 22 ++++++++++++++++---
 2 files changed, 23 insertions(+), 6 deletions(-)

diff --git a/include/boost/math/tools/cohen_acceleration.hpp b/include/boost/math/tools/cohen_acceleration.hpp
index 716245f72..2e76a75f8 100644
--- a/include/boost/math/tools/cohen_acceleration.hpp
+++ b/include/boost/math/tools/cohen_acceleration.hpp
@@ -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;
diff --git a/test/cohen_acceleration_test.cpp b/test/cohen_acceleration_test.cpp
index f300bdad6..ad70c3f82 100644
--- a/test/cohen_acceleration_test.cpp
+++ b/test/cohen_acceleration_test.cpp
@@ -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();
 }