Improve the performance of the Bessel functions, and update docs.

[SVN r59274]

include/boost/math/special_functions/bessel.hpp

@@ -21,6 +21,7 @@
 #include <boost/math/special_functions/detail/bessel_i0.hpp>
 #include <boost/math/special_functions/detail/bessel_i1.hpp>
 #include <boost/math/special_functions/detail/bessel_kn.hpp>
+#include <boost/math/special_functions/detail/iconv.hpp>
 #include <boost/math/special_functions/sin_pi.hpp>
 #include <boost/math/special_functions/cos_pi.hpp>
 #include <boost/math/special_functions/sinc.hpp>
@@ -145,8 +146,13 @@ T cyl_bessel_j_imp(T v, T x, const bessel_no_int_tag& t, const Policy& pol)
             "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)))
+   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);
    }
    
@@ -159,13 +165,10 @@ template <class T, class Policy>
 inline T cyl_bessel_j_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol)
 {
    BOOST_MATH_STD_USING  // ADL of std names.
-   typedef typename bessel_asymptotic_tag<T, Policy>::type tag_type;
-   if((fabs(v) < 200) && (floor(v) == v))
+   int ival = detail::iconv(v, pol);
+   if((abs(ival) < 200) && (0 == v - ival))
    {
-      if(fabs(x) > asymptotic_bessel_j_limit<T>(v, tag_type()))
-         return asymptotic_bessel_j_large_x_2(v, x);
-      else
-         return bessel_jn(iround(v, pol), x, pol);
+      return bessel_jn(ival/*iround(v, pol)*/, x, pol);
    }
    return cyl_bessel_j_imp(v, x, bessel_no_int_tag(), pol);
 }
@@ -174,16 +177,7 @@ template <class T, class Policy>
 inline T cyl_bessel_j_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;
-   if(fabs(x) > asymptotic_bessel_j_limit<T>(abs(v), tag_type()))
-   {
-      T r = asymptotic_bessel_j_large_x_2(static_cast<T>(abs(v)), x);
-      if((v < 0) && (v & 1))
-         r = -r;
-      return r;
-   }
-   else
-      return bessel_jn(v, x, pol);
+   return bessel_jn(v, x, pol);
 }
 
 template <class T, class Policy>

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

@@ -13,6 +13,7 @@
 #include <boost/math/special_functions/detail/bessel_j0.hpp>
 #include <boost/math/special_functions/detail/bessel_j1.hpp>
 #include <boost/math/special_functions/detail/bessel_jy.hpp>
+#include <boost/math/special_functions/detail/bessel_jy_asym.hpp>
 
 // Bessel function of the first kind of integer order
 // J_n(z) is the minimal solution
@@ -57,6 +58,10 @@ 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()))
+      return factor * asymptotic_bessel_j_large_x_2<T>(n, x);
+
     BOOST_ASSERT(n > 1);
     if (n < abs(x))                         // forward recurrence
     {

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

@@ -16,7 +16,6 @@
 #include <boost/math/special_functions/hypot.hpp>
 #include <boost/math/special_functions/sin_pi.hpp>
 #include <boost/math/special_functions/cos_pi.hpp>
-#include <boost/math/special_functions/detail/simple_complex.hpp>
 #include <boost/math/special_functions/detail/bessel_jy_asym.hpp>
 #include <boost/math/constants/constants.hpp>
 #include <boost/math/policies/error_handling.hpp>
@@ -30,6 +29,45 @@ namespace boost { namespace math {
 
 namespace detail {
 
+//
+// Simultaneous calculation of A&S 9.2.9 and 9.2.10
+// for use in A&S 9.2.5 and 9.2.6.
+// This series is quick to evaluate, but divergent unless
+// x is very large, in fact it's pretty hard to figure out
+// with any degree of precision when this series actually 
+// *will* converge!!  Consequently, we may just have to
+// try it and see...
+//
+template <class T, class Policy>
+bool hankel_PQ(T v, T x, T* p, T* q, const Policy& )
+{
+   BOOST_MATH_STD_USING
+   T tolerance = 2 * policies::get_epsilon<T, Policy>();
+   *p = 1;
+   *q = 0;
+   T k = 1;
+   T z8 = 8 * x;
+   T sq = 1;
+   T mu = 4 * v * v;
+   T term = 1;
+   bool ok = true;
+   do
+   {
+      term *= (mu - sq * sq) / (k * z8);
+      *q += term;
+      k += 1;
+      sq += 2;
+      T mult = (sq * sq - mu) / (k * z8);
+      ok = fabs(mult) < 0.5f;
+      term *= mult;
+      *p += term;
+      k += 1;
+      sq += 2;
+   }
+   while((fabs(term) > tolerance * *p) && ok);
+   return ok;
+}
+
 // Calculate Y(v, x) and Y(v+1, x) by Temme's method, see
 // Temme, Journal of Computational Physics, vol 21, 343 (1976)
 template <typename T, typename Policy>
@@ -76,7 +114,7 @@ int temme_jy(T v, T x, T* Y, T* Y1, const Policy& pol)
     T coef_mult = -x * x / 4;
 
     // series summation
-    tolerance = tools::epsilon<T>();
+    tolerance = policies::get_epsilon<T, Policy>();
     for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)
     {
         f = (k * f + p + q) / (k*k - v2);
@@ -115,7 +153,7 @@ int CF1_jy(T v, T x, T* fv, int* sign, const Policy& pol)
 
     // modified Lentz's method, see
     // Lentz, Applied Optics, vol 15, 668 (1976)
-    tolerance = 2 * tools::epsilon<T>();
+    tolerance = 2 * policies::get_epsilon<T, Policy>();;
     tiny = sqrt(tools::min_value<T>());
     C = f = tiny;                           // b0 = 0, replace with tiny
     D = 0;
@@ -140,59 +178,87 @@ int CF1_jy(T v, T x, T* fv, int* sign, const Policy& pol)
 
     return 0;
 }
-
-template <class T>
-struct complex_trait
-{
-   typedef typename mpl::if_<is_floating_point<T>,
-      std::complex<T>, sc::simple_complex<T> >::type type;
-};
-
-// Evaluate continued fraction p + iq = (J' + iY') / (J + iY), see
-// Press et al, Numerical Recipes in C, 2nd edition, 1992
+//
+// This algorithm was originally written by Xiaogang Zhang
+// using std::complex to perform the complex arithmetic.
+// However, that turns out to 10x or more slower than using
+// all real-valued arithmetic, so it's been rewritten using
+// real values only.
+//
 template <typename T, typename Policy>
 int CF2_jy(T v, T x, T* p, T* q, const Policy& pol)
 {
-    BOOST_MATH_STD_USING
+   BOOST_MATH_STD_USING
 
-    typedef typename complex_trait<T>::type complex_type;
+   T Cr, Ci, Dr, Di, fr, fi, a, br, bi, delta_r, delta_i, temp;
+   T tiny;
+   unsigned long k;
 
-    complex_type C, D, f, a, b, delta, one(1);
-    T tiny, zero(0);
-    unsigned long k;
+   // |x| >= |v|, CF2_jy converges rapidly
+   // |x| -> 0, CF2_jy fails to converge
+   BOOST_ASSERT(fabs(x) > 1);
 
-    // |x| >= |v|, CF2_jy converges rapidly
-    // |x| -> 0, CF2_jy fails to converge
-    BOOST_ASSERT(fabs(x) > 1);
+   // modified Lentz's method, complex numbers involved, see
+   // Lentz, Applied Optics, vol 15, 668 (1976)
+   T tolerance = 2 * policies::get_epsilon<T, Policy>();
+   tiny = sqrt(tools::min_value<T>());
+   Cr = fr = -0.5f / x;
+   Ci = fi = 1;
+   //Dr = Di = 0;
+   T v2 = v * v;
+   a = (0.25f - v2) / x; // Note complex this one time only!
+   br = 2 * x;
+   bi = 2;
+   temp = Cr * Cr + 1;
+   Ci = bi + a * Cr / temp;
+   Cr = br + a / temp;
+   Dr = br;
+   Di = bi;
+   //std::cout << "C = " << Cr << " " << Ci << std::endl;
+   //std::cout << "D = " << Dr << " " << Di << std::endl;
+   if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }
+   if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }
+   temp = Dr * Dr + Di * Di;
+   Dr = Dr / temp;
+   Di = -Di / temp;
+   delta_r = Cr * Dr - Ci * Di;
+   delta_i = Ci * Dr + Cr * Di;
+   temp = fr;
+   fr = temp * delta_r - fi * delta_i;
+   fi = temp * delta_i + fi * delta_r;
+   //std::cout << fr << " " << fi << std::endl;
+   for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++)
+   {
+      a = k - 0.5f;
+      a *= a;
+      a -= v2;
+      bi += 2;
+      temp = Cr * Cr + Ci * Ci;
+      Cr = br + a * Cr / temp;
+      Ci = bi - a * Ci / temp;
+      Dr = br + a * Dr;
+      Di = bi + a * Di;
+      //std::cout << "C = " << Cr << " " << Ci << std::endl;
+      //std::cout << "D = " << Dr << " " << Di << std::endl;
+      if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }
+      if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }
+      temp = Dr * Dr + Di * Di;
+      Dr = Dr / temp;
+      Di = -Di / temp;
+      delta_r = Cr * Dr - Ci * Di;
+      delta_i = Ci * Dr + Cr * Di;
+      temp = fr;
+      fr = temp * delta_r - fi * delta_i;
+      fi = temp * delta_i + fi * delta_r;
+      if (fabs(delta_r - 1) + fabs(delta_i) < tolerance)
+         break;
+      //std::cout << fr << " " << fi << std::endl;
+   }
+   policies::check_series_iterations("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol);
+   *p = fr;
+   *q = fi;
 
-    // modified Lentz's method, complex numbers involved, see
-    // Lentz, Applied Optics, vol 15, 668 (1976)
-    T tolerance = 2 * tools::epsilon<T>();
-    tiny = sqrt(tools::min_value<T>());
-    C = f = complex_type(-0.5f/x, 1);
-    D = 0;
-    for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)
-    {
-        a = (k - 0.5f)*(k - 0.5f) - v*v;
-        if (k == 1)
-        {
-            a *= complex_type(T(0), 1/x);
-        }
-        b = complex_type(2*x, T(2*k));
-        C = b + a / C;
-        D = b + a * D;
-        if (C == zero) { C = tiny; }
-        if (D == zero) { D = tiny; }
-        D = one / D;
-        delta = C * D;
-        f *= delta;
-        if (abs(delta - one) < tolerance) { break; }
-    }
-    policies::check_series_iterations("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol);
-    *p = real(f);
-    *q = imag(f);
-
-    return 0;
+   return 0;
 }
 
 enum
@@ -237,7 +303,22 @@ int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol)
 
     // x is positive until reflection
     W = T(2) / (x * pi<T>());               // Wronskian
-    if (x <= 2)                           // x in (0, 2]
+    if((x > 8) && (x < 1000) && hankel_PQ(v, x, &p, &q, pol))
+    {
+       //
+       // Hankel approximation: note that this method works best when x 
+       // is large, but in that case we end up calculating sines and cosines
+       // of large values, with horrendous resulting accuracy.  It is fast though
+       // when it works....
+       //
+       T chi = x - fmod(T(v / 2 + 0.25f), T(2)) * boost::math::constants::pi<T>();
+       T sc = sin(chi);
+       T cc = cos(chi);
+       chi = sqrt(2 / (boost::math::constants::pi<T>() * x));
+       Yv = chi * (p * sc + q * cc);
+       Jv = chi * (p * cc - q * sc);
+    }
+    else if (x <= 2)                           // x in (0, 2]
     {
         if(temme_jy(u, x, &Yu, &Yu1, pol))             // Temme series
         {

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

@@ -0,0 +1,41 @@
+//  Copyright (c) 2009 John Maddock
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_ICONV_HPP
+#define BOOST_MATH_ICONV_HPP
+
+#ifdef _MSC_VER
+#pragma once
+#endif
+
+#include <boost/math/special_functions/round.hpp>
+#include <boost/type_traits/is_convertible.hpp>
+
+namespace boost { namespace math { namespace detail{
+
+template <class T, class Policy>
+inline int iconv_imp(T v, Policy const&, mpl::true_ const&)
+{
+   return static_cast<int>(v);
+}
+
+template <class T, class Policy>
+inline int iconv_imp(T v, Policy const& pol, mpl::false_ const&)
+{
+   return boost::math::iround(v);
+}
+
+template <class T, class Policy>
+inline int iconv(T v, Policy const& pol)
+{
+   typedef typename boost::is_convertible<T, int>::type tag_type;
+   return iconv_imp(v, pol, tag_type());
+}
+
+
+}}} // namespaces
+
+#endif // BOOST_MATH_ICONV_HPP
+

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

@@ -1,172 +0,0 @@
-//  Copyright (c) 2007 John Maddock
-//  Use, modification and distribution are subject to the
-//  Boost Software License, Version 1.0. (See accompanying file
-//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
-
-#ifndef BOOST_MATH_SF_DETAIL_SIMPLE_COMPLEX_HPP
-#define BOOST_MATH_SF_DETAIL_SIMPLE_COMPLEX_HPP
-
-#ifdef _MSC_VER
-#pragma once
-#endif
-
-namespace boost{ namespace math{ namespace detail{ namespace sc{
-
-template <class T>
-class simple_complex
-{
-public:
-   simple_complex() : r(0), i(0) {}
-   simple_complex(T a) : r(a) {}
-   template <class U>
-   simple_complex(U a) : r(a) {}
-   simple_complex(T a, T b) : r(a), i(b) {}
-
-   simple_complex& operator += (const simple_complex& o)
-   {
-      r += o.r;
-      i += o.i;
-      return *this;
-   }
-   simple_complex& operator -= (const simple_complex& o)
-   {
-      r -= o.r;
-      i -= o.i;
-      return *this;
-   }
-   simple_complex& operator *= (const simple_complex& o)
-   { 
-      T lr = r * o.r - i * o.i;
-      T li = i * o.r + r * o.i;
-      r = lr;
-      i = li;
-      return *this;
-   }
-   simple_complex& operator /= (const simple_complex& o)
-   { 
-      BOOST_MATH_STD_USING
-      T lr;
-      T li;
-      if(fabs(o.r) > fabs(o.i))
-      {
-         T rat = o.i / o.r;
-         lr = r + i * rat;
-         li = i - r * rat;
-         rat = o.r + o.i * rat;
-         lr /= rat;
-         li /= rat;
-      }
-      else
-      {
-         T rat = o.r / o.i;
-         lr = i + r * rat;
-         li = i * rat - r;
-         rat = o.r * rat + o.i;
-         lr /= rat;
-         li /= rat;
-      }
-      r = lr;
-      i = li;
-      return *this;
-   }
-   bool operator == (const simple_complex& o)
-   {
-      return (r == o.r) && (i == o.i);
-   }
-   bool operator != (const simple_complex& o)
-   {
-      return !((r == o.r) && (i == o.i));
-   }
-   bool operator == (const T& o)
-   {
-      return (r == o) && (i == 0);
-   }
-   simple_complex& operator += (const T& o)
-   {
-      r += o;
-      return *this;
-   }
-   simple_complex& operator -= (const T& o)
-   {
-      r -= o;
-      return *this;
-   }
-   simple_complex& operator *= (const T& o)
-   { 
-      r *= o;
-      i *= o;
-      return *this;
-   }
-   simple_complex& operator /= (const T& o)
-   { 
-      r /= o;
-      i /= o;
-      return *this;
-   }
-   T real()const
-   {
-      return r;
-   }
-   T imag()const
-   {
-      return i;
-   }
-private:
-   T r, i;
-};
-
-template <class T>
-inline simple_complex<T> operator+(const simple_complex<T>& a, const simple_complex<T>& b)
-{
-   simple_complex<T> result(a);
-   result += b;
-   return result;
-}
-
-template <class T>
-inline simple_complex<T> operator-(const simple_complex<T>& a, const simple_complex<T>& b)
-{
-   simple_complex<T> result(a);
-   result -= b;
-   return result;
-}
-
-template <class T>
-inline simple_complex<T> operator*(const simple_complex<T>& a, const simple_complex<T>& b)
-{
-   simple_complex<T> result(a);
-   result *= b;
-   return result;
-}
-
-template <class T>
-inline simple_complex<T> operator/(const simple_complex<T>& a, const simple_complex<T>& b)
-{
-   simple_complex<T> result(a);
-   result /= b;
-   return result;
-}
-
-template <class T>
-inline T real(const simple_complex<T>& c)
-{
-   return c.real();
-}
-
-template <class T>
-inline T imag(const simple_complex<T>& c)
-{
-   return c.imag();
-}
-
-template <class T>
-inline T abs(const simple_complex<T>& c)
-{
-   return hypot(c.real(), c.imag());
-}
-
-}}}} // namespace
-
-#endif
-
-