]>
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_IK_HPP | |
7 | #define BOOST_MATH_BESSEL_IK_HPP | |
8 | ||
9 | #ifdef _MSC_VER | |
10 | #pragma once | |
11 | #endif | |
12 | ||
13 | #include <boost/math/special_functions/round.hpp> | |
14 | #include <boost/math/special_functions/gamma.hpp> | |
15 | #include <boost/math/special_functions/sin_pi.hpp> | |
16 | #include <boost/math/constants/constants.hpp> | |
17 | #include <boost/math/policies/error_handling.hpp> | |
18 | #include <boost/math/tools/config.hpp> | |
19 | ||
20 | // Modified Bessel functions of the first and second kind of fractional order | |
21 | ||
22 | namespace boost { namespace math { | |
23 | ||
24 | namespace detail { | |
25 | ||
26 | template <class T, class Policy> | |
27 | struct cyl_bessel_i_small_z | |
28 | { | |
29 | typedef T result_type; | |
30 | ||
31 | cyl_bessel_i_small_z(T v_, T z_) : k(0), v(v_), mult(z_*z_/4) | |
32 | { | |
33 | BOOST_MATH_STD_USING | |
34 | term = 1; | |
35 | } | |
36 | ||
37 | T operator()() | |
38 | { | |
39 | T result = term; | |
40 | ++k; | |
41 | term *= mult / k; | |
42 | term /= k + v; | |
43 | return result; | |
44 | } | |
45 | private: | |
46 | unsigned k; | |
47 | T v; | |
48 | T term; | |
49 | T mult; | |
50 | }; | |
51 | ||
52 | template <class T, class Policy> | |
53 | inline T bessel_i_small_z_series(T v, T x, const Policy& pol) | |
54 | { | |
55 | BOOST_MATH_STD_USING | |
56 | T prefix; | |
57 | if(v < max_factorial<T>::value) | |
58 | { | |
59 | prefix = pow(x / 2, v) / boost::math::tgamma(v + 1, pol); | |
60 | } | |
61 | else | |
62 | { | |
63 | prefix = v * log(x / 2) - boost::math::lgamma(v + 1, pol); | |
64 | prefix = exp(prefix); | |
65 | } | |
66 | if(prefix == 0) | |
67 | return prefix; | |
68 | ||
69 | cyl_bessel_i_small_z<T, Policy> s(v, x); | |
70 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
71 | #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) | |
72 | T zero = 0; | |
73 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero); | |
74 | #else | |
75 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
76 | #endif | |
77 | policies::check_series_iterations<T>("boost::math::bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol); | |
78 | return prefix * result; | |
79 | } | |
80 | ||
81 | // Calculate K(v, x) and K(v+1, x) by method analogous to | |
82 | // Temme, Journal of Computational Physics, vol 21, 343 (1976) | |
83 | template <typename T, typename Policy> | |
84 | int temme_ik(T v, T x, T* K, T* K1, const Policy& pol) | |
85 | { | |
86 | T f, h, p, q, coef, sum, sum1, tolerance; | |
87 | T a, b, c, d, sigma, gamma1, gamma2; | |
88 | unsigned long k; | |
89 | ||
90 | BOOST_MATH_STD_USING | |
91 | using namespace boost::math::tools; | |
92 | using namespace boost::math::constants; | |
93 | ||
94 | ||
95 | // |x| <= 2, Temme series converge rapidly | |
96 | // |x| > 2, the larger the |x|, the slower the convergence | |
97 | BOOST_ASSERT(abs(x) <= 2); | |
98 | BOOST_ASSERT(abs(v) <= 0.5f); | |
99 | ||
100 | T gp = boost::math::tgamma1pm1(v, pol); | |
101 | T gm = boost::math::tgamma1pm1(-v, pol); | |
102 | ||
103 | a = log(x / 2); | |
104 | b = exp(v * a); | |
105 | sigma = -a * v; | |
106 | c = abs(v) < tools::epsilon<T>() ? | |
107 | T(1) : T(boost::math::sin_pi(v) / (v * pi<T>())); | |
108 | d = abs(sigma) < tools::epsilon<T>() ? | |
109 | T(1) : T(sinh(sigma) / sigma); | |
110 | gamma1 = abs(v) < tools::epsilon<T>() ? | |
111 | T(-euler<T>()) : T((0.5f / v) * (gp - gm) * c); | |
112 | gamma2 = (2 + gp + gm) * c / 2; | |
113 | ||
114 | // initial values | |
115 | p = (gp + 1) / (2 * b); | |
116 | q = (1 + gm) * b / 2; | |
117 | f = (cosh(sigma) * gamma1 + d * (-a) * gamma2) / c; | |
118 | h = p; | |
119 | coef = 1; | |
120 | sum = coef * f; | |
121 | sum1 = coef * h; | |
122 | ||
123 | BOOST_MATH_INSTRUMENT_VARIABLE(p); | |
124 | BOOST_MATH_INSTRUMENT_VARIABLE(q); | |
125 | BOOST_MATH_INSTRUMENT_VARIABLE(f); | |
126 | BOOST_MATH_INSTRUMENT_VARIABLE(sigma); | |
127 | BOOST_MATH_INSTRUMENT_CODE(sinh(sigma)); | |
128 | BOOST_MATH_INSTRUMENT_VARIABLE(gamma1); | |
129 | BOOST_MATH_INSTRUMENT_VARIABLE(gamma2); | |
130 | BOOST_MATH_INSTRUMENT_VARIABLE(c); | |
131 | BOOST_MATH_INSTRUMENT_VARIABLE(d); | |
132 | BOOST_MATH_INSTRUMENT_VARIABLE(a); | |
133 | ||
134 | // series summation | |
135 | tolerance = tools::epsilon<T>(); | |
136 | for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++) | |
137 | { | |
138 | f = (k * f + p + q) / (k*k - v*v); | |
139 | p /= k - v; | |
140 | q /= k + v; | |
141 | h = p - k * f; | |
142 | coef *= x * x / (4 * k); | |
143 | sum += coef * f; | |
144 | sum1 += coef * h; | |
145 | if (abs(coef * f) < abs(sum) * tolerance) | |
146 | { | |
147 | break; | |
148 | } | |
149 | } | |
150 | policies::check_series_iterations<T>("boost::math::bessel_ik<%1%>(%1%,%1%) in temme_ik", k, pol); | |
151 | ||
152 | *K = sum; | |
153 | *K1 = 2 * sum1 / x; | |
154 | ||
155 | return 0; | |
156 | } | |
157 | ||
158 | // Evaluate continued fraction fv = I_(v+1) / I_v, derived from | |
159 | // Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73 | |
160 | template <typename T, typename Policy> | |
161 | int CF1_ik(T v, T x, T* fv, const Policy& pol) | |
162 | { | |
163 | T C, D, f, a, b, delta, tiny, tolerance; | |
164 | unsigned long k; | |
165 | ||
166 | BOOST_MATH_STD_USING | |
167 | ||
168 | // |x| <= |v|, CF1_ik converges rapidly | |
169 | // |x| > |v|, CF1_ik needs O(|x|) iterations to converge | |
170 | ||
171 | // modified Lentz's method, see | |
172 | // Lentz, Applied Optics, vol 15, 668 (1976) | |
173 | tolerance = 2 * tools::epsilon<T>(); | |
174 | BOOST_MATH_INSTRUMENT_VARIABLE(tolerance); | |
175 | tiny = sqrt(tools::min_value<T>()); | |
176 | BOOST_MATH_INSTRUMENT_VARIABLE(tiny); | |
177 | C = f = tiny; // b0 = 0, replace with tiny | |
178 | D = 0; | |
179 | for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++) | |
180 | { | |
181 | a = 1; | |
182 | b = 2 * (v + k) / x; | |
183 | C = b + a / C; | |
184 | D = b + a * D; | |
185 | if (C == 0) { C = tiny; } | |
186 | if (D == 0) { D = tiny; } | |
187 | D = 1 / D; | |
188 | delta = C * D; | |
189 | f *= delta; | |
190 | BOOST_MATH_INSTRUMENT_VARIABLE(delta-1); | |
191 | if (abs(delta - 1) <= tolerance) | |
192 | { | |
193 | break; | |
194 | } | |
195 | } | |
196 | BOOST_MATH_INSTRUMENT_VARIABLE(k); | |
197 | policies::check_series_iterations<T>("boost::math::bessel_ik<%1%>(%1%,%1%) in CF1_ik", k, pol); | |
198 | ||
199 | *fv = f; | |
200 | ||
201 | return 0; | |
202 | } | |
203 | ||
204 | // Calculate K(v, x) and K(v+1, x) by evaluating continued fraction | |
205 | // z1 / z0 = U(v+1.5, 2v+1, 2x) / U(v+0.5, 2v+1, 2x), see | |
206 | // Thompson and Barnett, Computer Physics Communications, vol 47, 245 (1987) | |
207 | template <typename T, typename Policy> | |
208 | int CF2_ik(T v, T x, T* Kv, T* Kv1, const Policy& pol) | |
209 | { | |
210 | BOOST_MATH_STD_USING | |
211 | using namespace boost::math::constants; | |
212 | ||
213 | T S, C, Q, D, f, a, b, q, delta, tolerance, current, prev; | |
214 | unsigned long k; | |
215 | ||
216 | // |x| >= |v|, CF2_ik converges rapidly | |
217 | // |x| -> 0, CF2_ik fails to converge | |
218 | ||
219 | BOOST_ASSERT(abs(x) > 1); | |
220 | ||
221 | // Steed's algorithm, see Thompson and Barnett, | |
222 | // Journal of Computational Physics, vol 64, 490 (1986) | |
223 | tolerance = tools::epsilon<T>(); | |
224 | a = v * v - 0.25f; | |
225 | b = 2 * (x + 1); // b1 | |
226 | D = 1 / b; // D1 = 1 / b1 | |
227 | f = delta = D; // f1 = delta1 = D1, coincidence | |
228 | prev = 0; // q0 | |
229 | current = 1; // q1 | |
230 | Q = C = -a; // Q1 = C1 because q1 = 1 | |
231 | S = 1 + Q * delta; // S1 | |
232 | BOOST_MATH_INSTRUMENT_VARIABLE(tolerance); | |
233 | BOOST_MATH_INSTRUMENT_VARIABLE(a); | |
234 | BOOST_MATH_INSTRUMENT_VARIABLE(b); | |
235 | BOOST_MATH_INSTRUMENT_VARIABLE(D); | |
236 | BOOST_MATH_INSTRUMENT_VARIABLE(f); | |
237 | ||
238 | for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++) // starting from 2 | |
239 | { | |
240 | // continued fraction f = z1 / z0 | |
241 | a -= 2 * (k - 1); | |
242 | b += 2; | |
243 | D = 1 / (b + a * D); | |
244 | delta *= b * D - 1; | |
245 | f += delta; | |
246 | ||
247 | // series summation S = 1 + \sum_{n=1}^{\infty} C_n * z_n / z_0 | |
248 | q = (prev - (b - 2) * current) / a; | |
249 | prev = current; | |
250 | current = q; // forward recurrence for q | |
251 | C *= -a / k; | |
252 | Q += C * q; | |
253 | S += Q * delta; | |
254 | // | |
255 | // Under some circumstances q can grow very small and C very | |
256 | // large, leading to under/overflow. This is particularly an | |
257 | // issue for types which have many digits precision but a narrow | |
258 | // exponent range. A typical example being a "double double" type. | |
259 | // To avoid this situation we can normalise q (and related prev/current) | |
260 | // and C. All other variables remain unchanged in value. A typical | |
261 | // test case occurs when x is close to 2, for example cyl_bessel_k(9.125, 2.125). | |
262 | // | |
263 | if(q < tools::epsilon<T>()) | |
264 | { | |
265 | C *= q; | |
266 | prev /= q; | |
267 | current /= q; | |
268 | q = 1; | |
269 | } | |
270 | ||
271 | // S converges slower than f | |
272 | BOOST_MATH_INSTRUMENT_VARIABLE(Q * delta); | |
273 | BOOST_MATH_INSTRUMENT_VARIABLE(abs(S) * tolerance); | |
274 | BOOST_MATH_INSTRUMENT_VARIABLE(S); | |
275 | if (abs(Q * delta) < abs(S) * tolerance) | |
276 | { | |
277 | break; | |
278 | } | |
279 | } | |
280 | policies::check_series_iterations<T>("boost::math::bessel_ik<%1%>(%1%,%1%) in CF2_ik", k, pol); | |
281 | ||
282 | if(x >= tools::log_max_value<T>()) | |
283 | *Kv = exp(0.5f * log(pi<T>() / (2 * x)) - x - log(S)); | |
284 | else | |
285 | *Kv = sqrt(pi<T>() / (2 * x)) * exp(-x) / S; | |
286 | *Kv1 = *Kv * (0.5f + v + x + (v * v - 0.25f) * f) / x; | |
287 | BOOST_MATH_INSTRUMENT_VARIABLE(*Kv); | |
288 | BOOST_MATH_INSTRUMENT_VARIABLE(*Kv1); | |
289 | ||
290 | return 0; | |
291 | } | |
292 | ||
293 | enum{ | |
294 | need_i = 1, | |
295 | need_k = 2 | |
296 | }; | |
297 | ||
298 | // Compute I(v, x) and K(v, x) simultaneously by Temme's method, see | |
299 | // Temme, Journal of Computational Physics, vol 19, 324 (1975) | |
300 | template <typename T, typename Policy> | |
301 | int bessel_ik(T v, T x, T* I, T* K, int kind, const Policy& pol) | |
302 | { | |
303 | // Kv1 = K_(v+1), fv = I_(v+1) / I_v | |
304 | // Ku1 = K_(u+1), fu = I_(u+1) / I_u | |
305 | T u, Iv, Kv, Kv1, Ku, Ku1, fv; | |
306 | T W, current, prev, next; | |
307 | bool reflect = false; | |
308 | unsigned n, k; | |
309 | int org_kind = kind; | |
310 | BOOST_MATH_INSTRUMENT_VARIABLE(v); | |
311 | BOOST_MATH_INSTRUMENT_VARIABLE(x); | |
312 | BOOST_MATH_INSTRUMENT_VARIABLE(kind); | |
313 | ||
314 | BOOST_MATH_STD_USING | |
315 | using namespace boost::math::tools; | |
316 | using namespace boost::math::constants; | |
317 | ||
318 | static const char* function = "boost::math::bessel_ik<%1%>(%1%,%1%)"; | |
319 | ||
320 | if (v < 0) | |
321 | { | |
322 | reflect = true; | |
323 | v = -v; // v is non-negative from here | |
324 | kind |= need_k; | |
325 | } | |
326 | n = iround(v, pol); | |
327 | u = v - n; // -1/2 <= u < 1/2 | |
328 | BOOST_MATH_INSTRUMENT_VARIABLE(n); | |
329 | BOOST_MATH_INSTRUMENT_VARIABLE(u); | |
330 | ||
331 | if (x < 0) | |
332 | { | |
333 | *I = *K = policies::raise_domain_error<T>(function, | |
334 | "Got x = %1% but real argument x must be non-negative, complex number result not supported.", x, pol); | |
335 | return 1; | |
336 | } | |
337 | if (x == 0) | |
338 | { | |
339 | Iv = (v == 0) ? static_cast<T>(1) : static_cast<T>(0); | |
340 | if(kind & need_k) | |
341 | { | |
342 | Kv = policies::raise_overflow_error<T>(function, 0, pol); | |
343 | } | |
344 | else | |
345 | { | |
346 | Kv = std::numeric_limits<T>::quiet_NaN(); // any value will do | |
347 | } | |
348 | ||
349 | if(reflect && (kind & need_i)) | |
350 | { | |
351 | T z = (u + n % 2); | |
352 | Iv = boost::math::sin_pi(z, pol) == 0 ? | |
353 | Iv : | |
354 | policies::raise_overflow_error<T>(function, 0, pol); // reflection formula | |
355 | } | |
356 | ||
357 | *I = Iv; | |
358 | *K = Kv; | |
359 | return 0; | |
360 | } | |
361 | ||
362 | // x is positive until reflection | |
363 | W = 1 / x; // Wronskian | |
364 | if (x <= 2) // x in (0, 2] | |
365 | { | |
366 | temme_ik(u, x, &Ku, &Ku1, pol); // Temme series | |
367 | } | |
368 | else // x in (2, \infty) | |
369 | { | |
370 | CF2_ik(u, x, &Ku, &Ku1, pol); // continued fraction CF2_ik | |
371 | } | |
372 | BOOST_MATH_INSTRUMENT_VARIABLE(Ku); | |
373 | BOOST_MATH_INSTRUMENT_VARIABLE(Ku1); | |
374 | prev = Ku; | |
375 | current = Ku1; | |
376 | T scale = 1; | |
377 | T scale_sign = 1; | |
378 | for (k = 1; k <= n; k++) // forward recurrence for K | |
379 | { | |
380 | T fact = 2 * (u + k) / x; | |
381 | if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current)) | |
382 | { | |
383 | prev /= current; | |
384 | scale /= current; | |
385 | scale_sign *= boost::math::sign(current); | |
386 | current = 1; | |
387 | } | |
388 | next = fact * current + prev; | |
389 | prev = current; | |
390 | current = next; | |
391 | } | |
392 | Kv = prev; | |
393 | Kv1 = current; | |
394 | BOOST_MATH_INSTRUMENT_VARIABLE(Kv); | |
395 | BOOST_MATH_INSTRUMENT_VARIABLE(Kv1); | |
396 | if(kind & need_i) | |
397 | { | |
398 | T lim = (4 * v * v + 10) / (8 * x); | |
399 | lim *= lim; | |
400 | lim *= lim; | |
401 | lim /= 24; | |
402 | if((lim < tools::epsilon<T>() * 10) && (x > 100)) | |
403 | { | |
404 | // x is huge compared to v, CF1 may be very slow | |
405 | // to converge so use asymptotic expansion for large | |
406 | // x case instead. Note that the asymptotic expansion | |
407 | // isn't very accurate - so it's deliberately very hard | |
408 | // to get here - probably we're going to overflow: | |
409 | Iv = asymptotic_bessel_i_large_x(v, x, pol); | |
410 | } | |
411 | else if((v > 0) && (x / v < 0.25)) | |
412 | { | |
413 | Iv = bessel_i_small_z_series(v, x, pol); | |
414 | } | |
415 | else | |
416 | { | |
417 | CF1_ik(v, x, &fv, pol); // continued fraction CF1_ik | |
418 | Iv = scale * W / (Kv * fv + Kv1); // Wronskian relation | |
419 | } | |
420 | } | |
421 | else | |
422 | Iv = std::numeric_limits<T>::quiet_NaN(); // any value will do | |
423 | ||
424 | if (reflect) | |
425 | { | |
426 | T z = (u + n % 2); | |
427 | T fact = (2 / pi<T>()) * (boost::math::sin_pi(z) * Kv); | |
428 | if(fact == 0) | |
429 | *I = Iv; | |
430 | else if(tools::max_value<T>() * scale < fact) | |
431 | *I = (org_kind & need_i) ? T(sign(fact) * scale_sign * policies::raise_overflow_error<T>(function, 0, pol)) : T(0); | |
432 | else | |
433 | *I = Iv + fact / scale; // reflection formula | |
434 | } | |
435 | else | |
436 | { | |
437 | *I = Iv; | |
438 | } | |
439 | if(tools::max_value<T>() * scale < Kv) | |
440 | *K = (org_kind & need_k) ? T(sign(Kv) * scale_sign * policies::raise_overflow_error<T>(function, 0, pol)) : T(0); | |
441 | else | |
442 | *K = Kv / scale; | |
443 | BOOST_MATH_INSTRUMENT_VARIABLE(*I); | |
444 | BOOST_MATH_INSTRUMENT_VARIABLE(*K); | |
445 | return 0; | |
446 | } | |
447 | ||
448 | }}} // namespaces | |
449 | ||
450 | #endif // BOOST_MATH_BESSEL_IK_HPP | |
451 |