diff --git a/include/boost/math/special_functions/detail/bessel_ik.hpp b/include/boost/math/special_functions/detail/bessel_ik.hpp index b0ae7cadb7..c5e7402936 100644 --- a/include/boost/math/special_functions/detail/bessel_ik.hpp +++ b/include/boost/math/special_functions/detail/bessel_ik.hpp @@ -208,7 +208,7 @@ BOOST_MATH_GPU_ENABLED int CF1_ik(T v, T x, T* fv, const Policy& pol) // z1 / z0 = U(v+1.5, 2v+1, 2x) / U(v+0.5, 2v+1, 2x), see // Thompson and Barnett, Computer Physics Communications, vol 47, 245 (1987) template -BOOST_MATH_GPU_ENABLED int CF2_ik(T v, T x, T* Kv, T* Kv1, const Policy& pol) +BOOST_MATH_GPU_ENABLED int CF2_ik(T v, T x, T* Kv, T* Kv1, T* Kv_scaled, T* Kv1_scaled, const Policy& pol) { BOOST_MATH_STD_USING using namespace boost::math::constants; @@ -282,11 +282,17 @@ BOOST_MATH_GPU_ENABLED int CF2_ik(T v, T x, T* Kv, T* Kv1, const Policy& pol) } policies::check_series_iterations("boost::math::bessel_ik<%1%>(%1%,%1%) in CF2_ik", k, pol); + const T ratio = (0.5f + v + x + (v * v - 0.25f) * f) / x; if(-x < tools::log_min_value()) *Kv = exp(0.5f * log(pi() / (2 * x)) - x - log(S)); else *Kv = sqrt(pi() / (2 * x)) * exp(-x) / S; - *Kv1 = *Kv * (0.5f + v + x + (v * v - 0.25f) * f) / x; + *Kv1 = *Kv * ratio; + if ((*Kv < tools::min_value()) || (*Kv1 < tools::min_value())) + { + *Kv_scaled = sqrt(pi() / 2) / sqrt(x) / S; + *Kv1_scaled = *Kv_scaled * ratio; + } BOOST_MATH_INSTRUMENT_VARIABLE(*Kv); BOOST_MATH_INSTRUMENT_VARIABLE(*Kv1); @@ -305,9 +311,10 @@ BOOST_MATH_GPU_ENABLED int bessel_ik(T v, T x, T* result_I, T* result_K, int kin { // Kv1 = K_(v+1), fv = I_(v+1) / I_v // Ku1 = K_(u+1), fu = I_(u+1) / I_u - T u, Iv, Kv, Kv1, Ku, Ku1, fv; + T u, Iv, Kv, Kv1, Ku, Ku1, Ku_scaled = 0, Ku1_scaled = 0, fv; T W, current, prev, next; bool reflect = false; + bool use_scaled_k = false; unsigned n, k; int org_kind = kind; BOOST_MATH_INSTRUMENT_VARIABLE(v); @@ -329,6 +336,7 @@ BOOST_MATH_GPU_ENABLED int bessel_ik(T v, T x, T* result_I, T* result_K, int kin T scale = 1; T scale_sign = 1; + int exp2_scale = 0; n = iround(v, pol); u = v - n; // -1/2 <= u < 1/2 @@ -356,12 +364,13 @@ BOOST_MATH_GPU_ENABLED int bessel_ik(T v, T x, T* result_I, T* result_K, int kin } else // x in (2, \infty) { - CF2_ik(u, x, &Ku, &Ku1, pol); // continued fraction CF2_ik + CF2_ik(u, x, &Ku, &Ku1, &Ku_scaled, &Ku1_scaled, pol); // continued fraction CF2_ik } BOOST_MATH_INSTRUMENT_VARIABLE(Ku); BOOST_MATH_INSTRUMENT_VARIABLE(Ku1); - prev = Ku; - current = Ku1; + use_scaled_k = ((kind & need_i) == 0) && (x > 2) && ((Ku < tools::min_value()) || (Ku1 < tools::min_value())); + prev = use_scaled_k ? Ku_scaled : Ku; + current = use_scaled_k ? Ku1_scaled : Ku1; for (k = 1; k <= n; k++) // forward recurrence for K { T fact = 2 * (u + k) / x; @@ -375,10 +384,22 @@ BOOST_MATH_GPU_ENABLED int bessel_ik(T v, T x, T* result_I, T* result_K, int kin : 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; + if (use_scaled_k) + { + // Rescale by a power of two so no rounding error is introduced + int e2 = 0; + (void)frexp(current, &e2); + prev = ldexp(prev, -e2); + current = ldexp(current, -e2); + exp2_scale += e2; + } + else + { + prev /= current; + scale /= current; + scale_sign *= ((boost::math::signbit)(current) ? -1 : 1); + current = 1; + } } next = fact * current + prev; prev = current; @@ -432,7 +453,40 @@ BOOST_MATH_GPU_ENABLED int bessel_ik(T v, T x, T* result_I, T* result_K, int kin { *result_I = Iv; } - if(tools::max_value() * scale < Kv) + if (use_scaled_k) + { + // Kv currently holds exp(x) * K_v(x) / 2^exp2_scale + int e2 = 0; + const T m = frexp(Kv, &e2); + exp2_scale += e2; + // Split exp(-x) into 2^j equal factors exp(-x / 2^j) that are each representable, + // and spread the binary exponent over them so no partial product leaves range. + T xr = x; + int parts = 1; + while (-xr < tools::log_min_value()) + { + xr /= 2; + parts *= 2; + } + const T factor = exp(-xr); + T result = m; + int e_rem = exp2_scale; + bool overflow = false; + for (int i = 0; i < parts; ++i) + { + const int e_i = e_rem / (parts - i); + e_rem -= e_i; + const T term = ldexp(factor, e_i); + if (result > tools::max_value() / term) + { + overflow = true; + break; + } + result *= term; + } + *result_K = overflow ? ((org_kind & need_k) ? policies::raise_overflow_error(function, nullptr, pol) : T(0)) : result; + } + else if(tools::max_value() * scale < Kv) *result_K = (org_kind & need_k) ? T(sign(Kv) * scale_sign * policies::raise_overflow_error(function, nullptr, pol)) : T(0); else *result_K = Kv / scale; diff --git a/include/boost/math/special_functions/detail/bessel_kn.hpp b/include/boost/math/special_functions/detail/bessel_kn.hpp index 41becc8aa9..a07bbe75c6 100644 --- a/include/boost/math/special_functions/detail/bessel_kn.hpp +++ b/include/boost/math/special_functions/detail/bessel_kn.hpp @@ -15,6 +15,7 @@ #include #include #include +#include #include #include @@ -60,6 +61,12 @@ BOOST_MATH_GPU_ENABLED T bessel_kn(int n, T x, const Policy& pol) { prev = bessel_k0(x); current = bessel_k1(x); + if ((prev < tools::min_value()) || (current < tools::min_value())) + { + T Iv, Kv; + bessel_ik(static_cast(n), x, &Iv, &Kv, need_k, pol); + return Kv; + } int k = 1; BOOST_MATH_ASSERT(k < n); T scale = 1; diff --git a/test/test_bessel_k.cpp b/test/test_bessel_k.cpp index 6053114585..e6af66f415 100644 --- a/test/test_bessel_k.cpp +++ b/test/test_bessel_k.cpp @@ -134,6 +134,30 @@ BOOST_AUTO_TEST_CASE( test_main ) test_bessel(0.1F, "float"); #endif test_bessel(0.1, "double"); + + { + // https://github.com/boostorg/math/issues/1229 + // Force evaluation at double precision so extended long double + // evaluation does not mask underflow in the starting K values. + typedef boost::math::policies::policy< + boost::math::policies::promote_double > no_promote_policy; + const no_promote_policy pol; + const double tol = 256 * std::numeric_limits::epsilon(); + + BOOST_CHECK_CLOSE_FRACTION(boost::math::cyl_bessel_k(1000, 747.0, pol), + 9.4914699277133192873957343540027064981e-66, tol); + BOOST_CHECK_CLOSE_FRACTION(boost::math::cyl_bessel_k(1000.25, 747.0, pol), + 1.2500811621640357229654772264173947351e-65, tol); + BOOST_CHECK_CLOSE_FRACTION(boost::math::cyl_bessel_k(-1000, 746.0, pol), + 5.0516775170486159420895532449059993901e-65, tol); + +#if !defined(BOOST_MATH_ENABLE_SYCL) + if (std::numeric_limits::has_denorm == std::denorm_present) + BOOST_CHECK_EQUAL(boost::math::cyl_bessel_k(100, 746.0, pol), + 8 * std::numeric_limits::denorm_min()); +#endif + } + #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS test_bessel(0.1L, "long double"); #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS