diff --git a/include/boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp b/include/boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp index 6f1461ffc9..07cb236db1 100644 --- a/include/boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp +++ b/include/boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp @@ -291,6 +291,14 @@ ak += 2; integer_part -= 2; } + if (ak - 1 == b) + { + // When ak - 1 == b are recursion coefficients dissappear to zero and + // we end up with a NaN result. Reduce the recursion steps by 1 to + // avoid this. We rely on |b| small and therefore no infinite recursion. + ak -= 1; + integer_part += 1; + } if (-integer_part > static_cast(policies::get_max_series_iterations())) return policies::raise_evaluation_error(function, "1F1 arguments sit in a range with a so negative that we have no evaluation method, got a = %1%", std::numeric_limits::quiet_NaN(), pol); diff --git a/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp b/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp index 3745f7e24d..8ed8c8086f 100644 --- a/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp +++ b/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp @@ -132,6 +132,7 @@ Real scaling_factor = exp(Real(log_scaling_factor)); Real term_m1; long long local_scaling = 0; + bool have_no_correct_bits = false; if ((aj.size() == 1) && (bj.size() == 0)) { @@ -238,10 +239,20 @@ break; if (abs_result * tol > abs(result)) { - // We have no correct bits in the result... just give up! - result = boost::math::policies::raise_evaluation_error("boost::math::hypergeometric_pFq<%1%>", "Cancellation is so severe that no bits in the reuslt are correct, last result was %1%", Real(result * exp(Real(log_scale))), pol); - return std::make_pair(result, result); + // Check if result is so small compared to abs_resuslt that there are no longer any + // correct bits... we require two consecutive passes here before aborting to + // avoid false positives when result transiently drops to near zero then rebounds. + if (have_no_correct_bits) + { + // We have no correct bits in the result... just give up! + result = boost::math::policies::raise_evaluation_error("boost::math::hypergeometric_pFq<%1%>", "Cancellation is so severe that no bits in the reuslt are correct, last result was %1%", Real(result * exp(Real(log_scale))), pol); + return std::make_pair(result, result); + } + else + have_no_correct_bits = true; } + else + have_no_correct_bits = false; term0 = term; } //std::cout << "result = " << result << std::endl; diff --git a/test/test_1F1.hpp b/test/test_1F1.hpp index 2e7805d442..5edf8cdc34 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, 91> hypergeometric_1F1_bugs = { { + static const std::array, 93> 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) }} , @@ -287,6 +287,11 @@ void test_spots6(T, const char* type_name) { { (T)std::ldexp((double)-9751199809536000, -45), (T)std::ldexp((double)-17654191685632000, -47), (T)std::ldexp((double)10587451850752000, -47), SC_(-1.9601415510439595625538337964298353914980331018955e+68) }}, { { (T)std::ldexp((double)-15233620754432000, -45), (T)std::ldexp((double)-12708283072512000, -46), (T)std::ldexp((double)10255461007360000, -46), SC_(-5.4344106361679075861859567858016187271235441673635e+125) }}, { { (T)std::ldexp((double)-11241354149888000, -45), (T)std::ldexp((double)-9580579905536000, -45), (T)std::ldexp((double)12224976846848000, -47), SC_(12046856548470067405870726490464935201150430438.035) }}, + // + // Bugs found while testing color maps: + // + { { SC_(0.078125000000000000), SC_(-0.039062500000000000), SC_(0.5), SC_(-0.3371910410915676603577770246237158427221) }}, + { { SC_(-19.492187500000000), SC_(0.50781250000000000), SC_(0.5), SC_(1.2551298228307647570646714060395253358015) }}, } }; static const std::array, 2> hypergeometric_1F1_big_bugs = { { #if DBL_MAX_EXP == LDBL_MAX_EXP