]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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> | |
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 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 |