boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

1F1 coverage.

Browse Source
This commit is contained in:
jzmaddock
2024-08-06 12:09:33 +01:00
parent 36d46448a9
commit 977feaf84a
2 changed files with 26 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files

27
include/boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp
View File

@@ -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<T>::quiet_NaN());
cache_offset = -cache_size;
return;
}
if ((term * bessel_cache[cache_size - 1] < tools::min_value<T>() / (tools::epsilon<T>() * tools::epsilon<T>())) || !(boost::math::isfinite)(term) || (!std::numeric_limits<T>::has_infinity && (fabs(term) > tools::max_value<T>())))
if ((fabs(term * bessel_cache[cache_size - 1]) < tools::min_value<T>() / (tools::epsilon<T>() * tools::epsilon<T>())) || !(boost::math::isfinite)(term) || (!std::numeric_limits<T>::has_infinity && (fabs(term) > tools::max_value<T>())))
{
term = -log(fabs(bessel_arg)) * b_minus_1_plus_n / 2;
log_scale = lltrunc(term);
@@ -88,15 +92,27 @@
if constexpr (std::numeric_limits<T>::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<T>::quiet_NaN());
}
}
else
if ((boost::math::isnan)(bessel_cache[cache_size - 1]) || (fabs(bessel_cache[cache_size - 1]) >= tools::max_value<T>()))
{
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<T>::quiet_NaN());
}
#else
if ((std::numeric_limits<T>::has_infinity && !(boost::math::isfinite)(bessel_cache[cache_size - 1]))
|| (!std::numeric_limits<T>::has_infinity && ((boost::math::isnan)(bessel_cache[cache_size - 1]) || (fabs(bessel_cache[cache_size - 1]) >= tools::max_value<T>()))))
{
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<T>::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;

3
test/test_1F1.hpp
View File

@@ -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>, 188> hypergeometric_1F1_bugs = { {
static const std::array<std::array<T, 4>, 189> hypergeometric_1F1_bugs = { {
{ { static_cast<T>(17955.561660766602), static_cast<T>(9.6968994205831605e-09), static_cast<T>(-82.406154185533524), SC_(6.98056008378736714088730927132364938220428678e-11) }},
{ { static_cast<T>(17955.561660766602), static_cast<T>(-9.6968994205831605e-09), static_cast<T>(-82.406154185533524), SC_(-6.98055306629610746072607353939306734740549551e-11) }},
{ { static_cast<T>(-17955.561660766602), static_cast<T>(-9.6968994205831605e-09), static_cast<T>(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<std::array<T, 4>, 2> hypergeometric_1F1_big_bugs = { {
#if DBL_MAX_EXP == LDBL_MAX_EXP