special_functions/detail/bessel_i1.hpp

0001 //  Copyright (c) 2017 John Maddock
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 // Modified Bessel function of the first kind of order zero
0007 // we use the approximating forms derived in:
0008 // "Rational Approximations for the Modified Bessel Function of the First Kind - I1(x) for Computations with Double Precision"
0009 // by Pavel Holoborodko,
0010 // see http://www.advanpix.com/2015/11/12/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i1-for-computations-with-double-precision/
0011 // The actual coefficients used are our own, and extend Pavel's work to precision's other than double.
0012
0013 #ifndef BOOST_MATH_BESSEL_I1_HPP
0014 #define BOOST_MATH_BESSEL_I1_HPP
0015
0016 #ifdef _MSC_VER
0017 #pragma once
0018 #endif
0019
0020 #include <boost/math/tools/rational.hpp>
0021 #include <boost/math/tools/big_constant.hpp>
0022 #include <boost/math/tools/assert.hpp>
0023
0024 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
0025 //
0026 // This is the only way we can avoid
0027 // warning: non-standard suffix on floating constant [-Wpedantic]
0028 // when building with -Wall -pedantic.  Neither __extension__
0029 // nor #pragma diagnostic ignored work :(
0030 //
0031 #pragma GCC system_header
0032 #endif
0033
0034 // Modified Bessel function of the first kind of order one
0035 // minimax rational approximations on intervals, see
0036 // Blair and Edwards, Chalk River Report AECL-4928, 1974
0037
0038 namespace boost { namespace math { namespace detail{
0039
0040 template <typename T>
0041 T bessel_i1(const T& x);
0042
0043 template <class T, class tag>
0044 struct bessel_i1_initializer
0045 {
0046    struct init
0047    {
0048       init()
0049       {
0050          do_init(tag());
0051       }
0052       static void do_init(const std::integral_constant<int, 64>&)
0053       {
0054          bessel_i1(T(1));
0055          bessel_i1(T(15));
0056          bessel_i1(T(80));
0057          bessel_i1(T(101));
0058       }
0059       static void do_init(const std::integral_constant<int, 113>&)
0060       {
0061          bessel_i1(T(1));
0062          bessel_i1(T(10));
0063          bessel_i1(T(14));
0064          bessel_i1(T(19));
0065          bessel_i1(T(34));
0066          bessel_i1(T(99));
0067          bessel_i1(T(101));
0068       }
0069       template <class U>
0070       static void do_init(const U&) {}
0071       void force_instantiate()const{}
0072    };
0073    static const init initializer;
0074    static void force_instantiate()
0075    {
0076       initializer.force_instantiate();
0077    }
0078 };
0079
0080 template <class T, class tag>
0081 const typename bessel_i1_initializer<T, tag>::init bessel_i1_initializer<T, tag>::initializer;
0082
0083 template <typename T, int N>
0084 T bessel_i1_imp(const T&, const std::integral_constant<int, N>&)
0085 {
0086    BOOST_MATH_ASSERT(0);
0087    return 0;
0088 }
0089
0090 template <typename T>
0091 T bessel_i1_imp(const T& x, const std::integral_constant<int, 24>&)
0092 {
0093    BOOST_MATH_STD_USING
0094       if(x < 7.75)
0095       {
0096          //Max error in interpolated form : 1.348e-08
0097          // Max Error found at float precision = Poly : 1.469121e-07
0098          static const float P[] = {
0099             8.333333221e-02f,
0100             6.944453712e-03f,
0101             3.472097211e-04f,
0102             1.158047174e-05f,
0103             2.739745142e-07f,
0104             5.135884609e-09f,
0105             5.262251502e-11f,
0106             1.331933703e-12f
0107          };
0108          T a = x * x / 4;
0109          T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) };
0110          return x * boost::math::tools::evaluate_polynomial(Q, a) / 2;
0111       }
0112       else
0113       {
0114          // Max error in interpolated form: 9.000e-08
0115          // Max Error found at float precision = Poly: 1.044345e-07
0116
0117          static const float P[] = {
0118             3.98942115977513013e-01f,
0119             -1.49581264836620262e-01f,
0120             -4.76475741878486795e-02f,
0121             -2.65157315524784407e-02f,
0122             -1.47148600683672014e-01f
0123          };
0124          T ex = exp(x / 2);
0125          T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0126          result *= ex;
0127          return result;
0128       }
0129 }
0130
0131 template <typename T>
0132 T bessel_i1_imp(const T& x, const std::integral_constant<int, 53>&)
0133 {
0134    BOOST_MATH_STD_USING
0135    if(x < 7.75)
0136    {
0137       // Bessel I0 over[10 ^ -16, 7.75]
0138       // Max error in interpolated form: 5.639e-17
0139       // Max Error found at double precision = Poly: 1.795559e-16
0140
0141       static const double P[] = {
0142          8.333333333333333803e-02,
0143          6.944444444444341983e-03,
0144          3.472222222225921045e-04,
0145          1.157407407354987232e-05,
0146          2.755731926254790268e-07,
0147          4.920949692800671435e-09,
0148          6.834657311305621830e-11,
0149          7.593969849687574339e-13,
0150          6.904822652741917551e-15,
0151          5.220157095351373194e-17,
0152          3.410720494727771276e-19,
0153          1.625212890947171108e-21,
0154          1.332898928162290861e-23
0155       };
0156       T a = x * x / 4;
0157       T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) };
0158       return x * boost::math::tools::evaluate_polynomial(Q, a) / 2;
0159    }
0160    else if(x < 500)
0161    {
0162       // Max error in interpolated form: 1.796e-16
0163       // Max Error found at double precision = Poly: 2.898731e-16
0164
0165       static const double P[] = {
0166          3.989422804014406054e-01,
0167          -1.496033551613111533e-01,
0168          -4.675104253598537322e-02,
0169          -4.090895951581637791e-02,
0170          -5.719036414430205390e-02,
0171          -1.528189554374492735e-01,
0172          3.458284470977172076e+00,
0173          -2.426181371595021021e+02,
0174          1.178785865993440669e+04,
0175          -4.404655582443487334e+05,
0176          1.277677779341446497e+07,
0177          -2.903390398236656519e+08,
0178          5.192386898222206474e+09,
0179          -7.313784438967834057e+10,
0180          8.087824484994859552e+11,
0181          -6.967602516005787001e+12,
0182          4.614040809616582764e+13,
0183          -2.298849639457172489e+14,
0184          8.325554073334618015e+14,
0185          -2.067285045778906105e+15,
0186          3.146401654361325073e+15,
0187          -2.213318202179221945e+15
0188       };
0189       return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0190    }
0191    else
0192    {
0193       // Max error in interpolated form: 1.320e-19
0194       // Max Error found at double precision = Poly: 7.065357e-17
0195       static const double P[] = {
0196          3.989422804014314820e-01,
0197          -1.496033551467584157e-01,
0198          -4.675105322571775911e-02,
0199          -4.090421597376992892e-02,
0200          -5.843630344778927582e-02
0201       };
0202       T ex = exp(x / 2);
0203       T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0204       result *= ex;
0205       return result;
0206    }
0207 }
0208
0209 template <typename T>
0210 T bessel_i1_imp(const T& x, const std::integral_constant<int, 64>&)
0211 {
0212    BOOST_MATH_STD_USING
0213       if(x < 7.75)
0214       {
0215          // Bessel I0 over[10 ^ -16, 7.75]
0216          // Max error in interpolated form: 8.086e-21
0217          // Max Error found at float80 precision = Poly: 7.225090e-20
0218          static const T P[] = {
0219             BOOST_MATH_BIG_CONSTANT(T, 64, 8.33333333333333333340071817e-02),
0220             BOOST_MATH_BIG_CONSTANT(T, 64, 6.94444444444444442462728070e-03),
0221             BOOST_MATH_BIG_CONSTANT(T, 64, 3.47222222222222318886683883e-04),
0222             BOOST_MATH_BIG_CONSTANT(T, 64, 1.15740740740738880709555060e-05),
0223             BOOST_MATH_BIG_CONSTANT(T, 64, 2.75573192240046222242685145e-07),
0224             BOOST_MATH_BIG_CONSTANT(T, 64, 4.92094986131253986838697503e-09),
0225             BOOST_MATH_BIG_CONSTANT(T, 64, 6.83465258979924922633502182e-11),
0226             BOOST_MATH_BIG_CONSTANT(T, 64, 7.59405830675154933645967137e-13),
0227             BOOST_MATH_BIG_CONSTANT(T, 64, 6.90369179710633344508897178e-15),
0228             BOOST_MATH_BIG_CONSTANT(T, 64, 5.23003610041709452814262671e-17),
0229             BOOST_MATH_BIG_CONSTANT(T, 64, 3.35291901027762552549170038e-19),
0230             BOOST_MATH_BIG_CONSTANT(T, 64, 1.83991379419781823063672109e-21),
0231             BOOST_MATH_BIG_CONSTANT(T, 64, 8.87732714140192556332037815e-24),
0232             BOOST_MATH_BIG_CONSTANT(T, 64, 3.32120654663773147206454247e-26),
0233             BOOST_MATH_BIG_CONSTANT(T, 64, 1.95294659305369207813486871e-28)
0234          };
0235          T a = x * x / 4;
0236          T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) };
0237          return x * boost::math::tools::evaluate_polynomial(Q, a) / 2;
0238       }
0239       else if(x < 20)
0240       {
0241          // Max error in interpolated form: 4.258e-20
0242          // Max Error found at float80 precision = Poly: 2.851105e-19
0243          // Maximum Deviation Found : 3.887e-20
0244          // Expected Error Term : 3.887e-20
0245          // Maximum Relative Change in Control Points : 1.681e-04
0246          static const T P[] = {
0247             BOOST_MATH_BIG_CONSTANT(T, 64, 3.98942260530218897338680e-01),
0248             BOOST_MATH_BIG_CONSTANT(T, 64, -1.49599542849073670179540e-01),
0249             BOOST_MATH_BIG_CONSTANT(T, 64, -4.70492865454119188276875e-02),
0250             BOOST_MATH_BIG_CONSTANT(T, 64, -3.12389893307392002405869e-02),
0251             BOOST_MATH_BIG_CONSTANT(T, 64, 1.49696126385202602071197e-01),
0252             BOOST_MATH_BIG_CONSTANT(T, 64, -3.84206507612717711565967e+01),
0253             BOOST_MATH_BIG_CONSTANT(T, 64, 2.14748094784412558689584e+03),
0254             BOOST_MATH_BIG_CONSTANT(T, 64, -7.70652726663596993005669e+04),
0255             BOOST_MATH_BIG_CONSTANT(T, 64, 2.01659736164815617174439e+06),
0256             BOOST_MATH_BIG_CONSTANT(T, 64, -4.04740659606466305607544e+07),
0257             BOOST_MATH_BIG_CONSTANT(T, 64, 6.38383394696382837263656e+08),
0258             BOOST_MATH_BIG_CONSTANT(T, 64, -8.00779638649147623107378e+09),
0259             BOOST_MATH_BIG_CONSTANT(T, 64, 8.02338237858684714480491e+10),
0260             BOOST_MATH_BIG_CONSTANT(T, 64, -6.41198553664947312995879e+11),
0261             BOOST_MATH_BIG_CONSTANT(T, 64, 4.05915186909564986897554e+12),
0262             BOOST_MATH_BIG_CONSTANT(T, 64, -2.00907636964168581116181e+13),
0263             BOOST_MATH_BIG_CONSTANT(T, 64, 7.60855263982359981275199e+13),
0264             BOOST_MATH_BIG_CONSTANT(T, 64, -2.12901817219239205393806e+14),
0265             BOOST_MATH_BIG_CONSTANT(T, 64, 4.14861794397709807823575e+14),
0266             BOOST_MATH_BIG_CONSTANT(T, 64, -5.02808138522587680348583e+14),
0267             BOOST_MATH_BIG_CONSTANT(T, 64, 2.85505477056514919387171e+14)
0268          };
0269          return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0270       }
0271       else if(x < 100)
0272       {
0273          // Bessel I0 over [15, 50]
0274          // Maximum Deviation Found:                     2.444e-20
0275          // Expected Error Term : 2.438e-20
0276          // Maximum Relative Change in Control Points : 2.101e-03
0277          // Max Error found at float80 precision = Poly : 6.029974e-20
0278
0279          static const T P[] = {
0280             BOOST_MATH_BIG_CONSTANT(T, 64, 3.98942280401431675205845e-01),
0281             BOOST_MATH_BIG_CONSTANT(T, 64, -1.49603355149968887210170e-01),
0282             BOOST_MATH_BIG_CONSTANT(T, 64, -4.67510486284376330257260e-02),
0283             BOOST_MATH_BIG_CONSTANT(T, 64, -4.09071458907089270559464e-02),
0284             BOOST_MATH_BIG_CONSTANT(T, 64, -5.75278280327696940044714e-02),
0285             BOOST_MATH_BIG_CONSTANT(T, 64, -1.10591299500956620739254e-01),
0286             BOOST_MATH_BIG_CONSTANT(T, 64, -2.77061766699949309115618e-01),
0287             BOOST_MATH_BIG_CONSTANT(T, 64, -5.42683771801837596371638e-01),
0288             BOOST_MATH_BIG_CONSTANT(T, 64, -9.17021412070404158464316e+00),
0289             BOOST_MATH_BIG_CONSTANT(T, 64, 1.04154379346763380543310e+02),
0290             BOOST_MATH_BIG_CONSTANT(T, 64, -1.43462345357478348323006e+03),
0291             BOOST_MATH_BIG_CONSTANT(T, 64, 9.98109660274422449523837e+03),
0292             BOOST_MATH_BIG_CONSTANT(T, 64, -3.74438822767781410362757e+04)
0293          };
0294          return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0295       }
0296       else
0297       {
0298          // Bessel I0 over[100, INF]
0299          // Max error in interpolated form: 2.456e-20
0300          // Max Error found at float80 precision = Poly: 5.446356e-20
0301          static const T P[] = {
0302             BOOST_MATH_BIG_CONSTANT(T, 64, 3.98942280401432677958445e-01),
0303             BOOST_MATH_BIG_CONSTANT(T, 64, -1.49603355150537411254359e-01),
0304             BOOST_MATH_BIG_CONSTANT(T, 64, -4.67510484842456251368526e-02),
0305             BOOST_MATH_BIG_CONSTANT(T, 64, -4.09071676503922479645155e-02),
0306             BOOST_MATH_BIG_CONSTANT(T, 64, -5.75256179814881566010606e-02),
0307             BOOST_MATH_BIG_CONSTANT(T, 64, -1.10754910257965227825040e-01),
0308             BOOST_MATH_BIG_CONSTANT(T, 64, -2.67858639515616079840294e-01),
0309             BOOST_MATH_BIG_CONSTANT(T, 64, -9.17266479586791298924367e-01)
0310          };
0311          T ex = exp(x / 2);
0312          T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0313          result *= ex;
0314          return result;
0315       }
0316 }
0317
0318 template <typename T>
0319 T bessel_i1_imp(const T& x, const std::integral_constant<int, 113>&)
0320 {
0321    BOOST_MATH_STD_USING
0322    if(x < 7.75)
0323    {
0324       // Bessel I0 over[10 ^ -34, 7.75]
0325       // Max error in interpolated form: 1.835e-35
0326       // Max Error found at float128 precision = Poly: 1.645036e-34
0327
0328       static const T P[] = {
0329          BOOST_MATH_BIG_CONSTANT(T, 113, 8.3333333333333333333333333333333331804098e-02),
0330          BOOST_MATH_BIG_CONSTANT(T, 113, 6.9444444444444444444444444444445418303082e-03),
0331          BOOST_MATH_BIG_CONSTANT(T, 113, 3.4722222222222222222222222222119082346591e-04),
0332          BOOST_MATH_BIG_CONSTANT(T, 113, 1.1574074074074074074074074078415867655987e-05),
0333          BOOST_MATH_BIG_CONSTANT(T, 113, 2.7557319223985890652557318255143448192453e-07),
0334          BOOST_MATH_BIG_CONSTANT(T, 113, 4.9209498614260519022423916850415000626427e-09),
0335          BOOST_MATH_BIG_CONSTANT(T, 113, 6.8346525853139609753354247043900442393686e-11),
0336          BOOST_MATH_BIG_CONSTANT(T, 113, 7.5940584281266233060080535940234144302217e-13),
0337          BOOST_MATH_BIG_CONSTANT(T, 113, 6.9036894801151120925605467963949641957095e-15),
0338          BOOST_MATH_BIG_CONSTANT(T, 113, 5.2300677879659941472662086395055636394839e-17),
0339          BOOST_MATH_BIG_CONSTANT(T, 113, 3.3526075563884539394691458717439115962233e-19),
0340          BOOST_MATH_BIG_CONSTANT(T, 113, 1.8420920639497841692288943167036233338434e-21),
0341          BOOST_MATH_BIG_CONSTANT(T, 113, 8.7718669711748690065381181691546032291365e-24),
0342          BOOST_MATH_BIG_CONSTANT(T, 113, 3.6549445715236427401845636880769861424730e-26),
0343          BOOST_MATH_BIG_CONSTANT(T, 113, 1.3437296196812697924703896979250126739676e-28),
0344          BOOST_MATH_BIG_CONSTANT(T, 113, 4.3912734588619073883015937023564978854893e-31),
0345          BOOST_MATH_BIG_CONSTANT(T, 113, 1.2839967682792395867255384448052781306897e-33),
0346          BOOST_MATH_BIG_CONSTANT(T, 113, 3.3790094235693528861015312806394354114982e-36),
0347          BOOST_MATH_BIG_CONSTANT(T, 113, 8.0423861671932104308662362292359563970482e-39),
0348          BOOST_MATH_BIG_CONSTANT(T, 113, 1.7493858979396446292135661268130281652945e-41),
0349          BOOST_MATH_BIG_CONSTANT(T, 113, 3.2786079392547776769387921361408303035537e-44),
0350          BOOST_MATH_BIG_CONSTANT(T, 113, 8.2335693685833531118863552173880047183822e-47)
0351       };
0352       T a = x * x / 4;
0353       T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) };
0354       return x * boost::math::tools::evaluate_polynomial(Q, a) / 2;
0355    }
0356    else if(x < 11)
0357    {
0358       // Max error in interpolated form: 8.574e-36
0359       // Maximum Deviation Found : 4.689e-36
0360       // Expected Error Term : 3.760e-36
0361       // Maximum Relative Change in Control Points : 5.204e-03
0362       // Max Error found at float128 precision = Poly : 2.882561e-34
0363
0364       static const T P[] = {
0365          BOOST_MATH_BIG_CONSTANT(T, 113, 8.333333333333333326889717360850080939e-02),
0366          BOOST_MATH_BIG_CONSTANT(T, 113, 6.944444444444444511272790848815114507e-03),
0367          BOOST_MATH_BIG_CONSTANT(T, 113, 3.472222222222221892451965054394153443e-04),
0368          BOOST_MATH_BIG_CONSTANT(T, 113, 1.157407407407408437378868534321538798e-05),
0369          BOOST_MATH_BIG_CONSTANT(T, 113, 2.755731922398566216824909767320161880e-07),
0370          BOOST_MATH_BIG_CONSTANT(T, 113, 4.920949861426434829568192525456800388e-09),
0371          BOOST_MATH_BIG_CONSTANT(T, 113, 6.834652585308926245465686943255486934e-11),
0372          BOOST_MATH_BIG_CONSTANT(T, 113, 7.594058428179852047689599244015979196e-13),
0373          BOOST_MATH_BIG_CONSTANT(T, 113, 6.903689479655006062822949671528763738e-15),
0374          BOOST_MATH_BIG_CONSTANT(T, 113, 5.230067791254403974475987777406992984e-17),
0375          BOOST_MATH_BIG_CONSTANT(T, 113, 3.352607536815161679702105115200693346e-19),
0376          BOOST_MATH_BIG_CONSTANT(T, 113, 1.842092161364672561828681848278567885e-21),
0377          BOOST_MATH_BIG_CONSTANT(T, 113, 8.771862912600611801856514076709932773e-24),
0378          BOOST_MATH_BIG_CONSTANT(T, 113, 3.654958704184380914803366733193713605e-26),
0379          BOOST_MATH_BIG_CONSTANT(T, 113, 1.343688672071130980471207297730607625e-28),
0380          BOOST_MATH_BIG_CONSTANT(T, 113, 4.392252844664709532905868749753463950e-31),
0381          BOOST_MATH_BIG_CONSTANT(T, 113, 1.282086786672692641959912811902298600e-33),
0382          BOOST_MATH_BIG_CONSTANT(T, 113, 3.408812012322547015191398229942864809e-36),
0383          BOOST_MATH_BIG_CONSTANT(T, 113, 7.681220437734066258673404589233009892e-39),
0384          BOOST_MATH_BIG_CONSTANT(T, 113, 2.072417451640733785626701738789290055e-41),
0385          BOOST_MATH_BIG_CONSTANT(T, 113, 1.352218520142636864158849446833681038e-44),
0386          BOOST_MATH_BIG_CONSTANT(T, 113, 1.407918492276267527897751358794783640e-46)
0387       };
0388       T a = x * x / 4;
0389       T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) };
0390       return x * boost::math::tools::evaluate_polynomial(Q, a) / 2;
0391    }
0392    else if(x < 15)
0393    {
0394       //Max error in interpolated form: 7.599e-36
0395       // Maximum Deviation Found : 1.766e-35
0396       // Expected Error Term : 1.021e-35
0397       // Maximum Relative Change in Control Points : 6.228e-03
0398       static const T P[] = {
0399          BOOST_MATH_BIG_CONSTANT(T, 113, 8.333333333333255774414858563409941233e-02),
0400          BOOST_MATH_BIG_CONSTANT(T, 113, 6.944444444444897867884955912228700291e-03),
0401          BOOST_MATH_BIG_CONSTANT(T, 113, 3.472222222220954970397343617150959467e-04),
0402          BOOST_MATH_BIG_CONSTANT(T, 113, 1.157407407409660682751155024932538578e-05),
0403          BOOST_MATH_BIG_CONSTANT(T, 113, 2.755731922369973706427272809014190998e-07),
0404          BOOST_MATH_BIG_CONSTANT(T, 113, 4.920949861702265600960449699129258153e-09),
0405          BOOST_MATH_BIG_CONSTANT(T, 113, 6.834652583208361401197752793379677147e-11),
0406          BOOST_MATH_BIG_CONSTANT(T, 113, 7.594058441128280500819776168239988143e-13),
0407          BOOST_MATH_BIG_CONSTANT(T, 113, 6.903689413939268702265479276217647209e-15),
0408          BOOST_MATH_BIG_CONSTANT(T, 113, 5.230068069012898202890718644753625569e-17),
0409          BOOST_MATH_BIG_CONSTANT(T, 113, 3.352606552027491657204243201021677257e-19),
0410          BOOST_MATH_BIG_CONSTANT(T, 113, 1.842095100698532984651921750204843362e-21),
0411          BOOST_MATH_BIG_CONSTANT(T, 113, 8.771789051329870174925649852681844169e-24),
0412          BOOST_MATH_BIG_CONSTANT(T, 113, 3.655114381199979536997025497438385062e-26),
0413          BOOST_MATH_BIG_CONSTANT(T, 113, 1.343415732516712339472538688374589373e-28),
0414          BOOST_MATH_BIG_CONSTANT(T, 113, 4.396177019032432392793591204647901390e-31),
0415          BOOST_MATH_BIG_CONSTANT(T, 113, 1.277563309255167951005939802771456315e-33),
0416          BOOST_MATH_BIG_CONSTANT(T, 113, 3.449201419305514579791370198046544736e-36),
0417          BOOST_MATH_BIG_CONSTANT(T, 113, 7.415430703400740634202379012388035255e-39),
0418          BOOST_MATH_BIG_CONSTANT(T, 113, 2.195458831864936225409005027914934499e-41),
0419          BOOST_MATH_BIG_CONSTANT(T, 113, 8.829726762743879793396637797534668039e-45),
0420          BOOST_MATH_BIG_CONSTANT(T, 113, 1.698302711685624490806751012380215488e-46),
0421          BOOST_MATH_BIG_CONSTANT(T, 113, -2.062520475425422618494185821587228317e-49),
0422          BOOST_MATH_BIG_CONSTANT(T, 113, 6.732372906742845717148185173723304360e-52)
0423       };
0424       T a = x * x / 4;
0425       T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) };
0426       return x * boost::math::tools::evaluate_polynomial(Q, a) / 2;
0427    }
0428    else if(x < 20)
0429    {
0430       // Max error in interpolated form: 8.864e-36
0431       // Max Error found at float128 precision = Poly: 8.522841e-35
0432       static const T P[] = {
0433          BOOST_MATH_BIG_CONSTANT(T, 113, 3.989422793693152031514179994954750043e-01),
0434          BOOST_MATH_BIG_CONSTANT(T, 113, -1.496029423752889591425633234009799670e-01),
0435          BOOST_MATH_BIG_CONSTANT(T, 113, -4.682975926820553021482820043377990241e-02),
0436          BOOST_MATH_BIG_CONSTANT(T, 113, -3.138871171577224532369979905856458929e-02),
0437          BOOST_MATH_BIG_CONSTANT(T, 113, -8.765350219426341341990447005798111212e-01),
0438          BOOST_MATH_BIG_CONSTANT(T, 113, 5.321389275507714530941178258122955540e+01),
0439          BOOST_MATH_BIG_CONSTANT(T, 113, -2.727748393898888756515271847678850411e+03),
0440          BOOST_MATH_BIG_CONSTANT(T, 113, 1.123040820686242586086564998713862335e+05),
0441          BOOST_MATH_BIG_CONSTANT(T, 113, -3.784112378374753535335272752884808068e+06),
0442          BOOST_MATH_BIG_CONSTANT(T, 113, 1.054920416060932189433079126269416563e+08),
0443          BOOST_MATH_BIG_CONSTANT(T, 113, -2.450129415468060676827180524327749553e+09),
0444          BOOST_MATH_BIG_CONSTANT(T, 113, 4.758831882046487398739784498047935515e+10),
0445          BOOST_MATH_BIG_CONSTANT(T, 113, -7.736936520262204842199620784338052937e+11),
0446          BOOST_MATH_BIG_CONSTANT(T, 113, 1.051128683324042629513978256179115439e+13),
0447          BOOST_MATH_BIG_CONSTANT(T, 113, -1.188008285959794869092624343537262342e+14),
0448          BOOST_MATH_BIG_CONSTANT(T, 113, 1.108530004906954627420484180793165669e+15),
0449          BOOST_MATH_BIG_CONSTANT(T, 113, -8.441516828490144766650287123765318484e+15),
0450          BOOST_MATH_BIG_CONSTANT(T, 113, 5.158251664797753450664499268756393535e+16),
0451          BOOST_MATH_BIG_CONSTANT(T, 113, -2.467314522709016832128790443932896401e+17),
0452          BOOST_MATH_BIG_CONSTANT(T, 113, 8.896222045367960462945885220710294075e+17),
0453          BOOST_MATH_BIG_CONSTANT(T, 113, -2.273382139594876997203657902425653079e+18),
0454          BOOST_MATH_BIG_CONSTANT(T, 113, 3.669871448568623680543943144842394531e+18),
0455          BOOST_MATH_BIG_CONSTANT(T, 113, -2.813923031370708069940575240509912588e+18)
0456       };
0457       return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0458    }
0459    else if(x < 35)
0460    {
0461       // Max error in interpolated form: 6.028e-35
0462       // Max Error found at float128 precision = Poly: 1.368313e-34
0463
0464       static const T P[] = {
0465          BOOST_MATH_BIG_CONSTANT(T, 113, 3.989422804012941975429616956496046931e-01),
0466          BOOST_MATH_BIG_CONSTANT(T, 113, -1.496033550576049830976679315420681402e-01),
0467          BOOST_MATH_BIG_CONSTANT(T, 113, -4.675107835141866009896710750800622147e-02),
0468          BOOST_MATH_BIG_CONSTANT(T, 113, -4.090104965125365961928716504473692957e-02),
0469          BOOST_MATH_BIG_CONSTANT(T, 113, -5.842241652296980863361375208605487570e-02),
0470          BOOST_MATH_BIG_CONSTANT(T, 113, -1.063604828033747303936724279018650633e-02),
0471          BOOST_MATH_BIG_CONSTANT(T, 113, -9.113375972811586130949401996332817152e+00),
0472          BOOST_MATH_BIG_CONSTANT(T, 113, 6.334748570425075872639817839399823709e+02),
0473          BOOST_MATH_BIG_CONSTANT(T, 113, -3.759150758768733692594821032784124765e+04),
0474          BOOST_MATH_BIG_CONSTANT(T, 113, 1.863672813448915255286274382558526321e+06),
0475          BOOST_MATH_BIG_CONSTANT(T, 113, -7.798248643371718775489178767529282534e+07),
0476          BOOST_MATH_BIG_CONSTANT(T, 113, 2.769963173932801026451013022000669267e+09),
0477          BOOST_MATH_BIG_CONSTANT(T, 113, -8.381780137198278741566746511015220011e+10),
0478          BOOST_MATH_BIG_CONSTANT(T, 113, 2.163891337116820832871382141011952931e+12),
0479          BOOST_MATH_BIG_CONSTANT(T, 113, -4.764325864671438675151635117936912390e+13),
0480          BOOST_MATH_BIG_CONSTANT(T, 113, 8.925668307403332887856809510525154955e+14),
0481          BOOST_MATH_BIG_CONSTANT(T, 113, -1.416692606589060039334938090985713641e+16),
0482          BOOST_MATH_BIG_CONSTANT(T, 113, 1.892398600219306424294729851605944429e+17),
0483          BOOST_MATH_BIG_CONSTANT(T, 113, -2.107232903741874160308537145391245060e+18),
0484          BOOST_MATH_BIG_CONSTANT(T, 113, 1.930223393531877588898224144054112045e+19),
0485          BOOST_MATH_BIG_CONSTANT(T, 113, -1.427759576167665663373350433236061007e+20),
0486          BOOST_MATH_BIG_CONSTANT(T, 113, 8.306019279465532835530812122374386654e+20),
0487          BOOST_MATH_BIG_CONSTANT(T, 113, -3.653753000392125229440044977239174472e+21),
0488          BOOST_MATH_BIG_CONSTANT(T, 113, 1.140760686989511568435076842569804906e+22),
0489          BOOST_MATH_BIG_CONSTANT(T, 113, -2.249149337812510200795436107962504749e+22),
0490          BOOST_MATH_BIG_CONSTANT(T, 113, 2.101619088427348382058085685849420866e+22)
0491       };
0492       return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0493    }
0494    else if(x < 100)
0495    {
0496       // Max error in interpolated form: 5.494e-35
0497       // Max Error found at float128 precision = Poly: 1.214651e-34
0498
0499       static const T P[] = {
0500          BOOST_MATH_BIG_CONSTANT(T, 113, 3.989422804014326779399307367861631577e-01),
0501          BOOST_MATH_BIG_CONSTANT(T, 113, -1.496033551505372542086590873271571919e-01),
0502          BOOST_MATH_BIG_CONSTANT(T, 113, -4.675104848454290286276466276677172664e-02),
0503          BOOST_MATH_BIG_CONSTANT(T, 113, -4.090716742397105403027549796269213215e-02),
0504          BOOST_MATH_BIG_CONSTANT(T, 113, -5.752570419098513588311026680089351230e-02),
0505          BOOST_MATH_BIG_CONSTANT(T, 113, -1.107369803696534592906420980901195808e-01),
0506          BOOST_MATH_BIG_CONSTANT(T, 113, -2.699214194000085622941721628134575121e-01),
0507          BOOST_MATH_BIG_CONSTANT(T, 113, -7.953006169077813678478720427604462133e-01),
0508          BOOST_MATH_BIG_CONSTANT(T, 113, -2.746618809476524091493444128605380593e+00),
0509          BOOST_MATH_BIG_CONSTANT(T, 113, -1.084446249943196826652788161656973391e+01),
0510          BOOST_MATH_BIG_CONSTANT(T, 113, -5.020325182518980633783194648285500554e+01),
0511          BOOST_MATH_BIG_CONSTANT(T, 113, -1.510195971266257573425196228564489134e+02),
0512          BOOST_MATH_BIG_CONSTANT(T, 113, -5.241661863814900938075696173192225056e+03),
0513          BOOST_MATH_BIG_CONSTANT(T, 113, 1.323374362891993686413568398575539777e+05),
0514          BOOST_MATH_BIG_CONSTANT(T, 113, -4.112838452096066633754042734723911040e+06),
0515          BOOST_MATH_BIG_CONSTANT(T, 113, 9.369270194978310081563767560113534023e+07),
0516          BOOST_MATH_BIG_CONSTANT(T, 113, -1.704295412488936504389347368131134993e+09),
0517          BOOST_MATH_BIG_CONSTANT(T, 113, 2.320829576277038198439987439508754886e+10),
0518          BOOST_MATH_BIG_CONSTANT(T, 113, -2.258818139077875493434420764260185306e+11),
0519          BOOST_MATH_BIG_CONSTANT(T, 113, 1.396791306321498426110315039064592443e+12),
0520          BOOST_MATH_BIG_CONSTANT(T, 113, -4.217617301585849875301440316301068439e+12)
0521       };
0522       return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0523    }
0524    else
0525    {
0526       // Bessel I0 over[100, INF]
0527       // Max error in interpolated form: 6.081e-35
0528       // Max Error found at float128 precision = Poly: 1.407151e-34
0529       static const T P[] = {
0530          BOOST_MATH_BIG_CONSTANT(T, 113, 3.9894228040143267793994605993438200208417e-01),
0531          BOOST_MATH_BIG_CONSTANT(T, 113, -1.4960335515053725422747977247811372936584e-01),
0532          BOOST_MATH_BIG_CONSTANT(T, 113, -4.6751048484542891946087411826356811991039e-02),
0533          BOOST_MATH_BIG_CONSTANT(T, 113, -4.0907167423975030452875828826630006305665e-02),
0534          BOOST_MATH_BIG_CONSTANT(T, 113, -5.7525704189964886494791082898669060345483e-02),
0535          BOOST_MATH_BIG_CONSTANT(T, 113, -1.1073698056568248642163476807108190176386e-01),
0536          BOOST_MATH_BIG_CONSTANT(T, 113, -2.6992139012879749064623499618582631684228e-01),
0537          BOOST_MATH_BIG_CONSTANT(T, 113, -7.9530409594026597988098934027440110587905e-01),
0538          BOOST_MATH_BIG_CONSTANT(T, 113, -2.7462844478733532517044536719240098183686e+00),
0539          BOOST_MATH_BIG_CONSTANT(T, 113, -1.0870711340681926669381449306654104739256e+01),
0540          BOOST_MATH_BIG_CONSTANT(T, 113, -4.8510175413216969245241059608553222505228e+01),
0541          BOOST_MATH_BIG_CONSTANT(T, 113, -2.4094682286011573747064907919522894740063e+02),
0542          BOOST_MATH_BIG_CONSTANT(T, 113, -1.3128845936764406865199641778959502795443e+03),
0543          BOOST_MATH_BIG_CONSTANT(T, 113, -8.1655901321962541203257516341266838487359e+03),
0544          BOOST_MATH_BIG_CONSTANT(T, 113, -3.8019591025686295090160445920753823994556e+04),
0545          BOOST_MATH_BIG_CONSTANT(T, 113, -6.7008089049178178697338128837158732831105e+05)
0546       };
0547       T ex = exp(x / 2);
0548       T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
0549       result *= ex;
0550       return result;
0551    }
0552 }
0553
0554 template <typename T>
0555 T bessel_i1_imp(const T& x, const std::integral_constant<int, 0>&)
0556 {
0557    if(boost::math::tools::digits<T>() <= 24)
0558       return bessel_i1_imp(x, std::integral_constant<int, 24>());
0559    else if(boost::math::tools::digits<T>() <= 53)
0560       return bessel_i1_imp(x, std::integral_constant<int, 53>());
0561    else if(boost::math::tools::digits<T>() <= 64)
0562       return bessel_i1_imp(x, std::integral_constant<int, 64>());
0563    else if(boost::math::tools::digits<T>() <= 113)
0564       return bessel_i1_imp(x, std::integral_constant<int, 113>());
0565    BOOST_MATH_ASSERT(0);
0566    return 0;
0567 }
0568
0569 template <typename T>
0570 inline T bessel_i1(const T& x)
0571 {
0572    typedef std::integral_constant<int,
0573       ((std::numeric_limits<T>::digits == 0) || (std::numeric_limits<T>::radix != 2)) ?
0574       0 :
0575       std::numeric_limits<T>::digits <= 24 ?
0576       24 :
0577       std::numeric_limits<T>::digits <= 53 ?
0578       53 :
0579       std::numeric_limits<T>::digits <= 64 ?
0580       64 :
0581       std::numeric_limits<T>::digits <= 113 ?
0582       113 : -1
0583    > tag_type;
0584
0585    bessel_i1_initializer<T, tag_type>::force_instantiate();
0586    return bessel_i1_imp(x, tag_type());
0587 }
0588
0589 }}} // namespaces
0590
0591 #endif // BOOST_MATH_BESSEL_I1_HPP
0592