From dfc29348653bbc976c9941a9fbc9f49d68cb954e Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 14 Jun 2024 18:34:24 +0100 Subject: [PATCH] Update Bessel functions at infinity. (#1144) * 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. --- .../boost/math/special_functions/bessel.hpp | 2 +- .../special_functions/detail/bessel_ik.hpp | 141 ++++++++++-------- .../detail/bessel_jy_asym.hpp | 4 + include/boost/math/special_functions/sinc.hpp | 7 +- .../boost/math/special_functions/sinhc.hpp | 32 +++- test/test_bessel_j.cpp | 4 +- test/test_bessel_j.hpp | 16 ++ test/test_bessel_k.hpp | 5 + test/test_bessel_y.hpp | 8 + test/test_sinc.cpp | 27 ++++ 10 files changed, 170 insertions(+), 76 deletions(-) diff --git a/include/boost/math/special_functions/bessel.hpp b/include/boost/math/special_functions/bessel.hpp index 82bed6e490..e9677d3c79 100644 --- a/include/boost/math/special_functions/bessel.hpp +++ b/include/boost/math/special_functions/bessel.hpp @@ -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: // diff --git a/include/boost/math/special_functions/detail/bessel_ik.hpp b/include/boost/math/special_functions/detail/bessel_ik.hpp index da2310866f..0c653b4753 100644 --- a/include/boost/math/special_functions/detail/bessel_ik.hpp +++ b/include/boost/math/special_functions/detail/bessel_ik.hpp @@ -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 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() * (1 - fact) > fabs(prev) - : false; - if(!will_overflow && ((tools::max_value() - 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; + // A&S 9.7.2 + Iv = std::numeric_limits::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() / sqrt(2 * x); } - Kv = prev; - Kv1 = current; - BOOST_MATH_INSTRUMENT_VARIABLE(Kv); - BOOST_MATH_INSTRUMENT_VARIABLE(Kv1); - if(kind & need_i) + else { - T lim = (4 * v * v + 10) / (8 * x); - lim *= lim; - lim *= lim; - lim /= 24; - if((lim < tools::epsilon() * 10) && (x > 100)) + 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] { - // 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); + temme_ik(u, x, &Ku, &Ku1, pol); // Temme series } - else if((v > 0) && (x / v < 0.25)) + else // x in (2, \infty) { - Iv = bessel_i_small_z_series(v, x, pol); + CF2_ik(u, x, &Ku, &Ku1, pol); // continued fraction CF2_ik } - else + 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() * (1 - fact) > fabs(prev) + : false; + if (!will_overflow && ((tools::max_value() - 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) { - CF1_ik(v, x, &fv, pol); // continued fraction CF1_ik - Iv = scale * W / (Kv * fv + Kv1); // Wronskian relation + T lim = (4 * v * v + 10) / (8 * x); + lim *= lim; + lim *= lim; + lim /= 24; + if ((lim < tools::epsilon() * 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::quiet_NaN(); // any value will do } - else - Iv = std::numeric_limits::quiet_NaN(); // any value will do - if (reflect) { T z = (u + n % 2); diff --git a/include/boost/math/special_functions/detail/bessel_jy_asym.hpp b/include/boost/math/special_functions/detail/bessel_jy_asym.hpp index 16b6543f90..cb09b202d5 100644 --- a/include/boost/math/special_functions/detail/bessel_jy_asym.hpp +++ b/include/boost/math/special_functions/detail/bessel_jy_asym.hpp @@ -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); diff --git a/include/boost/math/special_functions/sinc.hpp b/include/boost/math/special_functions/sinc.hpp index 7220d30b81..ff1b2e966b 100644 --- a/include/boost/math/special_functions/sinc.hpp +++ b/include/boost/math/special_functions/sinc.hpp @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -39,7 +40,11 @@ namespace boost { BOOST_MATH_STD_USING - if (abs(x) >= 3.3 * tools::forth_root_epsilon()) + if ((boost::math::isinf)(x)) + { + return 0; + } + else if (abs(x) >= 3.3 * tools::forth_root_epsilon()) { return(sin(x)/x); } diff --git a/include/boost/math/special_functions/sinhc.hpp b/include/boost/math/special_functions/sinhc.hpp index 6606ad3ea8..70f6c72e67 100644 --- a/include/boost/math/special_functions/sinhc.hpp +++ b/include/boost/math/special_functions/sinhc.hpp @@ -16,7 +16,9 @@ #endif #include +#include #include +#include #include #include #include @@ -32,8 +34,8 @@ namespace boost { // This is the "Hyperbolic Sinus Cardinal" of index Pi. - template - inline T sinhc_pi_imp(const T x) + template + 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("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 + inline T sinhc_pi_imp(const T x) + { + return sinhc_pi_imp(x, boost::math::policies::policy<>()); + } } // namespace detail - template - inline typename tools::promote_args::type sinhc_pi(T x) + template + inline typename tools::promote_args::type sinhc_pi(T x, const Policy& pol) { typedef typename tools::promote_args::type result_type; - return detail::sinhc_pi_imp(static_cast(x)); + return policies::checked_narrowing_cast(detail::sinhc_pi_imp(static_cast(x), pol), "sinhc(%1%)"); } - template - inline typename tools::promote_args::type sinhc_pi(T x, const Policy&) + template + inline typename tools::promote_args::type sinhc_pi(T x) { - return boost::math::sinhc_pi(x); + typedef typename tools::promote_args::type result_type; + return sinhc_pi(static_cast(x), policies::policy<>()); } template class U> diff --git a/test/test_bessel_j.cpp b/test/test_bessel_j.cpp index e2a7e83ccb..19a5f7426e 100644 --- a/test/test_bessel_j.cpp +++ b/test/test_bessel_j.cpp @@ -181,7 +181,7 @@ void expected_results() #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS - if((std::numeric_limits::digits != std::numeric_limits::digits) + BOOST_IF_CONSTEXPR ((std::numeric_limits::digits != std::numeric_limits::digits) && (std::numeric_limits::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::digits >= 90) + else BOOST_IF_CONSTEXPR (std::numeric_limits::digits >= 90) { add_expected_result( ".*", // compiler diff --git a/test/test_bessel_j.hpp b/test/test_bessel_j.hpp index 2d10bc6c94..82106213ea 100644 --- a/test/test_bessel_j.hpp +++ b/test/test_bessel_j.hpp @@ -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::has_infinity) + { + BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j(T(0), std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j(T(2), std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j(T(1.25), std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::sph_bessel(0, std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::sph_bessel(1, std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::sph_bessel(2, std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::sph_bessel(3, std::numeric_limits::infinity()), T(0)); + + BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j(T(0), -std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j(T(2), -std::numeric_limits::infinity()), T(0)); + } } diff --git a/test/test_bessel_k.hpp b/test/test_bessel_k.hpp index c0eb9aefdf..22df3218f0 100644 --- a/test/test_bessel_k.hpp +++ b/test/test_bessel_k.hpp @@ -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::infinity()); } + + BOOST_CHECK_EQUAL(boost::math::cyl_bessel_k(T(0), std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::cyl_bessel_k(T(1), std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::cyl_bessel_k(T(2), std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::cyl_bessel_k(T(2.25), std::numeric_limits::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); diff --git a/test/test_bessel_y.hpp b/test/test_bessel_y.hpp index 3e489d0952..28361a227c 100644 --- a/test/test_bessel_y.hpp +++ b/test/test_bessel_y.hpp @@ -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::infinity()); } + BOOST_CHECK_EQUAL(boost::math::cyl_neumann(T(0), std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::cyl_neumann(T(1), std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::cyl_neumann(T(2), std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::cyl_neumann(T(2.25), std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::sph_neumann(0, std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::sph_neumann(1, std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::sph_neumann(2, std::numeric_limits::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); diff --git a/test/test_sinc.cpp b/test/test_sinc.cpp index ad0eaeb940..e72c640cb3 100644 --- a/test/test_sinc.cpp +++ b/test/test_sinc.cpp @@ -7,6 +7,7 @@ #include "test_sinc.hpp" #include #include +#include // // 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::has_infinity) + { + BOOST_CHECK_EQUAL(boost::math::sinc_pi(std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::sinc_pi(-std::numeric_limits::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::infinity()), std::overflow_error); + BOOST_CHECK_THROW(boost::math::sinhc_pi(-std::numeric_limits::infinity()), std::overflow_error); + } + else + { + BOOST_CHECK_EQUAL(boost::math::sinhc_pi(std::numeric_limits::infinity()), std::numeric_limits::infinity()); + BOOST_CHECK_EQUAL(boost::math::sinhc_pi(-std::numeric_limits::infinity()), std::numeric_limits::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()), std::overflow_error); + BOOST_CHECK_THROW(boost::math::sinhc_pi(-boost::math::tools::max_value()), std::overflow_error); + } + else BOOST_IF_CONSTEXPR(std::numeric_limits::has_infinity) + { + BOOST_CHECK_EQUAL(boost::math::sinhc_pi(boost::math::tools::max_value()), std::numeric_limits::infinity()); + BOOST_CHECK_EQUAL(boost::math::sinhc_pi(-boost::math::tools::max_value()), std::numeric_limits::infinity()); + } }