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
27namespace 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>
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>
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>
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>
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>
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, nullptr, 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 //
469 // We can't allow next to completely cancel out or the subsequent logic breaks.
470 // Pretend that one bit did not cancel:
471 if (next == 0)
472 {
473 next = prev * tools::epsilon<T>() / 2;
474 }
475 prev = current;
476 current = next;
477 }
478 ratio = (s * init) / current; // scaling ratio
479 // can also call CF1_jy() to get fu, not much difference in precision
480 fu = prev / current;
481 }
482 else
483 {
484 //
485 // When v is large we may get overflow in this calculation
486 // leading to NaN's and other nasty surprises:
487 //
488 policies::check_series_iterations<T>(function, n, pol);
489 bool over = false;
490 for (k = n; k > 0; k--) // backward recurrence for J
491 {
492 T t = 2 * (u + k) / x;
493 if((t > 1) && (tools::max_value<T>() / t < current))
494 {
495 over = true;
496 break;
497 }
498 next = t * current - prev;
499 prev = current;
500 current = next;
501 }
502 if(!over)
503 {
504 ratio = (s * init) / current; // scaling ratio
505 // can also call CF1_jy() to get fu, not much difference in precision
506 fu = prev / current;
507 }
508 else
509 {
510 ratio = 0;
511 fu = 1;
512 }
513 }
514 CF2_jy(u, x, &p, &q, pol); // continued fraction CF2_jy
515 T t = u / x - fu; // t = J'/J
516 gamma = (p - t) / q;
517 //
518 // We can't allow gamma to cancel out to zero completely as it messes up
519 // the subsequent logic. So pretend that one bit didn't cancel out
520 // and set to a suitably small value. The only test case we've been able to
521 // find for this, is when v = 8.5 and x = 4*PI.
522 //
523 if(gamma == 0)
524 {
525 gamma = u * tools::epsilon<T>() / x;
526 }
527 BOOST_MATH_INSTRUMENT_VARIABLE(current);
528 BOOST_MATH_INSTRUMENT_VARIABLE(W);
529 BOOST_MATH_INSTRUMENT_VARIABLE(q);
530 BOOST_MATH_INSTRUMENT_VARIABLE(gamma);
531 BOOST_MATH_INSTRUMENT_VARIABLE(p);
532 BOOST_MATH_INSTRUMENT_VARIABLE(t);
533 Ju = sign(current) * sqrt(W / (q + gamma * (p - t)));
534 BOOST_MATH_INSTRUMENT_VARIABLE(Ju);
535
536 Jv = Ju * ratio; // normalization
537
538 Yu = gamma * Ju;
539 Yu1 = Yu * (u/x - p - q/gamma);
540
541 if(kind&need_y)
542 {
543 // compute Y:
544 prev = Yu;
545 current = Yu1;
546 policies::check_series_iterations<T>(function, n, pol);
547 for (k = 1; k <= n; k++) // forward recurrence for Y
548 {
549 T fact = 2 * (u + k) / x;
550 if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
551 {
552 prev /= current;
553 Yv_scale /= current;
554 current = 1;
555 }
556 next = fact * current - prev;
557 prev = current;
558 current = next;
559 }
560 Yv = prev;
561 }
562 else
563 Yv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
564 }
565
566 if (reflect)
567 {
568 if((sp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(sp * Yv)))
569 *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);
570 else
571 *J = cp * Jv - (sp == 0 ? T(0) : T((sp * Yv) / Yv_scale)); // reflection formula
572 if((cp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(cp * Yv)))
573 *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);
574 else
575 *Y = (sp != 0 ? sp * Jv : T(0)) + (cp == 0 ? T(0) : T((cp * Yv) / Yv_scale));
576 }
577 else
578 {
579 *J = Jv;
580 if(tools::max_value<T>() * fabs(Yv_scale) < fabs(Yv))
581 *Y = org_kind & need_y ? T(sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, nullptr, pol)) : T(0);
582 else
583 *Y = Yv / Yv_scale;
584 }
585
586 return 0;
587 }
588
589 } // namespace detail
590
591}} // namespaces
592
593#endif // BOOST_MATH_BESSEL_JY_HPP
594

source code of include/boost/math/special_functions/detail/bessel_jy.hpp