Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Minor 1F1 bug fixes. #778

Merged
merged 2 commits into from
Mar 13, 2022
Merged

Minor 1F1 bug fixes. #778

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
2519641
Two minor 1F1 bugfixes:
jzmaddock Mar 10, 2022
bd7fcf1
Correct conceptual usage of RealType.
jzmaddock Mar 11, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions include/boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<std::intmax_t>(policies::get_max_series_iterations<Policy>()))
return policies::raise_evaluation_error<T>(function, "1F1 arguments sit in a range with a so negative that we have no evaluation method, got a = %1%", std::numeric_limits<T>::quiet_NaN(), pol);
Expand Down
17 changes: 14 additions & 3 deletions include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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))
{
Expand Down Expand Up @@ -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;
Expand Down
7 changes: 6 additions & 1 deletion test/test_1F1.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,7 +162,7 @@ void test_spots5(T, const char* type_name)
template <class T>
void test_spots6(T, const char* type_name)
{
static const std::array<std::array<T, 4>, 91> hypergeometric_1F1_bugs = { {
static const std::array<std::array<T, 4>, 93> hypergeometric_1F1_bugs = { {
{ { static_cast<double>(17955.561660766602), static_cast<double>(9.6968994205831605e-09), static_cast<double>(-82.406154185533524), SC_(6.98056008378736714088730927132364938220428678e-11) }},
{ { static_cast<double>(17955.561660766602), static_cast<double>(-9.6968994205831605e-09), static_cast<double>(-82.406154185533524), SC_(-6.98055306629610746072607353939306734740549551e-11) }},
{ { static_cast<double>(-17955.561660766602), static_cast<double>(-9.6968994205831605e-09), static_cast<double>(82.406154185533524), SC_(-42897094853118832762870100.8669248353530950866) }} ,
Expand Down Expand Up @@ -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<std::array<T, 4>, 2> hypergeometric_1F1_big_bugs = { {
#if DBL_MAX_EXP == LDBL_MAX_EXP
Expand Down