boost_mirror/math
1
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2025-05-11 21:33:52 +00:00
Code Issues Packages Projects Releases Wiki Activity

Update Bessel functions at infinity. (#1144)

Browse Source 
* Update Bessel functions at infinity.
Also sinc functions, and update tests.
Fixes https://github.com/boostorg/math/issues/1143.

* Correct some test failures.

* Yikes, correct missing if.

* Temporary fix for multiprecision.
REMOVE THIS.
This commit is contained in:
jzmaddock 2024-06-14 18:34:24 +01:00 committed by GitHub
parent c5dc6c74f6
commit dfc2934865
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
10 changed files with 177 additions and 83 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
include/boost/math/special_functions/bessel.hpp
View File

@ -151,7 +151,7 @@ inline T sph_bessel_j_imp(unsigned n, T x, const Policy& pol)
// Special case, n == 0 resolves down to the sinus cardinal of x:
//
if(n == 0)
return boost::math::sinc_pi(x, pol);
return (boost::math::isinf)(x) ? T(0) : boost::math::sinc_pi(x, pol);
//
// Special case for x == 0:
//

153
include/boost/math/special_functions/detail/bessel_ik.hpp
View File

@ -322,83 +322,96 @@ int bessel_ik(T v, T x, T* result_I, T* result_K, int kind, const Policy& pol)
v = -v; // v is non-negative from here
kind |= need_k;
}
n = iround(v, pol);
u = v - n; // -1/2 <= u < 1/2
BOOST_MATH_INSTRUMENT_VARIABLE(n);
BOOST_MATH_INSTRUMENT_VARIABLE(u);
BOOST_MATH_ASSERT(x > 0); // Error handling for x <= 0 handled in cyl_bessel_i and cyl_bessel_k
// x is positive until reflection
W = 1 / x; // Wronskian
if (x <= 2) // x in (0, 2]
{
temme_ik(u, x, &Ku, &Ku1, pol); // Temme series
}
else // x in (2, \infty)
{
CF2_ik(u, x, &Ku, &Ku1, pol); // continued fraction CF2_ik
}
BOOST_MATH_INSTRUMENT_VARIABLE(Ku);
BOOST_MATH_INSTRUMENT_VARIABLE(Ku1);
prev = Ku;
current = Ku1;
T scale = 1;
T scale_sign = 1;
for (k = 1; k <= n; k++) // forward recurrence for K
if (((kind & need_i) == 0) && (fabs(4 * v * v - 25) / (8 * x) < tools::forth_root_epsilon<T>()))
{
T fact = 2 * (u + k) / x;
// Check for overflow: if (max - |prev|) / fact > max, then overflow
// (max - |prev|) / fact > max
// max * (1 - fact) > |prev|
// if fact < 1: safe to compute overflow check
// if fact >= 1: won't overflow
const bool will_overflow = (fact < 1)
? tools::max_value<T>() * (1 - fact) > fabs(prev)
: false;
if(!will_overflow && ((tools::max_value<T>() - fabs(prev)) / fact < fabs(current)))
{
prev /= current;
scale /= current;
scale_sign *= ((boost::math::signbit)(current) ? -1 : 1);
current = 1;
}
next = fact * current + prev;
prev = current;
current = next;
}
Kv = prev;
Kv1 = current;
BOOST_MATH_INSTRUMENT_VARIABLE(Kv);
BOOST_MATH_INSTRUMENT_VARIABLE(Kv1);
if(kind & need_i)
{
T lim = (4 * v * v + 10) / (8 * x);
lim *= lim;
lim *= lim;
lim /= 24;
if((lim < tools::epsilon<T>() * 10) && (x > 100))
{
// x is huge compared to v, CF1 may be very slow
// to converge so use asymptotic expansion for large
// x case instead. Note that the asymptotic expansion
// isn't very accurate - so it's deliberately very hard
// to get here - probably we're going to overflow:
Iv = asymptotic_bessel_i_large_x(v, x, pol);
}
else if((v > 0) && (x / v < 0.25))
{
Iv = bessel_i_small_z_series(v, x, pol);
}
else
{
CF1_ik(v, x, &fv, pol); // continued fraction CF1_ik
Iv = scale * W / (Kv * fv + Kv1); // Wronskian relation
}
// A&S 9.7.2
Iv = std::numeric_limits<T>::quiet_NaN(); // any value will do
T mu = 4 * v * v;
T eight_z = 8 * x;
Kv = 1 + (mu - 1) / eight_z + (mu - 1) * (mu - 9) / (2 * eight_z * eight_z) + (mu - 1) * (mu - 9) * (mu - 25) / (6 * eight_z * eight_z * eight_z);
Kv *= exp(-x) * constants::root_pi<T>() / sqrt(2 * x);
}
else
Iv = std::numeric_limits<T>::quiet_NaN(); // any value will do
{
n = iround(v, pol);
u = v - n; // -1/2 <= u < 1/2
BOOST_MATH_INSTRUMENT_VARIABLE(n);
BOOST_MATH_INSTRUMENT_VARIABLE(u);
BOOST_MATH_ASSERT(x > 0); // Error handling for x <= 0 handled in cyl_bessel_i and cyl_bessel_k
// x is positive until reflection
W = 1 / x; // Wronskian
if (x <= 2) // x in (0, 2]
{
temme_ik(u, x, &Ku, &Ku1, pol); // Temme series
}
else // x in (2, \infty)
{
CF2_ik(u, x, &Ku, &Ku1, pol); // continued fraction CF2_ik
}
BOOST_MATH_INSTRUMENT_VARIABLE(Ku);
BOOST_MATH_INSTRUMENT_VARIABLE(Ku1);
prev = Ku;
current = Ku1;
for (k = 1; k <= n; k++) // forward recurrence for K
{
T fact = 2 * (u + k) / x;
// Check for overflow: if (max - |prev|) / fact > max, then overflow
// (max - |prev|) / fact > max
// max * (1 - fact) > |prev|
// if fact < 1: safe to compute overflow check
// if fact >= 1: won't overflow
const bool will_overflow = (fact < 1)
? tools::max_value<T>() * (1 - fact) > fabs(prev)
: false;
if (!will_overflow && ((tools::max_value<T>() - fabs(prev)) / fact < fabs(current)))
{
prev /= current;
scale /= current;
scale_sign *= ((boost::math::signbit)(current) ? -1 : 1);
current = 1;
}
next = fact * current + prev;
prev = current;
current = next;
}
Kv = prev;
Kv1 = current;
BOOST_MATH_INSTRUMENT_VARIABLE(Kv);
BOOST_MATH_INSTRUMENT_VARIABLE(Kv1);
if (kind & need_i)
{
T lim = (4 * v * v + 10) / (8 * x);
lim *= lim;
lim *= lim;
lim /= 24;
if ((lim < tools::epsilon<T>() * 10) && (x > 100))
{
// x is huge compared to v, CF1 may be very slow
// to converge so use asymptotic expansion for large
// x case instead. Note that the asymptotic expansion
// isn't very accurate - so it's deliberately very hard
// to get here - probably we're going to overflow:
Iv = asymptotic_bessel_i_large_x(v, x, pol);
}
else if ((v > 0) && (x / v < 0.25))
{
Iv = bessel_i_small_z_series(v, x, pol);
}
else
{
CF1_ik(v, x, &fv, pol); // continued fraction CF1_ik
Iv = scale * W / (Kv * fv + Kv1); // Wronskian relation
}
}
else
Iv = std::numeric_limits<T>::quiet_NaN(); // any value will do
}
if (reflect)
{
T z = (u + n % 2);

4
include/boost/math/special_functions/detail/bessel_jy_asym.hpp
View File

@ -69,6 +69,8 @@ inline T asymptotic_bessel_y_large_x_2(T v, T x, const Policy& pol)
BOOST_MATH_STD_USING
// Get the phase and amplitude:
T ampl = asymptotic_bessel_amplitude(v, x);
if (0 == ampl)
return ampl;
T phase = asymptotic_bessel_phase_mx(v, x);
BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
BOOST_MATH_INSTRUMENT_VARIABLE(phase);
@ -97,6 +99,8 @@ inline T asymptotic_bessel_j_large_x_2(T v, T x, const Policy& pol)
BOOST_MATH_STD_USING
// Get the phase and amplitude:
T ampl = asymptotic_bessel_amplitude(v, x);
if (0 == ampl)
return ampl; // shortcut.
T phase = asymptotic_bessel_phase_mx(v, x);
BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
BOOST_MATH_INSTRUMENT_VARIABLE(phase);

7
include/boost/math/special_functions/sinc.hpp
View File

@ -19,6 +19,7 @@
#include <boost/math/tools/precision.hpp>
#include <boost/math/policies/policy.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
#include <limits>
#include <string>
#include <stdexcept>
@ -39,7 +40,11 @@ namespace boost
{
BOOST_MATH_STD_USING
if (abs(x) >= 3.3 * tools::forth_root_epsilon<T>())
if ((boost::math::isinf)(x))
{
return 0;
}
else if (abs(x) >= 3.3 * tools::forth_root_epsilon<T>())
{
return(sin(x)/x);
}

34
include/boost/math/special_functions/sinhc.hpp
View File

@ -16,7 +16,9 @@
#endif
#include <boost/math/tools/precision.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
#include <limits>
#include <string>
#include <stdexcept>
@ -32,8 +34,8 @@ namespace boost
{
// This is the "Hyperbolic Sinus Cardinal" of index Pi.
template<typename T>
inline T sinhc_pi_imp(const T x)
template<typename T, typename Policy>
inline T sinhc_pi_imp(const T x, const Policy&)
{
using ::std::abs;
using ::std::sinh;
@ -43,6 +45,10 @@ namespace boost
static T const taylor_2_bound = sqrt(taylor_0_bound);
static T const taylor_n_bound = sqrt(taylor_2_bound);
if((boost::math::isinf)(x))
{
return policies::raise_overflow_error<T>("sinhc(%1%)", nullptr, Policy());
}
if (abs(x) >= taylor_n_bound)
{
return(sinh(x)/x);
@ -69,20 +75,30 @@ namespace boost
return(result);
}
}
//
// Temporary fix to keep multiprecision happy.
// REMOVE ME!!
//
template<typename T>
inline T sinhc_pi_imp(const T x)
{
return sinhc_pi_imp(x, boost::math::policies::policy<>());
}
} // namespace detail
template <class T, class Policy>
inline typename tools::promote_args<T>::type sinhc_pi(T x, const Policy& pol)
{
typedef typename tools::promote_args<T>::type result_type;
return policies::checked_narrowing_cast<T, Policy>(detail::sinhc_pi_imp(static_cast<result_type>(x), pol), "sinhc(%1%)");
}
template <class T>
inline typename tools::promote_args<T>::type sinhc_pi(T x)
{
typedef typename tools::promote_args<T>::type result_type;
return detail::sinhc_pi_imp(static_cast<result_type>(x));
}
template <class T, class Policy>
inline typename tools::promote_args<T>::type sinhc_pi(T x, const Policy&)
{
return boost::math::sinhc_pi(x);
return sinhc_pi(static_cast<result_type>(x), policies::policy<>());
}
template<typename T, template<typename> class U>

4
test/test_bessel_j.cpp
View File

@ -181,7 +181,7 @@ void expected_results()
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
if((std::numeric_limits<double>::digits != std::numeric_limits<long double>::digits)
BOOST_IF_CONSTEXPR ((std::numeric_limits<double>::digits != std::numeric_limits<long double>::digits)
&& (std::numeric_limits<long double>::digits < 90))
{
// some errors spill over into type double as well:
@ -238,7 +238,7 @@ void expected_results()
".*(JN|j).*|.*Tricky.*", // test data group
".*", 33000, 20000); // test function
}
else if (std::numeric_limits<long double>::digits >= 90)
else BOOST_IF_CONSTEXPR (std::numeric_limits<long double>::digits >= 90)
{
add_expected_result(
".*", // compiler

16
test/test_bessel_j.hpp
View File

@ -280,5 +280,21 @@ void test_bessel(T, const char* name)
BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j(T(0), T(2.5)), boost::math::cyl_bessel_j(T(0), T(-2.5)));
BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j(T(1), T(2.5)), -boost::math::cyl_bessel_j(T(1), T(-2.5)));
BOOST_CHECK_CLOSE_FRACTION(boost::math::cyl_bessel_j(364, T(38.5)), SC_(1.793940496519190500748409872348034004417458734118663909894e-309), tolerance);
//
// Special cases at infinity:
//
BOOST_IF_CONSTEXPR (std::numeric_limits<T>::has_infinity)
{
BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j(T(0), std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j(T(2), std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j(T(1.25), std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::sph_bessel(0, std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::sph_bessel(1, std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::sph_bessel(2, std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::sph_bessel(3, std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j(T(0), -std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j(T(2), -std::numeric_limits<T>::infinity()), T(0));
}
}

5
test/test_bessel_k.hpp
View File

@ -184,6 +184,11 @@ void test_bessel(T, const char* name)
{
BOOST_CHECK_EQUAL(boost::math::cyl_bessel_k(16, ldexp(T(1), -1021)), std::numeric_limits<T>::infinity());
}
BOOST_CHECK_EQUAL(boost::math::cyl_bessel_k(T(0), std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::cyl_bessel_k(T(1), std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::cyl_bessel_k(T(2), std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::cyl_bessel_k(T(2.25), std::numeric_limits<T>::infinity()), T(0));
}
BOOST_CHECK_THROW(boost::math::cyl_bessel_k(T(1.25), T(0)), std::domain_error);
BOOST_CHECK_THROW(boost::math::cyl_bessel_k(T(-1.25), T(0)), std::domain_error);

8
test/test_bessel_y.hpp
View File

@ -232,7 +232,15 @@ void test_bessel(T, const char* name)
{
BOOST_CHECK_EQUAL(boost::math::cyl_neumann(T(121.25), T(0.25)), -std::numeric_limits<T>::infinity());
}
BOOST_CHECK_EQUAL(boost::math::cyl_neumann(T(0), std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::cyl_neumann(T(1), std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::cyl_neumann(T(2), std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::cyl_neumann(T(2.25), std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::sph_neumann(0, std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::sph_neumann(1, std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::sph_neumann(2, std::numeric_limits<T>::infinity()), T(0));
}
BOOST_CHECK_THROW(boost::math::cyl_neumann(T(0), T(-1)), std::domain_error);
BOOST_CHECK_THROW(boost::math::cyl_neumann(T(0.2), T(-1)), std::domain_error);
BOOST_CHECK_THROW(boost::math::cyl_neumann(T(2), T(0)), std::domain_error);

27
test/test_sinc.cpp
View File

@ -7,6 +7,7 @@
#include "test_sinc.hpp"
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/math/special_functions/next.hpp>
#include <boost/math/special_functions/sinhc.hpp>
//
// DESCRIPTION:
@ -80,6 +81,32 @@ void test_close_to_transition()
BOOST_CHECK_LE(boost::math::epsilon_difference(result, expected), 3);
val = boost::math::float_next(val);
}
BOOST_IF_CONSTEXPR(std::numeric_limits<T>::has_infinity)
{
BOOST_CHECK_EQUAL(boost::math::sinc_pi(std::numeric_limits<T>::infinity()), T(0));
BOOST_CHECK_EQUAL(boost::math::sinc_pi(-std::numeric_limits<T>::infinity()), T(0));
BOOST_IF_CONSTEXPR(boost::math::policies::policy<>::overflow_error_type::value == boost::math::policies::throw_on_error)
{
BOOST_CHECK_THROW(boost::math::sinhc_pi(std::numeric_limits<T>::infinity()), std::overflow_error);
BOOST_CHECK_THROW(boost::math::sinhc_pi(-std::numeric_limits<T>::infinity()), std::overflow_error);
}
else
{
BOOST_CHECK_EQUAL(boost::math::sinhc_pi(std::numeric_limits<T>::infinity()), std::numeric_limits<T>::infinity());
BOOST_CHECK_EQUAL(boost::math::sinhc_pi(-std::numeric_limits<T>::infinity()), std::numeric_limits<T>::infinity());
}
}
BOOST_IF_CONSTEXPR(boost::math::policies::policy<>::overflow_error_type::value == boost::math::policies::throw_on_error)
{
BOOST_CHECK_THROW(boost::math::sinhc_pi(boost::math::tools::max_value<T>()), std::overflow_error);
BOOST_CHECK_THROW(boost::math::sinhc_pi(-boost::math::tools::max_value<T>()), std::overflow_error);
}
else BOOST_IF_CONSTEXPR(std::numeric_limits<T>::has_infinity)
{
BOOST_CHECK_EQUAL(boost::math::sinhc_pi(boost::math::tools::max_value<T>()), std::numeric_limits<T>::infinity());
BOOST_CHECK_EQUAL(boost::math::sinhc_pi(-boost::math::tools::max_value<T>()), std::numeric_limits<T>::infinity());
}
}