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);
}