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
29namespace 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>
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>
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>
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>
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>
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 completely 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) * (Yv_scale != 0 ? sign(Yv_scale) : 1) * 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) * (Yv_scale != 0 ? sign(Yv_scale) : 1) * 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

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