// 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 #include #include #include #include #include // Modified Bessel functions of the first and second kind of fractional order namespace boost { namespace math { namespace detail { template struct cyl_bessel_i_small_z { typedef T result_type; cyl_bessel_i_small_z(T v_, T z_) : k(0), v(v_), mult(z_*z_/4) { BOOST_MATH_STD_USING term = 1; } T operator()() { T result = term; ++k; term *= mult / k; term /= k + v; return result; } private: unsigned k; T v; T term; T mult; }; template inline T bessel_i_small_z_series(T v, T x, const Policy& pol) { BOOST_MATH_STD_USING T prefix; if(v < max_factorial::value) { prefix = pow(x / 2, v) / boost::math::tgamma(v + 1, pol); } else { prefix = v * log(x / 2) - boost::math::lgamma(v + 1, pol); prefix = exp(prefix); } if(prefix == 0) return prefix; cyl_bessel_i_small_z s(v, x); boost::uintmax_t max_iter = policies::get_max_series_iterations(); #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) T zero = 0; T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon(), max_iter, zero); #else T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon(), max_iter); #endif policies::check_series_iterations("boost::math::bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol); return prefix * result; } // Calculate K(v, x) and K(v+1, x) by method analogous to // Temme, Journal of Computational Physics, vol 21, 343 (1976) template 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(1) : T(boost::math::sin_pi(v) / (v * pi())); d = abs(sigma) < tools::epsilon() ? T(1) : T(sinh(sigma) / sigma); gamma1 = abs(v) < tools::epsilon() ? T(-euler()) : 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; BOOST_MATH_INSTRUMENT_VARIABLE(p); BOOST_MATH_INSTRUMENT_VARIABLE(q); BOOST_MATH_INSTRUMENT_VARIABLE(f); BOOST_MATH_INSTRUMENT_VARIABLE(sigma); BOOST_MATH_INSTRUMENT_CODE(sinh(sigma)); BOOST_MATH_INSTRUMENT_VARIABLE(gamma1); BOOST_MATH_INSTRUMENT_VARIABLE(gamma2); BOOST_MATH_INSTRUMENT_VARIABLE(c); BOOST_MATH_INSTRUMENT_VARIABLE(d); BOOST_MATH_INSTRUMENT_VARIABLE(a); // series summation tolerance = tools::epsilon(); for (k = 1; k < policies::get_max_series_iterations(); 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 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(); BOOST_MATH_INSTRUMENT_VARIABLE(tolerance); tiny = sqrt(tools::min_value()); BOOST_MATH_INSTRUMENT_VARIABLE(tiny); C = f = tiny; // b0 = 0, replace with tiny D = 0; for (k = 1; k < policies::get_max_series_iterations(); 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 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(); 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(); 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; // // Under some circumstances q can grow very small and C very // large, leading to under/overflow. This is particularly an // issue for types which have many digits precision but a narrow // exponent range. A typical example being a "double double" type. // To avoid this situation we can normalise q (and related prev/current) // and C. All other variables remain unchanged in value. A typical // test case occurs when x is close to 2, for example cyl_bessel_k(9.125, 2.125). // if(q < tools::epsilon()) { C *= q; prev /= q; current /= q; q = 1; } // S converges slower than f BOOST_MATH_INSTRUMENT_VARIABLE(Q * delta); BOOST_MATH_INSTRUMENT_VARIABLE(abs(S) * tolerance); BOOST_MATH_INSTRUMENT_VARIABLE(S); if (abs(Q * delta) < abs(S) * tolerance) { break; } } policies::check_series_iterations("boost::math::bessel_ik<%1%>(%1%,%1%) in CF2_ik", k, pol); if(x >= tools::log_max_value()) *Kv = exp(0.5f * log(pi() / (2 * x)) - x - log(S)); else *Kv = sqrt(pi() / (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 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; int org_kind = kind; 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(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(1) : static_cast(0); if(kind & need_k) { Kv = policies::raise_overflow_error(function, 0, pol); } else { Kv = std::numeric_limits::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(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 } BOOST_MATH_INSTRUMENT_VARIABLE(Ku); BOOST_MATH_INSTRUMENT_VARIABLE(Ku1); prev = Ku; current = Ku1; T scale = 1; for (k = 1; k <= n; k++) // forward recurrence for K { T fact = 2 * (u + k) / x; if((tools::max_value() - fabs(prev)) / fact < fabs(current)) { prev /= current; scale /= current; current = 1; } next = fact * current + prev; prev = current; current = next; } Kv = prev; Kv1 = current; BOOST_MATH_INSTRUMENT_VARIABLE(Kv); BOOST_MATH_INSTRUMENT_VARIABLE(Kv1); if(kind & need_i) { T lim = (4 * v * v + 10) / (8 * x); lim *= lim; lim *= lim; lim /= 24; if((lim < tools::epsilon() * 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 if((v > 0) && (x / v < 0.25)) { Iv = bessel_i_small_z_series(v, x, pol); } else { CF1_ik(v, x, &fv, pol); // continued fraction CF1_ik Iv = scale * W / (Kv * fv + Kv1); // Wronskian relation } } else Iv = std::numeric_limits::quiet_NaN(); // any value will do if (reflect) { T z = (u + n % 2); T fact = (2 / pi()) * (boost::math::sin_pi(z) * Kv); if(fact == 0) *I = Iv; else if(tools::max_value() * scale < fact) *I = (org_kind & need_i) ? T(sign(fact) * sign(scale) * policies::raise_overflow_error(function, 0, pol)) : T(0); else *I = Iv + fact / scale; // reflection formula } else { *I = Iv; } if(tools::max_value() * scale < Kv) *K = (org_kind & need_k) ? T(sign(Kv) * sign(scale) * policies::raise_overflow_error(function, 0, pol)) : T(0); else *K = Kv / scale; BOOST_MATH_INSTRUMENT_VARIABLE(*I); BOOST_MATH_INSTRUMENT_VARIABLE(*K); return 0; } }}} // namespaces #endif // BOOST_MATH_BESSEL_IK_HPP