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_K1_HPP
0008 #define BOOST_MATH_BESSEL_K1_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/policies/error_handling.hpp>
0023 #include <boost/math/tools/assert.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/2016/01/05/rational-approximations-for-the-modified-bessel-function-of-the-second-kind-k1-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_k1(const T&);
0052 
0053    template <class T, class tag>
0054    struct bessel_k1_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_k1(T(0.5));
0065             bessel_k1(T(2));
0066             bessel_k1(T(6));
0067          }
0068          BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant<int, 64>&)
0069          {
0070             bessel_k1(T(0.5));
0071             bessel_k1(T(6));
0072          }
0073          template <class U>
0074          BOOST_MATH_GPU_ENABLED static void do_init(const U&) {}
0075          BOOST_MATH_GPU_ENABLED void force_instantiate()const {}
0076       };
0077       BOOST_MATH_STATIC const init initializer;
0078       BOOST_MATH_GPU_ENABLED static void force_instantiate()
0079       {
0080          #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0081          initializer.force_instantiate();
0082          #endif
0083       }
0084    };
0085 
0086    template <class T, class tag>
0087    const typename bessel_k1_initializer<T, tag>::init bessel_k1_initializer<T, tag>::initializer;
0088 
0089 
0090    template <typename T, int N>
0091    inline BOOST_MATH_GPU_ENABLED T bessel_k1_imp(const T&, const boost::math::integral_constant<int, N>&)
0092    {
0093       BOOST_MATH_ASSERT(0);
0094       return 0;
0095    }
0096 
0097    template <typename T>
0098    BOOST_MATH_GPU_ENABLED T bessel_k1_imp(const T& x, const boost::math::integral_constant<int, 24>&)
0099    {
0100       BOOST_MATH_STD_USING
0101       if(x <= 1)
0102       {
0103          // Maximum Deviation Found:                     3.090e-12
0104          // Expected Error Term : -3.053e-12
0105          // Maximum Relative Change in Control Points : 4.927e-02
0106          // Max Error found at float precision = Poly : 7.918347e-10
0107          BOOST_MATH_STATIC const T Y = 8.695471287e-02f;
0108          BOOST_MATH_STATIC const T P[] =
0109          {
0110             -3.621379531e-03f,
0111             7.131781976e-03f,
0112             -1.535278300e-05f
0113          };
0114          BOOST_MATH_STATIC const T Q[] =
0115          {
0116             1.000000000e+00f,
0117             -5.173102701e-02f,
0118             9.203530671e-04f
0119          };
0120 
0121          T a = x * x / 4;
0122          a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2;
0123 
0124          // Maximum Deviation Found:                     3.556e-08
0125          // Expected Error Term : -3.541e-08
0126          // Maximum Relative Change in Control Points : 8.203e-02
0127          BOOST_MATH_STATIC const T P2[] =
0128          {
0129             -3.079657469e-01f,
0130             -8.537108913e-02f,
0131             -4.640275408e-03f,
0132             -1.156442414e-04f
0133          };
0134 
0135          return tools::evaluate_polynomial(P2, T(x * x)) * x + 1 / x + log(x) * a;
0136       }
0137       else
0138       {
0139          // Maximum Deviation Found:                     3.369e-08
0140          // Expected Error Term : -3.227e-08
0141          // Maximum Relative Change in Control Points : 9.917e-02
0142          // Max Error found at float precision = Poly : 6.084411e-08
0143          BOOST_MATH_STATIC const T Y = 1.450342178f;
0144          BOOST_MATH_STATIC const T P[] =
0145          {
0146             -1.970280088e-01f,
0147             2.188747807e-02f,
0148             7.270394756e-01f,
0149             2.490678196e-01f
0150          };
0151          BOOST_MATH_STATIC const T Q[] =
0152          {
0153             1.000000000e+00f,
0154             2.274292882e+00f,
0155             9.904984851e-01f,
0156             4.585534549e-02f
0157          };
0158          if(x < tools::log_max_value<T>())
0159             return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
0160          else
0161          {
0162             T ex = exp(-x / 2);
0163             return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
0164          }
0165       }
0166    }
0167 
0168    template <typename T>
0169    BOOST_MATH_GPU_ENABLED T bessel_k1_imp(const T& x, const boost::math::integral_constant<int, 53>&)
0170    {
0171       BOOST_MATH_STD_USING
0172       if(x <= 1)
0173       {
0174          // Maximum Deviation Found:                     1.922e-17
0175          // Expected Error Term : 1.921e-17
0176          // Maximum Relative Change in Control Points : 5.287e-03
0177          // Max Error found at double precision = Poly : 2.004747e-17
0178          BOOST_MATH_STATIC const T Y = 8.69547128677368164e-02f;
0179          BOOST_MATH_STATIC const T P[] =
0180          {
0181             -3.62137953440350228e-03,
0182             7.11842087490330300e-03,
0183             1.00302560256614306e-05,
0184             1.77231085381040811e-06
0185          };
0186          BOOST_MATH_STATIC const T Q[] =
0187          {
0188             1.00000000000000000e+00,
0189             -4.80414794429043831e-02,
0190             9.85972641934416525e-04,
0191             -8.91196859397070326e-06
0192          };
0193 
0194          T a = x * x / 4;
0195          a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2;
0196 
0197          // Maximum Deviation Found:                     4.053e-17
0198          // Expected Error Term : -4.053e-17
0199          // Maximum Relative Change in Control Points : 3.103e-04
0200          // Max Error found at double precision = Poly : 1.246698e-16
0201 
0202          BOOST_MATH_STATIC const T P2[] =
0203          {
0204             -3.07965757829206184e-01,
0205             -7.80929703673074907e-02,
0206             -2.70619343754051620e-03,
0207             -2.49549522229072008e-05
0208          };
0209          BOOST_MATH_STATIC const T Q2[] = 
0210          {
0211             1.00000000000000000e+00,
0212             -2.36316836412163098e-02,
0213             2.64524577525962719e-04,
0214             -1.49749618004162787e-06
0215          };
0216 
0217          return tools::evaluate_rational(P2, Q2, T(x * x)) * x + 1 / x + log(x) * a;
0218       }
0219       else
0220       {
0221          // Maximum Deviation Found:                     8.883e-17
0222          // Expected Error Term : -1.641e-17
0223          // Maximum Relative Change in Control Points : 2.786e-01
0224          // Max Error found at double precision = Poly : 1.258798e-16
0225 
0226          BOOST_MATH_STATIC const T Y = 1.45034217834472656f;
0227          BOOST_MATH_STATIC const T P[] =
0228          {
0229             -1.97028041029226295e-01,
0230             -2.32408961548087617e+00,
0231             -7.98269784507699938e+00,
0232             -2.39968410774221632e+00,
0233             3.28314043780858713e+01,
0234             5.67713761158496058e+01,
0235             3.30907788466509823e+01,
0236             6.62582288933739787e+00,
0237             3.08851840645286691e-01
0238          };
0239          BOOST_MATH_STATIC const T Q[] =
0240          {
0241             1.00000000000000000e+00,
0242             1.41811409298826118e+01,
0243             7.35979466317556420e+01,
0244             1.77821793937080859e+02,
0245             2.11014501598705982e+02,
0246             1.19425262951064454e+02,
0247             2.88448064302447607e+01,
0248             2.27912927104139732e+00,
0249             2.50358186953478678e-02
0250          };
0251          if(x < tools::log_max_value<T>())
0252             return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
0253          else
0254          {
0255             T ex = exp(-x / 2);
0256             return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
0257          }
0258       }
0259    }
0260 
0261    template <typename T>
0262    BOOST_MATH_GPU_ENABLED T bessel_k1_imp(const T& x, const boost::math::integral_constant<int, 64>&)
0263    {
0264       BOOST_MATH_STD_USING
0265       if(x <= 1)
0266       {
0267          // Maximum Deviation Found:                     5.549e-23
0268          // Expected Error Term : -5.548e-23
0269          // Maximum Relative Change in Control Points : 2.002e-03
0270          // Max Error found at float80 precision = Poly : 9.352785e-22
0271          BOOST_MATH_STATIC const T Y = 8.695471286773681640625e-02f;
0272          BOOST_MATH_STATIC const T P[] =
0273          {
0274             BOOST_MATH_BIG_CONSTANT(T, 64, -3.621379534403483072861e-03),
0275             BOOST_MATH_BIG_CONSTANT(T, 64, 7.102135866103952705932e-03),
0276             BOOST_MATH_BIG_CONSTANT(T, 64, 4.167545240236717601167e-05),
0277             BOOST_MATH_BIG_CONSTANT(T, 64, 2.537484002571894870830e-06),
0278             BOOST_MATH_BIG_CONSTANT(T, 64, 6.603228256820000135990e-09)
0279          };
0280          BOOST_MATH_STATIC const T Q[] =
0281          {
0282             BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
0283             BOOST_MATH_BIG_CONSTANT(T, 64, -4.354457194045068370363e-02),
0284             BOOST_MATH_BIG_CONSTANT(T, 64, 8.709137201220209072820e-04),
0285             BOOST_MATH_BIG_CONSTANT(T, 64, -9.676151796359590545143e-06),
0286             BOOST_MATH_BIG_CONSTANT(T, 64, 5.162715192766245311659e-08)
0287          };
0288 
0289          T a = x * x / 4;
0290          a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2;
0291 
0292          // Maximum Deviation Found:                     1.995e-23
0293          // Expected Error Term : 1.995e-23
0294          // Maximum Relative Change in Control Points : 8.174e-04
0295          // Max Error found at float80 precision = Poly : 4.137325e-20
0296          BOOST_MATH_STATIC const T P2[] =
0297          {
0298             BOOST_MATH_BIG_CONSTANT(T, 64, -3.079657578292062244054e-01),
0299             BOOST_MATH_BIG_CONSTANT(T, 64, -7.963049154965966503231e-02),
0300             BOOST_MATH_BIG_CONSTANT(T, 64, -3.103277523735639924895e-03),
0301             BOOST_MATH_BIG_CONSTANT(T, 64, -4.023052834702215699504e-05),
0302             BOOST_MATH_BIG_CONSTANT(T, 64, -1.719459155018493821839e-07)
0303          };
0304          BOOST_MATH_STATIC const T Q2[] = 
0305          {
0306             BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
0307             BOOST_MATH_BIG_CONSTANT(T, 64, -1.863917670410152669768e-02),
0308             BOOST_MATH_BIG_CONSTANT(T, 64, 1.699367098849735298090e-04),
0309             BOOST_MATH_BIG_CONSTANT(T, 64, -9.309358790546076298429e-07),
0310             BOOST_MATH_BIG_CONSTANT(T, 64, 2.708893480271612711933e-09)
0311          };
0312 
0313          return tools::evaluate_rational(P2, Q2, T(x * x)) * x + 1 / x + log(x) * a;
0314       }
0315       else
0316       {
0317          // Maximum Deviation Found:                     9.785e-20
0318          // Expected Error Term : -3.302e-21
0319          // Maximum Relative Change in Control Points : 3.432e-01
0320          // Max Error found at float80 precision = Poly : 1.083755e-19
0321          BOOST_MATH_STATIC const T Y = 1.450342178344726562500e+00f;
0322          BOOST_MATH_STATIC const T P[] =
0323          {
0324             BOOST_MATH_BIG_CONSTANT(T, 64, -1.970280410292263112917e-01),
0325             BOOST_MATH_BIG_CONSTANT(T, 64, -4.058564803062959169322e+00),
0326             BOOST_MATH_BIG_CONSTANT(T, 64, -3.036658174194917777473e+01),
0327             BOOST_MATH_BIG_CONSTANT(T, 64, -9.576825392332820142173e+01),
0328             BOOST_MATH_BIG_CONSTANT(T, 64, -6.706969489248020941949e+01),
0329             BOOST_MATH_BIG_CONSTANT(T, 64, 3.264572499406168221382e+02),
0330             BOOST_MATH_BIG_CONSTANT(T, 64, 8.584972047303151034100e+02),
0331             BOOST_MATH_BIG_CONSTANT(T, 64, 8.422082733280017909550e+02),
0332             BOOST_MATH_BIG_CONSTANT(T, 64, 3.738005441471368178383e+02),
0333             BOOST_MATH_BIG_CONSTANT(T, 64, 7.016938390144121276609e+01),
0334             BOOST_MATH_BIG_CONSTANT(T, 64, 4.319614662598089438939e+00),
0335             BOOST_MATH_BIG_CONSTANT(T, 64, 3.710715864316521856193e-02)
0336          };
0337          BOOST_MATH_STATIC const T Q[] =
0338          {
0339             BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
0340             BOOST_MATH_BIG_CONSTANT(T, 64, 2.298433045824439052398e+01),
0341             BOOST_MATH_BIG_CONSTANT(T, 64, 2.082047745067709230037e+02),
0342             BOOST_MATH_BIG_CONSTANT(T, 64, 9.662367854250262046592e+02),
0343             BOOST_MATH_BIG_CONSTANT(T, 64, 2.504148628460454004686e+03),
0344             BOOST_MATH_BIG_CONSTANT(T, 64, 3.712730364911389908905e+03),
0345             BOOST_MATH_BIG_CONSTANT(T, 64, 3.108002081150068641112e+03),
0346             BOOST_MATH_BIG_CONSTANT(T, 64, 1.400149940532448553143e+03),
0347             BOOST_MATH_BIG_CONSTANT(T, 64, 3.083303048095846226299e+02),
0348             BOOST_MATH_BIG_CONSTANT(T, 64, 2.748706060530351833346e+01),
0349             BOOST_MATH_BIG_CONSTANT(T, 64, 6.321900849331506946977e-01),
0350          };
0351          if(x < tools::log_max_value<T>())
0352             return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
0353          else
0354          {
0355             T ex = exp(-x / 2);
0356             return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
0357          }
0358       }
0359    }
0360 
0361    template <typename T>
0362    BOOST_MATH_GPU_ENABLED T bessel_k1_imp(const T& x, const boost::math::integral_constant<int, 113>&)
0363    {
0364       BOOST_MATH_STD_USING
0365       if(x <= 1)
0366       {
0367          // Maximum Deviation Found:                     7.120e-35
0368          // Expected Error Term : -7.119e-35
0369          // Maximum Relative Change in Control Points : 1.207e-03
0370          // Max Error found at float128 precision = Poly : 7.143688e-35
0371          BOOST_MATH_STATIC const T Y = 8.695471286773681640625000000000000000e-02f;
0372          BOOST_MATH_STATIC const T P[] =
0373          {
0374             BOOST_MATH_BIG_CONSTANT(T, 113, -3.621379534403483072916666666666595475e-03),
0375             BOOST_MATH_BIG_CONSTANT(T, 113, 7.074117676930975433219826471336547627e-03),
0376             BOOST_MATH_BIG_CONSTANT(T, 113, 9.631337631362776369069668419033041661e-05),
0377             BOOST_MATH_BIG_CONSTANT(T, 113, 3.468935967870048731821071646104412775e-06),
0378             BOOST_MATH_BIG_CONSTANT(T, 113, 2.956705020559599861444492614737168261e-08),
0379             BOOST_MATH_BIG_CONSTANT(T, 113, 2.347140307321161346703214099534250263e-10),
0380             BOOST_MATH_BIG_CONSTANT(T, 113, 5.569608494081482873946791086435679661e-13)
0381          };
0382          BOOST_MATH_STATIC const T Q[] =
0383          {
0384             BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
0385             BOOST_MATH_BIG_CONSTANT(T, 113, -3.580768910152105375615558920428350204e-02),
0386             BOOST_MATH_BIG_CONSTANT(T, 113, 6.197467671701485365363068445534557369e-04),
0387             BOOST_MATH_BIG_CONSTANT(T, 113, -6.707466533308630411966030561446666237e-06),
0388             BOOST_MATH_BIG_CONSTANT(T, 113, 4.846687802282250112624373388491123527e-08),
0389             BOOST_MATH_BIG_CONSTANT(T, 113, -2.248493131151981569517383040323900343e-10),
0390             BOOST_MATH_BIG_CONSTANT(T, 113, 5.319279786372775264555728921709381080e-13)
0391          };
0392 
0393          T a = x * x / 4;
0394          a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2;
0395 
0396          // Maximum Deviation Found:                     4.473e-37
0397          // Expected Error Term : 4.473e-37
0398          // Maximum Relative Change in Control Points : 8.550e-04
0399          // Max Error found at float128 precision = Poly : 8.167701e-35
0400          BOOST_MATH_STATIC const T P2[] =
0401          {
0402             BOOST_MATH_BIG_CONSTANT(T, 113, -3.079657578292062244053600156878870690e-01),
0403             BOOST_MATH_BIG_CONSTANT(T, 113, -8.133183745732467770755578848987414875e-02),
0404             BOOST_MATH_BIG_CONSTANT(T, 113, -3.548968792764174773125420229299431951e-03),
0405             BOOST_MATH_BIG_CONSTANT(T, 113, -5.886125468718182876076972186152445490e-05),
0406             BOOST_MATH_BIG_CONSTANT(T, 113, -4.506712111733707245745396404449639865e-07),
0407             BOOST_MATH_BIG_CONSTANT(T, 113, -1.632502325880313239698965376754406011e-09),
0408             BOOST_MATH_BIG_CONSTANT(T, 113, -2.311973065898784812266544485665624227e-12)
0409          };
0410          BOOST_MATH_STATIC const T Q2[] = 
0411          {
0412             BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
0413             BOOST_MATH_BIG_CONSTANT(T, 113, -1.311471216733781016657962995723287450e-02),
0414             BOOST_MATH_BIG_CONSTANT(T, 113, 8.571876054797365417068164018709472969e-05),
0415             BOOST_MATH_BIG_CONSTANT(T, 113, -3.630181215268238731442496851497901293e-07),
0416             BOOST_MATH_BIG_CONSTANT(T, 113, 1.070176111227805048604885986867484807e-09),
0417             BOOST_MATH_BIG_CONSTANT(T, 113, -2.129046580769872602793220056461084761e-12),
0418             BOOST_MATH_BIG_CONSTANT(T, 113, 2.294906469421390890762001971790074432e-15)
0419          };
0420 
0421          return tools::evaluate_rational(P2, Q2, T(x * x)) * x + 1 / x + log(x) * a;
0422       }
0423       else if(x < 4)
0424       {
0425          // Max error in interpolated form: 5.307e-37
0426          // Max Error found at float128 precision = Poly: 7.087862e-35
0427          BOOST_MATH_STATIC const T Y = 1.5023040771484375f;
0428          BOOST_MATH_STATIC const T P[] =
0429          {
0430             BOOST_MATH_BIG_CONSTANT(T, 113, -2.489899398329369710528254347931380044e-01),
0431             BOOST_MATH_BIG_CONSTANT(T, 113, -6.819080211203854781858815596508456873e+00),
0432             BOOST_MATH_BIG_CONSTANT(T, 113, -7.599915699069767382647695624952723034e+01),
0433             BOOST_MATH_BIG_CONSTANT(T, 113, -4.450211910821295507926582231071300718e+02),
0434             BOOST_MATH_BIG_CONSTANT(T, 113, -1.451374687870925175794150513723956533e+03),
0435             BOOST_MATH_BIG_CONSTANT(T, 113, -2.405805746895098802803503988539098226e+03),
0436             BOOST_MATH_BIG_CONSTANT(T, 113, -5.638808326778389656403861103277220518e+02),
0437             BOOST_MATH_BIG_CONSTANT(T, 113, 5.513958744081268456191778822780865708e+03),
0438             BOOST_MATH_BIG_CONSTANT(T, 113, 1.121301640926540743072258116122834804e+04),
0439             BOOST_MATH_BIG_CONSTANT(T, 113, 1.080094900175649541266613109971296190e+04),
0440             BOOST_MATH_BIG_CONSTANT(T, 113, 5.896531083639613332407534434915552429e+03),
0441             BOOST_MATH_BIG_CONSTANT(T, 113, 1.856602122319645694042555107114028437e+03),
0442             BOOST_MATH_BIG_CONSTANT(T, 113, 3.237121918853145421414003823957537419e+02),
0443             BOOST_MATH_BIG_CONSTANT(T, 113, 2.842072954561323076230238664623893504e+01),
0444             BOOST_MATH_BIG_CONSTANT(T, 113, 1.039705646510167437971862966128055524e+00),
0445             BOOST_MATH_BIG_CONSTANT(T, 113, 1.008418100718254816100425022904039530e-02)
0446          };
0447          BOOST_MATH_STATIC const T Q[] =
0448          {
0449             BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
0450             BOOST_MATH_BIG_CONSTANT(T, 113, 2.927456835239137986889227412815459529e+01),
0451             BOOST_MATH_BIG_CONSTANT(T, 113, 3.598985593265577043711382994516531273e+02),
0452             BOOST_MATH_BIG_CONSTANT(T, 113, 2.449897377085510281395819892689690579e+03),
0453             BOOST_MATH_BIG_CONSTANT(T, 113, 1.025555887684561913263090023158085327e+04),
0454             BOOST_MATH_BIG_CONSTANT(T, 113, 2.774140447181062463181892531100679195e+04),
0455             BOOST_MATH_BIG_CONSTANT(T, 113, 4.962055507843204417243602332246120418e+04),
0456             BOOST_MATH_BIG_CONSTANT(T, 113, 5.908269326976180183216954452196772931e+04),
0457             BOOST_MATH_BIG_CONSTANT(T, 113, 4.655160454422016855911700790722577942e+04),
0458             BOOST_MATH_BIG_CONSTANT(T, 113, 2.383586885019548163464418964577684608e+04),
0459             BOOST_MATH_BIG_CONSTANT(T, 113, 7.679920375586960324298491662159976419e+03),
0460             BOOST_MATH_BIG_CONSTANT(T, 113, 1.478586421028842906987799049804565008e+03),
0461             BOOST_MATH_BIG_CONSTANT(T, 113, 1.565384974896746094224942654383537090e+02),
0462             BOOST_MATH_BIG_CONSTANT(T, 113, 7.902617937084010911005732488607114511e+00),
0463             BOOST_MATH_BIG_CONSTANT(T, 113, 1.429293010387921526110949911029094926e-01),
0464             BOOST_MATH_BIG_CONSTANT(T, 113, 3.880342607911083143560111853491047663e-04)
0465          };
0466          return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
0467       }
0468       else
0469       {
0470          // Maximum Deviation Found:                     4.359e-37
0471          // Expected Error Term : -6.565e-40
0472          // Maximum Relative Change in Control Points : 1.880e-01
0473          // Max Error found at float128 precision = Poly : 2.943572e-35
0474          BOOST_MATH_STATIC const T Y = 1.308816909790039062500000000000000000f;
0475          BOOST_MATH_STATIC const T P[] =
0476          {
0477             BOOST_MATH_BIG_CONSTANT(T, 113, -5.550277247453881129211735759447737350e-02),
0478             BOOST_MATH_BIG_CONSTANT(T, 113, -3.485883080219574328217554864956175929e+00),
0479             BOOST_MATH_BIG_CONSTANT(T, 113, -8.903760658131484239300875153154881958e+01),
0480             BOOST_MATH_BIG_CONSTANT(T, 113, -1.144813672213626237418235110712293337e+03),
0481             BOOST_MATH_BIG_CONSTANT(T, 113, -6.498400501156131446691826557494158173e+03),
0482             BOOST_MATH_BIG_CONSTANT(T, 113, 1.573531831870363502604119835922166116e+04),
0483             BOOST_MATH_BIG_CONSTANT(T, 113, 5.417416550054632009958262596048841154e+05),
0484             BOOST_MATH_BIG_CONSTANT(T, 113, 4.271266450613557412825896604269130661e+06),
0485             BOOST_MATH_BIG_CONSTANT(T, 113, 1.898386013314389952534433455681107783e+07),
0486             BOOST_MATH_BIG_CONSTANT(T, 113, 5.353798784656436259250791761023512750e+07),
0487             BOOST_MATH_BIG_CONSTANT(T, 113, 9.839619195427352438957774052763490067e+07),
0488             BOOST_MATH_BIG_CONSTANT(T, 113, 1.169246368651532232388152442538005637e+08),
0489             BOOST_MATH_BIG_CONSTANT(T, 113, 8.696368884166831199967845883371116431e+07),
0490             BOOST_MATH_BIG_CONSTANT(T, 113, 3.810226630422736458064005843327500169e+07),
0491             BOOST_MATH_BIG_CONSTANT(T, 113, 8.854996610560406127438950635716757614e+06),
0492             BOOST_MATH_BIG_CONSTANT(T, 113, 8.981057433937398731355768088809437625e+05),
0493             BOOST_MATH_BIG_CONSTANT(T, 113, 2.519440069856232098711793483639792952e+04)
0494          };
0495          BOOST_MATH_STATIC const T Q[] =
0496          {
0497             BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
0498             BOOST_MATH_BIG_CONSTANT(T, 113, 7.127348248283623146544565916604103560e+01),
0499             BOOST_MATH_BIG_CONSTANT(T, 113, 2.205092684176906740104488180754982065e+03),
0500             BOOST_MATH_BIG_CONSTANT(T, 113, 3.911249195069050636298346469740075758e+04),
0501             BOOST_MATH_BIG_CONSTANT(T, 113, 4.426103406579046249654548481377792614e+05),
0502             BOOST_MATH_BIG_CONSTANT(T, 113, 3.365861555422488771286500241966208541e+06),
0503             BOOST_MATH_BIG_CONSTANT(T, 113, 1.765377714160383676864913709252529840e+07),
0504             BOOST_MATH_BIG_CONSTANT(T, 113, 6.453822726931857253365138260720815246e+07),
0505             BOOST_MATH_BIG_CONSTANT(T, 113, 1.643207885048369990391975749439783892e+08),
0506             BOOST_MATH_BIG_CONSTANT(T, 113, 2.882540678243694621895816336640877878e+08),
0507             BOOST_MATH_BIG_CONSTANT(T, 113, 3.410120808992380266174106812005338148e+08),
0508             BOOST_MATH_BIG_CONSTANT(T, 113, 2.628138016559335882019310900426773027e+08),
0509             BOOST_MATH_BIG_CONSTANT(T, 113, 1.250794693811010646965360198541047961e+08),
0510             BOOST_MATH_BIG_CONSTANT(T, 113, 3.378723408195485594610593014072950078e+07),
0511             BOOST_MATH_BIG_CONSTANT(T, 113, 4.488253856312453816451380319061865560e+06),
0512             BOOST_MATH_BIG_CONSTANT(T, 113, 2.202167197882689873967723350537104582e+05),
0513             BOOST_MATH_BIG_CONSTANT(T, 113, 1.673233230356966539460728211412989843e+03)
0514          };
0515          if(x < tools::log_max_value<T>())
0516             return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
0517          else
0518          {
0519             T ex = exp(-x / 2);
0520             return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
0521          }
0522       }
0523     }
0524 
0525     template <typename T>
0526     BOOST_MATH_GPU_ENABLED T bessel_k1_imp(const T& x, const boost::math::integral_constant<int, 0>&)
0527     {
0528        if(boost::math::tools::digits<T>() <= 24)
0529           return bessel_k1_imp(x, boost::math::integral_constant<int, 24>());
0530        else if(boost::math::tools::digits<T>() <= 53)
0531           return bessel_k1_imp(x, boost::math::integral_constant<int, 53>());
0532        else if(boost::math::tools::digits<T>() <= 64)
0533           return bessel_k1_imp(x, boost::math::integral_constant<int, 64>());
0534        else if(boost::math::tools::digits<T>() <= 113)
0535           return bessel_k1_imp(x, boost::math::integral_constant<int, 113>());
0536        BOOST_MATH_ASSERT(0);
0537        return 0;
0538     }
0539 
0540    template <typename T>
0541    inline BOOST_MATH_GPU_ENABLED T bessel_k1(const T& x)
0542    {
0543       typedef boost::math::integral_constant<int,
0544          ((boost::math::numeric_limits<T>::digits == 0) || (boost::math::numeric_limits<T>::radix != 2)) ?
0545          0 :
0546          boost::math::numeric_limits<T>::digits <= 24 ?
0547          24 :
0548          boost::math::numeric_limits<T>::digits <= 53 ?
0549          53 :
0550          boost::math::numeric_limits<T>::digits <= 64 ?
0551          64 :
0552          boost::math::numeric_limits<T>::digits <= 113 ?
0553          113 : -1
0554       > tag_type;
0555 
0556       bessel_k1_initializer<T, tag_type>::force_instantiate();
0557       return bessel_k1_imp(x, tag_type());
0558    }
0559 
0560 }}} // namespaces
0561 
0562 #ifdef _MSC_VER
0563 #pragma warning(pop)
0564 #endif
0565 
0566 #endif // BOOST_MATH_BESSEL_K1_HPP
0567
[ 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!