Refactor Bessel function code:

* Remove unused dead code. * Improve and centralise asymptotic selection. * Reorder algorithm selection in bessel_jy. * Allow use of integer algorithms for arbitrary order - they're no slower than Steeds method which is also O(n). [SVN r82816]
2026-01-19 04:22:09 +00:00 · 2013-02-11 12:12:50 +00:00
parent 6eb440db06
commit 76cc581bce
6 changed files with 558 additions and 674 deletions
--- a/include/boost/math/special_functions/bessel.hpp
+++ b/include/boost/math/special_functions/bessel.hpp
@@ -100,22 +100,6 @@ T cyl_bessel_j_imp(T v, T x, const bessel_no_int_tag& t, const Policy& pol)
            function,
            "Got x = %1%, but we need x >= 0", x, pol);
   }
-   if(x == 0)
-      return (v == 0) ? 1 : (v > 0) ? 0 : 
-         policies::raise_domain_error<T>(
-            function, 
-            "Got v = %1%, but require v >= 0 or a negative integer: the result would be complex.", v, pol);
-   
-   
-   if((v >= 0) && ((x < 1) || (v > x * x / 4) || (x < 5)))
-   {
-      //
-      // This series will actually converge rapidly for all small
-      // x - say up to x < 20 - but the first few terms are large
-      // and divergent which leads to large errors :-(
-      //
-      return bessel_j_small_z_series(v, x, pol);
-   }
   
   T j, y;
   bessel_jy(v, x, &j, &y, need_j, pol);
@@ -127,9 +111,11 @@ inline T cyl_bessel_j_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& p
 {
   BOOST_MATH_STD_USING  // ADL of std names.
   int ival = detail::iconv(v, pol);
-   if((abs(ival) < 200) && (0 == v - ival))
+   // If v is an integer, use the integer recursion
+   // method, both that and Steeds method are O(v):
+   if((0 == v - ival))
   {
-      return bessel_jn(ival/*iround(v, pol)*/, x, pol);
+      return bessel_jn(ival, x, pol);
   }
   return cyl_bessel_j_imp(v, x, bessel_no_int_tag(), pol);
 }
@@ -296,14 +282,13 @@ template <class T, class Policy>
 inline T cyl_neumann_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol)
 {
   BOOST_MATH_STD_USING
-   typedef typename bessel_asymptotic_tag<T, Policy>::type tag_type;

   BOOST_MATH_INSTRUMENT_VARIABLE(v);
   BOOST_MATH_INSTRUMENT_VARIABLE(x);

   if(floor(v) == v)
   {
-      if((fabs(x) > asymptotic_bessel_y_limit<T>(tag_type())) && (4 * fabs(x) * sqrt(sqrt(sqrt(14 * tools::epsilon<T>() * (fabs(x) + fabs(v) / 2) / (5120 * fabs(x))))) > fabs(v)))
+      if(asymptotic_bessel_large_x_limit(v, x))
      {
         T r = asymptotic_bessel_y_large_x_2(static_cast<T>(abs(v)), x);
         if((v < 0) && (itrunc(v, pol) & 1))
@@ -327,12 +312,11 @@ template <class T, class Policy>
 inline T cyl_neumann_imp(int v, T x, const bessel_int_tag&, const Policy& pol)
 {
   BOOST_MATH_STD_USING
-   typedef typename bessel_asymptotic_tag<T, Policy>::type tag_type;

   BOOST_MATH_INSTRUMENT_VARIABLE(v);
   BOOST_MATH_INSTRUMENT_VARIABLE(x);

-   if((fabs(x) > asymptotic_bessel_y_limit<T>(tag_type())) && (4 * fabs(x) * sqrt(sqrt(sqrt(14 * tools::epsilon<T>() * (fabs(x) + abs(v) / 2) / (5120 * x)))) > abs(v)))
+   if(asymptotic_bessel_large_x_limit(T(v), x))
   {
      T r = asymptotic_bessel_y_large_x_2(static_cast<T>(abs(v)), x);
      if((v < 0) && (v & 1))
--- a/include/boost/math/special_functions/detail/bessel_jn.hpp
+++ b/include/boost/math/special_functions/detail/bessel_jn.hpp
@@ -64,8 +64,7 @@ T bessel_jn(int n, T x, const Policy& pol)
        return static_cast<T>(0);
    }

-    typedef typename bessel_asymptotic_tag<T, Policy>::type tag_type;
-    if(fabs(x) > asymptotic_bessel_j_limit<T>(n, tag_type()))
+    if(asymptotic_bessel_large_x_limit(T(n), x))
      return factor * asymptotic_bessel_j_large_x_2<T>(n, x);

    BOOST_ASSERT(n > 1);
@@ -74,6 +73,7 @@ T bessel_jn(int n, T x, const Policy& pol)
    {
        prev = bessel_j0(x);
        current = bessel_j1(x);
+        policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)", n, pol);
        for (int k = 1; k < n; k++)
        {
            T fact = 2 * k / x;
@@ -91,7 +91,7 @@ T bessel_jn(int n, T x, const Policy& pol)
            current = value;
        }
    }
-    else if(x < 1)
+    else if((x < 1) || (n > x * x / 4) || (x < 5))
    {
       return factor * bessel_j_small_z_series(T(n), x, pol);
    }
@@ -102,6 +102,8 @@ T bessel_jn(int n, T x, const Policy& pol)
        boost::math::detail::CF1_jy(static_cast<T>(n), x, &fn, &s, pol);
        prev = fn;
        current = 1;
+        // Check recursion won't go on too far:
+        policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)", n, pol);
        for (int k = n; k > 0; k--)
        {
            T fact = 2 * k / x;
--- a/include/boost/math/special_functions/detail/bessel_jy.hpp
+++ b/include/boost/math/special_functions/detail/bessel_jy.hpp
--- a/include/boost/math/special_functions/detail/bessel_jy_asym.hpp
+++ b/include/boost/math/special_functions/detail/bessel_jy_asym.hpp
@@ -20,61 +20,6 @@

 namespace boost{ namespace math{ namespace detail{

-template <class T>
-inline T asymptotic_bessel_j_large_x_P(T v, T x)
-{
-   // A&S 9.2.9
-   T s = 1;
-   T mu = 4 * v * v;
-   T ez2 = 8 * x;
-   ez2 *= ez2;
-   s -= (mu-1) * (mu-9) / (2 * ez2);
-   s += (mu-1) * (mu-9) * (mu-25) * (mu - 49) / (24 * ez2 * ez2);
-   return s;
-}
-
-template <class T>
-inline T asymptotic_bessel_j_large_x_Q(T v, T x)
-{
-   // A&S 9.2.10
-   T s = 0;
-   T mu = 4 * v * v;
-   T ez = 8*x;
-   s += (mu-1) / ez;
-   s -= (mu-1) * (mu-9) * (mu-25) / (6 * ez*ez*ez);
-   return s;
-}
-
-template <class T>
-inline T asymptotic_bessel_j_large_x(T v, T x)
-{
-   // 
-   // See http://functions.wolfram.com/BesselAiryStruveFunctions/BesselJ/06/02/02/0001/
-   //
-   // Also A&S 9.2.5
-   //
-   BOOST_MATH_STD_USING // ADL of std names
-   T chi = fabs(x) - constants::pi<T>() * (2 * v + 1) / 4;
-   return sqrt(2 / (constants::pi<T>() * x))
-      * (asymptotic_bessel_j_large_x_P(v, x) * cos(chi) 
-         - asymptotic_bessel_j_large_x_Q(v, x) * sin(chi));
-}
-
-template <class T>
-inline T asymptotic_bessel_y_large_x(T v, T x)
-{
-   // 
-   // See http://functions.wolfram.com/BesselAiryStruveFunctions/BesselJ/06/02/02/0001/
-   //
-   // Also A&S 9.2.5
-   //
-   BOOST_MATH_STD_USING // ADL of std names
-   T chi = fabs(x) - constants::pi<T>() * (2 * v + 1) / 4;
-   return sqrt(2 / (constants::pi<T>() * x))
-      * (asymptotic_bessel_j_large_x_P(v, x) * sin(chi) 
-         - asymptotic_bessel_j_large_x_Q(v, x) * cos(chi));
-}
-
 template <class T>
 inline T asymptotic_bessel_amplitude(T v, T x)
 {
@@ -174,89 +119,21 @@ inline T asymptotic_bessel_j_large_x_2(T v, T x)
   return sin_phase * ampl;
 }

-//
-// Various limits for the J and Y asymptotics
-// (the asympotic expansions are safe to use if
-// x is less than the limit given).
-// We assume that if we don't use these expansions then the
-// error will likely be >100eps, so the limits given are chosen
-// to lead to < 100eps truncation error.
-//
 template <class T>
-inline T asymptotic_bessel_y_limit(const mpl::int_<0>&)
+inline bool asymptotic_bessel_large_x_limit(const T& v, const T& x)
 {
-   // default case:
   BOOST_MATH_STD_USING
-   return 2.25 / pow(100 * tools::epsilon<T>() / T(0.001f), T(0.2f));
-}
-template <class T>
-inline T asymptotic_bessel_y_limit(const mpl::int_<53>&)
-{
-   // double case:
-   return 304 /*780*/;
-}
-template <class T>
-inline T asymptotic_bessel_y_limit(const mpl::int_<64>&)
-{
-   // 80-bit extended-double case:
-   return 1552 /*3500*/;
-}
-template <class T>
-inline T asymptotic_bessel_y_limit(const mpl::int_<113>&)
-{
-   // 128-bit long double case:
-   return 1245243 /*3128000*/;
-}
-
-template <class T, class Policy>
-struct bessel_asymptotic_tag
-{
-   typedef typename policies::precision<T, Policy>::type precision_type;
-   typedef typename mpl::if_<
-      mpl::or_<
-         mpl::equal_to<precision_type, mpl::int_<0> >,
-         mpl::greater<precision_type, mpl::int_<113> > >,
-      mpl::int_<0>,
-      typename mpl::if_<
-         mpl::greater<precision_type, mpl::int_<64> >,
-         mpl::int_<113>,
-         typename mpl::if_<
-            mpl::greater<precision_type, mpl::int_<53> >,
-            mpl::int_<64>,
-            mpl::int_<53>
-         >::type
-      >::type
-   >::type type;
-};
-
-template <class T>
-inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<0>&)
-{
-   // default case:
-   BOOST_MATH_STD_USING
-   T v2 = (std::max)(T(3), T(v * v));
-   return v2 / pow(100 * tools::epsilon<T>() / T(2e-5f), T(0.17f));
-}
-template <class T>
-inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<53>&)
-{
-   // double case:
-   T v2 = (std::max)(T(3), T(v * v));
-   return v2 * 33 /*73*/;
-}
-template <class T>
-inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<64>&)
-{
-   // 80-bit extended-double case:
-   T v2 = (std::max)(T(3), T(v * v));
-   return v2 * 121 /*266*/;
-}
-template <class T>
-inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<113>&)
-{
-   // 128-bit long double case:
-   T v2 = (std::max)(T(3), T(v * v));
-   return v2 * 39154 /*85700*/;
+   //
+   // Determines if x is large enough compared to v to take the asymptotic
+   // forms above.  From A&S 9.2.28 we require: 
+   //    v < x * eps^1/8
+   // and from A&S 9.2.29 we require:
+   //    v^12/10 < 1.5 * x * eps^1/10
+   // using the former seems to work OK in practice with broadly similar 
+   // error rates either side of the divide for v < 10000.
+   // At double precision eps^1/8 ~= 0.01.
+   //
+   return (std::max)(T(fabs(v)), T(1)) < x * sqrt(tools::forth_root_epsilon<T>());
 }

 template <class T, class Policy>
--- a/include/boost/math/special_functions/detail/bessel_yn.hpp
+++ b/include/boost/math/special_functions/detail/bessel_yn.hpp
@@ -75,6 +75,7 @@ T bessel_yn(int n, T x, const Policy& pol)
       current = bessel_y1(x, pol);
       int k = 1;
       BOOST_ASSERT(k < n);
+       policies::check_series_iterations<T>("boost::math::bessel_y_n<%1%>(%1%,%1%)", n, pol);
       do
       {
           T fact = 2 * k / x;
--- a/test/test_bessel_j.hpp
+++ b/test/test_bessel_j.hpp
@@ -186,10 +186,11 @@ void test_bessel(T, const char* name)
        {{ SC_(1), T(10667654)/(1024*1024), SC_(1.24591331097191900488116495350277530373473085499043086981229e-7) }},
    }};

-    static const boost::array<boost::array<typename table_type<T>::type, 3>, 16> jn_data = {{
+    static const boost::array<boost::array<typename table_type<T>::type, 3>, 17> jn_data = {{
        // This first one is a modified test case from https://svn.boost.org/trac/boost/ticket/2733
        {{ SC_(-1), SC_(1.25), SC_(-0.510623260319880467069474837274910375352924050139633057168856) }},
        {{ SC_(2), SC_(0), SC_(0) }},
+        {{ SC_(-2), SC_(0), SC_(0) }},
        {{ SC_(2), SC_(1e-02), SC_(1.249989583365885362413250958437642113452e-05) }},
        {{ SC_(5), SC_(10), SC_(-0.2340615281867936404436949416457777864635) }},
        {{ SC_(5), SC_(-10), SC_(0.2340615281867936404436949416457777864635) }},
@@ -217,8 +218,9 @@ void test_bessel(T, const char* name)
    do_test_cyl_bessel_j_int<T>(j1_tricky, name, "Bessel J1: Mathworld Data (tricky cases) (Integer Version)");
    do_test_cyl_bessel_j_int<T>(jn_data, name, "Bessel JN: Mathworld Data (Integer Version)");

-    static const boost::array<boost::array<T, 3>, 20> jv_data = {{
+    static const boost::array<boost::array<T, 3>, 21> jv_data = {{
        //SC_(-2.4), {{ SC_(0), std::numeric_limits<T>::infinity() }},
+        {{ T(22.5), T(0), SC_(0) }},
        {{ T(2457)/1024, T(1)/1024, SC_(3.80739920118603335646474073457326714709615200130620574875292e-9) }},
        {{ SC_(5.5), T(3217)/1024, SC_(0.0281933076257506091621579544064767140470089107926550720453038) }},
        {{ SC_(-5.5), T(3217)/1024, SC_(-2.55820064470647911823175836997490971806135336759164272675969) }},
@@ -263,5 +265,12 @@ void test_bessel(T, const char* name)

 #include "sph_bessel_data.ipp"
    do_test_sph_bessel_j<T>(sph_bessel_data, name, "Bessel j: Random Data");
+
+    //
+    // Special cases that are errors:
+    //
+    BOOST_CHECK_THROW(boost::math::cyl_bessel_j(T(-2.5), T(0)), std::domain_error);
+    BOOST_CHECK_THROW(boost::math::cyl_bessel_j(T(-2.5), T(-2)), std::domain_error);
+    BOOST_CHECK_THROW(boost::math::cyl_bessel_j(T(2.5), T(-2)), std::domain_error);
 }