Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-09-16 08:38:29 

0001 //  Copyright (c) 2006 Xiaogang Zhang
0002 //  Copyright (c) 2017 John Maddock
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_K0_HPP
0008 #define BOOST_MATH_BESSEL_K0_HPP
0009 
0010 #ifdef _MSC_VER
0011 #pragma once
0012 #pragma warning(push)
0013 #pragma warning(disable:4702) // Unreachable code (release mode only warning)
0014 #endif
0015 
0016 #include <boost/math/tools/config.hpp>
0017 #include <boost/math/tools/type_traits.hpp>
0018 #include <boost/math/tools/numeric_limits.hpp>
0019 #include <boost/math/tools/precision.hpp>
0020 #include <boost/math/tools/rational.hpp>
0021 #include <boost/math/tools/big_constant.hpp>
0022 #include <boost/math/tools/assert.hpp>
0023 #include <boost/math/policies/error_handling.hpp>
0024 
0025 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
0026 //
0027 // This is the only way we can avoid
0028 // warning: non-standard suffix on floating constant [-Wpedantic]
0029 // when building with -Wall -pedantic.  Neither __extension__
0030 // nor #pragma diagnostic ignored work :(
0031 //
0032 #pragma GCC system_header
0033 #endif
0034 
0035 // Modified Bessel function of the second kind of order zero
0036 // minimax rational approximations on intervals, see
0037 // Russon and Blair, Chalk River Report AECL-3461, 1969,
0038 // as revised by Pavel Holoborodko in "Rational Approximations 
0039 // for the Modified Bessel Function of the Second Kind - K0(x) 
0040 // for Computations with Double Precision", see 
0041 // http://www.advanpix.com/2015/11/25/rational-approximations-for-the-modified-bessel-function-of-the-second-kind-k0-for-computations-with-double-precision/
0042 //
0043 // The actual coefficients used are our own derivation (by JM)
0044 // since we extend to both greater and lesser precision than the
0045 // references above.  We can also improve performance WRT to
0046 // Holoborodko without loss of precision.
0047 
0048 namespace boost { namespace math { namespace detail{
0049 
0050 template <typename T>
0051 BOOST_MATH_GPU_ENABLED T bessel_k0(const T& x);
0052 
0053 template <class T, class tag>
0054 struct bessel_k0_initializer
0055 {
0056    struct init
0057    {
0058       BOOST_MATH_GPU_ENABLED init()
0059       {
0060          do_init(tag());
0061       }
0062       BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant<int, 113>&)
0063       {
0064          bessel_k0(T(0.5));
0065          bessel_k0(T(1.5));
0066       }
0067       BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant<int, 64>&)
0068       {
0069          bessel_k0(T(0.5));
0070          bessel_k0(T(1.5));
0071       }
0072       template <class U>
0073       BOOST_MATH_GPU_ENABLED static void do_init(const U&){}
0074       BOOST_MATH_GPU_ENABLED void force_instantiate()const{}
0075    };
0076    BOOST_MATH_STATIC const init initializer;
0077    BOOST_MATH_GPU_ENABLED static void force_instantiate()
0078    {
0079       #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0080       initializer.force_instantiate();
0081       #endif
0082    }
0083 };
0084 
0085 template <class T, class tag>
0086 const typename bessel_k0_initializer<T, tag>::init bessel_k0_initializer<T, tag>::initializer;
0087 
0088 
0089 template <typename T, int N>
0090 BOOST_MATH_GPU_ENABLED T bessel_k0_imp(const T&, const boost::math::integral_constant<int, N>&)
0091 {
0092    BOOST_MATH_ASSERT(0);
0093    return 0;
0094 }
0095 
0096 template <typename T>
0097 BOOST_MATH_GPU_ENABLED T bessel_k0_imp(const T& x, const boost::math::integral_constant<int, 24>&)
0098 {
0099    BOOST_MATH_STD_USING
0100    if(x <= 1)
0101    {
0102       // Maximum Deviation Found : 2.358e-09
0103       // Expected Error Term : -2.358e-09
0104       // Maximum Relative Change in Control Points : 9.552e-02
0105       // Max Error found at float precision = Poly : 4.448220e-08
0106       BOOST_MATH_STATIC const T Y = 1.137250900268554688f;
0107       BOOST_MATH_STATIC const T P[] = 
0108       {
0109          -1.372508979104259711e-01f,
0110          2.622545986273687617e-01f,
0111          5.047103728247919836e-03f
0112       };
0113       BOOST_MATH_STATIC const T Q[] = 
0114       {
0115          1.000000000000000000e+00f,
0116          -8.928694018000029415e-02f,
0117          2.985980684180969241e-03f
0118       };
0119       T a = x * x / 4;
0120       a = (tools::evaluate_rational(P, Q, a) + Y) * a + 1;
0121 
0122       // Maximum Deviation Found:                     1.346e-09
0123       // Expected Error Term : -1.343e-09
0124       // Maximum Relative Change in Control Points : 2.405e-02
0125       // Max Error found at float precision = Poly : 1.354814e-07
0126       BOOST_MATH_STATIC const T P2[] = {
0127          1.159315158e-01f,
0128          2.789828686e-01f,
0129          2.524902861e-02f,
0130          8.457241514e-04f,
0131          1.530051997e-05f
0132       };
0133       return tools::evaluate_polynomial(P2, T(x * x)) - log(x) * a;
0134    }
0135    else
0136    {
0137       // Maximum Deviation Found:                     1.587e-08
0138       // Expected Error Term : 1.531e-08
0139       // Maximum Relative Change in Control Points : 9.064e-02
0140       // Max Error found at float precision = Poly : 5.065020e-08
0141 
0142       BOOST_MATH_STATIC const T P[] =
0143       {
0144          2.533141220e-01f,
0145          5.221502603e-01f,
0146          6.380180669e-02f,
0147          -5.934976547e-02f
0148       };
0149       BOOST_MATH_STATIC const T Q[] =
0150       {
0151          1.000000000e+00f,
0152          2.679722431e+00f,
0153          1.561635813e+00f,
0154          1.573660661e-01f
0155       };
0156       if(x < tools::log_max_value<T>())
0157          return ((tools::evaluate_rational(P, Q, T(1 / x)) + 1) * exp(-x) / sqrt(x));
0158       else
0159       {
0160          T ex = exp(-x / 2);
0161          return ((tools::evaluate_rational(P, Q, T(1 / x)) + 1) * ex / sqrt(x)) * ex;
0162       }
0163    }
0164 }
0165 
0166 template <typename T>
0167 BOOST_MATH_GPU_ENABLED T bessel_k0_imp(const T& x, const boost::math::integral_constant<int, 53>&)
0168 {
0169    BOOST_MATH_STD_USING
0170    if(x <= 1)
0171    {
0172       // Maximum Deviation Found:                     6.077e-17
0173       // Expected Error Term : -6.077e-17
0174       // Maximum Relative Change in Control Points : 7.797e-02
0175       // Max Error found at double precision = Poly : 1.003156e-16
0176       BOOST_MATH_STATIC const T Y = 1.137250900268554688;
0177       BOOST_MATH_STATIC const T P[] =
0178       {
0179          -1.372509002685546267e-01,
0180          2.574916117833312855e-01,
0181          1.395474602146869316e-02,
0182          5.445476986653926759e-04,
0183          7.125159422136622118e-06
0184       };
0185       BOOST_MATH_STATIC const T Q[] =
0186       {
0187          1.000000000000000000e+00,
0188          -5.458333438017788530e-02,
0189          1.291052816975251298e-03,
0190          -1.367653946978586591e-05
0191       };
0192 
0193       T a = x * x / 4;
0194       a = (tools::evaluate_polynomial(P, a) / tools::evaluate_polynomial(Q, a) + Y) * a + 1;
0195 
0196       // Maximum Deviation Found:                     3.429e-18
0197       // Expected Error Term : 3.392e-18
0198       // Maximum Relative Change in Control Points : 2.041e-02
0199       // Max Error found at double precision = Poly : 2.513112e-16
0200       BOOST_MATH_STATIC const T P2[] =
0201       {
0202          1.159315156584124484e-01,
0203          2.789828789146031732e-01,
0204          2.524892993216121934e-02,
0205          8.460350907213637784e-04,
0206          1.491471924309617534e-05,
0207          1.627106892422088488e-07,
0208          1.208266102392756055e-09,
0209          6.611686391749704310e-12
0210       };
0211 
0212       return tools::evaluate_polynomial(P2, T(x * x)) - log(x) * a;
0213    }
0214    else
0215    {
0216       // Maximum Deviation Found:                     4.316e-17
0217       // Expected Error Term : 9.570e-18
0218       // Maximum Relative Change in Control Points : 2.757e-01
0219       // Max Error found at double precision = Poly : 1.001560e-16
0220 
0221       BOOST_MATH_STATIC const T Y = 1;
0222       BOOST_MATH_STATIC const T P[] =
0223       {
0224          2.533141373155002416e-01,
0225          3.628342133984595192e+00,
0226          1.868441889406606057e+01,
0227          4.306243981063412784e+01,
0228          4.424116209627428189e+01,
0229          1.562095339356220468e+01,
0230          -1.810138978229410898e+00,
0231          -1.414237994269995877e+00,
0232          -9.369168119754924625e-02
0233       };
0234       BOOST_MATH_STATIC const T Q[] =
0235       {
0236          1.000000000000000000e+00,
0237          1.494194694879908328e+01,
0238          8.265296455388554217e+01,
0239          2.162779506621866970e+02,
0240          2.845145155184222157e+02,
0241          1.851714491916334995e+02,
0242          5.486540717439723515e+01,
0243          6.118075837628957015e+00,
0244          1.586261269326235053e-01
0245       };
0246       if(x < tools::log_max_value<T>())
0247          return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
0248       else
0249       {
0250          T ex = exp(-x / 2);
0251          return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
0252       }
0253    }
0254 }
0255 
0256 template <typename T>
0257 BOOST_MATH_GPU_ENABLED T bessel_k0_imp(const T& x, const boost::math::integral_constant<int, 64>&)
0258 {
0259    BOOST_MATH_STD_USING
0260       if(x <= 1)
0261       {
0262          // Maximum Deviation Found:                     2.180e-22
0263          // Expected Error Term : 2.180e-22
0264          // Maximum Relative Change in Control Points : 2.943e-01
0265          // Max Error found at float80 precision = Poly : 3.923207e-20
0266          BOOST_MATH_STATIC const T Y = 1.137250900268554687500e+00;
0267          BOOST_MATH_STATIC const T P[] =
0268          {
0269             BOOST_MATH_BIG_CONSTANT(T, 64, -1.372509002685546875002e-01),
0270             BOOST_MATH_BIG_CONSTANT(T, 64, 2.566481981037407600436e-01),
0271             BOOST_MATH_BIG_CONSTANT(T, 64, 1.551881122448948854873e-02),
0272             BOOST_MATH_BIG_CONSTANT(T, 64, 6.646112454323276529650e-04),
0273             BOOST_MATH_BIG_CONSTANT(T, 64, 1.213747930378196492543e-05),
0274             BOOST_MATH_BIG_CONSTANT(T, 64, 9.423709328020389560844e-08)
0275          };
0276          BOOST_MATH_STATIC const T Q[] =
0277          {
0278             BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
0279             BOOST_MATH_BIG_CONSTANT(T, 64, -4.843828412587773008342e-02),
0280             BOOST_MATH_BIG_CONSTANT(T, 64, 1.088484822515098936140e-03),
0281             BOOST_MATH_BIG_CONSTANT(T, 64, -1.374724008530702784829e-05),
0282             BOOST_MATH_BIG_CONSTANT(T, 64, 8.452665455952581680339e-08)
0283          };
0284 
0285 
0286          T a = x * x / 4;
0287          a = (tools::evaluate_polynomial(P, a) / tools::evaluate_polynomial(Q, a) + Y) * a + 1;
0288 
0289          // Maximum Deviation Found:                     2.440e-21
0290          // Expected Error Term : -2.434e-21
0291          // Maximum Relative Change in Control Points : 2.459e-02
0292          // Max Error found at float80 precision = Poly : 1.482487e-19
0293          BOOST_MATH_STATIC const T P2[] =
0294          {
0295             BOOST_MATH_BIG_CONSTANT(T, 64, 1.159315156584124488110e-01),
0296             BOOST_MATH_BIG_CONSTANT(T, 64, 2.764832791416047889734e-01),
0297             BOOST_MATH_BIG_CONSTANT(T, 64, 1.926062887220923354112e-02),
0298             BOOST_MATH_BIG_CONSTANT(T, 64, 3.660777862036966089410e-04),
0299             BOOST_MATH_BIG_CONSTANT(T, 64, 2.094942446930673386849e-06)
0300          };
0301          BOOST_MATH_STATIC const T Q2[] =
0302          {
0303             BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
0304             BOOST_MATH_BIG_CONSTANT(T, 64, -2.156100313881251616320e-02),
0305             BOOST_MATH_BIG_CONSTANT(T, 64, 2.315993873344905957033e-04),
0306             BOOST_MATH_BIG_CONSTANT(T, 64, -1.529444499350703363451e-06),
0307             BOOST_MATH_BIG_CONSTANT(T, 64, 5.524988589917857531177e-09)
0308          };
0309          return tools::evaluate_rational(P2, Q2, T(x * x)) - log(x) * a;
0310       }
0311       else
0312       {
0313          // Maximum Deviation Found:                     4.291e-20
0314          // Expected Error Term : 2.236e-21
0315          // Maximum Relative Change in Control Points : 3.021e-01
0316          //Max Error found at float80 precision = Poly : 8.727378e-20
0317          BOOST_MATH_STATIC const T Y = 1;
0318          BOOST_MATH_STATIC const T P[] =
0319          {
0320             BOOST_MATH_BIG_CONSTANT(T, 64, 2.533141373155002512056e-01),
0321             BOOST_MATH_BIG_CONSTANT(T, 64, 5.417942070721928652715e+00),
0322             BOOST_MATH_BIG_CONSTANT(T, 64, 4.477464607463971754433e+01),
0323             BOOST_MATH_BIG_CONSTANT(T, 64, 1.838745728725943889876e+02),
0324             BOOST_MATH_BIG_CONSTANT(T, 64, 4.009736314927811202517e+02),
0325             BOOST_MATH_BIG_CONSTANT(T, 64, 4.557411293123609803452e+02),
0326             BOOST_MATH_BIG_CONSTANT(T, 64, 2.360222564015361268955e+02),
0327             BOOST_MATH_BIG_CONSTANT(T, 64, 2.385435333168505701022e+01),
0328             BOOST_MATH_BIG_CONSTANT(T, 64, -1.750195760942181592050e+01),
0329             BOOST_MATH_BIG_CONSTANT(T, 64, -4.059789241612946683713e+00),
0330             BOOST_MATH_BIG_CONSTANT(T, 64, -1.612783121537333908889e-01)
0331          };
0332          BOOST_MATH_STATIC const T Q[] =
0333          {
0334             BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
0335             BOOST_MATH_BIG_CONSTANT(T, 64, 2.200669254769325861404e+01),
0336             BOOST_MATH_BIG_CONSTANT(T, 64, 1.900177593527144126549e+02),
0337             BOOST_MATH_BIG_CONSTANT(T, 64, 8.361003989965786932682e+02),
0338             BOOST_MATH_BIG_CONSTANT(T, 64, 2.041319870804843395893e+03),
0339             BOOST_MATH_BIG_CONSTANT(T, 64, 2.828491555113790345068e+03),
0340             BOOST_MATH_BIG_CONSTANT(T, 64, 2.190342229261529076624e+03),
0341             BOOST_MATH_BIG_CONSTANT(T, 64, 9.003330795963812219852e+02),
0342             BOOST_MATH_BIG_CONSTANT(T, 64, 1.773371397243777891569e+02),
0343             BOOST_MATH_BIG_CONSTANT(T, 64, 1.368634935531158398439e+01),
0344             BOOST_MATH_BIG_CONSTANT(T, 64, 2.543310879400359967327e-01)
0345          };
0346          if(x < tools::log_max_value<T>())
0347             return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
0348          else
0349          {
0350             T ex = exp(-x / 2);
0351             return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
0352          }
0353       }
0354 }
0355 
0356 template <typename T>
0357 BOOST_MATH_GPU_ENABLED T bessel_k0_imp(const T& x, const boost::math::integral_constant<int, 113>&)
0358 {
0359    BOOST_MATH_STD_USING
0360       if(x <= 1)
0361       {
0362          // Maximum Deviation Found:                     5.682e-37
0363          // Expected Error Term : 5.682e-37
0364          // Maximum Relative Change in Control Points : 6.094e-04
0365          // Max Error found at float128 precision = Poly : 5.338213e-35
0366          BOOST_MATH_STATIC const T Y = 1.137250900268554687500000000000000000e+00f;
0367          BOOST_MATH_STATIC const T P[] =
0368          {
0369             BOOST_MATH_BIG_CONSTANT(T, 113, -1.372509002685546875000000000000000006e-01),
0370             BOOST_MATH_BIG_CONSTANT(T, 113, 2.556212905071072782462974351698081303e-01),
0371             BOOST_MATH_BIG_CONSTANT(T, 113, 1.742459135264203478530904179889103929e-02),
0372             BOOST_MATH_BIG_CONSTANT(T, 113, 8.077860530453688571555479526961318918e-04),
0373             BOOST_MATH_BIG_CONSTANT(T, 113, 1.868173911669241091399374307788635148e-05),
0374             BOOST_MATH_BIG_CONSTANT(T, 113, 2.496405768838992243478709145123306602e-07),
0375             BOOST_MATH_BIG_CONSTANT(T, 113, 1.752489221949580551692915881999762125e-09),
0376             BOOST_MATH_BIG_CONSTANT(T, 113, 5.243010555737173524710512824955368526e-12)
0377          };
0378          BOOST_MATH_STATIC const T Q[] =
0379          {
0380             BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
0381             BOOST_MATH_BIG_CONSTANT(T, 113, -4.095631064064621099785696980653193721e-02),
0382             BOOST_MATH_BIG_CONSTANT(T, 113, 8.313880983725212151967078809725835532e-04),
0383             BOOST_MATH_BIG_CONSTANT(T, 113, -1.095229912293480063501285562382835142e-05),
0384             BOOST_MATH_BIG_CONSTANT(T, 113, 1.022828799511943141130509410251996277e-07),
0385             BOOST_MATH_BIG_CONSTANT(T, 113, -6.860874007419812445494782795829046836e-10),
0386             BOOST_MATH_BIG_CONSTANT(T, 113, 3.107297802344970725756092082686799037e-12),
0387             BOOST_MATH_BIG_CONSTANT(T, 113, -7.460529579244623559164763757787600944e-15)
0388          };
0389          T a = x * x / 4;
0390          a = (tools::evaluate_rational(P, Q, a) + Y) * a + 1;
0391 
0392          // Maximum Deviation Found:                     5.173e-38
0393          // Expected Error Term : 5.105e-38
0394          // Maximum Relative Change in Control Points : 9.734e-03
0395          // Max Error found at float128 precision = Poly : 1.688806e-34
0396          BOOST_MATH_STATIC const T P2[] =
0397          {
0398             BOOST_MATH_BIG_CONSTANT(T, 113, 1.159315156584124488107200313757741370e-01),
0399             BOOST_MATH_BIG_CONSTANT(T, 113, 2.789828789146031122026800078439435369e-01),
0400             BOOST_MATH_BIG_CONSTANT(T, 113, 2.524892993216269451266750049024628432e-02),
0401             BOOST_MATH_BIG_CONSTANT(T, 113, 8.460350907082229957222453839935101823e-04),
0402             BOOST_MATH_BIG_CONSTANT(T, 113, 1.491471929926042875260452849503857976e-05),
0403             BOOST_MATH_BIG_CONSTANT(T, 113, 1.627105610481598430816014719558896866e-07),
0404             BOOST_MATH_BIG_CONSTANT(T, 113, 1.208426165007797264194914898538250281e-09),
0405             BOOST_MATH_BIG_CONSTANT(T, 113, 6.508697838747354949164182457073784117e-12),
0406             BOOST_MATH_BIG_CONSTANT(T, 113, 2.659784680639805301101014383907273109e-14),
0407             BOOST_MATH_BIG_CONSTANT(T, 113, 8.531090131964391104248859415958109654e-17),
0408             BOOST_MATH_BIG_CONSTANT(T, 113, 2.205195117066478034260323124669936314e-19),
0409             BOOST_MATH_BIG_CONSTANT(T, 113, 4.692219280289030165761119775783115426e-22),
0410             BOOST_MATH_BIG_CONSTANT(T, 113, 8.362350161092532344171965861545860747e-25),
0411             BOOST_MATH_BIG_CONSTANT(T, 113, 1.277990623924628999539014980773738258e-27)
0412          };
0413 
0414          return tools::evaluate_polynomial(P2, T(x * x)) - log(x) * a;
0415       }
0416       else
0417       {
0418          // Maximum Deviation Found:                     1.462e-34
0419          // Expected Error Term : 4.917e-40
0420          // Maximum Relative Change in Control Points : 3.385e-01
0421          // Max Error found at float128 precision = Poly : 1.567573e-34
0422          BOOST_MATH_STATIC const T Y = 1;
0423          BOOST_MATH_STATIC const T P[] =
0424          {
0425             BOOST_MATH_BIG_CONSTANT(T, 113, 2.533141373155002512078826424055226265e-01),
0426             BOOST_MATH_BIG_CONSTANT(T, 113, 2.001949740768235770078339977110749204e+01),
0427             BOOST_MATH_BIG_CONSTANT(T, 113, 6.991516715983883248363351472378349986e+02),
0428             BOOST_MATH_BIG_CONSTANT(T, 113, 1.429587951594593159075690819360687720e+04),
0429             BOOST_MATH_BIG_CONSTANT(T, 113, 1.911933815201948768044660065771258450e+05),
0430             BOOST_MATH_BIG_CONSTANT(T, 113, 1.769943016204926614862175317962439875e+06),
0431             BOOST_MATH_BIG_CONSTANT(T, 113, 1.170866154649560750500954150401105606e+07),
0432             BOOST_MATH_BIG_CONSTANT(T, 113, 5.634687099724383996792011977705727661e+07),
0433             BOOST_MATH_BIG_CONSTANT(T, 113, 1.989524036456492581597607246664394014e+08),
0434             BOOST_MATH_BIG_CONSTANT(T, 113, 5.160394785715328062088529400178080360e+08),
0435             BOOST_MATH_BIG_CONSTANT(T, 113, 9.778173054417826368076483100902201433e+08),
0436             BOOST_MATH_BIG_CONSTANT(T, 113, 1.335667778588806892764139643950439733e+09),
0437             BOOST_MATH_BIG_CONSTANT(T, 113, 1.283635100080306980206494425043706838e+09),
0438             BOOST_MATH_BIG_CONSTANT(T, 113, 8.300616188213640626577036321085025855e+08),
0439             BOOST_MATH_BIG_CONSTANT(T, 113, 3.277591957076162984986406540894621482e+08),
0440             BOOST_MATH_BIG_CONSTANT(T, 113, 5.564360536834214058158565361486115932e+07),
0441             BOOST_MATH_BIG_CONSTANT(T, 113, -1.043505161612403359098596828115690596e+07),
0442             BOOST_MATH_BIG_CONSTANT(T, 113, -7.217035248223503605127967970903027314e+06),
0443             BOOST_MATH_BIG_CONSTANT(T, 113, -1.422938158797326748375799596769964430e+06),
0444             BOOST_MATH_BIG_CONSTANT(T, 113, -1.229125746200586805278634786674745210e+05),
0445             BOOST_MATH_BIG_CONSTANT(T, 113, -4.201632288615609937883545928660649813e+03),
0446             BOOST_MATH_BIG_CONSTANT(T, 113, -3.690820607338480548346746717311811406e+01)
0447          };
0448          BOOST_MATH_STATIC const T Q[] =
0449          {
0450             BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
0451             BOOST_MATH_BIG_CONSTANT(T, 113, 7.964877874035741452203497983642653107e+01),
0452             BOOST_MATH_BIG_CONSTANT(T, 113, 2.808929943826193766839360018583294769e+03),
0453             BOOST_MATH_BIG_CONSTANT(T, 113, 5.814524004679994110944366890912384139e+04),
0454             BOOST_MATH_BIG_CONSTANT(T, 113, 7.897794522506725610540209610337355118e+05),
0455             BOOST_MATH_BIG_CONSTANT(T, 113, 7.456339470955813675629523617440433672e+06),
0456             BOOST_MATH_BIG_CONSTANT(T, 113, 5.057818717813969772198911392875127212e+07),
0457             BOOST_MATH_BIG_CONSTANT(T, 113, 2.513821619536852436424913886081133209e+08),
0458             BOOST_MATH_BIG_CONSTANT(T, 113, 9.255938846873380596038513316919990776e+08),
0459             BOOST_MATH_BIG_CONSTANT(T, 113, 2.537077551699028079347581816919572141e+09),
0460             BOOST_MATH_BIG_CONSTANT(T, 113, 5.176769339768120752974843214652367321e+09),
0461             BOOST_MATH_BIG_CONSTANT(T, 113, 7.828722317390455845253191337207432060e+09),
0462             BOOST_MATH_BIG_CONSTANT(T, 113, 8.698864296569996402006511705803675890e+09),
0463             BOOST_MATH_BIG_CONSTANT(T, 113, 7.007803261356636409943826918468544629e+09),
0464             BOOST_MATH_BIG_CONSTANT(T, 113, 4.016564631288740308993071395104715469e+09),
0465             BOOST_MATH_BIG_CONSTANT(T, 113, 1.595893010619754750655947035567624730e+09),
0466             BOOST_MATH_BIG_CONSTANT(T, 113, 4.241241839120481076862742189989406856e+08),
0467             BOOST_MATH_BIG_CONSTANT(T, 113, 7.168778094393076220871007550235840858e+07),
0468             BOOST_MATH_BIG_CONSTANT(T, 113, 7.156200301360388147635052029404211109e+06),
0469             BOOST_MATH_BIG_CONSTANT(T, 113, 3.752130382550379886741949463587008794e+05),
0470             BOOST_MATH_BIG_CONSTANT(T, 113, 8.370574966987293592457152146806662562e+03),
0471             BOOST_MATH_BIG_CONSTANT(T, 113, 4.871254714311063594080644835895740323e+01)
0472          };
0473          if(-x > tools::log_min_value<T>())
0474             return  ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
0475          else
0476          {
0477             T ex = exp(-x / 2);
0478             return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
0479          }
0480       }
0481 }
0482 
0483 template <typename T>
0484 BOOST_MATH_GPU_ENABLED T bessel_k0_imp(const T& x, const boost::math::integral_constant<int, 0>&)
0485 {
0486    if(boost::math::tools::digits<T>() <= 24)
0487       return bessel_k0_imp(x, boost::math::integral_constant<int, 24>());
0488    else if(boost::math::tools::digits<T>() <= 53)
0489       return bessel_k0_imp(x, boost::math::integral_constant<int, 53>());
0490    else if(boost::math::tools::digits<T>() <= 64)
0491       return bessel_k0_imp(x, boost::math::integral_constant<int, 64>());
0492    else if(boost::math::tools::digits<T>() <= 113)
0493       return bessel_k0_imp(x, boost::math::integral_constant<int, 113>());
0494    BOOST_MATH_ASSERT(0);
0495    return 0;
0496 }
0497 
0498 template <typename T>
0499 BOOST_MATH_GPU_ENABLED inline T bessel_k0(const T& x)
0500 {
0501    typedef boost::math::integral_constant<int,
0502       ((boost::math::numeric_limits<T>::digits == 0) || (boost::math::numeric_limits<T>::radix != 2)) ?
0503       0 :
0504       boost::math::numeric_limits<T>::digits <= 24 ?
0505       24 :
0506       boost::math::numeric_limits<T>::digits <= 53 ?
0507       53 :
0508       boost::math::numeric_limits<T>::digits <= 64 ?
0509       64 :
0510       boost::math::numeric_limits<T>::digits <= 113 ?
0511       113 : -1
0512    > tag_type;
0513 
0514    bessel_k0_initializer<T, tag_type>::force_instantiate();
0515    return bessel_k0_imp(x, tag_type());
0516 }
0517 
0518 }}} // namespaces
0519 
0520 #ifdef _MSC_VER
0521 #pragma warning(pop)
0522 #endif
0523 
0524 #endif // BOOST_MATH_BESSEL_K0_HPP
0525
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!