1 // Copyright (c) 2006 Xiaogang Zhang 2 // Use, modification and distribution are subject to the 3 // Boost Software License, Version 1.0. (See accompanying file 4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 5 6 #ifndef BOOST_MATH_BESSEL_JY_HPP 7 #define BOOST_MATH_BESSEL_JY_HPP 8 9 #ifdef _MSC_VER 10 #pragma once 11 #endif 12 13 #include <boost/math/tools/config.hpp> 14 #include <boost/math/special_functions/gamma.hpp> 15 #include <boost/math/special_functions/sign.hpp> 16 #include <boost/math/special_functions/hypot.hpp> 17 #include <boost/math/special_functions/sin_pi.hpp> 18 #include <boost/math/special_functions/cos_pi.hpp> 19 #include <boost/math/special_functions/detail/bessel_jy_asym.hpp> 20 #include <boost/math/special_functions/detail/bessel_jy_series.hpp> 21 #include <boost/math/constants/constants.hpp> 22 #include <boost/math/policies/error_handling.hpp> 23 #include <boost/mpl/if.hpp> 24 #include <boost/type_traits/is_floating_point.hpp> 25 #include <complex> 26 27 // Bessel functions of the first and second kind of fractional order 28 29 namespace boost { namespace math { 30 31 namespace detail { 32 33 // 34 // Simultaneous calculation of A&S 9.2.9 and 9.2.10 35 // for use in A&S 9.2.5 and 9.2.6. 36 // This series is quick to evaluate, but divergent unless 37 // x is very large, in fact it's pretty hard to figure out 38 // with any degree of precision when this series actually 39 // *will* converge!! Consequently, we may just have to 40 // try it and see... 41 // 42 template <class T, class Policy> hankel_PQ(T v,T x,T * p,T * q,const Policy &)43 bool hankel_PQ(T v, T x, T* p, T* q, const Policy& ) 44 { 45 BOOST_MATH_STD_USING 46 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 do 56 { 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, see 73 // Temme, Journal of Computational Physics, vol 21, 343 (1976) 74 template <typename T, typename Policy> temme_jy(T v,T x,T * Y,T * Y1,const Policy & pol)75 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_USING 82 using namespace boost::math::tools; 83 using namespace boost::math::constants; 84 85 BOOST_ASSERT(fabs(v) <= 0.5f); // precondition for using this routine 86 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 summation 118 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, see 142 // Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73 143 template <typename T, typename Policy> CF1_jy(T v,T x,T * fv,int * sign,const Policy & pol)144 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_USING 151 152 // |x| <= |v|, CF1_jy converges rapidly 153 // |x| > |v|, CF1_jy needs O(|x|) iterations to converge 154 155 // modified Lentz's method, see 156 // 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 tiny 160 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 denominator 179 180 return 0; 181 } 182 // 183 // This algorithm was originally written by Xiaogang Zhang 184 // using std::complex to perform the complex arithmetic. 185 // However, that turns out to 10x or more slower than using 186 // all real-valued arithmetic, so it's been rewritten using 187 // real values only. 188 // 189 template <typename T, typename Policy> CF2_jy(T v,T x,T * p,T * q,const Policy & pol)190 int CF2_jy(T v, T x, T* p, T* q, const Policy& pol) 191 { 192 BOOST_MATH_STD_USING 193 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 rapidly 199 // |x| -> 0, CF2_jy fails to converge 200 BOOST_ASSERT(fabs(x) > 1); 201 202 // modified Lentz's method, complex numbers involved, see 203 // 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 = 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 static const int need_j = 1; 260 static const int need_y = 2; 261 262 // Compute J(v, x) and Y(v, x) simultaneously by Steed's method, see 263 // Barnett et al, Computer Physics Communications, vol 8, 377 (1974) 264 template <typename T, typename Policy> bessel_jy(T v,T x,T * J,T * Y,int kind,const Policy & pol)265 int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol) 266 { 267 BOOST_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 static const char* function = "boost::math::bessel_jy<%1%>(%1%,%1%)"; 279 280 BOOST_MATH_STD_USING 281 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 here 288 } 289 if (v > static_cast<T>((std::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; 293 } 294 n = iround(v, pol); 295 u = v - n; // -1/2 <= u < 1/2 296 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 formula 304 } 305 306 if(x == 0) 307 { 308 if(v == 0) 309 *J = 1; 310 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 infinity 314 else 315 *J = std::numeric_limits<T>::quiet_NaN(); // any value will do, not using J. 316 317 if((kind & need_y) == 0) 318 *Y = std::numeric_limits<T>::quiet_NaN(); // any value will do, not using Y. 319 else if(v == 0) 320 *Y = -policies::raise_overflow_error<T>(function, 0, pol); 321 else 322 *Y = policies::raise_domain_error<T>(function, "Value of Bessel Y_v(x) is complex-infinity at %1%", x, pol); // complex infinity 323 return 1; 324 } 325 326 // x is positive until reflection 327 W = T(2) / (x * pi<T>()); // Wronskian 328 T Yv_scale = 1; 329 if(((kind & need_y) == 0) && ((x < 1) || (v > x * x / 4) || (x < 5))) 330 { 331 // 332 // This series will actually converge rapidly for all small 333 // x - say up to x < 20 - but the first few terms are large 334 // and divergent which leads to large errors :-( 335 // 336 Jv = bessel_j_small_z_series(v, x, pol); 337 Yv = std::numeric_limits<T>::quiet_NaN(); 338 } 339 else if((x < 1) && (u != 0) && (log(policies::get_epsilon<T, Policy>() / 2) > v * log((x/2) * (x/2) / v))) 340 { 341 // Evaluate using series representations. 342 // This is particularly important for x << v as in this 343 // area temme_jy may be slow to converge, if it converges at all. 344 // Requires x is not an integer. 345 if(kind&need_j) 346 Jv = bessel_j_small_z_series(v, x, pol); 347 else 348 Jv = std::numeric_limits<T>::quiet_NaN(); 349 if((org_kind&need_y && (!reflect || (cp != 0))) 350 || (org_kind & need_j && (reflect && (sp != 0)))) 351 { 352 // Only calculate if we need it, and if the reflection formula will actually use it: 353 Yv = bessel_y_small_z_series(v, x, &Yv_scale, pol); 354 } 355 else 356 Yv = std::numeric_limits<T>::quiet_NaN(); 357 } 358 else if((u == 0) && (x < policies::get_epsilon<T, Policy>())) 359 { 360 // Truncated series evaluation for small x and v an integer, 361 // much quicker in this area than temme_jy below. 362 if(kind&need_j) 363 Jv = bessel_j_small_z_series(v, x, pol); 364 else 365 Jv = std::numeric_limits<T>::quiet_NaN(); 366 if((org_kind&need_y && (!reflect || (cp != 0))) 367 || (org_kind & need_j && (reflect && (sp != 0)))) 368 { 369 // Only calculate if we need it, and if the reflection formula will actually use it: 370 Yv = bessel_yn_small_z(n, x, &Yv_scale, pol); 371 } 372 else 373 Yv = std::numeric_limits<T>::quiet_NaN(); 374 } 375 else if(asymptotic_bessel_large_x_limit(v, x)) 376 { 377 if(kind&need_y) 378 { 379 Yv = asymptotic_bessel_y_large_x_2(v, x); 380 } 381 else 382 Yv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it. 383 if(kind&need_j) 384 { 385 Jv = asymptotic_bessel_j_large_x_2(v, x); 386 } 387 else 388 Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it. 389 } 390 else if((x > 8) && hankel_PQ(v, x, &p, &q, pol)) 391 { 392 // 393 // Hankel approximation: note that this method works best when x 394 // is large, but in that case we end up calculating sines and cosines 395 // of large values, with horrendous resulting accuracy. It is fast though 396 // when it works.... 397 // 398 // Normally we calculate sin/cos(chi) where: 399 // 400 // chi = x - fmod(T(v / 2 + 0.25f), T(2)) * boost::math::constants::pi<T>(); 401 // 402 // But this introduces large errors, so use sin/cos addition formulae to 403 // improve accuracy: 404 // 405 T mod_v = fmod(T(v / 2 + 0.25f), T(2)); 406 T sx = sin(x); 407 T cx = cos(x); 408 T sv = sin_pi(mod_v); 409 T cv = cos_pi(mod_v); 410 411 T sc = sx * cv - sv * cx; // == sin(chi); 412 T cc = cx * cv + sx * sv; // == cos(chi); 413 T chi = boost::math::constants::root_two<T>() / (boost::math::constants::root_pi<T>() * sqrt(x)); //sqrt(2 / (boost::math::constants::pi<T>() * x)); 414 Yv = chi * (p * sc + q * cc); 415 Jv = chi * (p * cc - q * sc); 416 } 417 else if (x <= 2) // x in (0, 2] 418 { 419 if(temme_jy(u, x, &Yu, &Yu1, pol)) // Temme series 420 { 421 // domain error: 422 *J = *Y = Yu; 423 return 1; 424 } 425 prev = Yu; 426 current = Yu1; 427 T scale = 1; 428 policies::check_series_iterations<T>(function, n, pol); 429 for (k = 1; k <= n; k++) // forward recurrence for Y 430 { 431 T fact = 2 * (u + k) / x; 432 if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current)) 433 { 434 scale /= current; 435 prev /= current; 436 current = 1; 437 } 438 next = fact * current - prev; 439 prev = current; 440 current = next; 441 } 442 Yv = prev; 443 Yv1 = current; 444 if(kind&need_j) 445 { 446 CF1_jy(v, x, &fv, &s, pol); // continued fraction CF1_jy 447 Jv = scale * W / (Yv * fv - Yv1); // Wronskian relation 448 } 449 else 450 Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it. 451 Yv_scale = scale; 452 } 453 else // x in (2, \infty) 454 { 455 // Get Y(u, x): 456 457 T ratio; 458 CF1_jy(v, x, &fv, &s, pol); 459 // tiny initial value to prevent overflow 460 T init = sqrt(tools::min_value<T>()); 461 BOOST_MATH_INSTRUMENT_VARIABLE(init); 462 prev = fv * s * init; 463 current = s * init; 464 if(v < max_factorial<T>::value) 465 { 466 policies::check_series_iterations<T>(function, n, pol); 467 for (k = n; k > 0; k--) // backward recurrence for J 468 { 469 next = 2 * (u + k) * current / x - prev; 470 prev = current; 471 current = next; 472 } 473 ratio = (s * init) / current; // scaling ratio 474 // can also call CF1_jy() to get fu, not much difference in precision 475 fu = prev / current; 476 } 477 else 478 { 479 // 480 // When v is large we may get overflow in this calculation 481 // leading to NaN's and other nasty surprises: 482 // 483 policies::check_series_iterations<T>(function, n, pol); 484 bool over = false; 485 for (k = n; k > 0; k--) // backward recurrence for J 486 { 487 T t = 2 * (u + k) / x; 488 if((t > 1) && (tools::max_value<T>() / t < current)) 489 { 490 over = true; 491 break; 492 } 493 next = t * current - prev; 494 prev = current; 495 current = next; 496 } 497 if(!over) 498 { 499 ratio = (s * init) / current; // scaling ratio 500 // can also call CF1_jy() to get fu, not much difference in precision 501 fu = prev / current; 502 } 503 else 504 { 505 ratio = 0; 506 fu = 1; 507 } 508 } 509 CF2_jy(u, x, &p, &q, pol); // continued fraction CF2_jy 510 T t = u / x - fu; // t = J'/J 511 gamma = (p - t) / q; 512 // 513 // We can't allow gamma to cancel out to zero competely as it messes up 514 // the subsequent logic. So pretend that one bit didn't cancel out 515 // and set to a suitably small value. The only test case we've been able to 516 // find for this, is when v = 8.5 and x = 4*PI. 517 // 518 if(gamma == 0) 519 { 520 gamma = u * tools::epsilon<T>() / x; 521 } 522 BOOST_MATH_INSTRUMENT_VARIABLE(current); 523 BOOST_MATH_INSTRUMENT_VARIABLE(W); 524 BOOST_MATH_INSTRUMENT_VARIABLE(q); 525 BOOST_MATH_INSTRUMENT_VARIABLE(gamma); 526 BOOST_MATH_INSTRUMENT_VARIABLE(p); 527 BOOST_MATH_INSTRUMENT_VARIABLE(t); 528 Ju = sign(current) * sqrt(W / (q + gamma * (p - t))); 529 BOOST_MATH_INSTRUMENT_VARIABLE(Ju); 530 531 Jv = Ju * ratio; // normalization 532 533 Yu = gamma * Ju; 534 Yu1 = Yu * (u/x - p - q/gamma); 535 536 if(kind&need_y) 537 { 538 // compute Y: 539 prev = Yu; 540 current = Yu1; 541 policies::check_series_iterations<T>(function, n, pol); 542 for (k = 1; k <= n; k++) // forward recurrence for Y 543 { 544 T fact = 2 * (u + k) / x; 545 if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current)) 546 { 547 prev /= current; 548 Yv_scale /= current; 549 current = 1; 550 } 551 next = fact * current - prev; 552 prev = current; 553 current = next; 554 } 555 Yv = prev; 556 } 557 else 558 Yv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it. 559 } 560 561 if (reflect) 562 { 563 if((sp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(sp * Yv))) 564 *J = org_kind & need_j ? T(-sign(sp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0); 565 else 566 *J = cp * Jv - (sp == 0 ? T(0) : T((sp * Yv) / Yv_scale)); // reflection formula 567 if((cp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(cp * Yv))) 568 *Y = org_kind & need_y ? T(-sign(cp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0); 569 else 570 *Y = (sp != 0 ? sp * Jv : T(0)) + (cp == 0 ? T(0) : T((cp * Yv) / Yv_scale)); 571 } 572 else 573 { 574 *J = Jv; 575 if(tools::max_value<T>() * fabs(Yv_scale) < fabs(Yv)) 576 *Y = org_kind & need_y ? T(sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0); 577 else 578 *Y = Yv / Yv_scale; 579 } 580 581 return 0; 582 } 583 584 } // namespace detail 585 586 }} // namespaces 587 588 #endif // BOOST_MATH_BESSEL_JY_HPP 589 590