ceph/src/boost/boost/math/special_functions/detail/bessel_k1.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_K1_HPP
   8 #define BOOST_MATH_BESSEL_K1_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/math/tools/assert.hpp>
  20
  21 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
  22 //
  23 // This is the only way we can avoid
  24 // warning: non-standard suffix on floating constant [-Wpedantic]
  25 // when building with -Wall -pedantic.  Neither __extension__
  26 // nor #pragma diagnostic ignored work :(
  27 //
  28 #pragma GCC system_header
  29 #endif
  30
  31 // Modified Bessel function of the second kind of order zero
  32 // minimax rational approximations on intervals, see
  33 // Russon and Blair, Chalk River Report AECL-3461, 1969,
  34 // as revised by Pavel Holoborodko in "Rational Approximations
  35 // for the Modified Bessel Function of the Second Kind - K0(x)
  36 // for Computations with Double Precision", see
  37 // http://www.advanpix.com/2016/01/05/rational-approximations-for-the-modified-bessel-function-of-the-second-kind-k1-for-computations-with-double-precision/
  38 //
  39 // The actual coefficients used are our own derivation (by JM)
  40 // since we extend to both greater and lesser precision than the
  41 // references above.  We can also improve performance WRT to
  42 // Holoborodko without loss of precision.
  43
  44 namespace boost { namespace math { namespace detail{
  45
  46    template <typename T>
  47    T bessel_k1(const T&);
  48
  49    template <class T, class tag>
  50    struct bessel_k1_initializer
  51    {
  52       struct init
  53       {
  54          init()
  55          {
  56             do_init(tag());
  57          }
  58          static void do_init(const std::integral_constant<int, 113>&)
  59          {
  60             bessel_k1(T(0.5));
  61             bessel_k1(T(2));
  62             bessel_k1(T(6));
  63          }
  64          static void do_init(const std::integral_constant<int, 64>&)
  65          {
  66             bessel_k1(T(0.5));
  67             bessel_k1(T(6));
  68          }
  69          template <class U>
  70          static void do_init(const U&) {}
  71          void force_instantiate()const {}
  72       };
  73       static const init initializer;
  74       static void force_instantiate()
  75       {
  76          initializer.force_instantiate();
  77       }
  78    };
  79
  80    template <class T, class tag>
  81    const typename bessel_k1_initializer<T, tag>::init bessel_k1_initializer<T, tag>::initializer;
  82
  83
  84    template <typename T, int N>
  85    inline T bessel_k1_imp(const T&, const std::integral_constant<int, N>&)
  86    {
  87       BOOST_MATH_ASSERT(0);
  88       return 0;
  89    }
  90
  91    template <typename T>
  92    T bessel_k1_imp(const T& x, const std::integral_constant<int, 24>&)
  93    {
  94       BOOST_MATH_STD_USING
  95       if(x <= 1)
  96       {
  97          // Maximum Deviation Found:                     3.090e-12
  98          // Expected Error Term : -3.053e-12
  99          // Maximum Relative Change in Control Points : 4.927e-02
 100          // Max Error found at float precision = Poly : 7.918347e-10
 101          static const T Y = 8.695471287e-02f;
 102          static const T P[] =
 103          {
 104             -3.621379531e-03f,
 105             7.131781976e-03f,
 106             -1.535278300e-05f
 107          };
 108          static const T Q[] =
 109          {
 110             1.000000000e+00f,
 111             -5.173102701e-02f,
 112             9.203530671e-04f
 113          };
 114
 115          T a = x * x / 4;
 116          a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2;
 117
 118          // Maximum Deviation Found:                     3.556e-08
 119          // Expected Error Term : -3.541e-08
 120          // Maximum Relative Change in Control Points : 8.203e-02
 121          static const T P2[] =
 122          {
 123             -3.079657469e-01f,
 124             -8.537108913e-02f,
 125             -4.640275408e-03f,
 126             -1.156442414e-04f
 127          };
 128
 129          return tools::evaluate_polynomial(P2, T(x * x)) * x + 1 / x + log(x) * a;
 130       }
 131       else
 132       {
 133          // Maximum Deviation Found:                     3.369e-08
 134          // Expected Error Term : -3.227e-08
 135          // Maximum Relative Change in Control Points : 9.917e-02
 136          // Max Error found at float precision = Poly : 6.084411e-08
 137          static const T Y = 1.450342178f;
 138          static const T P[] =
 139          {
 140             -1.970280088e-01f,
 141             2.188747807e-02f,
 142             7.270394756e-01f,
 143             2.490678196e-01f
 144          };
 145          static const T Q[] =
 146          {
 147             1.000000000e+00f,
 148             2.274292882e+00f,
 149             9.904984851e-01f,
 150             4.585534549e-02f
 151          };
 152          if(x < tools::log_max_value<T>())
 153             return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
 154          else
 155          {
 156             T ex = exp(-x / 2);
 157             return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
 158          }
 159       }
 160    }
 161
 162    template <typename T>
 163    T bessel_k1_imp(const T& x, const std::integral_constant<int, 53>&)
 164    {
 165       BOOST_MATH_STD_USING
 166       if(x <= 1)
 167       {
 168          // Maximum Deviation Found:                     1.922e-17
 169          // Expected Error Term : 1.921e-17
 170          // Maximum Relative Change in Control Points : 5.287e-03
 171          // Max Error found at double precision = Poly : 2.004747e-17
 172          static const T Y = 8.69547128677368164e-02f;
 173          static const T P[] =
 174          {
 175             -3.62137953440350228e-03,
 176             7.11842087490330300e-03,
 177             1.00302560256614306e-05,
 178             1.77231085381040811e-06
 179          };
 180          static const T Q[] =
 181          {
 182             1.00000000000000000e+00,
 183             -4.80414794429043831e-02,
 184             9.85972641934416525e-04,
 185             -8.91196859397070326e-06
 186          };
 187
 188          T a = x * x / 4;
 189          a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2;
 190
 191          // Maximum Deviation Found:                     4.053e-17
 192          // Expected Error Term : -4.053e-17
 193          // Maximum Relative Change in Control Points : 3.103e-04
 194          // Max Error found at double precision = Poly : 1.246698e-16
 195
 196          static const T P2[] =
 197          {
 198             -3.07965757829206184e-01,
 199             -7.80929703673074907e-02,
 200             -2.70619343754051620e-03,
 201             -2.49549522229072008e-05
 202          };
 203          static const T Q2[] =
 204          {
 205             1.00000000000000000e+00,
 206             -2.36316836412163098e-02,
 207             2.64524577525962719e-04,
 208             -1.49749618004162787e-06
 209          };
 210
 211          return tools::evaluate_rational(P2, Q2, T(x * x)) * x + 1 / x + log(x) * a;
 212       }
 213       else
 214       {
 215          // Maximum Deviation Found:                     8.883e-17
 216          // Expected Error Term : -1.641e-17
 217          // Maximum Relative Change in Control Points : 2.786e-01
 218          // Max Error found at double precision = Poly : 1.258798e-16
 219
 220          static const T Y = 1.45034217834472656f;
 221          static const T P[] =
 222          {
 223             -1.97028041029226295e-01,
 224             -2.32408961548087617e+00,
 225             -7.98269784507699938e+00,
 226             -2.39968410774221632e+00,
 227             3.28314043780858713e+01,
 228             5.67713761158496058e+01,
 229             3.30907788466509823e+01,
 230             6.62582288933739787e+00,
 231             3.08851840645286691e-01
 232          };
 233          static const T Q[] =
 234          {
 235             1.00000000000000000e+00,
 236             1.41811409298826118e+01,
 237             7.35979466317556420e+01,
 238             1.77821793937080859e+02,
 239             2.11014501598705982e+02,
 240             1.19425262951064454e+02,
 241             2.88448064302447607e+01,
 242             2.27912927104139732e+00,
 243             2.50358186953478678e-02
 244          };
 245          if(x < tools::log_max_value<T>())
 246             return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
 247          else
 248          {
 249             T ex = exp(-x / 2);
 250             return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
 251          }
 252       }
 253    }
 254
 255    template <typename T>
 256    T bessel_k1_imp(const T& x, const std::integral_constant<int, 64>&)
 257    {
 258       BOOST_MATH_STD_USING
 259       if(x <= 1)
 260       {
 261          // Maximum Deviation Found:                     5.549e-23
 262          // Expected Error Term : -5.548e-23
 263          // Maximum Relative Change in Control Points : 2.002e-03
 264          // Max Error found at float80 precision = Poly : 9.352785e-22
 265          static const T Y = 8.695471286773681640625e-02f;
 266          static const T P[] =
 267          {
 268             BOOST_MATH_BIG_CONSTANT(T, 64, -3.621379534403483072861e-03),
 269             BOOST_MATH_BIG_CONSTANT(T, 64, 7.102135866103952705932e-03),
 270             BOOST_MATH_BIG_CONSTANT(T, 64, 4.167545240236717601167e-05),
 271             BOOST_MATH_BIG_CONSTANT(T, 64, 2.537484002571894870830e-06),
 272             BOOST_MATH_BIG_CONSTANT(T, 64, 6.603228256820000135990e-09)
 273          };
 274          static const T Q[] =
 275          {
 276             BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
 277             BOOST_MATH_BIG_CONSTANT(T, 64, -4.354457194045068370363e-02),
 278             BOOST_MATH_BIG_CONSTANT(T, 64, 8.709137201220209072820e-04),
 279             BOOST_MATH_BIG_CONSTANT(T, 64, -9.676151796359590545143e-06),
 280             BOOST_MATH_BIG_CONSTANT(T, 64, 5.162715192766245311659e-08)
 281          };
 282
 283          T a = x * x / 4;
 284          a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2;
 285
 286          // Maximum Deviation Found:                     1.995e-23
 287          // Expected Error Term : 1.995e-23
 288          // Maximum Relative Change in Control Points : 8.174e-04
 289          // Max Error found at float80 precision = Poly : 4.137325e-20
 290          static const T P2[] =
 291          {
 292             BOOST_MATH_BIG_CONSTANT(T, 64, -3.079657578292062244054e-01),
 293             BOOST_MATH_BIG_CONSTANT(T, 64, -7.963049154965966503231e-02),
 294             BOOST_MATH_BIG_CONSTANT(T, 64, -3.103277523735639924895e-03),
 295             BOOST_MATH_BIG_CONSTANT(T, 64, -4.023052834702215699504e-05),
 296             BOOST_MATH_BIG_CONSTANT(T, 64, -1.719459155018493821839e-07)
 297          };
 298          static const T Q2[] =
 299          {
 300             BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
 301             BOOST_MATH_BIG_CONSTANT(T, 64, -1.863917670410152669768e-02),
 302             BOOST_MATH_BIG_CONSTANT(T, 64, 1.699367098849735298090e-04),
 303             BOOST_MATH_BIG_CONSTANT(T, 64, -9.309358790546076298429e-07),
 304             BOOST_MATH_BIG_CONSTANT(T, 64, 2.708893480271612711933e-09)
 305          };
 306
 307          return tools::evaluate_rational(P2, Q2, T(x * x)) * x + 1 / x + log(x) * a;
 308       }
 309       else
 310       {
 311          // Maximum Deviation Found:                     9.785e-20
 312          // Expected Error Term : -3.302e-21
 313          // Maximum Relative Change in Control Points : 3.432e-01
 314          // Max Error found at float80 precision = Poly : 1.083755e-19
 315          static const T Y = 1.450342178344726562500e+00f;
 316          static const T P[] =
 317          {
 318             BOOST_MATH_BIG_CONSTANT(T, 64, -1.970280410292263112917e-01),
 319             BOOST_MATH_BIG_CONSTANT(T, 64, -4.058564803062959169322e+00),
 320             BOOST_MATH_BIG_CONSTANT(T, 64, -3.036658174194917777473e+01),
 321             BOOST_MATH_BIG_CONSTANT(T, 64, -9.576825392332820142173e+01),
 322             BOOST_MATH_BIG_CONSTANT(T, 64, -6.706969489248020941949e+01),
 323             BOOST_MATH_BIG_CONSTANT(T, 64, 3.264572499406168221382e+02),
 324             BOOST_MATH_BIG_CONSTANT(T, 64, 8.584972047303151034100e+02),
 325             BOOST_MATH_BIG_CONSTANT(T, 64, 8.422082733280017909550e+02),
 326             BOOST_MATH_BIG_CONSTANT(T, 64, 3.738005441471368178383e+02),
 327             BOOST_MATH_BIG_CONSTANT(T, 64, 7.016938390144121276609e+01),
 328             BOOST_MATH_BIG_CONSTANT(T, 64, 4.319614662598089438939e+00),
 329             BOOST_MATH_BIG_CONSTANT(T, 64, 3.710715864316521856193e-02)
 330          };
 331          static const T Q[] =
 332          {
 333             BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
 334             BOOST_MATH_BIG_CONSTANT(T, 64, 2.298433045824439052398e+01),
 335             BOOST_MATH_BIG_CONSTANT(T, 64, 2.082047745067709230037e+02),
 336             BOOST_MATH_BIG_CONSTANT(T, 64, 9.662367854250262046592e+02),
 337             BOOST_MATH_BIG_CONSTANT(T, 64, 2.504148628460454004686e+03),
 338             BOOST_MATH_BIG_CONSTANT(T, 64, 3.712730364911389908905e+03),
 339             BOOST_MATH_BIG_CONSTANT(T, 64, 3.108002081150068641112e+03),
 340             BOOST_MATH_BIG_CONSTANT(T, 64, 1.400149940532448553143e+03),
 341             BOOST_MATH_BIG_CONSTANT(T, 64, 3.083303048095846226299e+02),
 342             BOOST_MATH_BIG_CONSTANT(T, 64, 2.748706060530351833346e+01),
 343             BOOST_MATH_BIG_CONSTANT(T, 64, 6.321900849331506946977e-01),
 344          };
 345          if(x < tools::log_max_value<T>())
 346             return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
 347          else
 348          {
 349             T ex = exp(-x / 2);
 350             return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
 351          }
 352       }
 353    }
 354
 355    template <typename T>
 356    T bessel_k1_imp(const T& x, const std::integral_constant<int, 113>&)
 357    {
 358       BOOST_MATH_STD_USING
 359       if(x <= 1)
 360       {
 361          // Maximum Deviation Found:                     7.120e-35
 362          // Expected Error Term : -7.119e-35
 363          // Maximum Relative Change in Control Points : 1.207e-03
 364          // Max Error found at float128 precision = Poly : 7.143688e-35
 365          static const T Y = 8.695471286773681640625000000000000000e-02f;
 366          static const T P[] =
 367          {
 368             BOOST_MATH_BIG_CONSTANT(T, 113, -3.621379534403483072916666666666595475e-03),
 369             BOOST_MATH_BIG_CONSTANT(T, 113, 7.074117676930975433219826471336547627e-03),
 370             BOOST_MATH_BIG_CONSTANT(T, 113, 9.631337631362776369069668419033041661e-05),
 371             BOOST_MATH_BIG_CONSTANT(T, 113, 3.468935967870048731821071646104412775e-06),
 372             BOOST_MATH_BIG_CONSTANT(T, 113, 2.956705020559599861444492614737168261e-08),
 373             BOOST_MATH_BIG_CONSTANT(T, 113, 2.347140307321161346703214099534250263e-10),
 374             BOOST_MATH_BIG_CONSTANT(T, 113, 5.569608494081482873946791086435679661e-13)
 375          };
 376          static const T Q[] =
 377          {
 378             BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
 379             BOOST_MATH_BIG_CONSTANT(T, 113, -3.580768910152105375615558920428350204e-02),
 380             BOOST_MATH_BIG_CONSTANT(T, 113, 6.197467671701485365363068445534557369e-04),
 381             BOOST_MATH_BIG_CONSTANT(T, 113, -6.707466533308630411966030561446666237e-06),
 382             BOOST_MATH_BIG_CONSTANT(T, 113, 4.846687802282250112624373388491123527e-08),
 383             BOOST_MATH_BIG_CONSTANT(T, 113, -2.248493131151981569517383040323900343e-10),
 384             BOOST_MATH_BIG_CONSTANT(T, 113, 5.319279786372775264555728921709381080e-13)
 385          };
 386
 387          T a = x * x / 4;
 388          a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2;
 389
 390          // Maximum Deviation Found:                     4.473e-37
 391          // Expected Error Term : 4.473e-37
 392          // Maximum Relative Change in Control Points : 8.550e-04
 393          // Max Error found at float128 precision = Poly : 8.167701e-35
 394          static const T P2[] =
 395          {
 396             BOOST_MATH_BIG_CONSTANT(T, 113, -3.079657578292062244053600156878870690e-01),
 397             BOOST_MATH_BIG_CONSTANT(T, 113, -8.133183745732467770755578848987414875e-02),
 398             BOOST_MATH_BIG_CONSTANT(T, 113, -3.548968792764174773125420229299431951e-03),
 399             BOOST_MATH_BIG_CONSTANT(T, 113, -5.886125468718182876076972186152445490e-05),
 400             BOOST_MATH_BIG_CONSTANT(T, 113, -4.506712111733707245745396404449639865e-07),
 401             BOOST_MATH_BIG_CONSTANT(T, 113, -1.632502325880313239698965376754406011e-09),
 402             BOOST_MATH_BIG_CONSTANT(T, 113, -2.311973065898784812266544485665624227e-12)
 403          };
 404          static const T Q2[] =
 405          {
 406             BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
 407             BOOST_MATH_BIG_CONSTANT(T, 113, -1.311471216733781016657962995723287450e-02),
 408             BOOST_MATH_BIG_CONSTANT(T, 113, 8.571876054797365417068164018709472969e-05),
 409             BOOST_MATH_BIG_CONSTANT(T, 113, -3.630181215268238731442496851497901293e-07),
 410             BOOST_MATH_BIG_CONSTANT(T, 113, 1.070176111227805048604885986867484807e-09),
 411             BOOST_MATH_BIG_CONSTANT(T, 113, -2.129046580769872602793220056461084761e-12),
 412             BOOST_MATH_BIG_CONSTANT(T, 113, 2.294906469421390890762001971790074432e-15)
 413          };
 414
 415          return tools::evaluate_rational(P2, Q2, T(x * x)) * x + 1 / x + log(x) * a;
 416       }
 417       else if(x < 4)
 418       {
 419          // Max error in interpolated form: 5.307e-37
 420          // Max Error found at float128 precision = Poly: 7.087862e-35
 421          static const T Y = 1.5023040771484375f;
 422          static const T P[] =
 423          {
 424             BOOST_MATH_BIG_CONSTANT(T, 113, -2.489899398329369710528254347931380044e-01),
 425             BOOST_MATH_BIG_CONSTANT(T, 113, -6.819080211203854781858815596508456873e+00),
 426             BOOST_MATH_BIG_CONSTANT(T, 113, -7.599915699069767382647695624952723034e+01),
 427             BOOST_MATH_BIG_CONSTANT(T, 113, -4.450211910821295507926582231071300718e+02),
 428             BOOST_MATH_BIG_CONSTANT(T, 113, -1.451374687870925175794150513723956533e+03),
 429             BOOST_MATH_BIG_CONSTANT(T, 113, -2.405805746895098802803503988539098226e+03),
 430             BOOST_MATH_BIG_CONSTANT(T, 113, -5.638808326778389656403861103277220518e+02),
 431             BOOST_MATH_BIG_CONSTANT(T, 113, 5.513958744081268456191778822780865708e+03),
 432             BOOST_MATH_BIG_CONSTANT(T, 113, 1.121301640926540743072258116122834804e+04),
 433             BOOST_MATH_BIG_CONSTANT(T, 113, 1.080094900175649541266613109971296190e+04),
 434             BOOST_MATH_BIG_CONSTANT(T, 113, 5.896531083639613332407534434915552429e+03),
 435             BOOST_MATH_BIG_CONSTANT(T, 113, 1.856602122319645694042555107114028437e+03),
 436             BOOST_MATH_BIG_CONSTANT(T, 113, 3.237121918853145421414003823957537419e+02),
 437             BOOST_MATH_BIG_CONSTANT(T, 113, 2.842072954561323076230238664623893504e+01),
 438             BOOST_MATH_BIG_CONSTANT(T, 113, 1.039705646510167437971862966128055524e+00),
 439             BOOST_MATH_BIG_CONSTANT(T, 113, 1.008418100718254816100425022904039530e-02)
 440          };
 441          static const T Q[] =
 442          {
 443             BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
 444             BOOST_MATH_BIG_CONSTANT(T, 113, 2.927456835239137986889227412815459529e+01),
 445             BOOST_MATH_BIG_CONSTANT(T, 113, 3.598985593265577043711382994516531273e+02),
 446             BOOST_MATH_BIG_CONSTANT(T, 113, 2.449897377085510281395819892689690579e+03),
 447             BOOST_MATH_BIG_CONSTANT(T, 113, 1.025555887684561913263090023158085327e+04),
 448             BOOST_MATH_BIG_CONSTANT(T, 113, 2.774140447181062463181892531100679195e+04),
 449             BOOST_MATH_BIG_CONSTANT(T, 113, 4.962055507843204417243602332246120418e+04),
 450             BOOST_MATH_BIG_CONSTANT(T, 113, 5.908269326976180183216954452196772931e+04),
 451             BOOST_MATH_BIG_CONSTANT(T, 113, 4.655160454422016855911700790722577942e+04),
 452             BOOST_MATH_BIG_CONSTANT(T, 113, 2.383586885019548163464418964577684608e+04),
 453             BOOST_MATH_BIG_CONSTANT(T, 113, 7.679920375586960324298491662159976419e+03),
 454             BOOST_MATH_BIG_CONSTANT(T, 113, 1.478586421028842906987799049804565008e+03),
 455             BOOST_MATH_BIG_CONSTANT(T, 113, 1.565384974896746094224942654383537090e+02),
 456             BOOST_MATH_BIG_CONSTANT(T, 113, 7.902617937084010911005732488607114511e+00),
 457             BOOST_MATH_BIG_CONSTANT(T, 113, 1.429293010387921526110949911029094926e-01),
 458             BOOST_MATH_BIG_CONSTANT(T, 113, 3.880342607911083143560111853491047663e-04)
 459          };
 460          return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
 461       }
 462       else
 463       {
 464          // Maximum Deviation Found:                     4.359e-37
 465          // Expected Error Term : -6.565e-40
 466          // Maximum Relative Change in Control Points : 1.880e-01
 467          // Max Error found at float128 precision = Poly : 2.943572e-35
 468          static const T Y = 1.308816909790039062500000000000000000f;
 469          static const T P[] =
 470          {
 471             BOOST_MATH_BIG_CONSTANT(T, 113, -5.550277247453881129211735759447737350e-02),
 472             BOOST_MATH_BIG_CONSTANT(T, 113, -3.485883080219574328217554864956175929e+00),
 473             BOOST_MATH_BIG_CONSTANT(T, 113, -8.903760658131484239300875153154881958e+01),
 474             BOOST_MATH_BIG_CONSTANT(T, 113, -1.144813672213626237418235110712293337e+03),
 475             BOOST_MATH_BIG_CONSTANT(T, 113, -6.498400501156131446691826557494158173e+03),
 476             BOOST_MATH_BIG_CONSTANT(T, 113, 1.573531831870363502604119835922166116e+04),
 477             BOOST_MATH_BIG_CONSTANT(T, 113, 5.417416550054632009958262596048841154e+05),
 478             BOOST_MATH_BIG_CONSTANT(T, 113, 4.271266450613557412825896604269130661e+06),
 479             BOOST_MATH_BIG_CONSTANT(T, 113, 1.898386013314389952534433455681107783e+07),
 480             BOOST_MATH_BIG_CONSTANT(T, 113, 5.353798784656436259250791761023512750e+07),
 481             BOOST_MATH_BIG_CONSTANT(T, 113, 9.839619195427352438957774052763490067e+07),
 482             BOOST_MATH_BIG_CONSTANT(T, 113, 1.169246368651532232388152442538005637e+08),
 483             BOOST_MATH_BIG_CONSTANT(T, 113, 8.696368884166831199967845883371116431e+07),
 484             BOOST_MATH_BIG_CONSTANT(T, 113, 3.810226630422736458064005843327500169e+07),
 485             BOOST_MATH_BIG_CONSTANT(T, 113, 8.854996610560406127438950635716757614e+06),
 486             BOOST_MATH_BIG_CONSTANT(T, 113, 8.981057433937398731355768088809437625e+05),
 487             BOOST_MATH_BIG_CONSTANT(T, 113, 2.519440069856232098711793483639792952e+04)
 488          };
 489          static const T Q[] =
 490          {
 491             BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
 492             BOOST_MATH_BIG_CONSTANT(T, 113, 7.127348248283623146544565916604103560e+01),
 493             BOOST_MATH_BIG_CONSTANT(T, 113, 2.205092684176906740104488180754982065e+03),
 494             BOOST_MATH_BIG_CONSTANT(T, 113, 3.911249195069050636298346469740075758e+04),
 495             BOOST_MATH_BIG_CONSTANT(T, 113, 4.426103406579046249654548481377792614e+05),
 496             BOOST_MATH_BIG_CONSTANT(T, 113, 3.365861555422488771286500241966208541e+06),
 497             BOOST_MATH_BIG_CONSTANT(T, 113, 1.765377714160383676864913709252529840e+07),
 498             BOOST_MATH_BIG_CONSTANT(T, 113, 6.453822726931857253365138260720815246e+07),
 499             BOOST_MATH_BIG_CONSTANT(T, 113, 1.643207885048369990391975749439783892e+08),
 500             BOOST_MATH_BIG_CONSTANT(T, 113, 2.882540678243694621895816336640877878e+08),
 501             BOOST_MATH_BIG_CONSTANT(T, 113, 3.410120808992380266174106812005338148e+08),
 502             BOOST_MATH_BIG_CONSTANT(T, 113, 2.628138016559335882019310900426773027e+08),
 503             BOOST_MATH_BIG_CONSTANT(T, 113, 1.250794693811010646965360198541047961e+08),
 504             BOOST_MATH_BIG_CONSTANT(T, 113, 3.378723408195485594610593014072950078e+07),
 505             BOOST_MATH_BIG_CONSTANT(T, 113, 4.488253856312453816451380319061865560e+06),
 506             BOOST_MATH_BIG_CONSTANT(T, 113, 2.202167197882689873967723350537104582e+05),
 507             BOOST_MATH_BIG_CONSTANT(T, 113, 1.673233230356966539460728211412989843e+03)
 508          };
 509          if(x < tools::log_max_value<T>())
 510             return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
 511          else
 512          {
 513             T ex = exp(-x / 2);
 514             return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
 515          }
 516       }
 517     }
 518
 519     template <typename T>
 520     T bessel_k1_imp(const T& x, const std::integral_constant<int, 0>&)
 521     {
 522        if(boost::math::tools::digits<T>() <= 24)
 523           return bessel_k1_imp(x, std::integral_constant<int, 24>());
 524        else if(boost::math::tools::digits<T>() <= 53)
 525           return bessel_k1_imp(x, std::integral_constant<int, 53>());
 526        else if(boost::math::tools::digits<T>() <= 64)
 527           return bessel_k1_imp(x, std::integral_constant<int, 64>());
 528        else if(boost::math::tools::digits<T>() <= 113)
 529           return bessel_k1_imp(x, std::integral_constant<int, 113>());
 530        BOOST_MATH_ASSERT(0);
 531        return 0;
 532     }
 533
 534     template <typename T>
 535    inline T bessel_k1(const T& x)
 536    {
 537       typedef std::integral_constant<int,
 538          ((std::numeric_limits<T>::digits == 0) || (std::numeric_limits<T>::radix != 2)) ?
 539          0 :
 540          std::numeric_limits<T>::digits <= 24 ?
 541          24 :
 542          std::numeric_limits<T>::digits <= 53 ?
 543          53 :
 544          std::numeric_limits<T>::digits <= 64 ?
 545          64 :
 546          std::numeric_limits<T>::digits <= 113 ?
 547          113 : -1
 548       > tag_type;
 549
 550       bessel_k1_initializer<T, tag_type>::force_instantiate();
 551       return bessel_k1_imp(x, tag_type());
 552    }
 553
 554 }}} // namespaces
 555
 556 #ifdef _MSC_VER
 557 #pragma warning(pop)
 558 #endif
 559
 560 #endif // BOOST_MATH_BESSEL_K1_HPP
 561