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 <complex> 24 25 // Bessel functions of the first and second kind of fractional order 26 27 namespace boost { namespace math { 28 29 namespace detail { 30 31 // 32 // Simultaneous calculation of A&S 9.2.9 and 9.2.10 33 // for use in A&S 9.2.5 and 9.2.6. 34 // This series is quick to evaluate, but divergent unless 35 // x is very large, in fact it's pretty hard to figure out 36 // with any degree of precision when this series actually 37 // *will* converge!! Consequently, we may just have to 38 // try it and see... 39 // 40 template <class T, class Policy> hankel_PQ(T v,T x,T * p,T * q,const Policy &)41 bool hankel_PQ(T v, T x, T* p, T* q, const Policy& ) 42 { 43 BOOST_MATH_STD_USING 44 T tolerance = 2 * policies::get_epsilon<T, Policy>(); 45 *p = 1; 46 *q = 0; 47 T k = 1; 48 T z8 = 8 * x; 49 T sq = 1; 50 T mu = 4 * v * v; 51 T term = 1; 52 bool ok = true; 53 do 54 { 55 term *= (mu - sq * sq) / (k * z8); 56 *q += term; 57 k += 1; 58 sq += 2; 59 T mult = (sq * sq - mu) / (k * z8); 60 ok = fabs(mult) < 0.5f; 61 term *= mult; 62 *p += term; 63 k += 1; 64 sq += 2; 65 } 66 while((fabs(term) > tolerance * *p) && ok); 67 return ok; 68 } 69 70 // Calculate Y(v, x) and Y(v+1, x) by Temme's method, see 71 // Temme, Journal of Computational Physics, vol 21, 343 (1976) 72 template <typename T, typename Policy> temme_jy(T v,T x,T * Y,T * Y1,const Policy & pol)73 int temme_jy(T v, T x, T* Y, T* Y1, const Policy& pol) 74 { 75 T g, h, p, q, f, coef, sum, sum1, tolerance; 76 T a, d, e, sigma; 77 unsigned long k; 78 79 BOOST_MATH_STD_USING 80 using namespace boost::math::tools; 81 using namespace boost::math::constants; 82 83 BOOST_MATH_ASSERT(fabs(v) <= 0.5f); // precondition for using this routine 84 85 T gp = boost::math::tgamma1pm1(v, pol); 86 T gm = boost::math::tgamma1pm1(-v, pol); 87 T spv = boost::math::sin_pi(v, pol); 88 T spv2 = boost::math::sin_pi(v/2, pol); 89 T xp = pow(x/2, v); 90 91 a = log(x / 2); 92 sigma = -a * v; 93 d = abs(sigma) < tools::epsilon<T>() ? 94 T(1) : sinh(sigma) / sigma; 95 e = abs(v) < tools::epsilon<T>() ? T(v*pi<T>()*pi<T>() / 2) 96 : T(2 * spv2 * spv2 / v); 97 98 T g1 = (v == 0) ? T(-euler<T>()) : T((gp - gm) / ((1 + gp) * (1 + gm) * 2 * v)); 99 T g2 = (2 + gp + gm) / ((1 + gp) * (1 + gm) * 2); 100 T vspv = (fabs(v) < tools::epsilon<T>()) ? T(1/constants::pi<T>()) : T(v / spv); 101 f = (g1 * cosh(sigma) - g2 * a * d) * 2 * vspv; 102 103 p = vspv / (xp * (1 + gm)); 104 q = vspv * xp / (1 + gp); 105 106 g = f + e * q; 107 h = p; 108 coef = 1; 109 sum = coef * g; 110 sum1 = coef * h; 111 112 T v2 = v * v; 113 T coef_mult = -x * x / 4; 114 115 // series summation 116 tolerance = policies::get_epsilon<T, Policy>(); 117 for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++) 118 { 119 f = (k * f + p + q) / (k*k - v2); 120 p /= k - v; 121 q /= k + v; 122 g = f + e * q; 123 h = p - k * g; 124 coef *= coef_mult / k; 125 sum += coef * g; 126 sum1 += coef * h; 127 if (abs(coef * g) < abs(sum) * tolerance) 128 { 129 break; 130 } 131 } 132 policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in temme_jy", k, pol); 133 *Y = -sum; 134 *Y1 = -2 * sum1 / x; 135 136 return 0; 137 } 138 139 // Evaluate continued fraction fv = J_(v+1) / J_v, see 140 // Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73 141 template <typename T, typename Policy> CF1_jy(T v,T x,T * fv,int * sign,const Policy & pol)142 int CF1_jy(T v, T x, T* fv, int* sign, const Policy& pol) 143 { 144 T C, D, f, a, b, delta, tiny, tolerance; 145 unsigned long k; 146 int s = 1; 147 148 BOOST_MATH_STD_USING 149 150 // |x| <= |v|, CF1_jy converges rapidly 151 // |x| > |v|, CF1_jy needs O(|x|) iterations to converge 152 153 // modified Lentz's method, see 154 // Lentz, Applied Optics, vol 15, 668 (1976) 155 tolerance = 2 * policies::get_epsilon<T, Policy>(); 156 tiny = sqrt(tools::min_value<T>()); 157 C = f = tiny; // b0 = 0, replace with tiny 158 D = 0; 159 for (k = 1; k < policies::get_max_series_iterations<Policy>() * 100; k++) 160 { 161 a = -1; 162 b = 2 * (v + k) / x; 163 C = b + a / C; 164 D = b + a * D; 165 if (C == 0) { C = tiny; } 166 if (D == 0) { D = tiny; } 167 D = 1 / D; 168 delta = C * D; 169 f *= delta; 170 if (D < 0) { s = -s; } 171 if (abs(delta - 1) < tolerance) 172 { break; } 173 } 174 policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF1_jy", k / 100, pol); 175 *fv = -f; 176 *sign = s; // sign of denominator 177 178 return 0; 179 } 180 // 181 // This algorithm was originally written by Xiaogang Zhang 182 // using std::complex to perform the complex arithmetic. 183 // However, that turns out to 10x or more slower than using 184 // all real-valued arithmetic, so it's been rewritten using 185 // real values only. 186 // 187 template <typename T, typename Policy> CF2_jy(T v,T x,T * p,T * q,const Policy & pol)188 int CF2_jy(T v, T x, T* p, T* q, const Policy& pol) 189 { 190 BOOST_MATH_STD_USING 191 192 T Cr, Ci, Dr, Di, fr, fi, a, br, bi, delta_r, delta_i, temp; 193 T tiny; 194 unsigned long k; 195 196 // |x| >= |v|, CF2_jy converges rapidly 197 // |x| -> 0, CF2_jy fails to converge 198 BOOST_MATH_ASSERT(fabs(x) > 1); 199 200 // modified Lentz's method, complex numbers involved, see 201 // Lentz, Applied Optics, vol 15, 668 (1976) 202 T tolerance = 2 * policies::get_epsilon<T, Policy>(); 203 tiny = sqrt(tools::min_value<T>()); 204 Cr = fr = -0.5f / x; 205 Ci = fi = 1; 206 //Dr = Di = 0; 207 T v2 = v * v; 208 a = (0.25f - v2) / x; // Note complex this one time only! 209 br = 2 * x; 210 bi = 2; 211 temp = Cr * Cr + 1; 212 Ci = bi + a * Cr / temp; 213 Cr = br + a / temp; 214 Dr = br; 215 Di = bi; 216 if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; } 217 if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; } 218 temp = Dr * Dr + Di * Di; 219 Dr = Dr / temp; 220 Di = -Di / temp; 221 delta_r = Cr * Dr - Ci * Di; 222 delta_i = Ci * Dr + Cr * Di; 223 temp = fr; 224 fr = temp * delta_r - fi * delta_i; 225 fi = temp * delta_i + fi * delta_r; 226 for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++) 227 { 228 a = k - 0.5f; 229 a *= a; 230 a -= v2; 231 bi += 2; 232 temp = Cr * Cr + Ci * Ci; 233 Cr = br + a * Cr / temp; 234 Ci = bi - a * Ci / temp; 235 Dr = br + a * Dr; 236 Di = bi + a * Di; 237 if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; } 238 if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; } 239 temp = Dr * Dr + Di * Di; 240 Dr = Dr / temp; 241 Di = -Di / temp; 242 delta_r = Cr * Dr - Ci * Di; 243 delta_i = Ci * Dr + Cr * Di; 244 temp = fr; 245 fr = temp * delta_r - fi * delta_i; 246 fi = temp * delta_i + fi * delta_r; 247 if (fabs(delta_r - 1) + fabs(delta_i) < tolerance) 248 break; 249 } 250 policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol); 251 *p = fr; 252 *q = fi; 253 254 return 0; 255 } 256 257 static const int need_j = 1; 258 static const int need_y = 2; 259 260 // Compute J(v, x) and Y(v, x) simultaneously by Steed's method, see 261 // Barnett et al, Computer Physics Communications, vol 8, 377 (1974) 262 template <typename T, typename Policy> bessel_jy(T v,T x,T * J,T * Y,int kind,const Policy & pol)263 int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol) 264 { 265 BOOST_MATH_ASSERT(x >= 0); 266 267 T u, Jv, Ju, Yv, Yv1, Yu, Yu1(0), fv, fu; 268 T W, p, q, gamma, current, prev, next; 269 bool reflect = false; 270 unsigned n, k; 271 int s; 272 int org_kind = kind; 273 T cp = 0; 274 T sp = 0; 275 276 static const char* function = "boost::math::bessel_jy<%1%>(%1%,%1%)"; 277 278 BOOST_MATH_STD_USING 279 using namespace boost::math::tools; 280 using namespace boost::math::constants; 281 282 if (v < 0) 283 { 284 reflect = true; 285 v = -v; // v is non-negative from here 286 } 287 if (v > static_cast<T>((std::numeric_limits<int>::max)())) 288 { 289 *J = *Y = policies::raise_evaluation_error<T>(function, "Order of Bessel function is too large to evaluate: got %1%", v, pol); 290 return 1; 291 } 292 n = iround(v, pol); 293 u = v - n; // -1/2 <= u < 1/2 294 295 if(reflect) 296 { 297 T z = (u + n % 2); 298 cp = boost::math::cos_pi(z, pol); 299 sp = boost::math::sin_pi(z, pol); 300 if(u != 0) 301 kind = need_j|need_y; // need both for reflection formula 302 } 303 304 if(x == 0) 305 { 306 if(v == 0) 307 *J = 1; 308 else if((u == 0) || !reflect) 309 *J = 0; 310 else if(kind & need_j) 311 *J = policies::raise_domain_error<T>(function, "Value of Bessel J_v(x) is complex-infinity at %1%", x, pol); // complex infinity 312 else 313 *J = std::numeric_limits<T>::quiet_NaN(); // any value will do, not using J. 314 315 if((kind & need_y) == 0) 316 *Y = std::numeric_limits<T>::quiet_NaN(); // any value will do, not using Y. 317 else if(v == 0) 318 *Y = -policies::raise_overflow_error<T>(function, 0, pol); 319 else 320 *Y = policies::raise_domain_error<T>(function, "Value of Bessel Y_v(x) is complex-infinity at %1%", x, pol); // complex infinity 321 return 1; 322 } 323 324 // x is positive until reflection 325 W = T(2) / (x * pi<T>()); // Wronskian 326 T Yv_scale = 1; 327 if(((kind & need_y) == 0) && ((x < 1) || (v > x * x / 4) || (x < 5))) 328 { 329 // 330 // This series will actually converge rapidly for all small 331 // x - say up to x < 20 - but the first few terms are large 332 // and divergent which leads to large errors :-( 333 // 334 Jv = bessel_j_small_z_series(v, x, pol); 335 Yv = std::numeric_limits<T>::quiet_NaN(); 336 } 337 else if((x < 1) && (u != 0) && (log(policies::get_epsilon<T, Policy>() / 2) > v * log((x/2) * (x/2) / v))) 338 { 339 // Evaluate using series representations. 340 // This is particularly important for x << v as in this 341 // area temme_jy may be slow to converge, if it converges at all. 342 // Requires x is not an integer. 343 if(kind&need_j) 344 Jv = bessel_j_small_z_series(v, x, pol); 345 else 346 Jv = std::numeric_limits<T>::quiet_NaN(); 347 if((org_kind&need_y && (!reflect || (cp != 0))) 348 || (org_kind & need_j && (reflect && (sp != 0)))) 349 { 350 // Only calculate if we need it, and if the reflection formula will actually use it: 351 Yv = bessel_y_small_z_series(v, x, &Yv_scale, pol); 352 } 353 else 354 Yv = std::numeric_limits<T>::quiet_NaN(); 355 } 356 else if((u == 0) && (x < policies::get_epsilon<T, Policy>())) 357 { 358 // Truncated series evaluation for small x and v an integer, 359 // much quicker in this area than temme_jy below. 360 if(kind&need_j) 361 Jv = bessel_j_small_z_series(v, x, pol); 362 else 363 Jv = std::numeric_limits<T>::quiet_NaN(); 364 if((org_kind&need_y && (!reflect || (cp != 0))) 365 || (org_kind & need_j && (reflect && (sp != 0)))) 366 { 367 // Only calculate if we need it, and if the reflection formula will actually use it: 368 Yv = bessel_yn_small_z(n, x, &Yv_scale, pol); 369 } 370 else 371 Yv = std::numeric_limits<T>::quiet_NaN(); 372 } 373 else if(asymptotic_bessel_large_x_limit(v, x)) 374 { 375 if(kind&need_y) 376 { 377 Yv = asymptotic_bessel_y_large_x_2(v, x, pol); 378 } 379 else 380 Yv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it. 381 if(kind&need_j) 382 { 383 Jv = asymptotic_bessel_j_large_x_2(v, x, pol); 384 } 385 else 386 Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it. 387 } 388 else if((x > 8) && hankel_PQ(v, x, &p, &q, pol)) 389 { 390 // 391 // Hankel approximation: note that this method works best when x 392 // is large, but in that case we end up calculating sines and cosines 393 // of large values, with horrendous resulting accuracy. It is fast though 394 // when it works.... 395 // 396 // Normally we calculate sin/cos(chi) where: 397 // 398 // chi = x - fmod(T(v / 2 + 0.25f), T(2)) * boost::math::constants::pi<T>(); 399 // 400 // But this introduces large errors, so use sin/cos addition formulae to 401 // improve accuracy: 402 // 403 T mod_v = fmod(T(v / 2 + 0.25f), T(2)); 404 T sx = sin(x); 405 T cx = cos(x); 406 T sv = boost::math::sin_pi(mod_v, pol); 407 T cv = boost::math::cos_pi(mod_v, pol); 408 409 T sc = sx * cv - sv * cx; // == sin(chi); 410 T cc = cx * cv + sx * sv; // == cos(chi); 411 T chi = boost::math::constants::root_two<T>() / (boost::math::constants::root_pi<T>() * sqrt(x)); //sqrt(2 / (boost::math::constants::pi<T>() * x)); 412 Yv = chi * (p * sc + q * cc); 413 Jv = chi * (p * cc - q * sc); 414 } 415 else if (x <= 2) // x in (0, 2] 416 { 417 if(temme_jy(u, x, &Yu, &Yu1, pol)) // Temme series 418 { 419 // domain error: 420 *J = *Y = Yu; 421 return 1; 422 } 423 prev = Yu; 424 current = Yu1; 425 T scale = 1; 426 policies::check_series_iterations<T>(function, n, pol); 427 for (k = 1; k <= n; k++) // forward recurrence for Y 428 { 429 T fact = 2 * (u + k) / x; 430 if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current)) 431 { 432 scale /= current; 433 prev /= current; 434 current = 1; 435 } 436 next = fact * current - prev; 437 prev = current; 438 current = next; 439 } 440 Yv = prev; 441 Yv1 = current; 442 if(kind&need_j) 443 { 444 CF1_jy(v, x, &fv, &s, pol); // continued fraction CF1_jy 445 Jv = scale * W / (Yv * fv - Yv1); // Wronskian relation 446 } 447 else 448 Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it. 449 Yv_scale = scale; 450 } 451 else // x in (2, \infty) 452 { 453 // Get Y(u, x): 454 455 T ratio; 456 CF1_jy(v, x, &fv, &s, pol); 457 // tiny initial value to prevent overflow 458 T init = sqrt(tools::min_value<T>()); 459 BOOST_MATH_INSTRUMENT_VARIABLE(init); 460 prev = fv * s * init; 461 current = s * init; 462 if(v < max_factorial<T>::value) 463 { 464 policies::check_series_iterations<T>(function, n, pol); 465 for (k = n; k > 0; k--) // backward recurrence for J 466 { 467 next = 2 * (u + k) * current / x - prev; 468 prev = current; 469 current = next; 470 } 471 ratio = (s * init) / current; // scaling ratio 472 // can also call CF1_jy() to get fu, not much difference in precision 473 fu = prev / current; 474 } 475 else 476 { 477 // 478 // When v is large we may get overflow in this calculation 479 // leading to NaN's and other nasty surprises: 480 // 481 policies::check_series_iterations<T>(function, n, pol); 482 bool over = false; 483 for (k = n; k > 0; k--) // backward recurrence for J 484 { 485 T t = 2 * (u + k) / x; 486 if((t > 1) && (tools::max_value<T>() / t < current)) 487 { 488 over = true; 489 break; 490 } 491 next = t * current - prev; 492 prev = current; 493 current = next; 494 } 495 if(!over) 496 { 497 ratio = (s * init) / current; // scaling ratio 498 // can also call CF1_jy() to get fu, not much difference in precision 499 fu = prev / current; 500 } 501 else 502 { 503 ratio = 0; 504 fu = 1; 505 } 506 } 507 CF2_jy(u, x, &p, &q, pol); // continued fraction CF2_jy 508 T t = u / x - fu; // t = J'/J 509 gamma = (p - t) / q; 510 // 511 // We can't allow gamma to cancel out to zero completely as it messes up 512 // the subsequent logic. So pretend that one bit didn't cancel out 513 // and set to a suitably small value. The only test case we've been able to 514 // find for this, is when v = 8.5 and x = 4*PI. 515 // 516 if(gamma == 0) 517 { 518 gamma = u * tools::epsilon<T>() / x; 519 } 520 BOOST_MATH_INSTRUMENT_VARIABLE(current); 521 BOOST_MATH_INSTRUMENT_VARIABLE(W); 522 BOOST_MATH_INSTRUMENT_VARIABLE(q); 523 BOOST_MATH_INSTRUMENT_VARIABLE(gamma); 524 BOOST_MATH_INSTRUMENT_VARIABLE(p); 525 BOOST_MATH_INSTRUMENT_VARIABLE(t); 526 Ju = sign(current) * sqrt(W / (q + gamma * (p - t))); 527 BOOST_MATH_INSTRUMENT_VARIABLE(Ju); 528 529 Jv = Ju * ratio; // normalization 530 531 Yu = gamma * Ju; 532 Yu1 = Yu * (u/x - p - q/gamma); 533 534 if(kind&need_y) 535 { 536 // compute Y: 537 prev = Yu; 538 current = Yu1; 539 policies::check_series_iterations<T>(function, n, pol); 540 for (k = 1; k <= n; k++) // forward recurrence for Y 541 { 542 T fact = 2 * (u + k) / x; 543 if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current)) 544 { 545 prev /= current; 546 Yv_scale /= current; 547 current = 1; 548 } 549 next = fact * current - prev; 550 prev = current; 551 current = next; 552 } 553 Yv = prev; 554 } 555 else 556 Yv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it. 557 } 558 559 if (reflect) 560 { 561 if((sp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(sp * Yv))) 562 *J = org_kind & need_j ? T(-sign(sp) * sign(Yv) * (Yv_scale != 0 ? sign(Yv_scale) : 1) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0); 563 else 564 *J = cp * Jv - (sp == 0 ? T(0) : T((sp * Yv) / Yv_scale)); // reflection formula 565 if((cp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(cp * Yv))) 566 *Y = org_kind & need_y ? T(-sign(cp) * sign(Yv) * (Yv_scale != 0 ? sign(Yv_scale) : 1) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0); 567 else 568 *Y = (sp != 0 ? sp * Jv : T(0)) + (cp == 0 ? T(0) : T((cp * Yv) / Yv_scale)); 569 } 570 else 571 { 572 *J = Jv; 573 if(tools::max_value<T>() * fabs(Yv_scale) < fabs(Yv)) 574 *Y = org_kind & need_y ? T(sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0); 575 else 576 *Y = Yv / Yv_scale; 577 } 578 579 return 0; 580 } 581 582 } // namespace detail 583 584 }} // namespaces 585 586 #endif // BOOST_MATH_BESSEL_JY_HPP 587