diff --git a/include/boost/math/special_functions/detail/hypergeometric_asym.hpp b/include/boost/math/special_functions/detail/hypergeometric_asym.hpp index e603957180..41a565854d 100644 --- a/include/boost/math/special_functions/detail/hypergeometric_asym.hpp +++ b/include/boost/math/special_functions/detail/hypergeometric_asym.hpp @@ -44,8 +44,10 @@ } else { - e = z > static_cast((std::numeric_limits::max)()) ? (std::numeric_limits::max)() : lltrunc(z, pol); + e = z >= static_cast((std::numeric_limits::max)()) ? (std::numeric_limits::max)() : lltrunc(z, pol); log_scaling += e; + if (z - e > tools::log_max_value()) + return policies::raise_overflow_error(function, nullptr, pol); prefix = exp(z - e); } if ((fabs(a) < 10) && (fabs(b) < 10)) diff --git a/include/boost/math/special_functions/hypergeometric_1F1.hpp b/include/boost/math/special_functions/hypergeometric_1F1.hpp index 5b24c58cd9..ee857183a0 100644 --- a/include/boost/math/special_functions/hypergeometric_1F1.hpp +++ b/include/boost/math/special_functions/hypergeometric_1F1.hpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -637,7 +638,7 @@ namespace boost { namespace math { namespace detail { static const thread_local long long max_scaling = lltrunc(boost::math::tools::log_max_value()) - 2; static const thread_local T max_scale_factor = exp(T(max_scaling)); - while (log_scaling > max_scaling) + while (!(boost::math::isinf)(result) && (log_scaling > max_scaling)) { result *= max_scale_factor; log_scaling -= max_scaling; diff --git a/test/test_1F1.hpp b/test/test_1F1.hpp index b4c67faf1a..9333480245 100644 --- a/test/test_1F1.hpp +++ b/test/test_1F1.hpp @@ -420,6 +420,13 @@ void test_spots6(T, const char* type_name) T val = boost::math::hypergeometric_1F1(hypergeometric_1F1_big_bugs[i][0], hypergeometric_1F1_big_bugs[i][1], hypergeometric_1F1_big_bugs[i][2]); BOOST_CHECK((boost::math::isinf)(val)); } + // https://github.com/boostorg/math/issues/1288 + BOOST_CHECK_EQUAL(boost::math::hypergeometric_1F1(SC_(2.0), SC_(3.0), SC_(1e30)), + (std::numeric_limits::infinity)()); + // This is the largest log_hypergeometric_1F1 we can currently compute, beyond this the fix for the above + // bug causes premature overflow. However, there is no easy fix and Wolfram also fails in this case. It's **hard** :( + BOOST_CHECK_CLOSE_FRACTION(boost::math::log_hypergeometric_1F1(SC_(2.0), SC_(3.0), static_cast((std::numeric_limits::max)())), + SC_(9223372036854775765.02487480528339082), boost::math::tools::epsilon() * 10); } }