[Boost-commit] svn:boost r82816 - in trunk: boost/math/special_functions boost/math/special_functions/detail libs/math/test

Subject: [Boost-commit] svn:boost r82816 - in trunk: boost/math/special_functions boost/math/special_functions/detail libs/math/test
From: john_at_[hidden]
Date: 2013-02-11 07:12:53

Author: johnmaddock
Date: 2013-02-11 07:12:50 EST (Mon, 11 Feb 2013)
New Revision: 82816
URL: http://svn.boost.org/trac/boost/changeset/82816

Log:
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).
Text files modified: 
   trunk/boost/math/special_functions/bessel.hpp                |    28                                         
   trunk/boost/math/special_functions/detail/bessel_jn.hpp      |     8                                         
   trunk/boost/math/special_functions/detail/bessel_jy.hpp      |  1037 ++++++++++++++++++++------------------- 
   trunk/boost/math/special_functions/detail/bessel_jy_asym.hpp |   147 -----                                   
   trunk/boost/math/special_functions/detail/bessel_yn.hpp      |     1                                         
   trunk/libs/math/test/test_bessel_j.hpp                       |    13                                         
   6 files changed, 559 insertions(+), 675 deletions(-)

Modified: trunk/boost/math/special_functions/bessel.hpp
==============================================================================
--- trunk/boost/math/special_functions/bessel.hpp	(original)
+++ trunk/boost/math/special_functions/bessel.hpp	2013-02-11 07:12:50 EST (Mon, 11 Feb 2013)
@@ -100,22 +100,6 @@
             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 @@
 {
    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 @@
 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 @@
 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))

Modified: trunk/boost/math/special_functions/detail/bessel_jn.hpp
==============================================================================
--- trunk/boost/math/special_functions/detail/bessel_jn.hpp	(original)
+++ trunk/boost/math/special_functions/detail/bessel_jn.hpp	2013-02-11 07:12:50 EST (Mon, 11 Feb 2013)
@@ -64,8 +64,7 @@
         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 @@
     {
         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 @@
             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 @@
         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;

Modified: trunk/boost/math/special_functions/detail/bessel_jy.hpp
==============================================================================
--- trunk/boost/math/special_functions/detail/bessel_jy.hpp	(original)
+++ trunk/boost/math/special_functions/detail/bessel_jy.hpp	2013-02-11 07:12:50 EST (Mon, 11 Feb 2013)
@@ -28,448 +28,458 @@
 
 namespace boost { namespace math {
 
-namespace detail {
+   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>
-int temme_jy(T v, T x, T* Y, T* Y1, const Policy& pol)
-{
-    T g, h, p, q, f, coef, sum, sum1, tolerance;
-    T a, d, e, sigma;
-    unsigned long k;
-
-    BOOST_MATH_STD_USING
-    using namespace boost::math::tools;
-    using namespace boost::math::constants;
-
-    BOOST_ASSERT(fabs(v) <= 0.5f);  // precondition for using this routine
-
-    T gp = boost::math::tgamma1pm1(v, pol);
-    T gm = boost::math::tgamma1pm1(-v, pol);
-    T spv = boost::math::sin_pi(v, pol);
-    T spv2 = boost::math::sin_pi(v/2, pol);
-    T xp = pow(x/2, v);
-
-    a = log(x / 2);
-    sigma = -a * v;
-    d = abs(sigma) < tools::epsilon<T>() ?
-        T(1) : sinh(sigma) / sigma;
-    e = abs(v) < tools::epsilon<T>() ? T(v*pi<T>()*pi<T>() / 2)
-        : T(2 * spv2 * spv2 / v);
-
-    T g1 = (v == 0) ? T(-euler<T>()) : T((gp - gm) / ((1 + gp) * (1 + gm) * 2 * v));
-    T g2 = (2 + gp + gm) / ((1 + gp) * (1 + gm) * 2);
-    T vspv = (fabs(v) < tools::epsilon<T>()) ? T(1/constants::pi<T>()) : T(v / spv);
-    f = (g1 * cosh(sigma) - g2 * a * d) * 2 * vspv;
-
-    p = vspv / (xp * (1 + gm));
-    q = vspv * xp / (1 + gp);
-
-    g = f + e * q;
-    h = p;
-    coef = 1;
-    sum = coef * g;
-    sum1 = coef * h;
-
-    T v2 = v * v;
-    T coef_mult = -x * x / 4;
-
-    // series summation
-    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);
-        p /= k - v;
-        q /= k + v;
-        g = f + e * q;
-        h = p - k * g;
-        coef *= coef_mult / k;
-        sum += coef * g;
-        sum1 += coef * h;
-        if (abs(coef * g) < abs(sum) * tolerance) 
-        { 
-           break; 
-        }
-    }
-    policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in temme_jy", k, pol);
-    *Y = -sum;
-    *Y1 = -2 * sum1 / x;
-
-    return 0;
-}
-
-// Evaluate continued fraction fv = J_(v+1) / J_v, see
-// Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73
-template <typename T, typename Policy>
-int CF1_jy(T v, T x, T* fv, int* sign, const Policy& pol)
-{
-    T C, D, f, a, b, delta, tiny, tolerance;
-    unsigned long k;
-    int s = 1;
-
-    BOOST_MATH_STD_USING
-
-    // |x| <= |v|, CF1_jy converges rapidly
-    // |x| > |v|, CF1_jy needs O(|x|) iterations to converge
-
-    // modified Lentz's method, see
-    // Lentz, Applied Optics, vol 15, 668 (1976)
-    tolerance = 2 * policies::get_epsilon<T, Policy>();;
-    tiny = sqrt(tools::min_value<T>());
-    C = f = tiny;                           // b0 = 0, replace with tiny
-    D = 0;
-    for (k = 1; k < policies::get_max_series_iterations<Policy>() * 100; k++)
-    {
-        a = -1;
-        b = 2 * (v + k) / x;
-        C = b + a / C;
-        D = b + a * D;
-        if (C == 0) { C = tiny; }
-        if (D == 0) { D = tiny; }
-        D = 1 / D;
-        delta = C * D;
-        f *= delta;
-        if (D < 0) { s = -s; }
-        if (abs(delta - 1) < tolerance) 
-        { break; }
-    }
-    policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF1_jy", k / 100, pol);
-    *fv = -f;
-    *sign = s;                              // sign of denominator
-
-    return 0;
-}
-//
-// 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
-
-   T Cr, Ci, Dr, Di, fr, fi, a, br, bi, delta_r, delta_i, temp;
-   T tiny;
-   unsigned long k;
-
-   // |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;
-   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;
-   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;
-      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;
-   }
-   policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol);
-   *p = fr;
-   *q = fi;
-
-   return 0;
-}
-
-static const int need_j = 1;
-static const int need_y = 2;
-
-// Compute J(v, x) and Y(v, x) simultaneously by Steed's method, see
-// Barnett et al, Computer Physics Communications, vol 8, 377 (1974)
-template <typename T, typename Policy>
-int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol)
-{
-    BOOST_ASSERT(x >= 0);
-
-    T u, Jv, Ju, Yv, Yv1, Yu, Yu1(0), fv, fu;
-    T W, p, q, gamma, current, prev, next;
-    bool reflect = false;
-    unsigned n, k;
-    int s;
-    int org_kind = kind;
-    T cp = 0;
-    T sp = 0;
-
-    static const char* function = "boost::math::bessel_jy<%1%>(%1%,%1%)";
-
-    BOOST_MATH_STD_USING
-    using namespace boost::math::tools;
-    using namespace boost::math::constants;
-
-    if (v < 0)
-    {
-        reflect = true;
-        v = -v;                             // v is non-negative from here
-        kind = need_j|need_y;               // need both for reflection formula
-    }
-    n = iround(v, pol);
-    u = v - n;                              // -1/2 <= u < 1/2
-
-    if(reflect)
-    {
-        T z = (u + n % 2);
-        cp = boost::math::cos_pi(z, pol);
-        sp = boost::math::sin_pi(z, pol);
-    }
-
-    if (x == 0)
-    {
-       *J = *Y = policies::raise_overflow_error<T>(
-          function, 0, pol);
-       return 1;
-    }
-
-    // x is positive until reflection
-    W = T(2) / (x * pi<T>());               // Wronskian
-    T Yv_scale = 1;
-    if((x > 8) && 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....
-       //
-       // Normally we calculate sin/cos(chi) where:
-       //
-       // chi = x - fmod(T(v / 2 + 0.25f), T(2)) * boost::math::constants::pi<T>();
-       //
-       // But this introduces large errors, so use sin/cos addition formulae to
-       // improve accuracy:
-       //
-       T mod_v = fmod(T(v / 2 + 0.25f), T(2));
-       T sx = sin(x);
-       T cx = cos(x);
-       T sv = sin_pi(mod_v);
-       T cv = cos_pi(mod_v);
-
-       T sc = sx * cv - sv * cx; // == sin(chi);
-       T cc = cx * cv + sx * sv; // == cos(chi);
-       T chi = boost::math::constants::root_two<T>() / (boost::math::constants::root_pi<T>() * sqrt(x)); //sqrt(2 / (boost::math::constants::pi<T>() * x));
-       Yv = chi * (p * sc + q * cc);
-       Jv = chi * (p * cc - q * sc);
-    }
-    else if((x < 1) && (u != 0) && (log(policies::get_epsilon<T, Policy>() / 2) > v * log((x/2) * (x/2) / v)))
-    {
-       // Evaluate using series representations.
-       // This is particularly important for x << v as in this
-       // area temme_jy may be slow to converge, if it converges at all.
-       // Requires x is not an integer.
-       if(kind&need_j)
-          Jv = bessel_j_small_z_series(v, x, pol);
-       else
-          Jv = std::numeric_limits<T>::quiet_NaN();
-       if((org_kind&need_y && (!reflect || (cp != 0))) 
-          || (org_kind & need_j && (reflect && (sp != 0))))
-       {
-          // Only calculate if we need it, and if the reflection formula will actually use it:
-          Yv = bessel_y_small_z_series(v, x, &Yv_scale, pol);
-       }
-       else
-          Yv = std::numeric_limits<T>::quiet_NaN();
-    }
-    else if((u == 0) && (x < policies::get_epsilon<T, Policy>()))
-    {
-       // Truncated series evaluation for small x and v an integer,
-       // much quicker in this area than temme_jy below.
-       if(kind&need_j)
-          Jv = bessel_j_small_z_series(v, x, pol);
-       else
-          Jv = std::numeric_limits<T>::quiet_NaN();
-       if((org_kind&need_y && (!reflect || (cp != 0))) 
-          || (org_kind & need_j && (reflect && (sp != 0))))
-       {
-          // Only calculate if we need it, and if the reflection formula will actually use it:
-          Yv = bessel_yn_small_z(n, x, &Yv_scale, pol);
-       }
-       else
-          Yv = std::numeric_limits<T>::quiet_NaN();
-    }
-    else if (x <= 2)                           // x in (0, 2]
-    {
-        if(temme_jy(u, x, &Yu, &Yu1, pol))             // Temme series
-        {
-           // domain error:
-           *J = *Y = Yu;
-           return 1;
-        }
-        prev = Yu;
-        current = Yu1;
-        T scale = 1;
-        for (k = 1; k <= n; k++)            // forward recurrence for Y
-        {
-            T fact = 2 * (u + k) / x;
-            if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
-            {
-               scale /= current;
-               prev /= current;
-               current = 1;
-            }
-            next = fact * current - prev;
-            prev = current;
-            current = next;
-        }
-         Yv = prev;
-         Yv1 = current;
-         if(kind&need_j)
+      //
+      // 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
          {
-            CF1_jy(v, x, &fv, &s, pol);                 // continued fraction CF1_jy
-            Jv = scale * W / (Yv * fv - Yv1);           // Wronskian relation
+            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;
          }
-         else
-            Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
-         Yv_scale = scale;
-    }
-    else                                    // x in (2, \infty)
-    {
-        // Get Y(u, x):
-        // define tag type that will dispatch to right limits:
-        typedef typename bessel_asymptotic_tag<T, Policy>::type tag_type;
-
-        T lim, ratio;
-        switch(kind)
-        {
-        case need_j:
-           lim = asymptotic_bessel_j_limit<T>(v, tag_type());
-           break;
-        case need_y:
-           lim = asymptotic_bessel_y_limit<T>(tag_type());
-           break;
-        default:
-           lim = (std::max)(
-              asymptotic_bessel_j_limit<T>(v, tag_type()),
-              asymptotic_bessel_y_limit<T>(tag_type()));
-           break;
-        }
-        if(x > lim)
-        {
-           if(kind&need_y)
-           {
-              Yu = asymptotic_bessel_y_large_x_2(u, x);
-              Yu1 = asymptotic_bessel_y_large_x_2(T(u + 1), x);
-           }
-           else
-              Yu = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
-           if(kind&need_j)
-           {
-              Jv = asymptotic_bessel_j_large_x_2(v, x);
-           }
-           else
-              Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
-        }
-        else
-        {
-           CF1_jy(v, x, &fv, &s, pol);
-           // tiny initial value to prevent overflow
-           T init = sqrt(tools::min_value<T>());
-           prev = fv * s * init;
-           current = s * init;
-           if(v < max_factorial<T>::value)
-           {
-              for (k = n; k > 0; k--)             // backward recurrence for J
-              {
+         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>
+      int temme_jy(T v, T x, T* Y, T* Y1, const Policy& pol)
+      {
+         T g, h, p, q, f, coef, sum, sum1, tolerance;
+         T a, d, e, sigma;
+         unsigned long k;
+
+         BOOST_MATH_STD_USING
+            using namespace boost::math::tools;
+         using namespace boost::math::constants;
+
+         BOOST_ASSERT(fabs(v) <= 0.5f);  // precondition for using this routine
+
+         T gp = boost::math::tgamma1pm1(v, pol);
+         T gm = boost::math::tgamma1pm1(-v, pol);
+         T spv = boost::math::sin_pi(v, pol);
+         T spv2 = boost::math::sin_pi(v/2, pol);
+         T xp = pow(x/2, v);
+
+         a = log(x / 2);
+         sigma = -a * v;
+         d = abs(sigma) < tools::epsilon<T>() ?
+            T(1) : sinh(sigma) / sigma;
+         e = abs(v) < tools::epsilon<T>() ? T(v*pi<T>()*pi<T>() / 2)
+            : T(2 * spv2 * spv2 / v);
+
+         T g1 = (v == 0) ? T(-euler<T>()) : T((gp - gm) / ((1 + gp) * (1 + gm) * 2 * v));
+         T g2 = (2 + gp + gm) / ((1 + gp) * (1 + gm) * 2);
+         T vspv = (fabs(v) < tools::epsilon<T>()) ? T(1/constants::pi<T>()) : T(v / spv);
+         f = (g1 * cosh(sigma) - g2 * a * d) * 2 * vspv;
+
+         p = vspv / (xp * (1 + gm));
+         q = vspv * xp / (1 + gp);
+
+         g = f + e * q;
+         h = p;
+         coef = 1;
+         sum = coef * g;
+         sum1 = coef * h;
+
+         T v2 = v * v;
+         T coef_mult = -x * x / 4;
+
+         // series summation
+         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);
+            p /= k - v;
+            q /= k + v;
+            g = f + e * q;
+            h = p - k * g;
+            coef *= coef_mult / k;
+            sum += coef * g;
+            sum1 += coef * h;
+            if (abs(coef * g) < abs(sum) * tolerance) 
+            { 
+               break; 
+            }
+         }
+         policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in temme_jy", k, pol);
+         *Y = -sum;
+         *Y1 = -2 * sum1 / x;
+
+         return 0;
+      }
+
+      // Evaluate continued fraction fv = J_(v+1) / J_v, see
+      // Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73
+      template <typename T, typename Policy>
+      int CF1_jy(T v, T x, T* fv, int* sign, const Policy& pol)
+      {
+         T C, D, f, a, b, delta, tiny, tolerance;
+         unsigned long k;
+         int s = 1;
+
+         BOOST_MATH_STD_USING
+
+            // |x| <= |v|, CF1_jy converges rapidly
+            // |x| > |v|, CF1_jy needs O(|x|) iterations to converge
+
+            // modified Lentz's method, see
+            // Lentz, Applied Optics, vol 15, 668 (1976)
+            tolerance = 2 * policies::get_epsilon<T, Policy>();;
+         tiny = sqrt(tools::min_value<T>());
+         C = f = tiny;                           // b0 = 0, replace with tiny
+         D = 0;
+         for (k = 1; k < policies::get_max_series_iterations<Policy>() * 100; k++)
+         {
+            a = -1;
+            b = 2 * (v + k) / x;
+            C = b + a / C;
+            D = b + a * D;
+            if (C == 0) { C = tiny; }
+            if (D == 0) { D = tiny; }
+            D = 1 / D;
+            delta = C * D;
+            f *= delta;
+            if (D < 0) { s = -s; }
+            if (abs(delta - 1) < tolerance) 
+            { break; }
+         }
+         policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF1_jy", k / 100, pol);
+         *fv = -f;
+         *sign = s;                              // sign of denominator
+
+         return 0;
+      }
+      //
+      // 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
+
+            T Cr, Ci, Dr, Di, fr, fi, a, br, bi, delta_r, delta_i, temp;
+         T tiny;
+         unsigned long k;
+
+         // |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;
+         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;
+         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;
+            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;
+         }
+         policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol);
+         *p = fr;
+         *q = fi;
+
+         return 0;
+      }
+
+      static const int need_j = 1;
+      static const int need_y = 2;
+
+      // Compute J(v, x) and Y(v, x) simultaneously by Steed's method, see
+      // Barnett et al, Computer Physics Communications, vol 8, 377 (1974)
+      template <typename T, typename Policy>
+      int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol)
+      {
+         BOOST_ASSERT(x >= 0);
+
+         T u, Jv, Ju, Yv, Yv1, Yu, Yu1(0), fv, fu;
+         T W, p, q, gamma, current, prev, next;
+         bool reflect = false;
+         unsigned n, k;
+         int s;
+         int org_kind = kind;
+         T cp = 0;
+         T sp = 0;
+
+         static const char* function = "boost::math::bessel_jy<%1%>(%1%,%1%)";
+
+         BOOST_MATH_STD_USING
+            using namespace boost::math::tools;
+         using namespace boost::math::constants;
+
+         if (v < 0)
+         {
+            reflect = true;
+            v = -v;                             // v is non-negative from here
+         }
+         if(v > static_cast<T>((std::numeric_limits<int>::max)()))
+            policies::raise_evaluation_error<T>(function, "Order of Bessel function is too large to evaluate: got %1%", v, pol);
+         n = iround(v, pol);
+         u = v - n;                              // -1/2 <= u < 1/2
+
+         if(reflect)
+         {
+            T z = (u + n % 2);
+            cp = boost::math::cos_pi(z, pol);
+            sp = boost::math::sin_pi(z, pol);
+            if(u != 0)
+               kind = need_j|need_y;               // need both for reflection formula
+         }
+
+         if(x == 0)
+         {
+            if(v == 0)
+               *J = 1;
+            else if((u == 0) || !reflect)
+               *J = 0;
+            else if(kind & need_j)
+               *J = policies::raise_domain_error<T>(function, "Value of Bessel J_v(x) is complex-infinity at %1%", x, pol); // complex infinity
+            else
+               *J = std::numeric_limits<T>::quiet_NaN();  // any value will do, not using J.
+
+            if((kind & need_y) == 0)
+               *Y = std::numeric_limits<T>::quiet_NaN();  // any value will do, not using Y.
+            else if(v == 0)
+               *Y = -policies::raise_overflow_error<T>(function, 0, pol);
+            else
+               *Y = policies::raise_domain_error<T>(function, "Value of Bessel Y_v(x) is complex-infinity at %1%", x, pol); // complex infinity
+            return 1;
+         }
+
+         // x is positive until reflection
+         W = T(2) / (x * pi<T>());               // Wronskian
+         T Yv_scale = 1;
+         if(((kind & need_y) == 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 :-(
+            //
+            Jv = bessel_j_small_z_series(v, x, pol);
+            Yv = std::numeric_limits<T>::quiet_NaN();
+         }
+         else if((x < 1) && (u != 0) && (log(policies::get_epsilon<T, Policy>() / 2) > v * log((x/2) * (x/2) / v)))
+         {
+            // Evaluate using series representations.
+            // This is particularly important for x << v as in this
+            // area temme_jy may be slow to converge, if it converges at all.
+            // Requires x is not an integer.
+            if(kind&need_j)
+               Jv = bessel_j_small_z_series(v, x, pol);
+            else
+               Jv = std::numeric_limits<T>::quiet_NaN();
+            if((org_kind&need_y && (!reflect || (cp != 0))) 
+               || (org_kind & need_j && (reflect && (sp != 0))))
+            {
+               // Only calculate if we need it, and if the reflection formula will actually use it:
+               Yv = bessel_y_small_z_series(v, x, &Yv_scale, pol);
+            }
+            else
+               Yv = std::numeric_limits<T>::quiet_NaN();
+         }
+         else if((u == 0) && (x < policies::get_epsilon<T, Policy>()))
+         {
+            // Truncated series evaluation for small x and v an integer,
+            // much quicker in this area than temme_jy below.
+            if(kind&need_j)
+               Jv = bessel_j_small_z_series(v, x, pol);
+            else
+               Jv = std::numeric_limits<T>::quiet_NaN();
+            if((org_kind&need_y && (!reflect || (cp != 0))) 
+               || (org_kind & need_j && (reflect && (sp != 0))))
+            {
+               // Only calculate if we need it, and if the reflection formula will actually use it:
+               Yv = bessel_yn_small_z(n, x, &Yv_scale, pol);
+            }
+            else
+               Yv = std::numeric_limits<T>::quiet_NaN();
+         }
+         else if(asymptotic_bessel_large_x_limit(v, x))
+         {
+            if(kind&need_y)
+            {
+               Yv = asymptotic_bessel_y_large_x_2(v, x);
+            }
+            else
+               Yv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
+            if(kind&need_j)
+            {
+               Jv = asymptotic_bessel_j_large_x_2(v, x);
+            }
+            else
+               Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
+         }
+         else if((x > 8) && 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....
+            //
+            // Normally we calculate sin/cos(chi) where:
+            //
+            // chi = x - fmod(T(v / 2 + 0.25f), T(2)) * boost::math::constants::pi<T>();
+            //
+            // But this introduces large errors, so use sin/cos addition formulae to
+            // improve accuracy:
+            //
+            T mod_v = fmod(T(v / 2 + 0.25f), T(2));
+            T sx = sin(x);
+            T cx = cos(x);
+            T sv = sin_pi(mod_v);
+            T cv = cos_pi(mod_v);
+
+            T sc = sx * cv - sv * cx; // == sin(chi);
+            T cc = cx * cv + sx * sv; // == cos(chi);
+            T chi = boost::math::constants::root_two<T>() / (boost::math::constants::root_pi<T>() * sqrt(x)); //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
+            {
+               // domain error:
+               *J = *Y = Yu;
+               return 1;
+            }
+            prev = Yu;
+            current = Yu1;
+            T scale = 1;
+            policies::check_series_iterations<T>(function, n, pol);
+            for (k = 1; k <= n; k++)            // forward recurrence for Y
+            {
+               T fact = 2 * (u + k) / x;
+               if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
+               {
+                  scale /= current;
+                  prev /= current;
+                  current = 1;
+               }
+               next = fact * current - prev;
+               prev = current;
+               current = next;
+            }
+            Yv = prev;
+            Yv1 = current;
+            if(kind&need_j)
+            {
+               CF1_jy(v, x, &fv, &s, pol);                 // continued fraction CF1_jy
+               Jv = scale * W / (Yv * fv - Yv1);           // Wronskian relation
+            }
+            else
+               Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
+            Yv_scale = scale;
+         }
+         else                                    // x in (2, \infty)
+         {
+            // Get Y(u, x):
+
+            T ratio;
+            CF1_jy(v, x, &fv, &s, pol);
+            // tiny initial value to prevent overflow
+            T init = sqrt(tools::min_value<T>());
+            prev = fv * s * init;
+            current = s * init;
+            if(v < max_factorial<T>::value)
+            {
+               policies::check_series_iterations<T>(function, n, pol);
+               for (k = n; k > 0; k--)             // backward recurrence for J
+               {
                   next = 2 * (u + k) * current / x - prev;
                   prev = current;
                   current = next;
-              }
-              ratio = (s * init) / current;     // scaling ratio
-              // can also call CF1_jy() to get fu, not much difference in precision
-              fu = prev / current;
-           }
-           else
-           {
-              //
-              // When v is large we may get overflow in this calculation
-              // leading to NaN's and other nasty surprises:
-              //
-              bool over = false;
-              for (k = n; k > 0; k--)             // backward recurrence for J
-              {
+               }
+               ratio = (s * init) / current;     // scaling ratio
+               // can also call CF1_jy() to get fu, not much difference in precision
+               fu = prev / current;
+            }
+            else
+            {
+               //
+               // When v is large we may get overflow in this calculation
+               // leading to NaN's and other nasty surprises:
+               //
+               policies::check_series_iterations<T>(function, n, pol);
+               bool over = false;
+               for (k = n; k > 0; k--)             // backward recurrence for J
+               {
                   T t = 2 * (u + k) / x;
                   if((t > 1) && (tools::max_value<T>() / t < current))
                   {
@@ -479,87 +489,88 @@
                   next = t * current - prev;
                   prev = current;
                   current = next;
-              }
-              if(!over)
-              {
-                 ratio = (s * init) / current;     // scaling ratio
-                 // can also call CF1_jy() to get fu, not much difference in precision
-                 fu = prev / current;
-              }
-              else
-              {
-                 ratio = 0;
-                 fu = 1;
-              }
-           }
-           CF2_jy(u, x, &p, &q, pol);                  // continued fraction CF2_jy
-           T t = u / x - fu;                   // t = J'/J
-           gamma = (p - t) / q;
-           //
-           // We can't allow gamma to cancel out to zero competely as it messes up
-           // the subsequent logic.  So pretend that one bit didn't cancel out
-           // and set to a suitably small value.  The only test case we've been able to
-           // find for this, is when v = 8.5 and x = 4*PI.
-           //
-           if(gamma == 0)
-           {
-              gamma = u * tools::epsilon<T>() / x;
-           }
-           Ju = sign(current) * sqrt(W / (q + gamma * (p - t)));
-
-           Jv = Ju * ratio;                    // normalization
-
-           Yu = gamma * Ju;
-           Yu1 = Yu * (u/x - p - q/gamma);
-        }
-        if(kind&need_y)
-        {
-           // compute Y:
-           prev = Yu;
-           current = Yu1;
-           for (k = 1; k <= n; k++)            // forward recurrence for Y
-           {
-               T fact = 2 * (u + k) / x;
-               if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
+               }
+               if(!over)
                {
-                  prev /= current;
-                  Yv_scale /= current;
-                  current = 1;
+                  ratio = (s * init) / current;     // scaling ratio
+                  // can also call CF1_jy() to get fu, not much difference in precision
+                  fu = prev / current;
                }
-               next = fact * current - prev;
-               prev = current;
-               current = next;
-           }
-           Yv = prev;
-        }
-        else
-           Yv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
-    }
-
-    if (reflect)
-    {
-        if((sp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(sp * Yv)))
-           *J = org_kind & need_j ? T(-sign(sp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
-        else
-            *J = cp * Jv - (sp == 0 ? T(0) : T((sp * Yv) / Yv_scale));     // reflection formula
-        if((cp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(cp * Yv)))
-           *Y = org_kind & need_y ? T(-sign(cp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
-        else
-           *Y = sp * Jv + (cp == 0 ? T(0) : T((cp * Yv) / Yv_scale));
-    }
-    else
-    {
-        *J = Jv;
-        if(tools::max_value<T>() * fabs(Yv_scale) < fabs(Yv))
-           *Y = org_kind & need_y ? T(sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
-        else
-         *Y = Yv / Yv_scale;
-    }
+               else
+               {
+                  ratio = 0;
+                  fu = 1;
+               }
+            }
+            CF2_jy(u, x, &p, &q, pol);                  // continued fraction CF2_jy
+            T t = u / x - fu;                   // t = J'/J
+            gamma = (p - t) / q;
+            //
+            // We can't allow gamma to cancel out to zero competely as it messes up
+            // the subsequent logic.  So pretend that one bit didn't cancel out
+            // and set to a suitably small value.  The only test case we've been able to
+            // find for this, is when v = 8.5 and x = 4*PI.
+            //
+            if(gamma == 0)
+            {
+               gamma = u * tools::epsilon<T>() / x;
+            }
+            Ju = sign(current) * sqrt(W / (q + gamma * (p - t)));
+
+            Jv = Ju * ratio;                    // normalization
+
+            Yu = gamma * Ju;
+            Yu1 = Yu * (u/x - p - q/gamma);
+
+            if(kind&need_y)
+            {
+               // compute Y:
+               prev = Yu;
+               current = Yu1;
+               policies::check_series_iterations<T>(function, n, pol);
+               for (k = 1; k <= n; k++)            // forward recurrence for Y
+               {
+                  T fact = 2 * (u + k) / x;
+                  if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
+                  {
+                     prev /= current;
+                     Yv_scale /= current;
+                     current = 1;
+                  }
+                  next = fact * current - prev;
+                  prev = current;
+                  current = next;
+               }
+               Yv = prev;
+            }
+            else
+               Yv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
+         }
+
+         if (reflect)
+         {
+            if((sp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(sp * Yv)))
+               *J = org_kind & need_j ? T(-sign(sp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
+            else
+               *J = cp * Jv - (sp == 0 ? T(0) : T((sp * Yv) / Yv_scale));     // reflection formula
+            if((cp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(cp * Yv)))
+               *Y = org_kind & need_y ? T(-sign(cp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
+            else
+               *Y = sp * Jv + (cp == 0 ? T(0) : T((cp * Yv) / Yv_scale));
+         }
+         else
+         {
+            *J = Jv;
+            if(tools::max_value<T>() * fabs(Yv_scale) < fabs(Yv))
+               *Y = org_kind & need_y ? T(sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
+            else
+               *Y = Yv / Yv_scale;
+         }
 
-    return 0;
-}
+         return 0;
+      }
 
-} // namespace detail
+   } // namespace detail
 
 }} // namespaces
 

Modified: trunk/boost/math/special_functions/detail/bessel_jy_asym.hpp
==============================================================================
--- trunk/boost/math/special_functions/detail/bessel_jy_asym.hpp	(original)
+++ trunk/boost/math/special_functions/detail/bessel_jy_asym.hpp	2013-02-11 07:12:50 EST (Mon, 11 Feb 2013)
@@ -21,61 +21,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)
 {
    // Calculate the amplitude of J(v, x) and Y(v, x) for large
@@ -174,89 +119,21 @@
    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>&)
-{
-   // 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>&)
+inline bool asymptotic_bessel_large_x_limit(const T& v, const T& x)
 {
-   // 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>

Modified: trunk/boost/math/special_functions/detail/bessel_yn.hpp
==============================================================================
--- trunk/boost/math/special_functions/detail/bessel_yn.hpp	(original)
+++ trunk/boost/math/special_functions/detail/bessel_yn.hpp	2013-02-11 07:12:50 EST (Mon, 11 Feb 2013)
@@ -75,6 +75,7 @@
        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;

Modified: trunk/libs/math/test/test_bessel_j.hpp
==============================================================================
--- trunk/libs/math/test/test_bessel_j.hpp	(original)
+++ trunk/libs/math/test/test_bessel_j.hpp	2013-02-11 07:12:50 EST (Mon, 11 Feb 2013)
@@ -186,10 +186,11 @@
         {{ 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 @@
     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 @@
 
 #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);
 }