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

source code of boost/libs/math/include/boost/math/special_functions/detail/bessel_ik.hpp