special_functions/detail/bessel_ik.hpp

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