601 lines · plain
1// Copyright (c) 2006 Xiaogang Zhang2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_BESSEL_JY_HPP7#define BOOST_MATH_BESSEL_JY_HPP8 9#ifdef _MSC_VER10#pragma once11#endif12 13#include <boost/math/tools/config.hpp>14#include <boost/math/tools/numeric_limits.hpp>15#include <boost/math/tools/type_traits.hpp>16#include <boost/math/special_functions/gamma.hpp>17#include <boost/math/special_functions/sign.hpp>18#include <boost/math/special_functions/hypot.hpp>19#include <boost/math/special_functions/sin_pi.hpp>20#include <boost/math/special_functions/cos_pi.hpp>21#include <boost/math/special_functions/round.hpp>22#include <boost/math/special_functions/detail/bessel_jy_asym.hpp>23#include <boost/math/special_functions/detail/bessel_jy_series.hpp>24#include <boost/math/constants/constants.hpp>25#include <boost/math/policies/error_handling.hpp>26 27// Bessel functions of the first and second kind of fractional order28 29namespace boost { namespace math {30 31 namespace detail {32 33 //34 // Simultaneous calculation of A&S 9.2.9 and 9.2.1035 // for use in A&S 9.2.5 and 9.2.6.36 // This series is quick to evaluate, but divergent unless37 // x is very large, in fact it's pretty hard to figure out38 // with any degree of precision when this series actually39 // *will* converge!! Consequently, we may just have to40 // try it and see...41 //42 template <class T, class Policy>43 BOOST_MATH_GPU_ENABLED bool hankel_PQ(T v, T x, T* p, T* q, const Policy& )44 {45 BOOST_MATH_STD_USING46 T tolerance = 2 * policies::get_epsilon<T, Policy>();47 *p = 1;48 *q = 0;49 T k = 1;50 T z8 = 8 * x;51 T sq = 1;52 T mu = 4 * v * v;53 T term = 1;54 bool ok = true;55 do56 {57 term *= (mu - sq * sq) / (k * z8);58 *q += term;59 k += 1;60 sq += 2;61 T mult = (sq * sq - mu) / (k * z8);62 ok = fabs(mult) < 0.5f;63 term *= mult;64 *p += term;65 k += 1;66 sq += 2;67 }68 while((fabs(term) > tolerance * *p) && ok);69 return ok;70 }71 72 // Calculate Y(v, x) and Y(v+1, x) by Temme's method, see73 // Temme, Journal of Computational Physics, vol 21, 343 (1976)74 template <typename T, typename Policy>75 BOOST_MATH_GPU_ENABLED int temme_jy(T v, T x, T* Y, T* Y1, const Policy& pol)76 {77 T g, h, p, q, f, coef, sum, sum1, tolerance;78 T a, d, e, sigma;79 unsigned long k;80 81 BOOST_MATH_STD_USING82 using namespace boost::math::tools;83 using namespace boost::math::constants;84 85 BOOST_MATH_ASSERT(fabs(v) <= 0.5f); // precondition for using this routine86 87 T gp = boost::math::tgamma1pm1(v, pol);88 T gm = boost::math::tgamma1pm1(-v, pol);89 T spv = boost::math::sin_pi(v, pol);90 T spv2 = boost::math::sin_pi(v/2, pol);91 T xp = pow(x/2, v);92 93 a = log(x / 2);94 sigma = -a * v;95 d = abs(sigma) < tools::epsilon<T>() ?96 T(1) : sinh(sigma) / sigma;97 e = abs(v) < tools::epsilon<T>() ? T(v*pi<T>()*pi<T>() / 2)98 : T(2 * spv2 * spv2 / v);99 100 T g1 = (v == 0) ? T(-euler<T>()) : T((gp - gm) / ((1 + gp) * (1 + gm) * 2 * v));101 T g2 = (2 + gp + gm) / ((1 + gp) * (1 + gm) * 2);102 T vspv = (fabs(v) < tools::epsilon<T>()) ? T(1/constants::pi<T>()) : T(v / spv);103 f = (g1 * cosh(sigma) - g2 * a * d) * 2 * vspv;104 105 p = vspv / (xp * (1 + gm));106 q = vspv * xp / (1 + gp);107 108 g = f + e * q;109 h = p;110 coef = 1;111 sum = coef * g;112 sum1 = coef * h;113 114 T v2 = v * v;115 T coef_mult = -x * x / 4;116 117 // series summation118 tolerance = policies::get_epsilon<T, Policy>();119 for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)120 {121 f = (k * f + p + q) / (k*k - v2);122 p /= k - v;123 q /= k + v;124 g = f + e * q;125 h = p - k * g;126 coef *= coef_mult / k;127 sum += coef * g;128 sum1 += coef * h;129 if (abs(coef * g) < abs(sum) * tolerance)130 {131 break;132 }133 }134 policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in temme_jy", k, pol);135 *Y = -sum;136 *Y1 = -2 * sum1 / x;137 138 return 0;139 }140 141 // Evaluate continued fraction fv = J_(v+1) / J_v, see142 // Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73143 template <typename T, typename Policy>144 BOOST_MATH_GPU_ENABLED int CF1_jy(T v, T x, T* fv, int* sign, const Policy& pol)145 {146 T C, D, f, a, b, delta, tiny, tolerance;147 unsigned long k;148 int s = 1;149 150 BOOST_MATH_STD_USING151 152 // |x| <= |v|, CF1_jy converges rapidly153 // |x| > |v|, CF1_jy needs O(|x|) iterations to converge154 155 // modified Lentz's method, see156 // Lentz, Applied Optics, vol 15, 668 (1976)157 tolerance = 2 * policies::get_epsilon<T, Policy>();158 tiny = sqrt(tools::min_value<T>());159 C = f = tiny; // b0 = 0, replace with tiny160 D = 0;161 for (k = 1; k < policies::get_max_series_iterations<Policy>() * 100; k++)162 {163 a = -1;164 b = 2 * (v + k) / x;165 C = b + a / C;166 D = b + a * D;167 if (C == 0) { C = tiny; }168 if (D == 0) { D = tiny; }169 D = 1 / D;170 delta = C * D;171 f *= delta;172 if (D < 0) { s = -s; }173 if (abs(delta - 1) < tolerance)174 { break; }175 }176 policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF1_jy", k / 100, pol);177 *fv = -f;178 *sign = s; // sign of denominator179 180 return 0;181 }182 //183 // This algorithm was originally written by Xiaogang Zhang184 // using std::complex to perform the complex arithmetic.185 // However, that turns out to 10x or more slower than using186 // all real-valued arithmetic, so it's been rewritten using187 // real values only.188 //189 template <typename T, typename Policy>190 BOOST_MATH_GPU_ENABLED int CF2_jy(T v, T x, T* p, T* q, const Policy& pol)191 {192 BOOST_MATH_STD_USING193 194 T Cr, Ci, Dr, Di, fr, fi, a, br, bi, delta_r, delta_i, temp;195 T tiny;196 unsigned long k;197 198 // |x| >= |v|, CF2_jy converges rapidly199 // |x| -> 0, CF2_jy fails to converge200 BOOST_MATH_ASSERT(fabs(x) > 1);201 202 // modified Lentz's method, complex numbers involved, see203 // Lentz, Applied Optics, vol 15, 668 (1976)204 T tolerance = 2 * policies::get_epsilon<T, Policy>();205 tiny = sqrt(tools::min_value<T>());206 Cr = fr = -0.5f / x;207 Ci = fi = 1;208 //Dr = Di = 0;209 T v2 = v * v;210 a = (0.25f - v2) / x; // Note complex this one time only!211 br = 2 * x;212 bi = 2;213 temp = Cr * Cr + 1;214 Ci = bi + a * Cr / temp;215 Cr = br + a / temp;216 Dr = br;217 Di = bi;218 if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }219 if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }220 temp = Dr * Dr + Di * Di;221 Dr = Dr / temp;222 Di = -Di / temp;223 delta_r = Cr * Dr - Ci * Di;224 delta_i = Ci * Dr + Cr * Di;225 temp = fr;226 fr = temp * delta_r - fi * delta_i;227 fi = temp * delta_i + fi * delta_r;228 for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++)229 {230 a = static_cast<T>(k) - 0.5f;231 a *= a;232 a -= v2;233 bi += 2;234 temp = Cr * Cr + Ci * Ci;235 Cr = br + a * Cr / temp;236 Ci = bi - a * Ci / temp;237 Dr = br + a * Dr;238 Di = bi + a * Di;239 if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }240 if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }241 temp = Dr * Dr + Di * Di;242 Dr = Dr / temp;243 Di = -Di / temp;244 delta_r = Cr * Dr - Ci * Di;245 delta_i = Ci * Dr + Cr * Di;246 temp = fr;247 fr = temp * delta_r - fi * delta_i;248 fi = temp * delta_i + fi * delta_r;249 if (fabs(delta_r - 1) + fabs(delta_i) < tolerance)250 break;251 }252 policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol);253 *p = fr;254 *q = fi;255 256 return 0;257 }258 259 BOOST_MATH_STATIC const int need_j = 1;260 BOOST_MATH_STATIC const int need_y = 2;261 262 // Compute J(v, x) and Y(v, x) simultaneously by Steed's method, see263 // Barnett et al, Computer Physics Communications, vol 8, 377 (1974)264 template <typename T, typename Policy>265 BOOST_MATH_GPU_ENABLED int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol)266 {267 BOOST_MATH_ASSERT(x >= 0);268 269 T u, Jv, Ju, Yv, Yv1, Yu, Yu1(0), fv, fu;270 T W, p, q, gamma, current, prev, next;271 bool reflect = false;272 unsigned n, k;273 int s;274 int org_kind = kind;275 T cp = 0;276 T sp = 0;277 278 constexpr auto function = "boost::math::bessel_jy<%1%>(%1%,%1%)";279 280 BOOST_MATH_STD_USING281 using namespace boost::math::tools;282 using namespace boost::math::constants;283 284 if (v < 0)285 {286 reflect = true;287 v = -v; // v is non-negative from here288 }289 if (v > static_cast<T>((boost::math::numeric_limits<int>::max)()))290 {291 *J = *Y = policies::raise_evaluation_error<T>(function, "Order of Bessel function is too large to evaluate: got %1%", v, pol);292 return 1; // LCOV_EXCL_LINE previous line will throw.293 }294 n = static_cast<unsigned>(iround(v, pol));295 u = v - n; // -1/2 <= u < 1/2296 297 if(reflect)298 {299 T z = (u + n % 2);300 cp = boost::math::cos_pi(z, pol);301 sp = boost::math::sin_pi(z, pol);302 if(u != 0)303 kind = need_j|need_y; // need both for reflection formula304 }305 306 if(x == 0)307 {308 if (v == 0)309 *J = 1; // LCOV_EXCL_LINE multiprecision case only310 else if ((u == 0) || !reflect)311 *J = 0;312 else if(kind & need_j)313 *J = policies::raise_domain_error<T>(function, "Value of Bessel J_v(x) is complex-infinity at %1%", x, pol); // complex infinity314 else315 *J = boost::math::numeric_limits<T>::quiet_NaN(); // LCOV_EXCL_LINE, we should never get here, any value will do, not using J.316 317 if((kind & need_y) == 0)318 *Y = boost::math::numeric_limits<T>::quiet_NaN(); // any value will do, not using Y.319 else320 {321 // We shoud never get here:322 BOOST_MATH_ASSERT(x != 0); // LCOV_EXCL_LINE323 }324 return 1;325 }326 327 // x is positive until reflection328 W = T(2) / (x * pi<T>()); // Wronskian329 T Yv_scale = 1;330 if(((kind & need_y) == 0) && ((x < 1) || (v > x * x / 4) || (x < 5)))331 {332 //333 // This series will actually converge rapidly for all small334 // x - say up to x < 20 - but the first few terms are large335 // and divergent which leads to large errors :-(336 //337 Jv = bessel_j_small_z_series(v, x, pol);338 Yv = boost::math::numeric_limits<T>::quiet_NaN();339 }340 else if((x < 1) && (u != 0) && (log(policies::get_epsilon<T, Policy>() / 2) > v * log((x/2) * (x/2) / v)))341 {342 // Evaluate using series representations.343 // This is particularly important for x << v as in this344 // area temme_jy may be slow to converge, if it converges at all.345 // Requires x is not an integer.346 if(kind&need_j)347 Jv = bessel_j_small_z_series(v, x, pol);348 else349 Jv = boost::math::numeric_limits<T>::quiet_NaN();350 if((org_kind&need_y && (!reflect || (cp != 0)))351 || (org_kind & need_j && (reflect && (sp != 0))))352 {353 // Only calculate if we need it, and if the reflection formula will actually use it:354 Yv = bessel_y_small_z_series(v, x, &Yv_scale, pol);355 }356 else357 Yv = boost::math::numeric_limits<T>::quiet_NaN();358 }359 else if((u == 0) && (x < policies::get_epsilon<T, Policy>()))360 {361 // Truncated series evaluation for small x and v an integer,362 // much quicker in this area than temme_jy below.363 // This code is only used in the multiprecision case, otherwise364 // we go via bessel_jn.365 // LCOV_EXCL_START366 if(kind&need_j)367 Jv = bessel_j_small_z_series(v, x, pol);368 else369 Jv = boost::math::numeric_limits<T>::quiet_NaN();370 if((org_kind&need_y && (!reflect || (cp != 0)))371 || (org_kind & need_j && (reflect && (sp != 0))))372 {373 // Only calculate if we need it, and if the reflection formula will actually use it:374 Yv = bessel_yn_small_z(static_cast<int>(n), x, &Yv_scale, pol);375 }376 else377 Yv = boost::math::numeric_limits<T>::quiet_NaN();378 // LCOV_EXCL_STOP379 }380 else if(asymptotic_bessel_large_x_limit(v, x))381 {382 if(kind&need_y)383 {384 Yv = asymptotic_bessel_y_large_x_2(v, x, pol);385 }386 else387 Yv = boost::math::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.388 if(kind&need_j)389 {390 Jv = asymptotic_bessel_j_large_x_2(v, x, pol);391 }392 else393 Jv = boost::math::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.394 }395 else if((x > 8) && hankel_PQ(v, x, &p, &q, pol))396 {397 //398 // Hankel approximation: note that this method works best when x399 // is large, but in that case we end up calculating sines and cosines400 // of large values, with horrendous resulting accuracy. It is fast though401 // when it works....402 //403 // Normally we calculate sin/cos(chi) where:404 //405 // chi = x - fmod(T(v / 2 + 0.25f), T(2)) * boost::math::constants::pi<T>();406 //407 // But this introduces large errors, so use sin/cos addition formulae to408 // improve accuracy:409 //410 T mod_v = fmod(T(v / 2 + 0.25f), T(2));411 T sx = sin(x);412 T cx = cos(x);413 T sv = boost::math::sin_pi(mod_v, pol);414 T cv = boost::math::cos_pi(mod_v, pol);415 416 T sc = sx * cv - sv * cx; // == sin(chi);417 T cc = cx * cv + sx * sv; // == cos(chi);418 T chi = boost::math::constants::root_two<T>() / (boost::math::constants::root_pi<T>() * sqrt(x)); //sqrt(2 / (boost::math::constants::pi<T>() * x));419 Yv = chi * (p * sc + q * cc);420 Jv = chi * (p * cc - q * sc);421 }422 else if (x <= 2) // x in (0, 2]423 {424 if(temme_jy(u, x, &Yu, &Yu1, pol)) // Temme series425 {426 // domain error, this should really have already been handled.427 *J = *Y = Yu; // LCOV_EXCL_LINE428 return 1; // LCOV_EXCL_LINE429 }430 prev = Yu;431 current = Yu1;432 T scale = 1;433 policies::check_series_iterations<T>(function, n, pol);434 for (k = 1; k <= n; k++) // forward recurrence for Y435 {436 T fact = 2 * (u + k) / x;437 if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))438 {439 scale /= current;440 prev /= current;441 current = 1;442 }443 next = fact * current - prev;444 prev = current;445 current = next;446 }447 Yv = prev;448 Yv1 = current;449 if(kind&need_j)450 {451 CF1_jy(v, x, &fv, &s, pol); // continued fraction CF1_jy452 Jv = scale * W / (Yv * fv - Yv1); // Wronskian relation453 }454 else455 Jv = boost::math::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.456 Yv_scale = scale;457 }458 else // x in (2, \infty)459 {460 // Get Y(u, x):461 462 T ratio;463 CF1_jy(v, x, &fv, &s, pol);464 // tiny initial value to prevent overflow465 T init = sqrt(tools::min_value<T>());466 BOOST_MATH_INSTRUMENT_VARIABLE(init);467 prev = fv * s * init;468 current = s * init;469 if(v < max_factorial<T>::value)470 {471 policies::check_series_iterations<T>(function, n, pol);472 for (k = n; k > 0; k--) // backward recurrence for J473 {474 next = 2 * (u + k) * current / x - prev;475 //476 // We can't allow next to completely cancel out or the subsequent logic breaks.477 // Pretend that one bit did not cancel:478 if (next == 0)479 {480 next = prev * tools::epsilon<T>() / 2; // LCOV_EXCL_LINE requires specific hardware and rounding to trigger, does get tested on msvc481 }482 prev = current;483 current = next;484 }485 ratio = (s * init) / current; // scaling ratio486 // can also call CF1_jy() to get fu, not much difference in precision487 fu = prev / current;488 }489 else490 {491 //492 // When v is large we may get overflow in this calculation493 // leading to NaN's and other nasty surprises:494 //495 policies::check_series_iterations<T>(function, n, pol);496 bool over = false;497 for (k = n; k > 0; k--) // backward recurrence for J498 {499 T t = 2 * (u + k) / x;500 if((t > 1) && (tools::max_value<T>() / t < current))501 {502 over = true;503 break;504 }505 next = t * current - prev;506 prev = current;507 current = next;508 }509 if(!over)510 {511 ratio = (s * init) / current; // scaling ratio512 // can also call CF1_jy() to get fu, not much difference in precision513 fu = prev / current;514 }515 else516 {517 ratio = 0;518 fu = 1;519 }520 }521 CF2_jy(u, x, &p, &q, pol); // continued fraction CF2_jy522 T t = u / x - fu; // t = J'/J523 gamma = (p - t) / q;524 //525 // We can't allow gamma to cancel out to zero completely as it messes up526 // the subsequent logic. So pretend that one bit didn't cancel out527 // and set to a suitably small value. The only test case we've been able to528 // find for this, is when v = 8.5 and x = 4*PI.529 //530 if(gamma == 0)531 {532 gamma = u * tools::epsilon<T>() / x; // LCOV_EXCL_LINE requires specific hardware and rounding to trigger, does get tested on msvc533 }534 BOOST_MATH_INSTRUMENT_VARIABLE(current);535 BOOST_MATH_INSTRUMENT_VARIABLE(W);536 BOOST_MATH_INSTRUMENT_VARIABLE(q);537 BOOST_MATH_INSTRUMENT_VARIABLE(gamma);538 BOOST_MATH_INSTRUMENT_VARIABLE(p);539 BOOST_MATH_INSTRUMENT_VARIABLE(t);540 Ju = sign(current) * sqrt(W / (q + gamma * (p - t)));541 BOOST_MATH_INSTRUMENT_VARIABLE(Ju);542 543 Jv = Ju * ratio; // normalization544 545 Yu = gamma * Ju;546 Yu1 = Yu * (u/x - p - q/gamma);547 548 if(kind&need_y)549 {550 // compute Y:551 prev = Yu;552 current = Yu1;553 policies::check_series_iterations<T>(function, n, pol);554 for (k = 1; k <= n; k++) // forward recurrence for Y555 {556 T fact = 2 * (u + k) / x;557 if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))558 {559 prev /= current;560 Yv_scale /= current;561 current = 1;562 }563 next = fact * current - prev;564 prev = current;565 current = next;566 }567 Yv = prev;568 }569 else570 Yv = boost::math::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.571 }572 573 if (reflect)574 {575 if((sp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(sp * Yv)))576 *J = org_kind & need_j ? T(-sign(sp) * sign(Yv) * (Yv_scale != 0 ? sign(Yv_scale) : 1) * policies::raise_overflow_error<T>(function, nullptr, pol)) : T(0);577 else578 *J = cp * Jv - (sp == 0 ? T(0) : T((sp * Yv) / Yv_scale)); // reflection formula579 if((cp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(cp * Yv)))580 *Y = org_kind & need_y ? T(-sign(cp) * sign(Yv) * (Yv_scale != 0 ? sign(Yv_scale) : 1) * policies::raise_overflow_error<T>(function, nullptr, pol)) : T(0);581 else582 *Y = (sp != 0 ? sp * Jv : T(0)) + (cp == 0 ? T(0) : T((cp * Yv) / Yv_scale));583 }584 else585 {586 *J = Jv;587 if(tools::max_value<T>() * fabs(Yv_scale) < fabs(Yv))588 *Y = org_kind & need_y ? T(sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, nullptr, pol)) : T(0);589 else590 *Y = Yv / Yv_scale;591 }592 593 return 0;594 }595 596 } // namespace detail597 598}} // namespaces599 600#endif // BOOST_MATH_BESSEL_JY_HPP601