// Copyright (c) 2006 Xiaogang Zhang | |
// Use, modification and distribution are subject to the | |
// Boost Software License, Version 1.0. (See accompanying file | |
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
#ifndef BOOST_MATH_BESSEL_IK_HPP | |
#define BOOST_MATH_BESSEL_IK_HPP | |
#ifdef _MSC_VER | |
#pragma once | |
#endif | |
#include <boost/math/special_functions/round.hpp> | |
#include <boost/math/special_functions/gamma.hpp> | |
#include <boost/math/special_functions/sin_pi.hpp> | |
#include <boost/math/constants/constants.hpp> | |
#include <boost/math/policies/error_handling.hpp> | |
#include <boost/math/tools/config.hpp> | |
// Modified Bessel functions of the first and second kind of fractional order | |
namespace boost { namespace math { | |
namespace detail { | |
// Calculate K(v, x) and K(v+1, x) by method analogous to | |
// Temme, Journal of Computational Physics, vol 21, 343 (1976) | |
template <typename T, typename Policy> | |
int temme_ik(T v, T x, T* K, T* K1, const Policy& pol) | |
{ | |
T f, h, p, q, coef, sum, sum1, tolerance; | |
T a, b, c, d, sigma, gamma1, gamma2; | |
unsigned long k; | |
BOOST_MATH_STD_USING | |
using namespace boost::math::tools; | |
using namespace boost::math::constants; | |
// |x| <= 2, Temme series converge rapidly | |
// |x| > 2, the larger the |x|, the slower the convergence | |
BOOST_ASSERT(abs(x) <= 2); | |
BOOST_ASSERT(abs(v) <= 0.5f); | |
T gp = boost::math::tgamma1pm1(v, pol); | |
T gm = boost::math::tgamma1pm1(-v, pol); | |
a = log(x / 2); | |
b = exp(v * a); | |
sigma = -a * v; | |
c = abs(v) < tools::epsilon<T>() ? | |
T(1) : T(boost::math::sin_pi(v) / (v * pi<T>())); | |
d = abs(sigma) < tools::epsilon<T>() ? | |
T(1) : T(sinh(sigma) / sigma); | |
gamma1 = abs(v) < tools::epsilon<T>() ? | |
T(-euler<T>()) : T((0.5f / v) * (gp - gm) * c); | |
gamma2 = (2 + gp + gm) * c / 2; | |
// initial values | |
p = (gp + 1) / (2 * b); | |
q = (1 + gm) * b / 2; | |
f = (cosh(sigma) * gamma1 + d * (-a) * gamma2) / c; | |
h = p; | |
coef = 1; | |
sum = coef * f; | |
sum1 = coef * h; | |
// series summation | |
tolerance = tools::epsilon<T>(); | |
for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++) | |
{ | |
f = (k * f + p + q) / (k*k - v*v); | |
p /= k - v; | |
q /= k + v; | |
h = p - k * f; | |
coef *= x * x / (4 * k); | |
sum += coef * f; | |
sum1 += coef * h; | |
if (abs(coef * f) < abs(sum) * tolerance) | |
{ | |
break; | |
} | |
} | |
policies::check_series_iterations("boost::math::bessel_ik<%1%>(%1%,%1%) in temme_ik", k, pol); | |
*K = sum; | |
*K1 = 2 * sum1 / x; | |
return 0; | |
} | |
// Evaluate continued fraction fv = I_(v+1) / I_v, derived from | |
// Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73 | |
template <typename T, typename Policy> | |
int CF1_ik(T v, T x, T* fv, const Policy& pol) | |
{ | |
T C, D, f, a, b, delta, tiny, tolerance; | |
unsigned long k; | |
BOOST_MATH_STD_USING | |
// |x| <= |v|, CF1_ik converges rapidly | |
// |x| > |v|, CF1_ik needs O(|x|) iterations to converge | |
// modified Lentz's method, see | |
// Lentz, Applied Optics, vol 15, 668 (1976) | |
tolerance = 2 * tools::epsilon<T>(); | |
BOOST_MATH_INSTRUMENT_VARIABLE(tolerance); | |
tiny = sqrt(tools::min_value<T>()); | |
BOOST_MATH_INSTRUMENT_VARIABLE(tiny); | |
C = f = tiny; // b0 = 0, replace with tiny | |
D = 0; | |
for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++) | |
{ | |
a = 1; | |
b = 2 * (v + k) / x; | |
C = b + a / C; | |
D = b + a * D; | |
if (C == 0) { C = tiny; } | |
if (D == 0) { D = tiny; } | |
D = 1 / D; | |
delta = C * D; | |
f *= delta; | |
BOOST_MATH_INSTRUMENT_VARIABLE(delta-1); | |
if (abs(delta - 1) <= tolerance) | |
{ | |
break; | |
} | |
} | |
BOOST_MATH_INSTRUMENT_VARIABLE(k); | |
policies::check_series_iterations("boost::math::bessel_ik<%1%>(%1%,%1%) in CF1_ik", k, pol); | |
*fv = f; | |
return 0; | |
} | |
// Calculate K(v, x) and K(v+1, x) by evaluating continued fraction | |
// z1 / z0 = U(v+1.5, 2v+1, 2x) / U(v+0.5, 2v+1, 2x), see | |
// Thompson and Barnett, Computer Physics Communications, vol 47, 245 (1987) | |
template <typename T, typename Policy> | |
int CF2_ik(T v, T x, T* Kv, T* Kv1, const Policy& pol) | |
{ | |
BOOST_MATH_STD_USING | |
using namespace boost::math::constants; | |
T S, C, Q, D, f, a, b, q, delta, tolerance, current, prev; | |
unsigned long k; | |
// |x| >= |v|, CF2_ik converges rapidly | |
// |x| -> 0, CF2_ik fails to converge | |
BOOST_ASSERT(abs(x) > 1); | |
// Steed's algorithm, see Thompson and Barnett, | |
// Journal of Computational Physics, vol 64, 490 (1986) | |
tolerance = tools::epsilon<T>(); | |
a = v * v - 0.25f; | |
b = 2 * (x + 1); // b1 | |
D = 1 / b; // D1 = 1 / b1 | |
f = delta = D; // f1 = delta1 = D1, coincidence | |
prev = 0; // q0 | |
current = 1; // q1 | |
Q = C = -a; // Q1 = C1 because q1 = 1 | |
S = 1 + Q * delta; // S1 | |
BOOST_MATH_INSTRUMENT_VARIABLE(tolerance); | |
BOOST_MATH_INSTRUMENT_VARIABLE(a); | |
BOOST_MATH_INSTRUMENT_VARIABLE(b); | |
BOOST_MATH_INSTRUMENT_VARIABLE(D); | |
BOOST_MATH_INSTRUMENT_VARIABLE(f); | |
for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++) // starting from 2 | |
{ | |
// continued fraction f = z1 / z0 | |
a -= 2 * (k - 1); | |
b += 2; | |
D = 1 / (b + a * D); | |
delta *= b * D - 1; | |
f += delta; | |
// series summation S = 1 + \sum_{n=1}^{\infty} C_n * z_n / z_0 | |
q = (prev - (b - 2) * current) / a; | |
prev = current; | |
current = q; // forward recurrence for q | |
C *= -a / k; | |
Q += C * q; | |
S += Q * delta; | |
// S converges slower than f | |
BOOST_MATH_INSTRUMENT_VARIABLE(Q * delta); | |
BOOST_MATH_INSTRUMENT_VARIABLE(abs(S) * tolerance); | |
if (abs(Q * delta) < abs(S) * tolerance) | |
{ | |
break; | |
} | |
} | |
policies::check_series_iterations("boost::math::bessel_ik<%1%>(%1%,%1%) in CF2_ik", k, pol); | |
*Kv = sqrt(pi<T>() / (2 * x)) * exp(-x) / S; | |
*Kv1 = *Kv * (0.5f + v + x + (v * v - 0.25f) * f) / x; | |
BOOST_MATH_INSTRUMENT_VARIABLE(*Kv); | |
BOOST_MATH_INSTRUMENT_VARIABLE(*Kv1); | |
return 0; | |
} | |
enum{ | |
need_i = 1, | |
need_k = 2 | |
}; | |
// Compute I(v, x) and K(v, x) simultaneously by Temme's method, see | |
// Temme, Journal of Computational Physics, vol 19, 324 (1975) | |
template <typename T, typename Policy> | |
int bessel_ik(T v, T x, T* I, T* K, int kind, const Policy& pol) | |
{ | |
// Kv1 = K_(v+1), fv = I_(v+1) / I_v | |
// Ku1 = K_(u+1), fu = I_(u+1) / I_u | |
T u, Iv, Kv, Kv1, Ku, Ku1, fv; | |
T W, current, prev, next; | |
bool reflect = false; | |
unsigned n, k; | |
BOOST_MATH_INSTRUMENT_VARIABLE(v); | |
BOOST_MATH_INSTRUMENT_VARIABLE(x); | |
BOOST_MATH_INSTRUMENT_VARIABLE(kind); | |
BOOST_MATH_STD_USING | |
using namespace boost::math::tools; | |
using namespace boost::math::constants; | |
static const char* function = "boost::math::bessel_ik<%1%>(%1%,%1%)"; | |
if (v < 0) | |
{ | |
reflect = true; | |
v = -v; // v is non-negative from here | |
kind |= need_k; | |
} | |
n = iround(v, pol); | |
u = v - n; // -1/2 <= u < 1/2 | |
BOOST_MATH_INSTRUMENT_VARIABLE(n); | |
BOOST_MATH_INSTRUMENT_VARIABLE(u); | |
if (x < 0) | |
{ | |
*I = *K = policies::raise_domain_error<T>(function, | |
"Got x = %1% but real argument x must be non-negative, complex number result not supported.", x, pol); | |
return 1; | |
} | |
if (x == 0) | |
{ | |
Iv = (v == 0) ? static_cast<T>(1) : static_cast<T>(0); | |
if(kind & need_k) | |
{ | |
Kv = policies::raise_overflow_error<T>(function, 0, pol); | |
} | |
else | |
{ | |
Kv = std::numeric_limits<T>::quiet_NaN(); // any value will do | |
} | |
if(reflect && (kind & need_i)) | |
{ | |
T z = (u + n % 2); | |
Iv = boost::math::sin_pi(z, pol) == 0 ? | |
Iv : | |
policies::raise_overflow_error<T>(function, 0, pol); // reflection formula | |
} | |
*I = Iv; | |
*K = Kv; | |
return 0; | |
} | |
// x is positive until reflection | |
W = 1 / x; // Wronskian | |
if (x <= 2) // x in (0, 2] | |
{ | |
temme_ik(u, x, &Ku, &Ku1, pol); // Temme series | |
} | |
else // x in (2, \infty) | |
{ | |
CF2_ik(u, x, &Ku, &Ku1, pol); // continued fraction CF2_ik | |
} | |
prev = Ku; | |
current = Ku1; | |
for (k = 1; k <= n; k++) // forward recurrence for K | |
{ | |
next = 2 * (u + k) * current / x + prev; | |
prev = current; | |
current = next; | |
} | |
Kv = prev; | |
Kv1 = current; | |
if(kind & need_i) | |
{ | |
T lim = (4 * v * v + 10) / (8 * x); | |
lim *= lim; | |
lim *= lim; | |
lim /= 24; | |
if((lim < tools::epsilon<T>() * 10) && (x > 100)) | |
{ | |
// x is huge compared to v, CF1 may be very slow | |
// to converge so use asymptotic expansion for large | |
// x case instead. Note that the asymptotic expansion | |
// isn't very accurate - so it's deliberately very hard | |
// to get here - probably we're going to overflow: | |
Iv = asymptotic_bessel_i_large_x(v, x, pol); | |
} | |
else | |
{ | |
CF1_ik(v, x, &fv, pol); // continued fraction CF1_ik | |
Iv = W / (Kv * fv + Kv1); // Wronskian relation | |
} | |
} | |
else | |
Iv = std::numeric_limits<T>::quiet_NaN(); // any value will do | |
if (reflect) | |
{ | |
T z = (u + n % 2); | |
*I = Iv + (2 / pi<T>()) * boost::math::sin_pi(z) * Kv; // reflection formula | |
*K = Kv; | |
} | |
else | |
{ | |
*I = Iv; | |
*K = Kv; | |
} | |
BOOST_MATH_INSTRUMENT_VARIABLE(*I); | |
BOOST_MATH_INSTRUMENT_VARIABLE(*K); | |
return 0; | |
} | |
}}} // namespaces | |
#endif // BOOST_MATH_BESSEL_IK_HPP | |