boost/math/special_functions/detail/bessel_ik.hpp

   1 //  Copyright (c) 2006 Xiaogang Zhang
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_BESSEL_IK_HPP
   7 #define BOOST_MATH_BESSEL_IK_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/math/special_functions/round.hpp>
  14 #include <boost/math/special_functions/gamma.hpp>
  15 #include <boost/math/special_functions/sin_pi.hpp>
  16 #include <boost/math/constants/constants.hpp>
  17 #include <boost/math/policies/error_handling.hpp>
  18 #include <boost/math/tools/config.hpp>
  19
  20 // Modified Bessel functions of the first and second kind of fractional order
  21
  22 namespace boost { namespace math {
  23
  24 namespace detail {
  25
  26 template <class T, class Policy>
  27 struct cyl_bessel_i_small_z
  28 {
  29    typedef T result_type;
  30
  31    cyl_bessel_i_small_z(T v_, T z_) : k(0), v(v_), mult(z_*z_/4)
  32    {
  33       BOOST_MATH_STD_USING
  34       term = 1;
  35    }
  36
  37    T operator()()
  38    {
  39       T result = term;
  40       ++k;
  41       term *= mult / k;
  42       term /= k + v;
  43       return result;
  44    }
  45 private:
  46    unsigned k;
  47    T v;
  48    T term;
  49    T mult;
  50 };
  51
  52 template <class T, class Policy>
  53 inline T bessel_i_small_z_series(T v, T x, const Policy& pol)
  54 {
  55    BOOST_MATH_STD_USING
  56    T prefix;
  57    if(v < max_factorial<T>::value)
  58    {
  59       prefix = pow(x / 2, v) / boost::math::tgamma(v + 1, pol);
  60    }
  61    else
  62    {
  63       prefix = v * log(x / 2) - boost::math::lgamma(v + 1, pol);
  64       prefix = exp(prefix);
  65    }
  66    if(prefix == 0)
  67       return prefix;
  68
  69    cyl_bessel_i_small_z<T, Policy> s(v, x);
  70    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
  71 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
  72    T zero = 0;
  73    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
  74 #else
  75    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
  76 #endif
  77    policies::check_series_iterations<T>("boost::math::bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
  78    return prefix * result;
  79 }
  80
  81 // Calculate K(v, x) and K(v+1, x) by method analogous to
  82 // Temme, Journal of Computational Physics, vol 21, 343 (1976)
  83 template <typename T, typename Policy>
  84 int temme_ik(T v, T x, T* K, T* K1, const Policy& pol)
  85 {
  86     T f, h, p, q, coef, sum, sum1, tolerance;
  87     T a, b, c, d, sigma, gamma1, gamma2;
  88     unsigned long k;
  89
  90     BOOST_MATH_STD_USING
  91     using namespace boost::math::tools;
  92     using namespace boost::math::constants;
  93
  94
  95     // |x| <= 2, Temme series converge rapidly
  96     // |x| > 2, the larger the |x|, the slower the convergence
  97     BOOST_ASSERT(abs(x) <= 2);
  98     BOOST_ASSERT(abs(v) <= 0.5f);
  99
 100     T gp = boost::math::tgamma1pm1(v, pol);
 101     T gm = boost::math::tgamma1pm1(-v, pol);
 102
 103     a = log(x / 2);
 104     b = exp(v * a);
 105     sigma = -a * v;
 106     c = abs(v) < tools::epsilon<T>() ?
 107        T(1) : T(boost::math::sin_pi(v) / (v * pi<T>()));
 108     d = abs(sigma) < tools::epsilon<T>() ?
 109         T(1) : T(sinh(sigma) / sigma);
 110     gamma1 = abs(v) < tools::epsilon<T>() ?
 111         T(-euler<T>()) : T((0.5f / v) * (gp - gm) * c);
 112     gamma2 = (2 + gp + gm) * c / 2;
 113
 114     // initial values
 115     p = (gp + 1) / (2 * b);
 116     q = (1 + gm) * b / 2;
 117     f = (cosh(sigma) * gamma1 + d * (-a) * gamma2) / c;
 118     h = p;
 119     coef = 1;
 120     sum = coef * f;
 121     sum1 = coef * h;
 122
 123     BOOST_MATH_INSTRUMENT_VARIABLE(p);
 124     BOOST_MATH_INSTRUMENT_VARIABLE(q);
 125     BOOST_MATH_INSTRUMENT_VARIABLE(f);
 126     BOOST_MATH_INSTRUMENT_VARIABLE(sigma);
 127     BOOST_MATH_INSTRUMENT_CODE(sinh(sigma));
 128     BOOST_MATH_INSTRUMENT_VARIABLE(gamma1);
 129     BOOST_MATH_INSTRUMENT_VARIABLE(gamma2);
 130     BOOST_MATH_INSTRUMENT_VARIABLE(c);
 131     BOOST_MATH_INSTRUMENT_VARIABLE(d);
 132     BOOST_MATH_INSTRUMENT_VARIABLE(a);
 133
 134     // series summation
 135     tolerance = tools::epsilon<T>();
 136     for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)
 137     {
 138         f = (k * f + p + q) / (k*k - v*v);
 139         p /= k - v;
 140         q /= k + v;
 141         h = p - k * f;
 142         coef *= x * x / (4 * k);
 143         sum += coef * f;
 144         sum1 += coef * h;
 145         if (abs(coef * f) < abs(sum) * tolerance)
 146         {
 147            break;
 148         }
 149     }
 150     policies::check_series_iterations<T>("boost::math::bessel_ik<%1%>(%1%,%1%) in temme_ik", k, pol);
 151
 152     *K = sum;
 153     *K1 = 2 * sum1 / x;
 154
 155     return 0;
 156 }
 157
 158 // Evaluate continued fraction fv = I_(v+1) / I_v, derived from
 159 // Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73
 160 template <typename T, typename Policy>
 161 int CF1_ik(T v, T x, T* fv, const Policy& pol)
 162 {
 163     T C, D, f, a, b, delta, tiny, tolerance;
 164     unsigned long k;
 165
 166     BOOST_MATH_STD_USING
 167
 168     // |x| <= |v|, CF1_ik converges rapidly
 169     // |x| > |v|, CF1_ik needs O(|x|) iterations to converge
 170
 171     // modified Lentz's method, see
 172     // Lentz, Applied Optics, vol 15, 668 (1976)
 173     tolerance = 2 * tools::epsilon<T>();
 174     BOOST_MATH_INSTRUMENT_VARIABLE(tolerance);
 175     tiny = sqrt(tools::min_value<T>());
 176     BOOST_MATH_INSTRUMENT_VARIABLE(tiny);
 177     C = f = tiny;                           // b0 = 0, replace with tiny
 178     D = 0;
 179     for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)
 180     {
 181         a = 1;
 182         b = 2 * (v + k) / x;
 183         C = b + a / C;
 184         D = b + a * D;
 185         if (C == 0) { C = tiny; }
 186         if (D == 0) { D = tiny; }
 187         D = 1 / D;
 188         delta = C * D;
 189         f *= delta;
 190         BOOST_MATH_INSTRUMENT_VARIABLE(delta-1);
 191         if (abs(delta - 1) <= tolerance)
 192         {
 193            break;
 194         }
 195     }
 196     BOOST_MATH_INSTRUMENT_VARIABLE(k);
 197     policies::check_series_iterations<T>("boost::math::bessel_ik<%1%>(%1%,%1%) in CF1_ik", k, pol);
 198
 199     *fv = f;
 200
 201     return 0;
 202 }
 203
 204 // Calculate K(v, x) and K(v+1, x) by evaluating continued fraction
 205 // z1 / z0 = U(v+1.5, 2v+1, 2x) / U(v+0.5, 2v+1, 2x), see
 206 // Thompson and Barnett, Computer Physics Communications, vol 47, 245 (1987)
 207 template <typename T, typename Policy>
 208 int CF2_ik(T v, T x, T* Kv, T* Kv1, const Policy& pol)
 209 {
 210     BOOST_MATH_STD_USING
 211     using namespace boost::math::constants;
 212
 213     T S, C, Q, D, f, a, b, q, delta, tolerance, current, prev;
 214     unsigned long k;
 215
 216     // |x| >= |v|, CF2_ik converges rapidly
 217     // |x| -> 0, CF2_ik fails to converge
 218
 219     BOOST_ASSERT(abs(x) > 1);
 220
 221     // Steed's algorithm, see Thompson and Barnett,
 222     // Journal of Computational Physics, vol 64, 490 (1986)
 223     tolerance = tools::epsilon<T>();
 224     a = v * v - 0.25f;
 225     b = 2 * (x + 1);                              // b1
 226     D = 1 / b;                                    // D1 = 1 / b1
 227     f = delta = D;                                // f1 = delta1 = D1, coincidence
 228     prev = 0;                                     // q0
 229     current = 1;                                  // q1
 230     Q = C = -a;                                   // Q1 = C1 because q1 = 1
 231     S = 1 + Q * delta;                            // S1
 232     BOOST_MATH_INSTRUMENT_VARIABLE(tolerance);
 233     BOOST_MATH_INSTRUMENT_VARIABLE(a);
 234     BOOST_MATH_INSTRUMENT_VARIABLE(b);
 235     BOOST_MATH_INSTRUMENT_VARIABLE(D);
 236     BOOST_MATH_INSTRUMENT_VARIABLE(f);
 237
 238     for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++)     // starting from 2
 239     {
 240         // continued fraction f = z1 / z0
 241         a -= 2 * (k - 1);
 242         b += 2;
 243         D = 1 / (b + a * D);
 244         delta *= b * D - 1;
 245         f += delta;
 246
 247         // series summation S = 1 + \sum_{n=1}^{\infty} C_n * z_n / z_0
 248         q = (prev - (b - 2) * current) / a;
 249         prev = current;
 250         current = q;                        // forward recurrence for q
 251         C *= -a / k;
 252         Q += C * q;
 253         S += Q * delta;
 254         //
 255         // Under some circumstances q can grow very small and C very
 256         // large, leading to under/overflow.  This is particularly an
 257         // issue for types which have many digits precision but a narrow
 258         // exponent range.  A typical example being a "double double" type.
 259         // To avoid this situation we can normalise q (and related prev/current)
 260         // and C.  All other variables remain unchanged in value.  A typical
 261         // test case occurs when x is close to 2, for example cyl_bessel_k(9.125, 2.125).
 262         //
 263         if(q < tools::epsilon<T>())
 264         {
 265            C *= q;
 266            prev /= q;
 267            current /= q;
 268            q = 1;
 269         }
 270
 271         // S converges slower than f
 272         BOOST_MATH_INSTRUMENT_VARIABLE(Q * delta);
 273         BOOST_MATH_INSTRUMENT_VARIABLE(abs(S) * tolerance);
 274         if (abs(Q * delta) < abs(S) * tolerance)
 275         {
 276            break;
 277         }
 278     }
 279     policies::check_series_iterations<T>("boost::math::bessel_ik<%1%>(%1%,%1%) in CF2_ik", k, pol);
 280
 281     *Kv = sqrt(pi<T>() / (2 * x)) * exp(-x) / S;
 282     *Kv1 = *Kv * (0.5f + v + x + (v * v - 0.25f) * f) / x;
 283     BOOST_MATH_INSTRUMENT_VARIABLE(*Kv);
 284     BOOST_MATH_INSTRUMENT_VARIABLE(*Kv1);
 285
 286     return 0;
 287 }
 288
 289 enum{
 290    need_i = 1,
 291    need_k = 2
 292 };
 293
 294 // Compute I(v, x) and K(v, x) simultaneously by Temme's method, see
 295 // Temme, Journal of Computational Physics, vol 19, 324 (1975)
 296 template <typename T, typename Policy>
 297 int bessel_ik(T v, T x, T* I, T* K, int kind, const Policy& pol)
 298 {
 299     // Kv1 = K_(v+1), fv = I_(v+1) / I_v
 300     // Ku1 = K_(u+1), fu = I_(u+1) / I_u
 301     T u, Iv, Kv, Kv1, Ku, Ku1, fv;
 302     T W, current, prev, next;
 303     bool reflect = false;
 304     unsigned n, k;
 305     int org_kind = kind;
 306     BOOST_MATH_INSTRUMENT_VARIABLE(v);
 307     BOOST_MATH_INSTRUMENT_VARIABLE(x);
 308     BOOST_MATH_INSTRUMENT_VARIABLE(kind);
 309
 310     BOOST_MATH_STD_USING
 311     using namespace boost::math::tools;
 312     using namespace boost::math::constants;
 313
 314     static const char* function = "boost::math::bessel_ik<%1%>(%1%,%1%)";
 315
 316     if (v < 0)
 317     {
 318         reflect = true;
 319         v = -v;                             // v is non-negative from here
 320         kind |= need_k;
 321     }
 322     n = iround(v, pol);
 323     u = v - n;                              // -1/2 <= u < 1/2
 324     BOOST_MATH_INSTRUMENT_VARIABLE(n);
 325     BOOST_MATH_INSTRUMENT_VARIABLE(u);
 326
 327     if (x < 0)
 328     {
 329        *I = *K = policies::raise_domain_error<T>(function,
 330             "Got x = %1% but real argument x must be non-negative, complex number result not supported.", x, pol);
 331         return 1;
 332     }
 333     if (x == 0)
 334     {
 335        Iv = (v == 0) ? static_cast<T>(1) : static_cast<T>(0);
 336        if(kind & need_k)
 337        {
 338          Kv = policies::raise_overflow_error<T>(function, 0, pol);
 339        }
 340        else
 341        {
 342           Kv = std::numeric_limits<T>::quiet_NaN(); // any value will do
 343        }
 344
 345        if(reflect && (kind & need_i))
 346        {
 347            T z = (u + n % 2);
 348            Iv = boost::math::sin_pi(z, pol) == 0 ?
 349                Iv :
 350                policies::raise_overflow_error<T>(function, 0, pol);   // reflection formula
 351        }
 352
 353        *I = Iv;
 354        *K = Kv;
 355        return 0;
 356     }
 357
 358     // x is positive until reflection
 359     W = 1 / x;                                 // Wronskian
 360     if (x <= 2)                                // x in (0, 2]
 361     {
 362         temme_ik(u, x, &Ku, &Ku1, pol);             // Temme series
 363     }
 364     else                                       // x in (2, \infty)
 365     {
 366         CF2_ik(u, x, &Ku, &Ku1, pol);               // continued fraction CF2_ik
 367     }
 368     BOOST_MATH_INSTRUMENT_VARIABLE(Ku);
 369     BOOST_MATH_INSTRUMENT_VARIABLE(Ku1);
 370     prev = Ku;
 371     current = Ku1;
 372     T scale = 1;
 373     for (k = 1; k <= n; k++)                   // forward recurrence for K
 374     {
 375         T fact = 2 * (u + k) / x;
 376         if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
 377         {
 378            prev /= current;
 379            scale /= current;
 380            current = 1;
 381         }
 382         next = fact * current + prev;
 383         prev = current;
 384         current = next;
 385     }
 386     Kv = prev;
 387     Kv1 = current;
 388     BOOST_MATH_INSTRUMENT_VARIABLE(Kv);
 389     BOOST_MATH_INSTRUMENT_VARIABLE(Kv1);
 390     if(kind & need_i)
 391     {
 392        T lim = (4 * v * v + 10) / (8 * x);
 393        lim *= lim;
 394        lim *= lim;
 395        lim /= 24;
 396        if((lim < tools::epsilon<T>() * 10) && (x > 100))
 397        {
 398           // x is huge compared to v, CF1 may be very slow
 399           // to converge so use asymptotic expansion for large
 400           // x case instead.  Note that the asymptotic expansion
 401           // isn't very accurate - so it's deliberately very hard
 402           // to get here - probably we're going to overflow:
 403           Iv = asymptotic_bessel_i_large_x(v, x, pol);
 404        }
 405        else if((v > 0) && (x / v < 0.25))
 406        {
 407           Iv = bessel_i_small_z_series(v, x, pol);
 408        }
 409        else
 410        {
 411           CF1_ik(v, x, &fv, pol);                         // continued fraction CF1_ik
 412           Iv = scale * W / (Kv * fv + Kv1);                  // Wronskian relation
 413        }
 414     }
 415     else
 416        Iv = std::numeric_limits<T>::quiet_NaN(); // any value will do
 417
 418     if (reflect)
 419     {
 420         T z = (u + n % 2);
 421         T fact = (2 / pi<T>()) * (boost::math::sin_pi(z) * Kv);
 422         if(fact == 0)
 423            *I = Iv;
 424         else if(tools::max_value<T>() * scale < fact)
 425            *I = (org_kind & need_i) ? T(sign(fact) * sign(scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
 426         else
 427          *I = Iv + fact / scale;   // reflection formula
 428     }
 429     else
 430     {
 431         *I = Iv;
 432     }
 433     if(tools::max_value<T>() * scale < Kv)
 434       *K = (org_kind & need_k) ? T(sign(Kv) * sign(scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
 435     else
 436       *K = Kv / scale;
 437     BOOST_MATH_INSTRUMENT_VARIABLE(*I);
 438     BOOST_MATH_INSTRUMENT_VARIABLE(*K);
 439     return 0;
 440 }
 441
 442 }}} // namespaces
 443
 444 #endif // BOOST_MATH_BESSEL_IK_HPP
 445