From 76cc581bcef943c7d273eb8a14503354feabd5f6 Mon Sep 17 00:00:00 2001 From: John Maddock Date: Mon, 11 Feb 2013 12:12:50 +0000 Subject: [PATCH] Refactor Bessel function code: * Remove unused dead code. * Improve and centralise asymptotic selection. * Reorder algorithm selection in bessel_jy. * Allow use of integer algorithms for arbitrary order - they're no slower than Steeds method which is also O(n). [SVN r82816] --- .../boost/math/special_functions/bessel.hpp | 28 +- .../special_functions/detail/bessel_jn.hpp | 8 +- .../special_functions/detail/bessel_jy.hpp | 1035 +++++++++-------- .../detail/bessel_jy_asym.hpp | 147 +-- .../special_functions/detail/bessel_yn.hpp | 1 + test/test_bessel_j.hpp | 13 +- 6 files changed, 558 insertions(+), 674 deletions(-) diff --git a/include/boost/math/special_functions/bessel.hpp b/include/boost/math/special_functions/bessel.hpp index ae11c4e84..b0a43f38f 100644 --- a/include/boost/math/special_functions/bessel.hpp +++ b/include/boost/math/special_functions/bessel.hpp @@ -100,22 +100,6 @@ T cyl_bessel_j_imp(T v, T x, const bessel_no_int_tag& t, const Policy& pol) function, "Got x = %1%, but we need x >= 0", x, pol); } - if(x == 0) - return (v == 0) ? 1 : (v > 0) ? 0 : - policies::raise_domain_error( - function, - "Got v = %1%, but require v >= 0 or a negative integer: the result would be complex.", v, pol); - - - if((v >= 0) && ((x < 1) || (v > x * x / 4) || (x < 5))) - { - // - // This series will actually converge rapidly for all small - // x - say up to x < 20 - but the first few terms are large - // and divergent which leads to large errors :-( - // - return bessel_j_small_z_series(v, x, pol); - } T j, y; bessel_jy(v, x, &j, &y, need_j, pol); @@ -127,9 +111,11 @@ inline T cyl_bessel_j_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& p { BOOST_MATH_STD_USING // ADL of std names. int ival = detail::iconv(v, pol); - if((abs(ival) < 200) && (0 == v - ival)) + // If v is an integer, use the integer recursion + // method, both that and Steeds method are O(v): + if((0 == v - ival)) { - return bessel_jn(ival/*iround(v, pol)*/, x, pol); + return bessel_jn(ival, x, pol); } return cyl_bessel_j_imp(v, x, bessel_no_int_tag(), pol); } @@ -296,14 +282,13 @@ template 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::type tag_type; BOOST_MATH_INSTRUMENT_VARIABLE(v); BOOST_MATH_INSTRUMENT_VARIABLE(x); if(floor(v) == v) { - if((fabs(x) > asymptotic_bessel_y_limit(tag_type())) && (4 * fabs(x) * sqrt(sqrt(sqrt(14 * tools::epsilon() * (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(abs(v)), x); if((v < 0) && (itrunc(v, pol) & 1)) @@ -327,12 +312,11 @@ template 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::type tag_type; BOOST_MATH_INSTRUMENT_VARIABLE(v); BOOST_MATH_INSTRUMENT_VARIABLE(x); - if((fabs(x) > asymptotic_bessel_y_limit(tag_type())) && (4 * fabs(x) * sqrt(sqrt(sqrt(14 * tools::epsilon() * (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(abs(v)), x); if((v < 0) && (v & 1)) diff --git a/include/boost/math/special_functions/detail/bessel_jn.hpp b/include/boost/math/special_functions/detail/bessel_jn.hpp index 10fdfb4ac..3f15f9cd8 100644 --- a/include/boost/math/special_functions/detail/bessel_jn.hpp +++ b/include/boost/math/special_functions/detail/bessel_jn.hpp @@ -64,8 +64,7 @@ T bessel_jn(int n, T x, const Policy& pol) return static_cast(0); } - typedef typename bessel_asymptotic_tag::type tag_type; - if(fabs(x) > asymptotic_bessel_j_limit(n, tag_type())) + if(asymptotic_bessel_large_x_limit(T(n), x)) return factor * asymptotic_bessel_j_large_x_2(n, x); BOOST_ASSERT(n > 1); @@ -74,6 +73,7 @@ T bessel_jn(int n, T x, const Policy& pol) { prev = bessel_j0(x); current = bessel_j1(x); + policies::check_series_iterations("boost::math::bessel_j_n<%1%>(%1%,%1%)", n, pol); for (int k = 1; k < n; k++) { T fact = 2 * k / x; @@ -91,7 +91,7 @@ T bessel_jn(int n, T x, const Policy& pol) current = value; } } - else if(x < 1) + else if((x < 1) || (n > x * x / 4) || (x < 5)) { return factor * bessel_j_small_z_series(T(n), x, pol); } @@ -102,6 +102,8 @@ T bessel_jn(int n, T x, const Policy& pol) boost::math::detail::CF1_jy(static_cast(n), x, &fn, &s, pol); prev = fn; current = 1; + // Check recursion won't go on too far: + policies::check_series_iterations("boost::math::bessel_j_n<%1%>(%1%,%1%)", n, pol); for (int k = n; k > 0; k--) { T fact = 2 * k / x; diff --git a/include/boost/math/special_functions/detail/bessel_jy.hpp b/include/boost/math/special_functions/detail/bessel_jy.hpp index aa1f367e7..cbc19900b 100644 --- a/include/boost/math/special_functions/detail/bessel_jy.hpp +++ b/include/boost/math/special_functions/detail/bessel_jy.hpp @@ -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 -bool hankel_PQ(T v, T x, T* p, T* q, const Policy& ) -{ - BOOST_MATH_STD_USING - T tolerance = 2 * policies::get_epsilon(); - *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 -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(1) : sinh(sigma) / sigma; - e = abs(v) < tools::epsilon() ? T(v*pi()*pi() / 2) - : T(2 * spv2 * spv2 / v); - - T g1 = (v == 0) ? T(-euler()) : 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(1/constants::pi()) : 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(); - for (k = 1; k < policies::get_max_series_iterations(); 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("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 -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();; - tiny = sqrt(tools::min_value()); - C = f = tiny; // b0 = 0, replace with tiny - D = 0; - for (k = 1; k < policies::get_max_series_iterations() * 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("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 -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(); - tiny = sqrt(tools::min_value()); - 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(); 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("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 -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( - function, 0, pol); - return 1; - } - - // x is positive until reflection - W = T(2) / (x * pi()); // 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(); - // - // 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() / (boost::math::constants::root_pi() * sqrt(x)); //sqrt(2 / (boost::math::constants::pi() * x)); - Yv = chi * (p * sc + q * cc); - Jv = chi * (p * cc - q * sc); - } - else if((x < 1) && (u != 0) && (log(policies::get_epsilon() / 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::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::quiet_NaN(); - } - else if((u == 0) && (x < policies::get_epsilon())) - { - // 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::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::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() - 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 + bool hankel_PQ(T v, T x, T* p, T* q, const Policy& ) + { + BOOST_MATH_STD_USING + T tolerance = 2 * policies::get_epsilon(); + *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::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::type tag_type; + while((fabs(term) > tolerance * *p) && ok); + return ok; + } - T lim, ratio; - switch(kind) - { - case need_j: - lim = asymptotic_bessel_j_limit(v, tag_type()); - break; - case need_y: - lim = asymptotic_bessel_y_limit(tag_type()); - break; - default: - lim = (std::max)( - asymptotic_bessel_j_limit(v, tag_type()), - asymptotic_bessel_y_limit(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::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::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()); - prev = fv * s * init; - current = s * init; - if(v < max_factorial::value) - { - for (k = n; k > 0; k--) // backward recurrence for J - { + // Calculate Y(v, x) and Y(v+1, x) by Temme's method, see + // Temme, Journal of Computational Physics, vol 21, 343 (1976) + template + 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(1) : sinh(sigma) / sigma; + e = abs(v) < tools::epsilon() ? T(v*pi()*pi() / 2) + : T(2 * spv2 * spv2 / v); + + T g1 = (v == 0) ? T(-euler()) : 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(1/constants::pi()) : 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(); + for (k = 1; k < policies::get_max_series_iterations(); 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("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 + 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();; + tiny = sqrt(tools::min_value()); + C = f = tiny; // b0 = 0, replace with tiny + D = 0; + for (k = 1; k < policies::get_max_series_iterations() * 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("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 + 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(); + tiny = sqrt(tools::min_value()); + 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(); 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("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 + 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((std::numeric_limits::max)())) + policies::raise_evaluation_error(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(function, "Value of Bessel J_v(x) is complex-infinity at %1%", x, pol); // complex infinity + else + *J = std::numeric_limits::quiet_NaN(); // any value will do, not using J. + + if((kind & need_y) == 0) + *Y = std::numeric_limits::quiet_NaN(); // any value will do, not using Y. + else if(v == 0) + *Y = -policies::raise_overflow_error(function, 0, pol); + else + *Y = policies::raise_domain_error(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()); // 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::quiet_NaN(); + } + else if((x < 1) && (u != 0) && (log(policies::get_epsilon() / 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::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::quiet_NaN(); + } + else if((u == 0) && (x < policies::get_epsilon())) + { + // 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::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::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::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::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(); + // + // 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() / (boost::math::constants::root_pi() * sqrt(x)); //sqrt(2 / (boost::math::constants::pi() * 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(function, n, pol); + for (k = 1; k <= n; k++) // forward recurrence for Y + { + T fact = 2 * (u + k) / x; + if((tools::max_value() - 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::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()); + prev = fv * s * init; + current = s * init; + if(v < max_factorial::value) + { + policies::check_series_iterations(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(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 < current)) { @@ -479,87 +489,88 @@ int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol) 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() / 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() - 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::quiet_NaN(); // any value will do, we're not using it. - } + 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() / x; + } + Ju = sign(current) * sqrt(W / (q + gamma * (p - t))); - if (reflect) - { - if((sp != 0) && (tools::max_value() * fabs(Yv_scale) < fabs(sp * Yv))) - *J = org_kind & need_j ? T(-sign(sp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error(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() * fabs(Yv_scale) < fabs(cp * Yv))) - *Y = org_kind & need_y ? T(-sign(cp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error(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() * fabs(Yv_scale) < fabs(Yv)) - *Y = org_kind & need_y ? T(sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error(function, 0, pol)) : T(0); - else - *Y = Yv / Yv_scale; - } + Jv = Ju * ratio; // normalization - return 0; -} + Yu = gamma * Ju; + Yu1 = Yu * (u/x - p - q/gamma); -} // namespace detail + if(kind&need_y) + { + // compute Y: + prev = Yu; + current = Yu1; + policies::check_series_iterations(function, n, pol); + for (k = 1; k <= n; k++) // forward recurrence for Y + { + T fact = 2 * (u + k) / x; + if((tools::max_value() - 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::quiet_NaN(); // any value will do, we're not using it. + } + + if (reflect) + { + if((sp != 0) && (tools::max_value() * fabs(Yv_scale) < fabs(sp * Yv))) + *J = org_kind & need_j ? T(-sign(sp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error(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() * fabs(Yv_scale) < fabs(cp * Yv))) + *Y = org_kind & need_y ? T(-sign(cp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error(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() * fabs(Yv_scale) < fabs(Yv)) + *Y = org_kind & need_y ? T(sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error(function, 0, pol)) : T(0); + else + *Y = Yv / Yv_scale; + } + + return 0; + } + + } // namespace detail }} // namespaces diff --git a/include/boost/math/special_functions/detail/bessel_jy_asym.hpp b/include/boost/math/special_functions/detail/bessel_jy_asym.hpp index ee44adaa8..81f6238e5 100644 --- a/include/boost/math/special_functions/detail/bessel_jy_asym.hpp +++ b/include/boost/math/special_functions/detail/bessel_jy_asym.hpp @@ -20,61 +20,6 @@ namespace boost{ namespace math{ namespace detail{ -template -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 -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 -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() * (2 * v + 1) / 4; - return sqrt(2 / (constants::pi() * x)) - * (asymptotic_bessel_j_large_x_P(v, x) * cos(chi) - - asymptotic_bessel_j_large_x_Q(v, x) * sin(chi)); -} - -template -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() * (2 * v + 1) / 4; - return sqrt(2 / (constants::pi() * x)) - * (asymptotic_bessel_j_large_x_P(v, x) * sin(chi) - - asymptotic_bessel_j_large_x_Q(v, x) * cos(chi)); -} - template inline T asymptotic_bessel_amplitude(T v, T x) { @@ -174,89 +119,21 @@ inline T asymptotic_bessel_j_large_x_2(T v, T x) return sin_phase * ampl; } -// -// Various limits for the J and Y asymptotics -// (the asympotic expansions are safe to use if -// x is less than the limit given). -// We assume that if we don't use these expansions then the -// error will likely be >100eps, so the limits given are chosen -// to lead to < 100eps truncation error. -// template -inline T asymptotic_bessel_y_limit(const mpl::int_<0>&) +inline bool asymptotic_bessel_large_x_limit(const T& v, const T& x) { - // default case: BOOST_MATH_STD_USING - return 2.25 / pow(100 * tools::epsilon() / T(0.001f), T(0.2f)); -} -template -inline T asymptotic_bessel_y_limit(const mpl::int_<53>&) -{ - // double case: - return 304 /*780*/; -} -template -inline T asymptotic_bessel_y_limit(const mpl::int_<64>&) -{ - // 80-bit extended-double case: - return 1552 /*3500*/; -} -template -inline T asymptotic_bessel_y_limit(const mpl::int_<113>&) -{ - // 128-bit long double case: - return 1245243 /*3128000*/; -} - -template -struct bessel_asymptotic_tag -{ - typedef typename policies::precision::type precision_type; - typedef typename mpl::if_< - mpl::or_< - mpl::equal_to >, - mpl::greater > >, - mpl::int_<0>, - typename mpl::if_< - mpl::greater >, - mpl::int_<113>, - typename mpl::if_< - mpl::greater >, - mpl::int_<64>, - mpl::int_<53> - >::type - >::type - >::type type; -}; - -template -inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<0>&) -{ - // default case: - BOOST_MATH_STD_USING - T v2 = (std::max)(T(3), T(v * v)); - return v2 / pow(100 * tools::epsilon() / T(2e-5f), T(0.17f)); -} -template -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 -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 -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()); } template diff --git a/include/boost/math/special_functions/detail/bessel_yn.hpp b/include/boost/math/special_functions/detail/bessel_yn.hpp index d3ee5a58c..0509062bb 100644 --- a/include/boost/math/special_functions/detail/bessel_yn.hpp +++ b/include/boost/math/special_functions/detail/bessel_yn.hpp @@ -75,6 +75,7 @@ T bessel_yn(int n, T x, const Policy& pol) current = bessel_y1(x, pol); int k = 1; BOOST_ASSERT(k < n); + policies::check_series_iterations("boost::math::bessel_y_n<%1%>(%1%,%1%)", n, pol); do { T fact = 2 * k / x; diff --git a/test/test_bessel_j.hpp b/test/test_bessel_j.hpp index 43a578d36..f89991c66 100644 --- a/test/test_bessel_j.hpp +++ b/test/test_bessel_j.hpp @@ -186,10 +186,11 @@ void test_bessel(T, const char* name) {{ SC_(1), T(10667654)/(1024*1024), SC_(1.24591331097191900488116495350277530373473085499043086981229e-7) }}, }}; - static const boost::array::type, 3>, 16> jn_data = {{ + static const boost::array::type, 3>, 17> jn_data = {{ // This first one is a modified test case from https://svn.boost.org/trac/boost/ticket/2733 {{ SC_(-1), SC_(1.25), SC_(-0.510623260319880467069474837274910375352924050139633057168856) }}, {{ SC_(2), SC_(0), SC_(0) }}, + {{ SC_(-2), SC_(0), SC_(0) }}, {{ SC_(2), SC_(1e-02), SC_(1.249989583365885362413250958437642113452e-05) }}, {{ SC_(5), SC_(10), SC_(-0.2340615281867936404436949416457777864635) }}, {{ SC_(5), SC_(-10), SC_(0.2340615281867936404436949416457777864635) }}, @@ -217,8 +218,9 @@ void test_bessel(T, const char* name) do_test_cyl_bessel_j_int(j1_tricky, name, "Bessel J1: Mathworld Data (tricky cases) (Integer Version)"); do_test_cyl_bessel_j_int(jn_data, name, "Bessel JN: Mathworld Data (Integer Version)"); - static const boost::array, 20> jv_data = {{ + static const boost::array, 21> jv_data = {{ //SC_(-2.4), {{ SC_(0), std::numeric_limits::infinity() }}, + {{ T(22.5), T(0), SC_(0) }}, {{ T(2457)/1024, T(1)/1024, SC_(3.80739920118603335646474073457326714709615200130620574875292e-9) }}, {{ SC_(5.5), T(3217)/1024, SC_(0.0281933076257506091621579544064767140470089107926550720453038) }}, {{ SC_(-5.5), T(3217)/1024, SC_(-2.55820064470647911823175836997490971806135336759164272675969) }}, @@ -263,5 +265,12 @@ void test_bessel(T, const char* name) #include "sph_bessel_data.ipp" do_test_sph_bessel_j(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); }