1 //  Copyright (c) 2006 Xiaogang Zhang
2 //  Use, modification and distribution are subject to the
3 //  Boost Software License, Version 1.0. (See accompanying file
4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5 
6 #ifndef BOOST_MATH_BESSEL_I0_HPP
7 #define BOOST_MATH_BESSEL_I0_HPP
8 
9 #ifdef _MSC_VER
10 #pragma once
11 #endif
12 
13 #include <boost/math/tools/rational.hpp>
14 #include <boost/math/tools/big_constant.hpp>
15 #include <boost/assert.hpp>
16 
17 // Modified Bessel function of the first kind of order zero
18 // minimax rational approximations on intervals, see
19 // Blair and Edwards, Chalk River Report AECL-4928, 1974
20 
21 namespace boost { namespace math { namespace detail{
22 
23 template <typename T>
24 T bessel_i0(T x);
25 
26 template <class T>
27 struct bessel_i0_initializer
28 {
29    struct init
30    {
initboost::math::detail::bessel_i0_initializer::init31       init()
32       {
33          do_init();
34       }
do_initboost::math::detail::bessel_i0_initializer::init35       static void do_init()
36       {
37          bessel_i0(T(1));
38       }
force_instantiateboost::math::detail::bessel_i0_initializer::init39       void force_instantiate()const{}
40    };
41    static const init initializer;
force_instantiateboost::math::detail::bessel_i0_initializer42    static void force_instantiate()
43    {
44       initializer.force_instantiate();
45    }
46 };
47 
48 template <class T>
49 const typename bessel_i0_initializer<T>::init bessel_i0_initializer<T>::initializer;
50 
51 template <typename T>
bessel_i0(T x)52 T bessel_i0(T x)
53 {
54     bessel_i0_initializer<T>::force_instantiate();
55 
56     static const T P1[] = {
57         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2335582639474375249e+15)),
58         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.5050369673018427753e+14)),
59         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.2940087627407749166e+13)),
60         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -8.4925101247114157499e+11)),
61         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1912746104985237192e+10)),
62         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0313066708737980747e+08)),
63         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.9545626019847898221e+05)),
64         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.4125195876041896775e+03)),
65         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -7.0935347449210549190e+00)),
66         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5453977791786851041e-02)),
67         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.5172644670688975051e-05)),
68         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.0517226450451067446e-08)),
69         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.6843448573468483278e-11)),
70         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5982226675653184646e-14)),
71         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.2487866627945699800e-18)),
72     };
73     static const T Q1[] = {
74         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2335582639474375245e+15)),
75         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.8858692566751002988e+12)),
76         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.2207067397808979846e+10)),
77         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0377081058062166144e+07)),
78         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.8527560179962773045e+03)),
79         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
80     };
81     static const T P2[] = {
82         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2210262233306573296e-04)),
83         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3067392038106924055e-02)),
84         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4700805721174453923e-01)),
85         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.5674518371240761397e+00)),
86         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.3517945679239481621e+01)),
87         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.1611322818701131207e+01)),
88         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.6090021968656180000e+00)),
89     };
90     static const T Q2[] = {
91         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.5194330231005480228e-04)),
92         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.2547697594819615062e-02)),
93         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1151759188741312645e+00)),
94         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3982595353892851542e+01)),
95         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.0228002066743340583e+01)),
96         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.5539563258012929600e+01)),
97         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.1446690275135491500e+01)),
98         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
99     };
100     T value, factor, r;
101 
102     BOOST_MATH_STD_USING
103     using namespace boost::math::tools;
104 
105     BOOST_ASSERT(x >= 0); // negative x is handled before we get here
106     if (x == 0)
107     {
108         return static_cast<T>(1);
109     }
110     if (x <= 15)                        // x in (0, 15]
111     {
112         T y = x * x;
113         value = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);
114     }
115     else                                // x in (15, \infty)
116     {
117         T y = 1 / x - T(1) / 15;
118         r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);
119         factor = exp(x) / sqrt(x);
120         value = factor * r;
121     }
122 
123     return value;
124 }
125 
126 }}} // namespaces
127 
128 #endif // BOOST_MATH_BESSEL_I0_HPP
129 
130