448 lines · plain
1// Copyright (c) 2006 Xiaogang Zhang2// Copyright (c) 2024 Matt Borland3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_BESSEL_IK_HPP8#define BOOST_MATH_BESSEL_IK_HPP9 10#ifdef _MSC_VER11#pragma once12#endif13 14#include <boost/math/tools/config.hpp>15#include <boost/math/tools/cstdint.hpp>16#include <boost/math/tools/numeric_limits.hpp>17#include <boost/math/tools/type_traits.hpp>18#include <boost/math/tools/series.hpp>19#include <boost/math/special_functions/sign.hpp>20#include <boost/math/special_functions/round.hpp>21#include <boost/math/special_functions/gamma.hpp>22#include <boost/math/special_functions/sin_pi.hpp>23#include <boost/math/constants/constants.hpp>24#include <boost/math/policies/error_handling.hpp>25 26// Modified Bessel functions of the first and second kind of fractional order27 28namespace boost { namespace math {29 30namespace detail {31 32template <class T, class Policy>33struct cyl_bessel_i_small_z34{35 typedef T result_type;36 37 BOOST_MATH_GPU_ENABLED cyl_bessel_i_small_z(T v_, T z_) : k(0), v(v_), mult(z_*z_/4)38 {39 BOOST_MATH_STD_USING40 term = 1;41 }42 43 BOOST_MATH_GPU_ENABLED T operator()()44 {45 T result = term;46 ++k;47 term *= mult / k;48 term /= k + v;49 return result;50 }51private:52 unsigned k;53 T v;54 T term;55 T mult;56};57 58template <class T, class Policy>59BOOST_MATH_GPU_ENABLED inline T bessel_i_small_z_series(T v, T x, const Policy& pol)60{61 BOOST_MATH_STD_USING62 T prefix;63 if(v < max_factorial<T>::value)64 {65 prefix = pow(x / 2, v) / boost::math::tgamma(v + 1, pol);66 }67 else68 {69 prefix = v * log(x / 2) - boost::math::lgamma(v + 1, pol);70 prefix = exp(prefix);71 }72 if(prefix == 0)73 return prefix;74 75 cyl_bessel_i_small_z<T, Policy> s(v, x);76 boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();77 78 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);79 80 policies::check_series_iterations<T>("boost::math::bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);81 return prefix * result;82}83 84// Calculate K(v, x) and K(v+1, x) by method analogous to85// Temme, Journal of Computational Physics, vol 21, 343 (1976)86template <typename T, typename Policy>87BOOST_MATH_GPU_ENABLED int temme_ik(T v, T x, T* result_K, T* K1, const Policy& pol)88{89 T f, h, p, q, coef, sum, sum1, tolerance;90 T a, b, c, d, sigma, gamma1, gamma2;91 unsigned long k;92 93 BOOST_MATH_STD_USING94 using namespace boost::math::tools;95 using namespace boost::math::constants;96 97 98 // |x| <= 2, Temme series converge rapidly99 // |x| > 2, the larger the |x|, the slower the convergence100 BOOST_MATH_ASSERT(abs(x) <= 2);101 BOOST_MATH_ASSERT(abs(v) <= 0.5f);102 103 T gp = boost::math::tgamma1pm1(v, pol);104 T gm = boost::math::tgamma1pm1(-v, pol);105 106 a = log(x / 2);107 b = exp(v * a);108 sigma = -a * v;109 c = abs(v) < tools::epsilon<T>() ?110 T(1) : T(boost::math::sin_pi(v, pol) / (v * pi<T>()));111 d = abs(sigma) < tools::epsilon<T>() ?112 T(1) : T(sinh(sigma) / sigma);113 gamma1 = abs(v) < tools::epsilon<T>() ?114 T(-euler<T>()) : T((0.5f / v) * (gp - gm) * c);115 gamma2 = (2 + gp + gm) * c / 2;116 117 // initial values118 p = (gp + 1) / (2 * b);119 q = (1 + gm) * b / 2;120 f = (cosh(sigma) * gamma1 + d * (-a) * gamma2) / c;121 h = p;122 coef = 1;123 sum = coef * f;124 sum1 = coef * h;125 126 BOOST_MATH_INSTRUMENT_VARIABLE(p);127 BOOST_MATH_INSTRUMENT_VARIABLE(q);128 BOOST_MATH_INSTRUMENT_VARIABLE(f);129 BOOST_MATH_INSTRUMENT_VARIABLE(sigma);130 BOOST_MATH_INSTRUMENT_CODE(sinh(sigma));131 BOOST_MATH_INSTRUMENT_VARIABLE(gamma1);132 BOOST_MATH_INSTRUMENT_VARIABLE(gamma2);133 BOOST_MATH_INSTRUMENT_VARIABLE(c);134 BOOST_MATH_INSTRUMENT_VARIABLE(d);135 BOOST_MATH_INSTRUMENT_VARIABLE(a);136 137 // series summation138 tolerance = tools::epsilon<T>();139 for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)140 {141 f = (k * f + p + q) / (k*k - v*v);142 p /= k - v;143 q /= k + v;144 h = p - k * f;145 coef *= x * x / (4 * k);146 sum += coef * f;147 sum1 += coef * h;148 if (abs(coef * f) < abs(sum) * tolerance)149 {150 break;151 }152 }153 policies::check_series_iterations<T>("boost::math::bessel_ik<%1%>(%1%,%1%) in temme_ik", k, pol);154 155 *result_K = sum;156 *K1 = 2 * sum1 / x;157 158 return 0;159}160 161// Evaluate continued fraction fv = I_(v+1) / I_v, derived from162// Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73163template <typename T, typename Policy>164BOOST_MATH_GPU_ENABLED int CF1_ik(T v, T x, T* fv, const Policy& pol)165{166 T C, D, f, a, b, delta, tiny, tolerance;167 unsigned long k;168 169 BOOST_MATH_STD_USING170 171 // |x| <= |v|, CF1_ik converges rapidly172 // |x| > |v|, CF1_ik needs O(|x|) iterations to converge173 174 // modified Lentz's method, see175 // Lentz, Applied Optics, vol 15, 668 (1976)176 tolerance = 2 * tools::epsilon<T>();177 BOOST_MATH_INSTRUMENT_VARIABLE(tolerance);178 tiny = sqrt(tools::min_value<T>());179 BOOST_MATH_INSTRUMENT_VARIABLE(tiny);180 C = f = tiny; // b0 = 0, replace with tiny181 D = 0;182 for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)183 {184 a = 1;185 b = 2 * (v + k) / x;186 C = b + a / C;187 D = b + a * D;188 if (C == 0) { C = tiny; }189 if (D == 0) { D = tiny; }190 D = 1 / D;191 delta = C * D;192 f *= delta;193 BOOST_MATH_INSTRUMENT_VARIABLE(delta-1);194 if (abs(delta - 1) <= tolerance)195 {196 break;197 }198 }199 BOOST_MATH_INSTRUMENT_VARIABLE(k);200 policies::check_series_iterations<T>("boost::math::bessel_ik<%1%>(%1%,%1%) in CF1_ik", k, pol);201 202 *fv = f;203 204 return 0;205}206 207// Calculate K(v, x) and K(v+1, x) by evaluating continued fraction208// z1 / z0 = U(v+1.5, 2v+1, 2x) / U(v+0.5, 2v+1, 2x), see209// Thompson and Barnett, Computer Physics Communications, vol 47, 245 (1987)210template <typename T, typename Policy>211BOOST_MATH_GPU_ENABLED int CF2_ik(T v, T x, T* Kv, T* Kv1, const Policy& pol)212{213 BOOST_MATH_STD_USING214 using namespace boost::math::constants;215 216 T S, C, Q, D, f, a, b, q, delta, tolerance, current, prev;217 unsigned long k;218 219 // |x| >= |v|, CF2_ik converges rapidly220 // |x| -> 0, CF2_ik fails to converge221 222 BOOST_MATH_ASSERT(abs(x) > 1);223 224 // Steed's algorithm, see Thompson and Barnett,225 // Journal of Computational Physics, vol 64, 490 (1986)226 tolerance = tools::epsilon<T>();227 a = v * v - 0.25f;228 b = 2 * (x + 1); // b1229 D = 1 / b; // D1 = 1 / b1230 f = delta = D; // f1 = delta1 = D1, coincidence231 prev = 0; // q0232 current = 1; // q1233 Q = C = -a; // Q1 = C1 because q1 = 1234 S = 1 + Q * delta; // S1235 BOOST_MATH_INSTRUMENT_VARIABLE(tolerance);236 BOOST_MATH_INSTRUMENT_VARIABLE(a);237 BOOST_MATH_INSTRUMENT_VARIABLE(b);238 BOOST_MATH_INSTRUMENT_VARIABLE(D);239 BOOST_MATH_INSTRUMENT_VARIABLE(f);240 241 for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++) // starting from 2242 {243 // continued fraction f = z1 / z0244 a -= 2 * (k - 1);245 b += 2;246 D = 1 / (b + a * D);247 delta *= b * D - 1;248 f += delta;249 250 // series summation S = 1 + \sum_{n=1}^{\infty} C_n * z_n / z_0251 q = (prev - (b - 2) * current) / a;252 prev = current;253 current = q; // forward recurrence for q254 C *= -a / k;255 Q += C * q;256 S += Q * delta;257 //258 // Under some circumstances q can grow very small and C very259 // large, leading to under/overflow. This is particularly an260 // issue for types which have many digits precision but a narrow261 // exponent range. A typical example being a "double double" type.262 // To avoid this situation we can normalise q (and related prev/current)263 // and C. All other variables remain unchanged in value. A typical264 // test case occurs when x is close to 2, for example cyl_bessel_k(9.125, 2.125).265 //266 if(q < tools::epsilon<T>())267 {268 C *= q;269 prev /= q;270 current /= q;271 q = 1;272 }273 274 // S converges slower than f275 BOOST_MATH_INSTRUMENT_VARIABLE(Q * delta);276 BOOST_MATH_INSTRUMENT_VARIABLE(abs(S) * tolerance);277 BOOST_MATH_INSTRUMENT_VARIABLE(S);278 if (abs(Q * delta) < abs(S) * tolerance)279 {280 break;281 }282 }283 policies::check_series_iterations<T>("boost::math::bessel_ik<%1%>(%1%,%1%) in CF2_ik", k, pol);284 285 if(-x < tools::log_min_value<T>())286 *Kv = exp(0.5f * log(pi<T>() / (2 * x)) - x - log(S));287 else288 *Kv = sqrt(pi<T>() / (2 * x)) * exp(-x) / S;289 *Kv1 = *Kv * (0.5f + v + x + (v * v - 0.25f) * f) / x;290 BOOST_MATH_INSTRUMENT_VARIABLE(*Kv);291 BOOST_MATH_INSTRUMENT_VARIABLE(*Kv1);292 293 return 0;294}295 296enum{297 need_i = 1,298 need_k = 2299};300 301// Compute I(v, x) and K(v, x) simultaneously by Temme's method, see302// Temme, Journal of Computational Physics, vol 19, 324 (1975)303template <typename T, typename Policy>304BOOST_MATH_GPU_ENABLED int bessel_ik(T v, T x, T* result_I, T* result_K, int kind, const Policy& pol)305{306 // Kv1 = K_(v+1), fv = I_(v+1) / I_v307 // Ku1 = K_(u+1), fu = I_(u+1) / I_u308 T u, Iv, Kv, Kv1, Ku, Ku1, fv;309 T W, current, prev, next;310 bool reflect = false;311 unsigned n, k;312 int org_kind = kind;313 BOOST_MATH_INSTRUMENT_VARIABLE(v);314 BOOST_MATH_INSTRUMENT_VARIABLE(x);315 BOOST_MATH_INSTRUMENT_VARIABLE(kind);316 317 BOOST_MATH_STD_USING318 using namespace boost::math::tools;319 using namespace boost::math::constants;320 321 constexpr auto function = "boost::math::bessel_ik<%1%>(%1%,%1%)";322 323 if (v < 0)324 {325 reflect = true;326 v = -v; // v is non-negative from here327 kind |= need_k;328 }329 330 T scale = 1;331 T scale_sign = 1;332 333 n = iround(v, pol);334 u = v - n; // -1/2 <= u < 1/2335 BOOST_MATH_INSTRUMENT_VARIABLE(n);336 BOOST_MATH_INSTRUMENT_VARIABLE(u);337 338 if (((kind & need_i) == 0) && (fabs(4 * v * v - 25) / (8 * x) < tools::forth_root_epsilon<T>()))339 {340 // A&S 9.7.2341 Iv = boost::math::numeric_limits<T>::quiet_NaN(); // any value will do342 T mu = 4 * v * v;343 T eight_z = 8 * x;344 Kv = 1 + (mu - 1) / eight_z + (mu - 1) * (mu - 9) / (2 * eight_z * eight_z) + (mu - 1) * (mu - 9) * (mu - 25) / (6 * eight_z * eight_z * eight_z);345 Kv *= exp(-x) * constants::root_pi<T>() / sqrt(2 * x);346 }347 else348 {349 BOOST_MATH_ASSERT(x > 0); // Error handling for x <= 0 handled in cyl_bessel_i and cyl_bessel_k350 351 // x is positive until reflection352 W = 1 / x; // Wronskian353 if (x <= 2) // x in (0, 2]354 {355 temme_ik(u, x, &Ku, &Ku1, pol); // Temme series356 }357 else // x in (2, \infty)358 {359 CF2_ik(u, x, &Ku, &Ku1, pol); // continued fraction CF2_ik360 }361 BOOST_MATH_INSTRUMENT_VARIABLE(Ku);362 BOOST_MATH_INSTRUMENT_VARIABLE(Ku1);363 prev = Ku;364 current = Ku1;365 for (k = 1; k <= n; k++) // forward recurrence for K366 {367 T fact = 2 * (u + k) / x;368 // Check for overflow: if (max - |prev|) / fact > max, then overflow369 // (max - |prev|) / fact > max370 // max * (1 - fact) > |prev|371 // if fact < 1: safe to compute overflow check372 // if fact >= 1: won't overflow373 const bool will_overflow = (fact < 1)374 ? tools::max_value<T>() * (1 - fact) > fabs(prev)375 : false;376 if (!will_overflow && ((tools::max_value<T>() - fabs(prev)) / fact < fabs(current)))377 {378 prev /= current;379 scale /= current;380 scale_sign *= ((boost::math::signbit)(current) ? -1 : 1);381 current = 1;382 }383 next = fact * current + prev;384 prev = current;385 current = next;386 }387 Kv = prev;388 Kv1 = current;389 BOOST_MATH_INSTRUMENT_VARIABLE(Kv);390 BOOST_MATH_INSTRUMENT_VARIABLE(Kv1);391 if (kind & need_i)392 {393 T lim = (4 * v * v + 10) / (8 * x);394 lim *= lim;395 lim *= lim;396 lim /= 24;397 if ((lim < tools::epsilon<T>() * 10) && (x > 100))398 {399 // x is huge compared to v, CF1 may be very slow400 // to converge so use asymptotic expansion for large401 // x case instead. Note that the asymptotic expansion402 // isn't very accurate - so it's deliberately very hard403 // to get here - probably we're going to overflow:404 Iv = asymptotic_bessel_i_large_x(v, x, pol);405 }406 else if ((v > 0) && (x / v < 0.25))407 {408 Iv = bessel_i_small_z_series(v, x, pol);409 }410 else411 {412 CF1_ik(v, x, &fv, pol); // continued fraction CF1_ik413 Iv = scale * W / (Kv * fv + Kv1); // Wronskian relation414 }415 }416 else417 Iv = boost::math::numeric_limits<T>::quiet_NaN(); // any value will do418 }419 if (reflect && (kind & need_i))420 {421 BOOST_MATH_ASSERT(fabs(v - n - u) < tools::forth_root_epsilon<T>());422 T z = (u + n % 2);423 T fact = (2 / pi<T>()) * (boost::math::sin_pi(z, pol) * Kv);424 if(fact == 0)425 *result_I = Iv;426 else if(tools::max_value<T>() * scale < fact)427 *result_I = (org_kind & need_i) ? T(sign(fact) * scale_sign * policies::raise_overflow_error<T>(function, nullptr, pol)) : T(0);428 else429 *result_I = Iv + fact / scale; // reflection formula430 }431 else432 {433 *result_I = Iv;434 }435 if(tools::max_value<T>() * scale < Kv)436 *result_K = (org_kind & need_k) ? T(sign(Kv) * scale_sign * policies::raise_overflow_error<T>(function, nullptr, pol)) : T(0);437 else438 *result_K = Kv / scale;439 BOOST_MATH_INSTRUMENT_VARIABLE(*result_I);440 BOOST_MATH_INSTRUMENT_VARIABLE(*result_K);441 return 0;442}443 444}}} // namespaces445 446#endif // BOOST_MATH_BESSEL_IK_HPP447 448