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
cyl_bessel_i_small_zboost::math::detail::cyl_bessel_i_small_z31 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
operator ()boost::math::detail::cyl_bessel_i_small_z37 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>
bessel_i_small_z_series(T v,T x,const Policy & pol)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>
temme_ik(T v,T x,T * K,T * K1,const Policy & pol)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>
CF1_ik(T v,T x,T * fv,const Policy & pol)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>
CF2_ik(T v,T x,T * Kv,T * Kv1,const Policy & pol)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>
bessel_ik(T v,T x,T * I,T * K,int kind,const Policy & pol)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 for (k = 1; k <= n; k++) // forward recurrence for K
378 {
379 T fact = 2 * (u + k) / x;
380 if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
381 {
382 prev /= current;
383 scale /= current;
384 current = 1;
385 }
386 next = fact * current + prev;
387 prev = current;
388 current = next;
389 }
390 Kv = prev;
391 Kv1 = current;
392 BOOST_MATH_INSTRUMENT_VARIABLE(Kv);
393 BOOST_MATH_INSTRUMENT_VARIABLE(Kv1);
394 if(kind & need_i)
395 {
396 T lim = (4 * v * v + 10) / (8 * x);
397 lim *= lim;
398 lim *= lim;
399 lim /= 24;
400 if((lim < tools::epsilon<T>() * 10) && (x > 100))
401 {
402 // x is huge compared to v, CF1 may be very slow
403 // to converge so use asymptotic expansion for large
404 // x case instead. Note that the asymptotic expansion
405 // isn't very accurate - so it's deliberately very hard
406 // to get here - probably we're going to overflow:
407 Iv = asymptotic_bessel_i_large_x(v, x, pol);
408 }
409 else if((v > 0) && (x / v < 0.25))
410 {
411 Iv = bessel_i_small_z_series(v, x, pol);
412 }
413 else
414 {
415 CF1_ik(v, x, &fv, pol); // continued fraction CF1_ik
416 Iv = scale * W / (Kv * fv + Kv1); // Wronskian relation
417 }
418 }
419 else
420 Iv = std::numeric_limits<T>::quiet_NaN(); // any value will do
421
422 if (reflect)
423 {
424 T z = (u + n % 2);
425 T fact = (2 / pi<T>()) * (boost::math::sin_pi(z) * Kv);
426 if(fact == 0)
427 *I = Iv;
428 else if(tools::max_value<T>() * scale < fact)
429 *I = (org_kind & need_i) ? T(sign(fact) * sign(scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
430 else
431 *I = Iv + fact / scale; // reflection formula
432 }
433 else
434 {
435 *I = Iv;
436 }
437 if(tools::max_value<T>() * scale < Kv)
438 *K = (org_kind & need_k) ? T(sign(Kv) * sign(scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
439 else
440 *K = Kv / scale;
441 BOOST_MATH_INSTRUMENT_VARIABLE(*I);
442 BOOST_MATH_INSTRUMENT_VARIABLE(*K);
443 return 0;
444 }
445
446 }}} // namespaces
447
448 #endif // BOOST_MATH_BESSEL_IK_HPP
449
450