special_functions/detail/bessel_ik.hpp

0001 //  Copyright (c) 2006 Xiaogang Zhang
0002 //  Copyright (c) 2024 Matt Borland
0003 //  Use, modification and distribution are subject to the
0004 //  Boost Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_MATH_BESSEL_IK_HPP
0008 #define BOOST_MATH_BESSEL_IK_HPP
0009
0010 #ifdef _MSC_VER
0011 #pragma once
0012 #endif
0013
0014 #include <boost/math/tools/config.hpp>
0015 #include <boost/math/tools/cstdint.hpp>
0016 #include <boost/math/tools/numeric_limits.hpp>
0017 #include <boost/math/tools/type_traits.hpp>
0018 #include <boost/math/tools/series.hpp>
0019 #include <boost/math/special_functions/sign.hpp>
0020 #include <boost/math/special_functions/round.hpp>
0021 #include <boost/math/special_functions/gamma.hpp>
0022 #include <boost/math/special_functions/sin_pi.hpp>
0023 #include <boost/math/constants/constants.hpp>
0024 #include <boost/math/policies/error_handling.hpp>
0025
0026 // Modified Bessel functions of the first and second kind of fractional order
0027
0028 namespace boost { namespace math {
0029
0030 namespace detail {
0031
0032 template <class T, class Policy>
0033 struct cyl_bessel_i_small_z
0034 {
0035    typedef T result_type;
0036
0037    BOOST_MATH_GPU_ENABLED cyl_bessel_i_small_z(T v_, T z_) : k(0), v(v_), mult(z_*z_/4)
0038    {
0039       BOOST_MATH_STD_USING
0040       term = 1;
0041    }
0042
0043    BOOST_MATH_GPU_ENABLED T operator()()
0044    {
0045       T result = term;
0046       ++k;
0047       term *= mult / k;
0048       term /= k + v;
0049       return result;
0050    }
0051 private:
0052    unsigned k;
0053    T v;
0054    T term;
0055    T mult;
0056 };
0057
0058 template <class T, class Policy>
0059 BOOST_MATH_GPU_ENABLED inline T bessel_i_small_z_series(T v, T x, const Policy& pol)
0060 {
0061    BOOST_MATH_STD_USING
0062    T prefix;
0063    if(v < max_factorial<T>::value)
0064    {
0065       prefix = pow(x / 2, v) / boost::math::tgamma(v + 1, pol);
0066    }
0067    else
0068    {
0069       prefix = v * log(x / 2) - boost::math::lgamma(v + 1, pol);
0070       prefix = exp(prefix);
0071    }
0072    if(prefix == 0)
0073       return prefix;
0074
0075    cyl_bessel_i_small_z<T, Policy> s(v, x);
0076    boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0077
0078    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0079
0080    policies::check_series_iterations<T>("boost::math::bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
0081    return prefix * result;
0082 }
0083
0084 // Calculate K(v, x) and K(v+1, x) by method analogous to
0085 // Temme, Journal of Computational Physics, vol 21, 343 (1976)
0086 template <typename T, typename Policy>
0087 BOOST_MATH_GPU_ENABLED int temme_ik(T v, T x, T* result_K, T* K1, const Policy& pol)
0088 {
0089     T f, h, p, q, coef, sum, sum1, tolerance;
0090     T a, b, c, d, sigma, gamma1, gamma2;
0091     unsigned long k;
0092
0093     BOOST_MATH_STD_USING
0094     using namespace boost::math::tools;
0095     using namespace boost::math::constants;
0096
0097
0098     // |x| <= 2, Temme series converge rapidly
0099     // |x| > 2, the larger the |x|, the slower the convergence
0100     BOOST_MATH_ASSERT(abs(x) <= 2);
0101     BOOST_MATH_ASSERT(abs(v) <= 0.5f);
0102
0103     T gp = boost::math::tgamma1pm1(v, pol);
0104     T gm = boost::math::tgamma1pm1(-v, pol);
0105
0106     a = log(x / 2);
0107     b = exp(v * a);
0108     sigma = -a * v;
0109     c = abs(v) < tools::epsilon<T>() ?
0110        T(1) : T(boost::math::sin_pi(v, pol) / (v * pi<T>()));
0111     d = abs(sigma) < tools::epsilon<T>() ?
0112         T(1) : T(sinh(sigma) / sigma);
0113     gamma1 = abs(v) < tools::epsilon<T>() ?
0114         T(-euler<T>()) : T((0.5f / v) * (gp - gm) * c);
0115     gamma2 = (2 + gp + gm) * c / 2;
0116
0117     // initial values
0118     p = (gp + 1) / (2 * b);
0119     q = (1 + gm) * b / 2;
0120     f = (cosh(sigma) * gamma1 + d * (-a) * gamma2) / c;
0121     h = p;
0122     coef = 1;
0123     sum = coef * f;
0124     sum1 = coef * h;
0125
0126     BOOST_MATH_INSTRUMENT_VARIABLE(p);
0127     BOOST_MATH_INSTRUMENT_VARIABLE(q);
0128     BOOST_MATH_INSTRUMENT_VARIABLE(f);
0129     BOOST_MATH_INSTRUMENT_VARIABLE(sigma);
0130     BOOST_MATH_INSTRUMENT_CODE(sinh(sigma));
0131     BOOST_MATH_INSTRUMENT_VARIABLE(gamma1);
0132     BOOST_MATH_INSTRUMENT_VARIABLE(gamma2);
0133     BOOST_MATH_INSTRUMENT_VARIABLE(c);
0134     BOOST_MATH_INSTRUMENT_VARIABLE(d);
0135     BOOST_MATH_INSTRUMENT_VARIABLE(a);
0136
0137     // series summation
0138     tolerance = tools::epsilon<T>();
0139     for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)
0140     {
0141         f = (k * f + p + q) / (k*k - v*v);
0142         p /= k - v;
0143         q /= k + v;
0144         h = p - k * f;
0145         coef *= x * x / (4 * k);
0146         sum += coef * f;
0147         sum1 += coef * h;
0148         if (abs(coef * f) < abs(sum) * tolerance)
0149         {
0150            break;
0151         }
0152     }
0153     policies::check_series_iterations<T>("boost::math::bessel_ik<%1%>(%1%,%1%) in temme_ik", k, pol);
0154
0155     *result_K = sum;
0156     *K1 = 2 * sum1 / x;
0157
0158     return 0;
0159 }
0160
0161 // Evaluate continued fraction fv = I_(v+1) / I_v, derived from
0162 // Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73
0163 template <typename T, typename Policy>
0164 BOOST_MATH_GPU_ENABLED int CF1_ik(T v, T x, T* fv, const Policy& pol)
0165 {
0166     T C, D, f, a, b, delta, tiny, tolerance;
0167     unsigned long k;
0168
0169     BOOST_MATH_STD_USING
0170
0171     // |x| <= |v|, CF1_ik converges rapidly
0172     // |x| > |v|, CF1_ik needs O(|x|) iterations to converge
0173
0174     // modified Lentz's method, see
0175     // Lentz, Applied Optics, vol 15, 668 (1976)
0176     tolerance = 2 * tools::epsilon<T>();
0177     BOOST_MATH_INSTRUMENT_VARIABLE(tolerance);
0178     tiny = sqrt(tools::min_value<T>());
0179     BOOST_MATH_INSTRUMENT_VARIABLE(tiny);
0180     C = f = tiny;                           // b0 = 0, replace with tiny
0181     D = 0;
0182     for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)
0183     {
0184         a = 1;
0185         b = 2 * (v + k) / x;
0186         C = b + a / C;
0187         D = b + a * D;
0188         if (C == 0) { C = tiny; }
0189         if (D == 0) { D = tiny; }
0190         D = 1 / D;
0191         delta = C * D;
0192         f *= delta;
0193         BOOST_MATH_INSTRUMENT_VARIABLE(delta-1);
0194         if (abs(delta - 1) <= tolerance)
0195         {
0196            break;
0197         }
0198     }
0199     BOOST_MATH_INSTRUMENT_VARIABLE(k);
0200     policies::check_series_iterations<T>("boost::math::bessel_ik<%1%>(%1%,%1%) in CF1_ik", k, pol);
0201
0202     *fv = f;
0203
0204     return 0;
0205 }
0206
0207 // Calculate K(v, x) and K(v+1, x) by evaluating continued fraction
0208 // z1 / z0 = U(v+1.5, 2v+1, 2x) / U(v+0.5, 2v+1, 2x), see
0209 // Thompson and Barnett, Computer Physics Communications, vol 47, 245 (1987)
0210 template <typename T, typename Policy>
0211 BOOST_MATH_GPU_ENABLED int CF2_ik(T v, T x, T* Kv, T* Kv1, const Policy& pol)
0212 {
0213     BOOST_MATH_STD_USING
0214     using namespace boost::math::constants;
0215
0216     T S, C, Q, D, f, a, b, q, delta, tolerance, current, prev;
0217     unsigned long k;
0218
0219     // |x| >= |v|, CF2_ik converges rapidly
0220     // |x| -> 0, CF2_ik fails to converge
0221
0222     BOOST_MATH_ASSERT(abs(x) > 1);
0223
0224     // Steed's algorithm, see Thompson and Barnett,
0225     // Journal of Computational Physics, vol 64, 490 (1986)
0226     tolerance = tools::epsilon<T>();
0227     a = v * v - 0.25f;
0228     b = 2 * (x + 1);                              // b1
0229     D = 1 / b;                                    // D1 = 1 / b1
0230     f = delta = D;                                // f1 = delta1 = D1, coincidence
0231     prev = 0;                                     // q0
0232     current = 1;                                  // q1
0233     Q = C = -a;                                   // Q1 = C1 because q1 = 1
0234     S = 1 + Q * delta;                            // S1
0235     BOOST_MATH_INSTRUMENT_VARIABLE(tolerance);
0236     BOOST_MATH_INSTRUMENT_VARIABLE(a);
0237     BOOST_MATH_INSTRUMENT_VARIABLE(b);
0238     BOOST_MATH_INSTRUMENT_VARIABLE(D);
0239     BOOST_MATH_INSTRUMENT_VARIABLE(f);
0240
0241     for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++)     // starting from 2
0242     {
0243         // continued fraction f = z1 / z0
0244         a -= 2 * (k - 1);
0245         b += 2;
0246         D = 1 / (b + a * D);
0247         delta *= b * D - 1;
0248         f += delta;
0249
0250         // series summation S = 1 + \sum_{n=1}^{\infty} C_n * z_n / z_0
0251         q = (prev - (b - 2) * current) / a;
0252         prev = current;
0253         current = q;                        // forward recurrence for q
0254         C *= -a / k;
0255         Q += C * q;
0256         S += Q * delta;
0257         //
0258         // Under some circumstances q can grow very small and C very
0259         // large, leading to under/overflow.  This is particularly an
0260         // issue for types which have many digits precision but a narrow
0261         // exponent range.  A typical example being a "double double" type.
0262         // To avoid this situation we can normalise q (and related prev/current)
0263         // and C.  All other variables remain unchanged in value.  A typical
0264         // test case occurs when x is close to 2, for example cyl_bessel_k(9.125, 2.125).
0265         //
0266         if(q < tools::epsilon<T>())
0267         {
0268            C *= q;
0269            prev /= q;
0270            current /= q;
0271            q = 1;
0272         }
0273
0274         // S converges slower than f
0275         BOOST_MATH_INSTRUMENT_VARIABLE(Q * delta);
0276         BOOST_MATH_INSTRUMENT_VARIABLE(abs(S) * tolerance);
0277         BOOST_MATH_INSTRUMENT_VARIABLE(S);
0278         if (abs(Q * delta) < abs(S) * tolerance)
0279         {
0280            break;
0281         }
0282     }
0283     policies::check_series_iterations<T>("boost::math::bessel_ik<%1%>(%1%,%1%) in CF2_ik", k, pol);
0284
0285     if(x >= tools::log_max_value<T>())
0286        *Kv = exp(0.5f * log(pi<T>() / (2 * x)) - x - log(S));
0287     else
0288       *Kv = sqrt(pi<T>() / (2 * x)) * exp(-x) / S;
0289     *Kv1 = *Kv * (0.5f + v + x + (v * v - 0.25f) * f) / x;
0290     BOOST_MATH_INSTRUMENT_VARIABLE(*Kv);
0291     BOOST_MATH_INSTRUMENT_VARIABLE(*Kv1);
0292
0293     return 0;
0294 }
0295
0296 enum{
0297    need_i = 1,
0298    need_k = 2
0299 };
0300
0301 // Compute I(v, x) and K(v, x) simultaneously by Temme's method, see
0302 // Temme, Journal of Computational Physics, vol 19, 324 (1975)
0303 template <typename T, typename Policy>
0304 BOOST_MATH_GPU_ENABLED int bessel_ik(T v, T x, T* result_I, T* result_K, int kind, const Policy& pol)
0305 {
0306     // Kv1 = K_(v+1), fv = I_(v+1) / I_v
0307     // Ku1 = K_(u+1), fu = I_(u+1) / I_u
0308     T u, Iv, Kv, Kv1, Ku, Ku1, fv;
0309     T W, current, prev, next;
0310     bool reflect = false;
0311     unsigned n, k;
0312     int org_kind = kind;
0313     BOOST_MATH_INSTRUMENT_VARIABLE(v);
0314     BOOST_MATH_INSTRUMENT_VARIABLE(x);
0315     BOOST_MATH_INSTRUMENT_VARIABLE(kind);
0316
0317     BOOST_MATH_STD_USING
0318     using namespace boost::math::tools;
0319     using namespace boost::math::constants;
0320
0321     constexpr auto function = "boost::math::bessel_ik<%1%>(%1%,%1%)";
0322
0323     if (v < 0)
0324     {
0325         reflect = true;
0326         v = -v;                             // v is non-negative from here
0327         kind |= need_k;
0328     }
0329
0330     T scale = 1;
0331     T scale_sign = 1;
0332
0333     if (((kind & need_i) == 0) && (fabs(4 * v * v - 25) / (8 * x) < tools::forth_root_epsilon<T>()))
0334     {
0335        // A&S 9.7.2
0336        Iv = boost::math::numeric_limits<T>::quiet_NaN(); // any value will do
0337        T mu = 4 * v * v;
0338        T eight_z = 8 * x;
0339        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);
0340        Kv *= exp(-x) * constants::root_pi<T>() / sqrt(2 * x);
0341     }
0342     else
0343     {
0344        n = iround(v, pol);
0345        u = v - n;                              // -1/2 <= u < 1/2
0346        BOOST_MATH_INSTRUMENT_VARIABLE(n);
0347        BOOST_MATH_INSTRUMENT_VARIABLE(u);
0348
0349        BOOST_MATH_ASSERT(x > 0); // Error handling for x <= 0 handled in cyl_bessel_i and cyl_bessel_k
0350
0351        // x is positive until reflection
0352        W = 1 / x;                                 // Wronskian
0353        if (x <= 2)                                // x in (0, 2]
0354        {
0355           temme_ik(u, x, &Ku, &Ku1, pol);             // Temme series
0356        }
0357        else                                       // x in (2, \infty)
0358        {
0359           CF2_ik(u, x, &Ku, &Ku1, pol);               // continued fraction CF2_ik
0360        }
0361        BOOST_MATH_INSTRUMENT_VARIABLE(Ku);
0362        BOOST_MATH_INSTRUMENT_VARIABLE(Ku1);
0363        prev = Ku;
0364        current = Ku1;
0365        for (k = 1; k <= n; k++)                   // forward recurrence for K
0366        {
0367           T fact = 2 * (u + k) / x;
0368           // Check for overflow: if (max - |prev|) / fact > max, then overflow
0369           // (max - |prev|) / fact > max
0370           // max * (1 - fact) > |prev|
0371           // if fact < 1: safe to compute overflow check
0372           // if fact >= 1:  won't overflow
0373           const bool will_overflow = (fact < 1)
0374              ? tools::max_value<T>() * (1 - fact) > fabs(prev)
0375              : false;
0376           if (!will_overflow && ((tools::max_value<T>() - fabs(prev)) / fact < fabs(current)))
0377           {
0378              prev /= current;
0379              scale /= current;
0380              scale_sign *= ((boost::math::signbit)(current) ? -1 : 1);
0381              current = 1;
0382           }
0383           next = fact * current + prev;
0384           prev = current;
0385           current = next;
0386        }
0387        Kv = prev;
0388        Kv1 = current;
0389        BOOST_MATH_INSTRUMENT_VARIABLE(Kv);
0390        BOOST_MATH_INSTRUMENT_VARIABLE(Kv1);
0391        if (kind & need_i)
0392        {
0393           T lim = (4 * v * v + 10) / (8 * x);
0394           lim *= lim;
0395           lim *= lim;
0396           lim /= 24;
0397           if ((lim < tools::epsilon<T>() * 10) && (x > 100))
0398           {
0399              // x is huge compared to v, CF1 may be very slow
0400              // to converge so use asymptotic expansion for large
0401              // x case instead.  Note that the asymptotic expansion
0402              // isn't very accurate - so it's deliberately very hard
0403              // to get here - probably we're going to overflow:
0404              Iv = asymptotic_bessel_i_large_x(v, x, pol);
0405           }
0406           else if ((v > 0) && (x / v < 0.25))
0407           {
0408              Iv = bessel_i_small_z_series(v, x, pol);
0409           }
0410           else
0411           {
0412              CF1_ik(v, x, &fv, pol);                         // continued fraction CF1_ik
0413              Iv = scale * W / (Kv * fv + Kv1);                  // Wronskian relation
0414           }
0415        }
0416        else
0417           Iv = boost::math::numeric_limits<T>::quiet_NaN(); // any value will do
0418     }
0419     if (reflect)
0420     {
0421         T z = (u + n % 2);
0422         T fact = (2 / pi<T>()) * (boost::math::sin_pi(z, pol) * Kv);
0423         if(fact == 0)
0424            *result_I = Iv;
0425         else if(tools::max_value<T>() * scale < fact)
0426            *result_I = (org_kind & need_i) ? T(sign(fact) * scale_sign * policies::raise_overflow_error<T>(function, nullptr, pol)) : T(0);
0427         else
0428          *result_I = Iv + fact / scale;   // reflection formula
0429     }
0430     else
0431     {
0432         *result_I = Iv;
0433     }
0434     if(tools::max_value<T>() * scale < Kv)
0435        *result_K = (org_kind & need_k) ? T(sign(Kv) * scale_sign * policies::raise_overflow_error<T>(function, nullptr, pol)) : T(0);
0436     else
0437       *result_K = Kv / scale;
0438     BOOST_MATH_INSTRUMENT_VARIABLE(*result_I);
0439     BOOST_MATH_INSTRUMENT_VARIABLE(*result_K);
0440     return 0;
0441 }
0442
0443 }}} // namespaces
0444
0445 #endif // BOOST_MATH_BESSEL_IK_HPP
0446