Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:39:59 

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_I0_HPP
0008 #define BOOST_MATH_BESSEL_I0_HPP
0009 
0010 #ifdef _MSC_VER
0011 #pragma once
0012 #endif
0013 
0014 #include <boost/math/tools/rational.hpp>
0015 #include <boost/math/tools/big_constant.hpp>
0016 #include <boost/math/tools/assert.hpp>
0017 
0018 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
0019 //
0020 // This is the only way we can avoid
0021 // warning: non-standard suffix on floating constant [-Wpedantic]
0022 // when building with -Wall -pedantic.  Neither __extension__
0023 // nor #pragma diagnostic ignored work :(
0024 //
0025 #pragma GCC system_header
0026 #endif
0027 
0028 // Modified Bessel function of the first kind of order zero
0029 // we use the approximating forms derived in:
0030 // "Rational Approximations for the Modified Bessel Function of the First Kind - I0(x) for Computations with Double Precision"
0031 // by Pavel Holoborodko, 
0032 // see http://www.advanpix.com/2015/11/11/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i0-computations-double-precision
0033 // The actual coefficients used are our own, and extend Pavel's work to precision's other than double.
0034 
0035 namespace boost { namespace math { namespace detail{
0036 
0037 template <typename T>
0038 T bessel_i0(const T& x);
0039 
0040 template <class T, class tag>
0041 struct bessel_i0_initializer
0042 {
0043    struct init
0044    {
0045       init()
0046       {
0047          do_init(tag());
0048       }
0049       static void do_init(const std::integral_constant<int, 64>&)
0050       {
0051          bessel_i0(T(1));
0052          bessel_i0(T(8));
0053          bessel_i0(T(12));
0054          bessel_i0(T(40));
0055          bessel_i0(T(101));
0056       }
0057       static void do_init(const std::integral_constant<int, 113>&)
0058       {
0059          bessel_i0(T(1));
0060          bessel_i0(T(10));
0061          bessel_i0(T(20));
0062          bessel_i0(T(40));
0063          bessel_i0(T(101));
0064       }
0065       template <class U>
0066       static void do_init(const U&) {}
0067       void force_instantiate()const {}
0068    };
0069    static const init initializer;
0070    static void force_instantiate()
0071    {
0072       initializer.force_instantiate();
0073    }
0074 };
0075 
0076 template <class T, class tag>
0077 const typename bessel_i0_initializer<T, tag>::init bessel_i0_initializer<T, tag>::initializer;
0078 
0079 template <typename T, int N>
0080 T bessel_i0_imp(const T&, const std::integral_constant<int, N>&)
0081 {
0082    BOOST_MATH_ASSERT(0);
0083    return 0;
0084 }
0085 
0086 template <typename T>
0087 T bessel_i0_imp(const T& x, const std::integral_constant<int, 24>&)
0088 {
0089    BOOST_MATH_STD_USING
0090    if(x < 7.75)
0091    {
0092       // Max error in interpolated form: 3.929e-08
0093       // Max Error found at float precision = Poly: 1.991226e-07
0094       static const float P[] = {
0095          1.00000003928615375e+00f,
0096          2.49999576572179639e-01f,
0097          2.77785268558399407e-02f,
0098          1.73560257755821695e-03f,
0099          6.96166518788906424e-05f,
0100          1.89645733877137904e-06f,
0101          4.29455004657565361e-08f,
0102          3.90565476357034480e-10f,
0103          1.48095934745267240e-11f
0104       };
0105       T a = x * x / 4;
0106       return a * boost::math::tools::evaluate_polynomial(P, a) + 1;
0107    }
0108    else if(x < 50)
0109    {
0110       // Max error in interpolated form: 5.195e-08
0111       // Max Error found at float precision = Poly: 8.502534e-08
0112       static const float P[] = {
0113          3.98942651588301770e-01f,
0114          4.98327234176892844e-02f,
0115          2.91866904423115499e-02f,
0116          1.35614940793742178e-02f,
0117          1.31409251787866793e-01f
0118       };
0119       return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0120    }
0121    else
0122    {
0123       // Max error in interpolated form: 1.782e-09
0124       // Max Error found at float precision = Poly: 6.473568e-08
0125       static const float P[] = {
0126          3.98942391532752700e-01f,
0127          4.98455950638200020e-02f,
0128          2.94835666900682535e-02f
0129       };
0130       T ex = exp(x / 2);
0131       T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0132       result *= ex;
0133       return result;
0134    }
0135 }
0136 
0137 template <typename T>
0138 T bessel_i0_imp(const T& x, const std::integral_constant<int, 53>&)
0139 {
0140    BOOST_MATH_STD_USING
0141    if(x < 7.75)
0142    {
0143       // Bessel I0 over[10 ^ -16, 7.75]
0144       // Max error in interpolated form : 3.042e-18
0145       // Max Error found at double precision = Poly : 5.106609e-16 Cheb : 5.239199e-16
0146       static const double P[] = {
0147          1.00000000000000000e+00,
0148          2.49999999999999909e-01,
0149          2.77777777777782257e-02,
0150          1.73611111111023792e-03,
0151          6.94444444453352521e-05,
0152          1.92901234513219920e-06,
0153          3.93675991102510739e-08,
0154          6.15118672704439289e-10,
0155          7.59407002058973446e-12,
0156          7.59389793369836367e-14,
0157          6.27767773636292611e-16,
0158          4.34709704153272287e-18,
0159          2.63417742690109154e-20,
0160          1.13943037744822825e-22,
0161          9.07926920085624812e-25
0162       };
0163       T a = x * x / 4;
0164       return a * boost::math::tools::evaluate_polynomial(P, a) + 1;
0165    }
0166    else if(x < 500)
0167    {
0168       // Max error in interpolated form : 1.685e-16
0169       // Max Error found at double precision = Poly : 2.575063e-16 Cheb : 2.247615e+00
0170       static const double P[] = {
0171          3.98942280401425088e-01,
0172          4.98677850604961985e-02,
0173          2.80506233928312623e-02,
0174          2.92211225166047873e-02,
0175          4.44207299493659561e-02,
0176          1.30970574605856719e-01,
0177          -3.35052280231727022e+00,
0178          2.33025711583514727e+02,
0179          -1.13366350697172355e+04,
0180          4.24057674317867331e+05,
0181          -1.23157028595698731e+07,
0182          2.80231938155267516e+08,
0183          -5.01883999713777929e+09,
0184          7.08029243015109113e+10,
0185          -7.84261082124811106e+11,
0186          6.76825737854096565e+12,
0187          -4.49034849696138065e+13,
0188          2.24155239966958995e+14,
0189          -8.13426467865659318e+14,
0190          2.02391097391687777e+15,
0191          -3.08675715295370878e+15,
0192          2.17587543863819074e+15
0193       };
0194       return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0195    }
0196    else
0197    {
0198       // Max error in interpolated form : 2.437e-18
0199       // Max Error found at double precision = Poly : 1.216719e-16
0200       static const double P[] = {
0201          3.98942280401432905e-01,
0202          4.98677850491434560e-02,
0203          2.80506308916506102e-02,
0204          2.92179096853915176e-02,
0205          4.53371208762579442e-02
0206       };
0207       T ex = exp(x / 2);
0208       T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0209       result *= ex;
0210       return result;
0211    }
0212 }
0213 
0214 template <typename T>
0215 T bessel_i0_imp(const T& x, const std::integral_constant<int, 64>&)
0216 {
0217    BOOST_MATH_STD_USING
0218    if(x < 7.75)
0219    {
0220       // Bessel I0 over[10 ^ -16, 7.75]
0221       // Max error in interpolated form : 3.899e-20
0222       // Max Error found at float80 precision = Poly : 1.770840e-19
0223       static const T P[] = {
0224          BOOST_MATH_BIG_CONSTANT(T, 64, 9.99999999999999999961011629e-01),
0225          BOOST_MATH_BIG_CONSTANT(T, 64, 2.50000000000000001321873912e-01),
0226          BOOST_MATH_BIG_CONSTANT(T, 64, 2.77777777777777703400424216e-02),
0227          BOOST_MATH_BIG_CONSTANT(T, 64, 1.73611111111112764793802701e-03),
0228          BOOST_MATH_BIG_CONSTANT(T, 64, 6.94444444444251461247253525e-05),
0229          BOOST_MATH_BIG_CONSTANT(T, 64, 1.92901234569262206386118739e-06),
0230          BOOST_MATH_BIG_CONSTANT(T, 64, 3.93675988851131457141005209e-08),
0231          BOOST_MATH_BIG_CONSTANT(T, 64, 6.15118734688297476454205352e-10),
0232          BOOST_MATH_BIG_CONSTANT(T, 64, 7.59405797058091016449222685e-12),
0233          BOOST_MATH_BIG_CONSTANT(T, 64, 7.59406599631719800679835140e-14),
0234          BOOST_MATH_BIG_CONSTANT(T, 64, 6.27598961062070013516660425e-16),
0235          BOOST_MATH_BIG_CONSTANT(T, 64, 4.35920318970387940278362992e-18),
0236          BOOST_MATH_BIG_CONSTANT(T, 64, 2.57372492687715452949437981e-20),
0237          BOOST_MATH_BIG_CONSTANT(T, 64, 1.33908663475949906992942204e-22),
0238          BOOST_MATH_BIG_CONSTANT(T, 64, 5.15976668870980234582896010e-25),
0239          BOOST_MATH_BIG_CONSTANT(T, 64, 3.46240478946376069211156548e-27)
0240       };
0241       T a = x * x / 4;
0242       return a * boost::math::tools::evaluate_polynomial(P, a) + 1;
0243    }
0244    else if(x < 10)
0245    {
0246       // Maximum Deviation Found:                     6.906e-21
0247       // Expected Error Term : -6.903e-21
0248       // Maximum Relative Change in Control Points : 1.631e-04
0249       // Max Error found at float80 precision = Poly : 7.811948e-21
0250       static const T Y = 4.051098823547363281250e-01f;
0251       static const T P[] = {
0252          BOOST_MATH_BIG_CONSTANT(T, 64, -6.158081780620616479492e-03),
0253          BOOST_MATH_BIG_CONSTANT(T, 64, 4.883635969834048766148e-02),
0254          BOOST_MATH_BIG_CONSTANT(T, 64, 7.892782002476195771920e-02),
0255          BOOST_MATH_BIG_CONSTANT(T, 64, -1.478784996478070170327e+00),
0256          BOOST_MATH_BIG_CONSTANT(T, 64, 2.988611837308006851257e+01),
0257          BOOST_MATH_BIG_CONSTANT(T, 64, -4.140133766747436806179e+02),
0258          BOOST_MATH_BIG_CONSTANT(T, 64, 4.117316447921276453271e+03),
0259          BOOST_MATH_BIG_CONSTANT(T, 64, -2.942353667455141676001e+04),
0260          BOOST_MATH_BIG_CONSTANT(T, 64, 1.493482682461387081534e+05),
0261          BOOST_MATH_BIG_CONSTANT(T, 64, -5.228100538921466124653e+05),
0262          BOOST_MATH_BIG_CONSTANT(T, 64, 1.195279248600467989454e+06),
0263          BOOST_MATH_BIG_CONSTANT(T, 64, -1.601530760654337045917e+06),
0264          BOOST_MATH_BIG_CONSTANT(T, 64, 9.504921137873298402679e+05)
0265       };
0266       return exp(x) * (boost::math::tools::evaluate_polynomial(P, T(1 / x)) + Y) / sqrt(x);
0267    }
0268    else if(x < 15)
0269    {
0270       // Maximum Deviation Found:                     4.083e-21
0271       // Expected Error Term : -4.025e-21
0272       // Maximum Relative Change in Control Points : 1.304e-03
0273       // Max Error found at float80 precision = Poly : 2.303527e-20
0274       static const T Y = 4.033188819885253906250e-01f;
0275       static const T P[] = {
0276          BOOST_MATH_BIG_CONSTANT(T, 64, -4.376373876116109401062e-03),
0277          BOOST_MATH_BIG_CONSTANT(T, 64, 4.982899138682911273321e-02),
0278          BOOST_MATH_BIG_CONSTANT(T, 64, 3.109477529533515397644e-02),
0279          BOOST_MATH_BIG_CONSTANT(T, 64, -1.163760580110576407673e-01),
0280          BOOST_MATH_BIG_CONSTANT(T, 64, 4.776501832837367371883e+00),
0281          BOOST_MATH_BIG_CONSTANT(T, 64, -1.101478069227776656318e+02),
0282          BOOST_MATH_BIG_CONSTANT(T, 64, 1.892071912448960299773e+03),
0283          BOOST_MATH_BIG_CONSTANT(T, 64, -2.417739279982328117483e+04),
0284          BOOST_MATH_BIG_CONSTANT(T, 64, 2.296963447724067390552e+05),
0285          BOOST_MATH_BIG_CONSTANT(T, 64, -1.598589306710589358747e+06),
0286          BOOST_MATH_BIG_CONSTANT(T, 64, 7.903662411851774878322e+06),
0287          BOOST_MATH_BIG_CONSTANT(T, 64, -2.622677059040339516093e+07),
0288          BOOST_MATH_BIG_CONSTANT(T, 64, 5.227776578828667629347e+07),
0289          BOOST_MATH_BIG_CONSTANT(T, 64, -4.727797957441040896878e+07)
0290       };
0291       return exp(x) * (boost::math::tools::evaluate_polynomial(P, T(1 / x)) + Y) / sqrt(x);
0292    }
0293    else if(x < 50)
0294    {
0295       // Max error in interpolated form: 1.035e-21
0296       // Max Error found at float80 precision = Poly: 1.885872e-21
0297       static const T Y = 4.011702537536621093750e-01f;
0298       static const T P[] = {
0299          BOOST_MATH_BIG_CONSTANT(T, 64, -2.227973351806078464328e-03),
0300          BOOST_MATH_BIG_CONSTANT(T, 64, 4.986778486088017419036e-02),
0301          BOOST_MATH_BIG_CONSTANT(T, 64, 2.805066823812285310011e-02),
0302          BOOST_MATH_BIG_CONSTANT(T, 64, 2.921443721160964964623e-02),
0303          BOOST_MATH_BIG_CONSTANT(T, 64, 4.517504941996594744052e-02),
0304          BOOST_MATH_BIG_CONSTANT(T, 64, 6.316922639868793684401e-02),
0305          BOOST_MATH_BIG_CONSTANT(T, 64, 1.535891099168810015433e+00),
0306          BOOST_MATH_BIG_CONSTANT(T, 64, -4.706078229522448308087e+01),
0307          BOOST_MATH_BIG_CONSTANT(T, 64, 1.351015763079160914632e+03),
0308          BOOST_MATH_BIG_CONSTANT(T, 64, -2.948809013999277355098e+04),
0309          BOOST_MATH_BIG_CONSTANT(T, 64, 4.967598958582595361757e+05),
0310          BOOST_MATH_BIG_CONSTANT(T, 64, -6.346924657995383019558e+06),
0311          BOOST_MATH_BIG_CONSTANT(T, 64, 5.998794574259956613472e+07),
0312          BOOST_MATH_BIG_CONSTANT(T, 64, -4.016371355801690142095e+08),
0313          BOOST_MATH_BIG_CONSTANT(T, 64, 1.768791455631826490838e+09),
0314          BOOST_MATH_BIG_CONSTANT(T, 64, -4.441995678177349895640e+09),
0315          BOOST_MATH_BIG_CONSTANT(T, 64, 4.482292669974971387738e+09)
0316       };
0317       return exp(x) * (boost::math::tools::evaluate_polynomial(P, T(1 / x)) + Y) / sqrt(x);
0318    }
0319    else
0320    {
0321       // Bessel I0 over[50, INF]
0322       // Max error in interpolated form : 5.587e-20
0323       // Max Error found at float80 precision = Poly : 8.776852e-20
0324       static const T P[] = {
0325          BOOST_MATH_BIG_CONSTANT(T, 64, 3.98942280401432677955074061e-01),
0326          BOOST_MATH_BIG_CONSTANT(T, 64, 4.98677850501789875615574058e-02),
0327          BOOST_MATH_BIG_CONSTANT(T, 64, 2.80506290908675604202206833e-02),
0328          BOOST_MATH_BIG_CONSTANT(T, 64, 2.92194052159035901631494784e-02),
0329          BOOST_MATH_BIG_CONSTANT(T, 64, 4.47422430732256364094681137e-02),
0330          BOOST_MATH_BIG_CONSTANT(T, 64, 9.05971614435738691235525172e-02),
0331          BOOST_MATH_BIG_CONSTANT(T, 64, 2.29180522595459823234266708e-01),
0332          BOOST_MATH_BIG_CONSTANT(T, 64, 6.15122547776140254569073131e-01),
0333          BOOST_MATH_BIG_CONSTANT(T, 64, 7.48491812136365376477357324e+00),
0334          BOOST_MATH_BIG_CONSTANT(T, 64, -2.45569740166506688169730713e+02),
0335          BOOST_MATH_BIG_CONSTANT(T, 64, 9.66857566379480730407063170e+03),
0336          BOOST_MATH_BIG_CONSTANT(T, 64, -2.71924083955641197750323901e+05),
0337          BOOST_MATH_BIG_CONSTANT(T, 64, 5.74276685704579268845870586e+06),
0338          BOOST_MATH_BIG_CONSTANT(T, 64, -8.89753803265734681907148778e+07),
0339          BOOST_MATH_BIG_CONSTANT(T, 64, 9.82590905134996782086242180e+08),
0340          BOOST_MATH_BIG_CONSTANT(T, 64, -7.30623197145529889358596301e+09),
0341          BOOST_MATH_BIG_CONSTANT(T, 64, 3.27310000726207055200805893e+10),
0342          BOOST_MATH_BIG_CONSTANT(T, 64, -6.64365417189215599168817064e+10)
0343       };
0344       T ex = exp(x / 2);
0345       T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0346       result *= ex;
0347       return result;
0348    }
0349 }
0350 
0351 template <typename T>
0352 T bessel_i0_imp(const T& x, const std::integral_constant<int, 113>&)
0353 {
0354    BOOST_MATH_STD_USING
0355    if(x < 7.75)
0356    {
0357       // Bessel I0 over[10 ^ -34, 7.75]
0358       // Max error in interpolated form : 1.274e-34
0359       // Max Error found at float128 precision = Poly : 3.096091e-34
0360       static const T P[] = {
0361          BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000001273856e+00),
0362          BOOST_MATH_BIG_CONSTANT(T, 113, 2.4999999999999999999999999999999107477496e-01),
0363          BOOST_MATH_BIG_CONSTANT(T, 113, 2.7777777777777777777777777777881795230918e-02),
0364          BOOST_MATH_BIG_CONSTANT(T, 113, 1.7361111111111111111111111106290091648808e-03),
0365          BOOST_MATH_BIG_CONSTANT(T, 113, 6.9444444444444444444444445629960334523101e-05),
0366          BOOST_MATH_BIG_CONSTANT(T, 113, 1.9290123456790123456790105563456483249753e-06),
0367          BOOST_MATH_BIG_CONSTANT(T, 113, 3.9367598891408415217940836339080514004844e-08),
0368          BOOST_MATH_BIG_CONSTANT(T, 113, 6.1511873267825648777900014857992724731476e-10),
0369          BOOST_MATH_BIG_CONSTANT(T, 113, 7.5940584281266233066162999610732449709209e-12),
0370          BOOST_MATH_BIG_CONSTANT(T, 113, 7.5940584281266232783124723601470051895304e-14),
0371          BOOST_MATH_BIG_CONSTANT(T, 113, 6.2760813455591936763439337059117957836078e-16),
0372          BOOST_MATH_BIG_CONSTANT(T, 113, 4.3583898233049738471136482147779094353096e-18),
0373          BOOST_MATH_BIG_CONSTANT(T, 113, 2.5789288895299965395422423848480340736308e-20),
0374          BOOST_MATH_BIG_CONSTANT(T, 113, 1.3157800456718804437960453545507623434606e-22),
0375          BOOST_MATH_BIG_CONSTANT(T, 113, 5.8479113149412360748032684260932041506493e-25),
0376          BOOST_MATH_BIG_CONSTANT(T, 113, 2.2843403488398038539283241944594140493394e-27),
0377          BOOST_MATH_BIG_CONSTANT(T, 113, 7.9042925594356556196790242908697582021825e-30),
0378          BOOST_MATH_BIG_CONSTANT(T, 113, 2.4395919891312152120710245152115597111101e-32),
0379          BOOST_MATH_BIG_CONSTANT(T, 113, 6.7580986145276689333214547502373003196707e-35),
0380          BOOST_MATH_BIG_CONSTANT(T, 113, 1.6886514018062348877723837017198859723889e-37),
0381          BOOST_MATH_BIG_CONSTANT(T, 113, 3.8540558465757554512570197585002702777999e-40),
0382          BOOST_MATH_BIG_CONSTANT(T, 113, 7.4684706070226893763741850944911705726436e-43),
0383          BOOST_MATH_BIG_CONSTANT(T, 113, 2.0210715309399646335858150349406935414314e-45)
0384       };
0385       T a = x * x / 4;
0386       return a * boost::math::tools::evaluate_polynomial(P, a) + 1;
0387    }
0388    else if(x < 15)
0389    {
0390       // Bessel I0 over[7.75, 15]
0391       // Max error in interpolated form : 7.534e-35
0392       // Max Error found at float128 precision = Poly : 6.123912e-34
0393       static const T P[] = {
0394          BOOST_MATH_BIG_CONSTANT(T, 113, 9.9999999999999999992388573069504617493518e-01),
0395          BOOST_MATH_BIG_CONSTANT(T, 113, 2.5000000000000000007304739268173096975340e-01),
0396          BOOST_MATH_BIG_CONSTANT(T, 113, 2.7777777777777777744261405400543564492074e-02),
0397          BOOST_MATH_BIG_CONSTANT(T, 113, 1.7361111111111111209006987259719750726867e-03),
0398          BOOST_MATH_BIG_CONSTANT(T, 113, 6.9444444444444442399703186871329381908321e-05),
0399          BOOST_MATH_BIG_CONSTANT(T, 113, 1.9290123456790126709286741580242189785431e-06),
0400          BOOST_MATH_BIG_CONSTANT(T, 113, 3.9367598891408374246503061422528266924389e-08),
0401          BOOST_MATH_BIG_CONSTANT(T, 113, 6.1511873267826068395343047827801353170966e-10),
0402          BOOST_MATH_BIG_CONSTANT(T, 113, 7.5940584281262673459688011737168286944521e-12),
0403          BOOST_MATH_BIG_CONSTANT(T, 113, 7.5940584281291583769928563167645746144508e-14),
0404          BOOST_MATH_BIG_CONSTANT(T, 113, 6.2760813455438840231126529638737436950274e-16),
0405          BOOST_MATH_BIG_CONSTANT(T, 113, 4.3583898233839583885132809584770578894948e-18),
0406          BOOST_MATH_BIG_CONSTANT(T, 113, 2.5789288891798658971960571838369339742994e-20),
0407          BOOST_MATH_BIG_CONSTANT(T, 113, 1.3157800470129311623308216856009970266088e-22),
0408          BOOST_MATH_BIG_CONSTANT(T, 113, 5.8479112701534604520063520412207286692581e-25),
0409          BOOST_MATH_BIG_CONSTANT(T, 113, 2.2843404822552330714586265081801727491890e-27),
0410          BOOST_MATH_BIG_CONSTANT(T, 113, 7.9042888166225242675881424439818162458179e-30),
0411          BOOST_MATH_BIG_CONSTANT(T, 113, 2.4396027771820721384198604723320045236973e-32),
0412          BOOST_MATH_BIG_CONSTANT(T, 113, 6.7577659910606076328136207973456511895030e-35),
0413          BOOST_MATH_BIG_CONSTANT(T, 113, 1.6896548123724136624716224328803899914646e-37),
0414          BOOST_MATH_BIG_CONSTANT(T, 113, 3.8285850162160539150210466453921758781984e-40),
0415          BOOST_MATH_BIG_CONSTANT(T, 113, 7.9419071894227736216423562425429524883562e-43),
0416          BOOST_MATH_BIG_CONSTANT(T, 113, 1.4720374049498608905571855665134539425038e-45),
0417          BOOST_MATH_BIG_CONSTANT(T, 113, 2.7763533278527958112907118930154738930378e-48),
0418          BOOST_MATH_BIG_CONSTANT(T, 113, 3.1213839473168678646697528580511702663617e-51),
0419          BOOST_MATH_BIG_CONSTANT(T, 113, 1.0648035313124146852372607519737686740964e-53),
0420          -BOOST_MATH_BIG_CONSTANT(T, 113, 5.1255595184052024349371058585102280860878e-57),
0421          BOOST_MATH_BIG_CONSTANT(T, 113, 3.4652470895944157957727948355523715335882e-59)
0422       };
0423       T a = x * x / 4;
0424       return a * boost::math::tools::evaluate_polynomial(P, a) + 1;
0425    }
0426    else if(x < 30)
0427    {
0428       // Max error in interpolated form : 1.808e-34
0429       // Max Error found at float128 precision = Poly : 2.399403e-34
0430       static const T P[] = {
0431          BOOST_MATH_BIG_CONSTANT(T, 113, 3.9894228040870793650581242239624530714032e-01),
0432          BOOST_MATH_BIG_CONSTANT(T, 113, 4.9867780576714783790784348982178607842250e-02),
0433          BOOST_MATH_BIG_CONSTANT(T, 113, 2.8051948347934462928487999569249907599510e-02),
0434          BOOST_MATH_BIG_CONSTANT(T, 113, 2.8971143420388958551176254291160976367263e-02),
0435          BOOST_MATH_BIG_CONSTANT(T, 113, 7.8197359701715582763961322341827341098897e-02),
0436          BOOST_MATH_BIG_CONSTANT(T, 113, -3.3430484862908317377522273217643346601271e+00),
0437          BOOST_MATH_BIG_CONSTANT(T, 113, 2.7884507603213662610604413960838990199224e+02),
0438          BOOST_MATH_BIG_CONSTANT(T, 113, -1.8304926482356755790062999202373909300514e+04),
0439          BOOST_MATH_BIG_CONSTANT(T, 113, 9.8867173178574875515293357145875120137676e+05),
0440          BOOST_MATH_BIG_CONSTANT(T, 113, -4.4261178812193528551544261731796888257644e+07),
0441          BOOST_MATH_BIG_CONSTANT(T, 113, 1.6453010340778116475788083817762403540097e+09),
0442          BOOST_MATH_BIG_CONSTANT(T, 113, -5.0432401330113978669454035365747869477960e+10),
0443          BOOST_MATH_BIG_CONSTANT(T, 113, 1.2462165331309799059332310595587606836357e+12),
0444          BOOST_MATH_BIG_CONSTANT(T, 113, -2.3299800389951335932792950236410844978273e+13),
0445          BOOST_MATH_BIG_CONSTANT(T, 113, 2.5748218240248714177527965706790413406639e+14),
0446          BOOST_MATH_BIG_CONSTANT(T, 113, 1.8330014378766930869945511450377736037385e+15),
0447          BOOST_MATH_BIG_CONSTANT(T, 113, -1.8494610073827453236940544799030787866218e+17),
0448          BOOST_MATH_BIG_CONSTANT(T, 113, 5.7244661371420647691301043350229977856476e+18),
0449          BOOST_MATH_BIG_CONSTANT(T, 113, -1.2386378807889388140099109087465781254321e+20),
0450          BOOST_MATH_BIG_CONSTANT(T, 113, 2.1104000573102013529518477353943384110982e+21),
0451          BOOST_MATH_BIG_CONSTANT(T, 113, -2.9426541092239879262282594572224300191016e+22),
0452          BOOST_MATH_BIG_CONSTANT(T, 113, 3.4061439136301913488512592402635688101020e+23),
0453          BOOST_MATH_BIG_CONSTANT(T, 113, -3.2836554760521986358980180942859101564671e+24),
0454          BOOST_MATH_BIG_CONSTANT(T, 113, 2.6270285589905206294944214795661236766988e+25),
0455          BOOST_MATH_BIG_CONSTANT(T, 113, -1.7278631455211972017740134341610659484259e+26),
0456          BOOST_MATH_BIG_CONSTANT(T, 113, 9.1971734473772196124736986948034978906801e+26),
0457          BOOST_MATH_BIG_CONSTANT(T, 113, -3.8669270707172568763908838463689093500098e+27),
0458          BOOST_MATH_BIG_CONSTANT(T, 113, 1.2368879358870281916900125550129211146626e+28),
0459          BOOST_MATH_BIG_CONSTANT(T, 113, -2.8296235063297831758204519071113999839858e+28),
0460          BOOST_MATH_BIG_CONSTANT(T, 113, 4.1253861666023020670144616019148954773662e+28),
0461          BOOST_MATH_BIG_CONSTANT(T, 113, -2.8809536950051955163648980306847791014734e+28) };
0462       return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0463    }
0464    else if(x < 100)
0465    {
0466       // Bessel I0 over[30, 100]
0467       // Max error in interpolated form : 1.487e-34
0468       // Max Error found at float128 precision = Poly : 1.929924e-34
0469       static const T P[] = {
0470          BOOST_MATH_BIG_CONSTANT(T, 113, 3.9894228040143267793996798658172135362278e-01),
0471          BOOST_MATH_BIG_CONSTANT(T, 113, 4.9867785050179084714910130342157246539820e-02),
0472          BOOST_MATH_BIG_CONSTANT(T, 113, 2.8050629090725751585266360464766768437048e-02),
0473          BOOST_MATH_BIG_CONSTANT(T, 113, 2.9219405302833158254515212437025679637597e-02),
0474          BOOST_MATH_BIG_CONSTANT(T, 113, 4.4742214371598631578107310396249912330627e-02),
0475          BOOST_MATH_BIG_CONSTANT(T, 113, 9.0602983776478659136184969363625092585520e-02),
0476          BOOST_MATH_BIG_CONSTANT(T, 113, 2.2839507231977478205885469900971893734770e-01),
0477          BOOST_MATH_BIG_CONSTANT(T, 113, 6.8925739165733823730525449511456529001868e-01),
0478          BOOST_MATH_BIG_CONSTANT(T, 113, 2.4238082222874015159424842335385854632223e+00),
0479          BOOST_MATH_BIG_CONSTANT(T, 113, 9.6759648427182491050716309699208988458050e+00),
0480          BOOST_MATH_BIG_CONSTANT(T, 113, 4.7292246491169360014875196108746167872215e+01),
0481          BOOST_MATH_BIG_CONSTANT(T, 113, 3.1001411442786230340015781205680362993575e+01),
0482          BOOST_MATH_BIG_CONSTANT(T, 113, 9.8277628835804873490331739499978938078848e+03),
0483          BOOST_MATH_BIG_CONSTANT(T, 113, -3.1208326312801432038715638596517882759639e+05),
0484          BOOST_MATH_BIG_CONSTANT(T, 113, 9.4813611580683862051838126076298945680803e+06),
0485          BOOST_MATH_BIG_CONSTANT(T, 113, -2.1278197693321821164135890132925119054391e+08),
0486          BOOST_MATH_BIG_CONSTANT(T, 113, 3.3190303792682886967459489059860595063574e+09),
0487          BOOST_MATH_BIG_CONSTANT(T, 113, -2.1580767338646580750893606158043485767644e+10),
0488          BOOST_MATH_BIG_CONSTANT(T, 113, -5.0256008808415702780816006134784995506549e+11),
0489          BOOST_MATH_BIG_CONSTANT(T, 113, 1.9044186472918017896554580836514681614475e+13),
0490          BOOST_MATH_BIG_CONSTANT(T, 113, -3.2521078890073151875661384381880225635135e+14),
0491          BOOST_MATH_BIG_CONSTANT(T, 113, 3.3620352486836976842181057590770636605454e+15),
0492          BOOST_MATH_BIG_CONSTANT(T, 113, -2.0375525734060401555856465179734887312420e+16),
0493          BOOST_MATH_BIG_CONSTANT(T, 113, 5.6392664899881014534361728644608549445131e+16)
0494       };
0495       return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0496    }
0497    else
0498    {
0499       // Bessel I0 over[100, INF]
0500       // Max error in interpolated form : 5.459e-35
0501       // Max Error found at float128 precision = Poly : 1.472240e-34
0502       static const T P[] = {
0503          BOOST_MATH_BIG_CONSTANT(T, 113, 3.9894228040143267793994605993438166526772e-01),
0504          BOOST_MATH_BIG_CONSTANT(T, 113, 4.9867785050179084742493257495245185241487e-02),
0505          BOOST_MATH_BIG_CONSTANT(T, 113, 2.8050629090725735167652437695397756897920e-02),
0506          BOOST_MATH_BIG_CONSTANT(T, 113, 2.9219405302839307466358297347675795965363e-02),
0507          BOOST_MATH_BIG_CONSTANT(T, 113, 4.4742214369972689474366968442268908028204e-02),
0508          BOOST_MATH_BIG_CONSTANT(T, 113, 9.0602984099194778006610058410222616383078e-02),
0509          BOOST_MATH_BIG_CONSTANT(T, 113, 2.2839502241666629677015839125593079416327e-01),
0510          BOOST_MATH_BIG_CONSTANT(T, 113, 6.8926354981801627920292655818232972385750e-01),
0511          BOOST_MATH_BIG_CONSTANT(T, 113, 2.4231921590621824187100989532173995000655e+00),
0512          BOOST_MATH_BIG_CONSTANT(T, 113, 9.7264260959693775207585700654645245723497e+00),
0513          BOOST_MATH_BIG_CONSTANT(T, 113, 4.3890136225398811195878046856373030127018e+01),
0514          BOOST_MATH_BIG_CONSTANT(T, 113, 2.1999720924619285464910452647408431234369e+02),
0515          BOOST_MATH_BIG_CONSTANT(T, 113, 1.2076909538525038580501368530598517194748e+03),
0516          BOOST_MATH_BIG_CONSTANT(T, 113, 7.5684635141332367730007149159063086133399e+03),
0517          BOOST_MATH_BIG_CONSTANT(T, 113, 3.5178192543258299267923025833141286569141e+04),
0518          BOOST_MATH_BIG_CONSTANT(T, 113, 6.2966297919851965784482163987240461837728e+05) };
0519       T ex = exp(x / 2);
0520       T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0521       result *= ex;
0522       return result;
0523    }
0524 }
0525 
0526 template <typename T>
0527 T bessel_i0_imp(const T& x, const std::integral_constant<int, 0>&)
0528 {
0529    if(boost::math::tools::digits<T>() <= 24)
0530       return bessel_i0_imp(x, std::integral_constant<int, 24>());
0531    else if(boost::math::tools::digits<T>() <= 53)
0532       return bessel_i0_imp(x, std::integral_constant<int, 53>());
0533    else if(boost::math::tools::digits<T>() <= 64)
0534       return bessel_i0_imp(x, std::integral_constant<int, 64>());
0535    else if(boost::math::tools::digits<T>() <= 113)
0536       return bessel_i0_imp(x, std::integral_constant<int, 113>());
0537    BOOST_MATH_ASSERT(0);
0538    return 0;
0539 }
0540 
0541 template <typename T>
0542 inline T bessel_i0(const T& x)
0543 {
0544    typedef std::integral_constant<int,
0545       ((std::numeric_limits<T>::digits == 0) || (std::numeric_limits<T>::radix != 2)) ?
0546       0 :
0547       std::numeric_limits<T>::digits <= 24 ?
0548       24 :
0549       std::numeric_limits<T>::digits <= 53 ?
0550       53 :
0551       std::numeric_limits<T>::digits <= 64 ?
0552       64 :
0553       std::numeric_limits<T>::digits <= 113 ?
0554       113 : -1
0555    > tag_type;
0556 
0557    bessel_i0_initializer<T, tag_type>::force_instantiate();
0558    return bessel_i0_imp(x, tag_type());
0559 }
0560 
0561 }}} // namespaces
0562 
0563 #endif // BOOST_MATH_BESSEL_I0_HPP
0564
[ 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!