boost/math/special_functions/detail/bessel_k0.hpp

   1 //  Copyright (c) 2006 Xiaogang Zhang
   2 //  Copyright (c) 2017 John Maddock
   3 //  Use, modification and distribution are subject to the
   4 //  Boost Software License, Version 1.0. (See accompanying file
   5 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #ifndef BOOST_MATH_BESSEL_K0_HPP
   8 #define BOOST_MATH_BESSEL_K0_HPP
   9
  10 #ifdef _MSC_VER
  11 #pragma once
  12 #pragma warning(push)
  13 #pragma warning(disable:4702) // Unreachable code (release mode only warning)
  14 #endif
  15
  16 #include <boost/math/tools/rational.hpp>
  17 #include <boost/math/tools/big_constant.hpp>
  18 #include <boost/math/policies/error_handling.hpp>
  19 #include <boost/assert.hpp>
  20
  21 // Modified Bessel function of the second kind of order zero
  22 // minimax rational approximations on intervals, see
  23 // Russon and Blair, Chalk River Report AECL-3461, 1969,
  24 // as revised by Pavel Holoborodko in "Rational Approximations
  25 // for the Modified Bessel Function of the Second Kind - K0(x)
  26 // for Computations with Double Precision", see
  27 // http://www.advanpix.com/2015/11/25/rational-approximations-for-the-modified-bessel-function-of-the-second-kind-k0-for-computations-with-double-precision/
  28 //
  29 // The actual coefficients used are our own derivation (by JM)
  30 // since we extend to both greater and lesser precision than the
  31 // references above.  We can also improve performance WRT to
  32 // Holoborodko without loss of precision.
  33
  34 namespace boost { namespace math { namespace detail{
  35
  36 template <typename T>
  37 T bessel_k0(const T& x);
  38
  39 template <class T, class tag>
  40 struct bessel_k0_initializer
  41 {
  42    struct init
  43    {
  44       init()
  45       {
  46          do_init(tag());
  47       }
  48       static void do_init(const mpl::int_<113>&)
  49       {
  50          bessel_k0(T(0.5));
  51          bessel_k0(T(1.5));
  52       }
  53       static void do_init(const mpl::int_<64>&)
  54       {
  55          bessel_k0(T(0.5));
  56          bessel_k0(T(1.5));
  57       }
  58       template <class U>
  59       static void do_init(const U&){}
  60       void force_instantiate()const{}
  61    };
  62    static const init initializer;
  63    static void force_instantiate()
  64    {
  65       initializer.force_instantiate();
  66    }
  67 };
  68
  69 template <class T, class tag>
  70 const typename bessel_k0_initializer<T, tag>::init bessel_k0_initializer<T, tag>::initializer;
  71
  72
  73 template <typename T, int N>
  74 T bessel_k0_imp(const T& x, const mpl::int_<N>&)
  75 {
  76    BOOST_ASSERT(0);
  77    return 0;
  78 }
  79
  80 template <typename T>
  81 T bessel_k0_imp(const T& x, const mpl::int_<24>&)
  82 {
  83    BOOST_MATH_STD_USING
  84    if(x <= 1)
  85    {
  86       // Maximum Deviation Found : 2.358e-09
  87       // Expected Error Term : -2.358e-09
  88       // Maximum Relative Change in Control Points : 9.552e-02
  89       // Max Error found at float precision = Poly : 4.448220e-08
  90       static const T Y = 1.137250900268554688f;
  91       static const T P[] =
  92       {
  93          -1.372508979104259711e-01f,
  94          2.622545986273687617e-01f,
  95          5.047103728247919836e-03f
  96       };
  97       static const T Q[] =
  98       {
  99          1.000000000000000000e+00f,
 100          -8.928694018000029415e-02f,
 101          2.985980684180969241e-03f
 102       };
 103       T a = x * x / 4;
 104       a = (tools::evaluate_rational(P, Q, a) + Y) * a + 1;
 105
 106       // Maximum Deviation Found:                     1.346e-09
 107       // Expected Error Term : -1.343e-09
 108       // Maximum Relative Change in Control Points : 2.405e-02
 109       // Max Error found at float precision = Poly : 1.354814e-07
 110       static const T P2[] = {
 111          1.159315158e-01f,
 112          2.789828686e-01f,
 113          2.524902861e-02f,
 114          8.457241514e-04f,
 115          1.530051997e-05f
 116       };
 117       return tools::evaluate_polynomial(P2, T(x * x)) - log(x) * a;
 118    }
 119    else
 120    {
 121       // Maximum Deviation Found:                     1.587e-08
 122       // Expected Error Term : 1.531e-08
 123       // Maximum Relative Change in Control Points : 9.064e-02
 124       // Max Error found at float precision = Poly : 5.065020e-08
 125
 126       static const T P[] =
 127       {
 128          2.533141220e-01,
 129          5.221502603e-01,
 130          6.380180669e-02,
 131          -5.934976547e-02
 132       };
 133       static const T Q[] =
 134       {
 135          1.000000000e+00,
 136          2.679722431e+00,
 137          1.561635813e+00,
 138          1.573660661e-01
 139       };
 140       if(x < tools::log_max_value<T>())
 141          return ((tools::evaluate_rational(P, Q, T(1 / x)) + 1) * exp(-x) / sqrt(x));
 142       else
 143       {
 144          T ex = exp(-x / 2);
 145          return ((tools::evaluate_rational(P, Q, T(1 / x)) + 1) * ex / sqrt(x)) * ex;
 146       }
 147    }
 148 }
 149
 150 template <typename T>
 151 T bessel_k0_imp(const T& x, const mpl::int_<53>&)
 152 {
 153    BOOST_MATH_STD_USING
 154    if(x <= 1)
 155    {
 156       // Maximum Deviation Found:                     6.077e-17
 157       // Expected Error Term : -6.077e-17
 158       // Maximum Relative Change in Control Points : 7.797e-02
 159       // Max Error found at double precision = Poly : 1.003156e-16
 160       static const T Y = 1.137250900268554688;
 161       static const T P[] =
 162       {
 163          -1.372509002685546267e-01,
 164          2.574916117833312855e-01,
 165          1.395474602146869316e-02,
 166          5.445476986653926759e-04,
 167          7.125159422136622118e-06
 168       };
 169       static const T Q[] =
 170       {
 171          1.000000000000000000e+00,
 172          -5.458333438017788530e-02,
 173          1.291052816975251298e-03,
 174          -1.367653946978586591e-05
 175       };
 176
 177       T a = x * x / 4;
 178       a = (tools::evaluate_polynomial(P, a) / tools::evaluate_polynomial(Q, a) + Y) * a + 1;
 179
 180       // Maximum Deviation Found:                     3.429e-18
 181       // Expected Error Term : 3.392e-18
 182       // Maximum Relative Change in Control Points : 2.041e-02
 183       // Max Error found at double precision = Poly : 2.513112e-16
 184       static const T P2[] =
 185       {
 186          1.159315156584124484e-01,
 187          2.789828789146031732e-01,
 188          2.524892993216121934e-02,
 189          8.460350907213637784e-04,
 190          1.491471924309617534e-05,
 191          1.627106892422088488e-07,
 192          1.208266102392756055e-09,
 193          6.611686391749704310e-12
 194       };
 195
 196       return tools::evaluate_polynomial(P2, T(x * x)) - log(x) * a;
 197    }
 198    else
 199    {
 200       // Maximum Deviation Found:                     4.316e-17
 201       // Expected Error Term : 9.570e-18
 202       // Maximum Relative Change in Control Points : 2.757e-01
 203       // Max Error found at double precision = Poly : 1.001560e-16
 204
 205       static const T Y = 1;
 206       static const T P[] =
 207       {
 208          2.533141373155002416e-01,
 209          3.628342133984595192e+00,
 210          1.868441889406606057e+01,
 211          4.306243981063412784e+01,
 212          4.424116209627428189e+01,
 213          1.562095339356220468e+01,
 214          -1.810138978229410898e+00,
 215          -1.414237994269995877e+00,
 216          -9.369168119754924625e-02
 217       };
 218       static const T Q[] =
 219       {
 220          1.000000000000000000e+00,
 221          1.494194694879908328e+01,
 222          8.265296455388554217e+01,
 223          2.162779506621866970e+02,
 224          2.845145155184222157e+02,
 225          1.851714491916334995e+02,
 226          5.486540717439723515e+01,
 227          6.118075837628957015e+00,
 228          1.586261269326235053e-01
 229       };
 230       if(x < tools::log_max_value<T>())
 231          return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
 232       else
 233       {
 234          T ex = exp(-x / 2);
 235          return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
 236       }
 237    }
 238 }
 239
 240 template <typename T>
 241 T bessel_k0_imp(const T& x, const mpl::int_<64>&)
 242 {
 243    BOOST_MATH_STD_USING
 244       if(x <= 1)
 245       {
 246          // Maximum Deviation Found:                     2.180e-22
 247          // Expected Error Term : 2.180e-22
 248          // Maximum Relative Change in Control Points : 2.943e-01
 249          // Max Error found at float80 precision = Poly : 3.923207e-20
 250          static const T Y = 1.137250900268554687500e+00;
 251          static const T P[] =
 252          {
 253             BOOST_MATH_BIG_CONSTANT(T, 64, -1.372509002685546875002e-01),
 254             BOOST_MATH_BIG_CONSTANT(T, 64, 2.566481981037407600436e-01),
 255             BOOST_MATH_BIG_CONSTANT(T, 64, 1.551881122448948854873e-02),
 256             BOOST_MATH_BIG_CONSTANT(T, 64, 6.646112454323276529650e-04),
 257             BOOST_MATH_BIG_CONSTANT(T, 64, 1.213747930378196492543e-05),
 258             BOOST_MATH_BIG_CONSTANT(T, 64, 9.423709328020389560844e-08)
 259          };
 260          static const T Q[] =
 261          {
 262             BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
 263             BOOST_MATH_BIG_CONSTANT(T, 64, -4.843828412587773008342e-02),
 264             BOOST_MATH_BIG_CONSTANT(T, 64, 1.088484822515098936140e-03),
 265             BOOST_MATH_BIG_CONSTANT(T, 64, -1.374724008530702784829e-05),
 266             BOOST_MATH_BIG_CONSTANT(T, 64, 8.452665455952581680339e-08)
 267          };
 268
 269
 270          T a = x * x / 4;
 271          a = (tools::evaluate_polynomial(P, a) / tools::evaluate_polynomial(Q, a) + Y) * a + 1;
 272
 273          // Maximum Deviation Found:                     2.440e-21
 274          // Expected Error Term : -2.434e-21
 275          // Maximum Relative Change in Control Points : 2.459e-02
 276          // Max Error found at float80 precision = Poly : 1.482487e-19
 277          static const T P2[] =
 278          {
 279             BOOST_MATH_BIG_CONSTANT(T, 64, 1.159315156584124488110e-01),
 280             BOOST_MATH_BIG_CONSTANT(T, 64, 2.764832791416047889734e-01),
 281             BOOST_MATH_BIG_CONSTANT(T, 64, 1.926062887220923354112e-02),
 282             BOOST_MATH_BIG_CONSTANT(T, 64, 3.660777862036966089410e-04),
 283             BOOST_MATH_BIG_CONSTANT(T, 64, 2.094942446930673386849e-06)
 284          };
 285          static const T Q2[] =
 286          {
 287             BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
 288             BOOST_MATH_BIG_CONSTANT(T, 64, -2.156100313881251616320e-02),
 289             BOOST_MATH_BIG_CONSTANT(T, 64, 2.315993873344905957033e-04),
 290             BOOST_MATH_BIG_CONSTANT(T, 64, -1.529444499350703363451e-06),
 291             BOOST_MATH_BIG_CONSTANT(T, 64, 5.524988589917857531177e-09)
 292          };
 293          return tools::evaluate_rational(P2, Q2, T(x * x)) - log(x) * a;
 294       }
 295       else
 296       {
 297          // Maximum Deviation Found:                     4.291e-20
 298          // Expected Error Term : 2.236e-21
 299          // Maximum Relative Change in Control Points : 3.021e-01
 300          //Max Error found at float80 precision = Poly : 8.727378e-20
 301          static const T Y = 1;
 302          static const T P[] =
 303          {
 304             BOOST_MATH_BIG_CONSTANT(T, 64, 2.533141373155002512056e-01),
 305             BOOST_MATH_BIG_CONSTANT(T, 64, 5.417942070721928652715e+00),
 306             BOOST_MATH_BIG_CONSTANT(T, 64, 4.477464607463971754433e+01),
 307             BOOST_MATH_BIG_CONSTANT(T, 64, 1.838745728725943889876e+02),
 308             BOOST_MATH_BIG_CONSTANT(T, 64, 4.009736314927811202517e+02),
 309             BOOST_MATH_BIG_CONSTANT(T, 64, 4.557411293123609803452e+02),
 310             BOOST_MATH_BIG_CONSTANT(T, 64, 2.360222564015361268955e+02),
 311             BOOST_MATH_BIG_CONSTANT(T, 64, 2.385435333168505701022e+01),
 312             BOOST_MATH_BIG_CONSTANT(T, 64, -1.750195760942181592050e+01),
 313             BOOST_MATH_BIG_CONSTANT(T, 64, -4.059789241612946683713e+00),
 314             BOOST_MATH_BIG_CONSTANT(T, 64, -1.612783121537333908889e-01)
 315          };
 316          static const T Q[] =
 317          {
 318             BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
 319             BOOST_MATH_BIG_CONSTANT(T, 64, 2.200669254769325861404e+01),
 320             BOOST_MATH_BIG_CONSTANT(T, 64, 1.900177593527144126549e+02),
 321             BOOST_MATH_BIG_CONSTANT(T, 64, 8.361003989965786932682e+02),
 322             BOOST_MATH_BIG_CONSTANT(T, 64, 2.041319870804843395893e+03),
 323             BOOST_MATH_BIG_CONSTANT(T, 64, 2.828491555113790345068e+03),
 324             BOOST_MATH_BIG_CONSTANT(T, 64, 2.190342229261529076624e+03),
 325             BOOST_MATH_BIG_CONSTANT(T, 64, 9.003330795963812219852e+02),
 326             BOOST_MATH_BIG_CONSTANT(T, 64, 1.773371397243777891569e+02),
 327             BOOST_MATH_BIG_CONSTANT(T, 64, 1.368634935531158398439e+01),
 328             BOOST_MATH_BIG_CONSTANT(T, 64, 2.543310879400359967327e-01)
 329          };
 330          if(x < tools::log_max_value<T>())
 331             return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
 332          else
 333          {
 334             T ex = exp(-x / 2);
 335             return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
 336          }
 337       }
 338 }
 339
 340 template <typename T>
 341 T bessel_k0_imp(const T& x, const mpl::int_<113>&)
 342 {
 343    BOOST_MATH_STD_USING
 344       if(x <= 1)
 345       {
 346          // Maximum Deviation Found:                     5.682e-37
 347          // Expected Error Term : 5.682e-37
 348          // Maximum Relative Change in Control Points : 6.094e-04
 349          // Max Error found at float128 precision = Poly : 5.338213e-35
 350          static const T Y = 1.137250900268554687500000000000000000e+00f;
 351          static const T P[] =
 352          {
 353             BOOST_MATH_BIG_CONSTANT(T, 113, -1.372509002685546875000000000000000006e-01),
 354             BOOST_MATH_BIG_CONSTANT(T, 113, 2.556212905071072782462974351698081303e-01),
 355             BOOST_MATH_BIG_CONSTANT(T, 113, 1.742459135264203478530904179889103929e-02),
 356             BOOST_MATH_BIG_CONSTANT(T, 113, 8.077860530453688571555479526961318918e-04),
 357             BOOST_MATH_BIG_CONSTANT(T, 113, 1.868173911669241091399374307788635148e-05),
 358             BOOST_MATH_BIG_CONSTANT(T, 113, 2.496405768838992243478709145123306602e-07),
 359             BOOST_MATH_BIG_CONSTANT(T, 113, 1.752489221949580551692915881999762125e-09),
 360             BOOST_MATH_BIG_CONSTANT(T, 113, 5.243010555737173524710512824955368526e-12)
 361          };
 362          static const T Q[] =
 363          {
 364             BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
 365             BOOST_MATH_BIG_CONSTANT(T, 113, -4.095631064064621099785696980653193721e-02),
 366             BOOST_MATH_BIG_CONSTANT(T, 113, 8.313880983725212151967078809725835532e-04),
 367             BOOST_MATH_BIG_CONSTANT(T, 113, -1.095229912293480063501285562382835142e-05),
 368             BOOST_MATH_BIG_CONSTANT(T, 113, 1.022828799511943141130509410251996277e-07),
 369             BOOST_MATH_BIG_CONSTANT(T, 113, -6.860874007419812445494782795829046836e-10),
 370             BOOST_MATH_BIG_CONSTANT(T, 113, 3.107297802344970725756092082686799037e-12),
 371             BOOST_MATH_BIG_CONSTANT(T, 113, -7.460529579244623559164763757787600944e-15)
 372          };
 373          T a = x * x / 4;
 374          a = (tools::evaluate_rational(P, Q, a) + Y) * a + 1;
 375
 376          // Maximum Deviation Found:                     5.173e-38
 377          // Expected Error Term : 5.105e-38
 378          // Maximum Relative Change in Control Points : 9.734e-03
 379          // Max Error found at float128 precision = Poly : 1.688806e-34
 380          static const T P2[] =
 381          {
 382             BOOST_MATH_BIG_CONSTANT(T, 113, 1.159315156584124488107200313757741370e-01),
 383             BOOST_MATH_BIG_CONSTANT(T, 113, 2.789828789146031122026800078439435369e-01),
 384             BOOST_MATH_BIG_CONSTANT(T, 113, 2.524892993216269451266750049024628432e-02),
 385             BOOST_MATH_BIG_CONSTANT(T, 113, 8.460350907082229957222453839935101823e-04),
 386             BOOST_MATH_BIG_CONSTANT(T, 113, 1.491471929926042875260452849503857976e-05),
 387             BOOST_MATH_BIG_CONSTANT(T, 113, 1.627105610481598430816014719558896866e-07),
 388             BOOST_MATH_BIG_CONSTANT(T, 113, 1.208426165007797264194914898538250281e-09),
 389             BOOST_MATH_BIG_CONSTANT(T, 113, 6.508697838747354949164182457073784117e-12),
 390             BOOST_MATH_BIG_CONSTANT(T, 113, 2.659784680639805301101014383907273109e-14),
 391             BOOST_MATH_BIG_CONSTANT(T, 113, 8.531090131964391104248859415958109654e-17),
 392             BOOST_MATH_BIG_CONSTANT(T, 113, 2.205195117066478034260323124669936314e-19),
 393             BOOST_MATH_BIG_CONSTANT(T, 113, 4.692219280289030165761119775783115426e-22),
 394             BOOST_MATH_BIG_CONSTANT(T, 113, 8.362350161092532344171965861545860747e-25),
 395             BOOST_MATH_BIG_CONSTANT(T, 113, 1.277990623924628999539014980773738258e-27)
 396          };
 397
 398          return tools::evaluate_polynomial(P2, T(x * x)) - log(x) * a;
 399       }
 400       else
 401       {
 402          // Maximum Deviation Found:                     1.462e-34
 403          // Expected Error Term : 4.917e-40
 404          // Maximum Relative Change in Control Points : 3.385e-01
 405          // Max Error found at float128 precision = Poly : 1.567573e-34
 406          static const T Y = 1;
 407          static const T P[] =
 408          {
 409             BOOST_MATH_BIG_CONSTANT(T, 113, 2.533141373155002512078826424055226265e-01),
 410             BOOST_MATH_BIG_CONSTANT(T, 113, 2.001949740768235770078339977110749204e+01),
 411             BOOST_MATH_BIG_CONSTANT(T, 113, 6.991516715983883248363351472378349986e+02),
 412             BOOST_MATH_BIG_CONSTANT(T, 113, 1.429587951594593159075690819360687720e+04),
 413             BOOST_MATH_BIG_CONSTANT(T, 113, 1.911933815201948768044660065771258450e+05),
 414             BOOST_MATH_BIG_CONSTANT(T, 113, 1.769943016204926614862175317962439875e+06),
 415             BOOST_MATH_BIG_CONSTANT(T, 113, 1.170866154649560750500954150401105606e+07),
 416             BOOST_MATH_BIG_CONSTANT(T, 113, 5.634687099724383996792011977705727661e+07),
 417             BOOST_MATH_BIG_CONSTANT(T, 113, 1.989524036456492581597607246664394014e+08),
 418             BOOST_MATH_BIG_CONSTANT(T, 113, 5.160394785715328062088529400178080360e+08),
 419             BOOST_MATH_BIG_CONSTANT(T, 113, 9.778173054417826368076483100902201433e+08),
 420             BOOST_MATH_BIG_CONSTANT(T, 113, 1.335667778588806892764139643950439733e+09),
 421             BOOST_MATH_BIG_CONSTANT(T, 113, 1.283635100080306980206494425043706838e+09),
 422             BOOST_MATH_BIG_CONSTANT(T, 113, 8.300616188213640626577036321085025855e+08),
 423             BOOST_MATH_BIG_CONSTANT(T, 113, 3.277591957076162984986406540894621482e+08),
 424             BOOST_MATH_BIG_CONSTANT(T, 113, 5.564360536834214058158565361486115932e+07),
 425             BOOST_MATH_BIG_CONSTANT(T, 113, -1.043505161612403359098596828115690596e+07),
 426             BOOST_MATH_BIG_CONSTANT(T, 113, -7.217035248223503605127967970903027314e+06),
 427             BOOST_MATH_BIG_CONSTANT(T, 113, -1.422938158797326748375799596769964430e+06),
 428             BOOST_MATH_BIG_CONSTANT(T, 113, -1.229125746200586805278634786674745210e+05),
 429             BOOST_MATH_BIG_CONSTANT(T, 113, -4.201632288615609937883545928660649813e+03),
 430             BOOST_MATH_BIG_CONSTANT(T, 113, -3.690820607338480548346746717311811406e+01)
 431          };
 432          static const T Q[] =
 433          {
 434             BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
 435             BOOST_MATH_BIG_CONSTANT(T, 113, 7.964877874035741452203497983642653107e+01),
 436             BOOST_MATH_BIG_CONSTANT(T, 113, 2.808929943826193766839360018583294769e+03),
 437             BOOST_MATH_BIG_CONSTANT(T, 113, 5.814524004679994110944366890912384139e+04),
 438             BOOST_MATH_BIG_CONSTANT(T, 113, 7.897794522506725610540209610337355118e+05),
 439             BOOST_MATH_BIG_CONSTANT(T, 113, 7.456339470955813675629523617440433672e+06),
 440             BOOST_MATH_BIG_CONSTANT(T, 113, 5.057818717813969772198911392875127212e+07),
 441             BOOST_MATH_BIG_CONSTANT(T, 113, 2.513821619536852436424913886081133209e+08),
 442             BOOST_MATH_BIG_CONSTANT(T, 113, 9.255938846873380596038513316919990776e+08),
 443             BOOST_MATH_BIG_CONSTANT(T, 113, 2.537077551699028079347581816919572141e+09),
 444             BOOST_MATH_BIG_CONSTANT(T, 113, 5.176769339768120752974843214652367321e+09),
 445             BOOST_MATH_BIG_CONSTANT(T, 113, 7.828722317390455845253191337207432060e+09),
 446             BOOST_MATH_BIG_CONSTANT(T, 113, 8.698864296569996402006511705803675890e+09),
 447             BOOST_MATH_BIG_CONSTANT(T, 113, 7.007803261356636409943826918468544629e+09),
 448             BOOST_MATH_BIG_CONSTANT(T, 113, 4.016564631288740308993071395104715469e+09),
 449             BOOST_MATH_BIG_CONSTANT(T, 113, 1.595893010619754750655947035567624730e+09),
 450             BOOST_MATH_BIG_CONSTANT(T, 113, 4.241241839120481076862742189989406856e+08),
 451             BOOST_MATH_BIG_CONSTANT(T, 113, 7.168778094393076220871007550235840858e+07),
 452             BOOST_MATH_BIG_CONSTANT(T, 113, 7.156200301360388147635052029404211109e+06),
 453             BOOST_MATH_BIG_CONSTANT(T, 113, 3.752130382550379886741949463587008794e+05),
 454             BOOST_MATH_BIG_CONSTANT(T, 113, 8.370574966987293592457152146806662562e+03),
 455             BOOST_MATH_BIG_CONSTANT(T, 113, 4.871254714311063594080644835895740323e+01)
 456          };
 457          if(x < tools::log_max_value<T>())
 458             return  ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
 459          else
 460          {
 461             T ex = exp(-x / 2);
 462             return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
 463          }
 464       }
 465 }
 466
 467 template <typename T>
 468 T bessel_k0_imp(const T& x, const mpl::int_<0>&)
 469 {
 470    if(boost::math::tools::digits<T>() <= 24)
 471       return bessel_k0_imp(x, mpl::int_<24>());
 472    else if(boost::math::tools::digits<T>() <= 53)
 473       return bessel_k0_imp(x, mpl::int_<53>());
 474    else if(boost::math::tools::digits<T>() <= 64)
 475       return bessel_k0_imp(x, mpl::int_<64>());
 476    else if(boost::math::tools::digits<T>() <= 113)
 477       return bessel_k0_imp(x, mpl::int_<113>());
 478    BOOST_ASSERT(0);
 479    return 0;
 480 }
 481
 482 template <typename T>
 483 inline T bessel_k0(const T& x)
 484 {
 485    typedef mpl::int_<
 486       std::numeric_limits<T>::digits == 0 ?
 487       0 :
 488       std::numeric_limits<T>::digits <= 24 ?
 489       24 :
 490       std::numeric_limits<T>::digits <= 53 ?
 491       53 :
 492       std::numeric_limits<T>::digits <= 64 ?
 493       64 :
 494       std::numeric_limits<T>::digits <= 113 ?
 495       113 : -1
 496    > tag_type;
 497
 498    bessel_k0_initializer<T, tag_type>::force_instantiate();
 499    return bessel_k0_imp(x, tag_type());
 500 }
 501
 502 }}} // namespaces
 503
 504 #ifdef _MSC_VER
 505 #pragma warning(pop)
 506 #endif
 507
 508 #endif // BOOST_MATH_BESSEL_K0_HPP
 509