Fix underflow in BesselK

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 <typename T, typename Policy>
-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<T>("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<T>())
        *Kv = exp(0.5f * log(pi<T>() / (2 * x)) - x - log(S));
     else
       *Kv = sqrt(pi<T>() / (2 * x)) * exp(-x) / S;
-    *Kv1 = *Kv * (0.5f + v + x + (v * v - 0.25f) * f) / x;
+    *Kv1 = *Kv * ratio;
+    if ((*Kv == 0) || (*Kv1 == 0))
+    {
+       *Kv_scaled = sqrt(pi<T>() / 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;
+    T log_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 == 0) || (Ku1 == 0));
+       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;
@@ -376,7 +385,10 @@ BOOST_MATH_GPU_ENABLED int bessel_ik(T v, T x, T* result_I, T* result_K, int kin
           if (!will_overflow && ((tools::max_value<T>() - fabs(prev)) / fact < fabs(current)))
           {
              prev /= current;
-             scale /= current;
+             if (use_scaled_k)
+                log_scale += log(current);
+             else
+                scale /= current;
              scale_sign *= ((boost::math::signbit)(current) ? -1 : 1);
              current = 1;
           }
@@ -432,7 +444,15 @@ 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<T>() * scale < Kv)
+    if (use_scaled_k)
+    {
+       const T log_Kv = log(Kv) + log_scale - x;
+       if (log_Kv > tools::log_max_value<T>())
+          *result_K = (org_kind & need_k) ? policies::raise_overflow_error<T>(function, nullptr, pol) : T(0);
+       else
+          *result_K = exp(log_Kv);
+    }
+    else if(tools::max_value<T>() * scale < Kv)
        *result_K = (org_kind & need_k) ? T(sign(Kv) * scale_sign * policies::raise_overflow_error<T>(function, nullptr, pol)) : T(0);
     else
       *result_K = Kv / scale;

include/boost/math/special_functions/detail/bessel_kn.hpp

@@ -15,6 +15,7 @@
 #include <boost/math/policies/error_handling.hpp>
 #include <boost/math/special_functions/detail/bessel_k0.hpp>
 #include <boost/math/special_functions/detail/bessel_k1.hpp>
+#include <boost/math/special_functions/detail/bessel_ik.hpp>
 #include <boost/math/special_functions/sign.hpp>
 #include <boost/math/policies/error_handling.hpp>
 
@@ -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 == 0) || (current == 0))
+       {
+          T Iv, Kv;
+          bessel_ik(static_cast<T>(n), x, &Iv, &Kv, need_k, pol);
+          return Kv;
+       }
        int k = 1;
        BOOST_MATH_ASSERT(k < n);
        T scale = 1;

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<false> > no_promote_policy;
+      const no_promote_policy pol;
+      const double tol = 256 * std::numeric_limits<double>::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<double>::has_denorm == std::denorm_present)
+         BOOST_CHECK_EQUAL(boost::math::cyl_bessel_k(100, 746.0, pol),
+            8 * std::numeric_limits<double>::denorm_min());
+#endif
+   }
+
 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
    test_bessel(0.1L, "long double");
 #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS