Update Bessel functions at infinity. (#1144)

* Update Bessel functions at infinity. Also sinc functions, and update tests. Fixes #1143. * Correct some test failures. * Yikes, correct missing if. * Temporary fix for multiprecision. REMOVE THIS.
boostorg · Jun 14, 2024 · dfc2934 · dfc2934
1 parent c5dc6c7
commit dfc2934
Show file tree

Hide file tree

Showing 10 changed files with 170 additions and 76 deletions.
diff --git 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
@@ -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;
+       // 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);
     }
-    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<T>() * 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<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)
        {
-          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<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
     }
-    else
-       Iv = std::numeric_limits<T>::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
@@ -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
@@ -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);
             }

diff --git a/include/boost/math/special_functions/sinhc.hpp b/include/boost/math/special_functions/sinhc.hpp
@@ -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>
-       inline typename tools::promote_args<T>::type sinhc_pi(T x)
+       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 detail::sinhc_pi_imp(static_cast<result_type>(x));
+          return policies::checked_narrowing_cast<T, Policy>(detail::sinhc_pi_imp(static_cast<result_type>(x), pol), "sinhc(%1%)");
        }
 
-       template <class T, class Policy>
-       inline typename tools::promote_args<T>::type sinhc_pi(T x, const Policy&)
+       template <class T>
+       inline typename tools::promote_args<T>::type sinhc_pi(T x)
        {
-          return boost::math::sinhc_pi(x);
+          typedef typename tools::promote_args<T>::type result_type;
+          return sinhc_pi(static_cast<result_type>(x), policies::policy<>());
        }
 
         template<typename T, template<typename> class U>

diff --git 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<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

diff --git 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<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));
+    }
 }
 
diff --git 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<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);

diff --git 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<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);

diff --git a/test/test_sinc.cpp b/test/test_sinc.cpp
@@ -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());
+   }
 }