diff --git a/include/boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp b/include/boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp index 860217b05..07d9dabf8 100644 --- a/include/boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp +++ b/include/boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp @@ -74,8 +74,12 @@ { // We get very limited precision due to rapid denormalisation/underflow of the Bessel values, raise an exception and try something else: policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Underflow in Bessel functions", bessel_cache[cache_size - 1], pol); + // Exceptions are off if we get here, just fill the cache with NaN's and we'll let this method fail and fallback later: + std::fill(bessel_cache.begin(), bessel_cache.end(), std::numeric_limits::quiet_NaN()); + cache_offset = -cache_size; + return; } - if ((term * bessel_cache[cache_size - 1] < tools::min_value() / (tools::epsilon() * tools::epsilon())) || !(boost::math::isfinite)(term) || (!std::numeric_limits::has_infinity && (fabs(term) > tools::max_value()))) + if ((fabs(term * bessel_cache[cache_size - 1]) < tools::min_value() / (tools::epsilon() * tools::epsilon())) || !(boost::math::isfinite)(term) || (!std::numeric_limits::has_infinity && (fabs(term) > tools::max_value()))) { term = -log(fabs(bessel_arg)) * b_minus_1_plus_n / 2; log_scale = lltrunc(term); @@ -88,15 +92,27 @@ if constexpr (std::numeric_limits::has_infinity) { if (!(boost::math::isfinite)(bessel_cache[cache_size - 1])) + { policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Expected finite Bessel function result but got %1%", bessel_cache[cache_size - 1], pol); + // Exceptions are off if we get here, just fill the cache with NaN's and we'll let this method fail and fallback later: + std::fill(bessel_cache.begin(), bessel_cache.end(), std::numeric_limits::quiet_NaN()); + } } else if ((boost::math::isnan)(bessel_cache[cache_size - 1]) || (fabs(bessel_cache[cache_size - 1]) >= tools::max_value())) + { policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Expected finite Bessel function result but got %1%", bessel_cache[cache_size - 1], pol); + // Exceptions are off if we get here, just fill the cache with NaN's and we'll let this method fail and fallback later: + std::fill(bessel_cache.begin(), bessel_cache.end(), std::numeric_limits::quiet_NaN()); + } #else if ((std::numeric_limits::has_infinity && !(boost::math::isfinite)(bessel_cache[cache_size - 1])) || (!std::numeric_limits::has_infinity && ((boost::math::isnan)(bessel_cache[cache_size - 1]) || (fabs(bessel_cache[cache_size - 1]) >= tools::max_value())))) + { policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Expected finite Bessel function result but got %1%", bessel_cache[cache_size - 1], pol); + // Exceptions are off if we get here, just fill the cache with NaN's and we'll let this method fail and fallback later: + std::fill(bessel_cache.begin(), bessel_cache.end(), std::numeric_limits::quiet_NaN()); + } #endif cache_offset = -cache_size; refill_cache(); @@ -108,8 +124,13 @@ // very small (or zero) when b == 2a: // BOOST_MATH_STD_USING + // + // Except in the multiprecision case, we have probably illiminated anything + // would need more than the default 64 Bessel Functions. Anything more + // than that risks becoming a divergent series anyway... + // if(n - 2 - cache_offset >= cache_size) - refill_cache(); + refill_cache(); // LCOV_EXCL_LINE T result = A_minus_2 * term * bessel_cache[n - 2 - cache_offset]; term *= mult; ++n; @@ -122,7 +143,7 @@ if (A_minus_2 != 0) { if (n - 2 - cache_offset >= cache_size) - refill_cache(); + refill_cache(); // LCOV_EXCL_LINE result += A_minus_2 * term * bessel_cache[n - 2 - cache_offset]; } term *= mult; diff --git a/test/test_1F1.hpp b/test/test_1F1.hpp index 1d4b9b307..b4c67faf1 100644 --- a/test/test_1F1.hpp +++ b/test/test_1F1.hpp @@ -162,7 +162,7 @@ void test_spots5(T, const char* type_name) template void test_spots6(T, const char* type_name) { - static const std::array, 188> hypergeometric_1F1_bugs = { { + static const std::array, 189> hypergeometric_1F1_bugs = { { { { static_cast(17955.561660766602), static_cast(9.6968994205831605e-09), static_cast(-82.406154185533524), SC_(6.98056008378736714088730927132364938220428678e-11) }}, { { static_cast(17955.561660766602), static_cast(-9.6968994205831605e-09), static_cast(-82.406154185533524), SC_(-6.98055306629610746072607353939306734740549551e-11) }}, { { static_cast(-17955.561660766602), static_cast(-9.6968994205831605e-09), static_cast(82.406154185533524), SC_(-42897094853118832762870100.8669248353530950866) }} , @@ -399,6 +399,7 @@ void test_spots6(T, const char* type_name) // Uncovered code: { { SC_(3.125), SC_(-4.25), SC_(0.25), SC_(0.8377654489787288706170314748630804839990046835429477943509538691)}}, { { SC_(3.125), SC_(-4.25), SC_(6.25), SC_(-1.0732328736644883882525050418994056179896661423277961620744e8)}}, + { { SC_(-1792.0615542826642), SC_(1252.5126953125000), SC_(15.873352050781250), SC_(9.5429538115592332178745238256843491134483360764267909811922e-11)}}, } }; static const std::array, 2> hypergeometric_1F1_big_bugs = { { #if DBL_MAX_EXP == LDBL_MAX_EXP