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_J1_HPP
7 #define BOOST_MATH_BESSEL_J1_HPP
8
9 #ifdef _MSC_VER
10 #pragma once
11 #endif
12
13 #include <boost/math/constants/constants.hpp>
14 #include <boost/math/tools/rational.hpp>
15 #include <boost/math/tools/big_constant.hpp>
16 #include <boost/assert.hpp>
17
18 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
19 //
20 // This is the only way we can avoid
21 // warning: non-standard suffix on floating constant [-Wpedantic]
22 // when building with -Wall -pedantic. Neither __extension__
23 // nor #pragma diagnostic ignored work :(
24 //
25 #pragma GCC system_header
26 #endif
27
28 // Bessel function of the first kind of order one
29 // x <= 8, minimax rational approximations on root-bracketing intervals
30 // x > 8, Hankel asymptotic expansion in Hart, Computer Approximations, 1968
31
32 namespace boost { namespace math{ namespace detail{
33
34 template <typename T>
35 T bessel_j1(T x);
36
37 template <class T>
38 struct bessel_j1_initializer
39 {
40 struct init
41 {
initboost::math::detail::bessel_j1_initializer::init42 init()
43 {
44 do_init();
45 }
do_initboost::math::detail::bessel_j1_initializer::init46 static void do_init()
47 {
48 bessel_j1(T(1));
49 }
force_instantiateboost::math::detail::bessel_j1_initializer::init50 void force_instantiate()const{}
51 };
52 static const init initializer;
force_instantiateboost::math::detail::bessel_j1_initializer53 static void force_instantiate()
54 {
55 initializer.force_instantiate();
56 }
57 };
58
59 template <class T>
60 const typename bessel_j1_initializer<T>::init bessel_j1_initializer<T>::initializer;
61
62 template <typename T>
bessel_j1(T x)63 T bessel_j1(T x)
64 {
65 bessel_j1_initializer<T>::force_instantiate();
66
67 static const T P1[] = {
68 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4258509801366645672e+11)),
69 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6781041261492395835e+09)),
70 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1548696764841276794e+08)),
71 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 9.8062904098958257677e+05)),
72 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4615792982775076130e+03)),
73 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0650724020080236441e+01)),
74 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0767857011487300348e-02))
75 };
76 static const T Q1[] = {
77 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.1868604460820175290e+12)),
78 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.2091902282580133541e+10)),
79 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.0228375140097033958e+08)),
80 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.9117614494174794095e+05)),
81 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0742272239517380498e+03)),
82 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
83 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0))
84 };
85 static const T P2[] = {
86 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.7527881995806511112e+16)),
87 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.6608531731299018674e+15)),
88 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.6658018905416665164e+13)),
89 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.5580665670910619166e+11)),
90 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.8113931269860667829e+09)),
91 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.0793266148011179143e+06)),
92 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -7.5023342220781607561e+03)),
93 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.6179191852758252278e+00))
94 };
95 static const T Q2[] = {
96 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7253905888447681194e+18)),
97 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7128800897135812012e+16)),
98 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.4899346165481429307e+13)),
99 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.7622777286244082666e+11)),
100 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.4872502899596389593e+08)),
101 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1267125065029138050e+06)),
102 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3886978985861357615e+03)),
103 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))
104 };
105 static const T PC[] = {
106 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278571e+06)),
107 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9422465050776411957e+06)),
108 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.6033732483649391093e+06)),
109 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5235293511811373833e+06)),
110 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0982405543459346727e+05)),
111 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.6116166443246101165e+03)),
112 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0))
113 };
114 static const T QC[] = {
115 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278568e+06)),
116 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9341243899345856590e+06)),
117 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.5853394797230870728e+06)),
118 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5118095066341608816e+06)),
119 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0726385991103820119e+05)),
120 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4550094401904961825e+03)),
121 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))
122 };
123 static const T PS[] = {
124 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.3220913409857223519e+04)),
125 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.5145160675335701966e+04)),
126 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6178836581270835179e+04)),
127 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8494262873223866797e+04)),
128 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7063754290207680021e+03)),
129 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.5265133846636032186e+01)),
130 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0))
131 };
132 static const T QS[] = {
133 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.0871281941028743574e+05)),
134 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8194580422439972989e+06)),
135 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.4194606696037208929e+06)),
136 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0029443582266975117e+05)),
137 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7890229745772202641e+04)),
138 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.6383677696049909675e+02)),
139 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))
140 };
141 static const T x1 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.8317059702075123156e+00)),
142 x2 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.0155866698156187535e+00)),
143 x11 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 9.810e+02)),
144 x12 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.2527979248768438556e-04)),
145 x21 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7960e+03)),
146 x22 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.8330184381246462950e-05));
147
148 T value, factor, r, rc, rs, w;
149
150 BOOST_MATH_STD_USING
151 using namespace boost::math::tools;
152 using namespace boost::math::constants;
153
154 w = abs(x);
155 if (x == 0)
156 {
157 return static_cast<T>(0);
158 }
159 if (w <= 4) // w in (0, 4]
160 {
161 T y = x * x;
162 BOOST_ASSERT(sizeof(P1) == sizeof(Q1));
163 r = evaluate_rational(P1, Q1, y);
164 factor = w * (w + x1) * ((w - x11/256) - x12);
165 value = factor * r;
166 }
167 else if (w <= 8) // w in (4, 8]
168 {
169 T y = x * x;
170 BOOST_ASSERT(sizeof(P2) == sizeof(Q2));
171 r = evaluate_rational(P2, Q2, y);
172 factor = w * (w + x2) * ((w - x21/256) - x22);
173 value = factor * r;
174 }
175 else // w in (8, \infty)
176 {
177 T y = 8 / w;
178 T y2 = y * y;
179 BOOST_ASSERT(sizeof(PC) == sizeof(QC));
180 BOOST_ASSERT(sizeof(PS) == sizeof(QS));
181 rc = evaluate_rational(PC, QC, y2);
182 rs = evaluate_rational(PS, QS, y2);
183 factor = 1 / (sqrt(w) * constants::root_pi<T>());
184 //
185 // What follows is really just:
186 //
187 // T z = w - 0.75f * pi<T>();
188 // value = factor * (rc * cos(z) - y * rs * sin(z));
189 //
190 // but using the sin/cos addition rules plus constants
191 // for the values of sin/cos of 3PI/4 which then cancel
192 // out with corresponding terms in "factor".
193 //
194 T sx = sin(x);
195 T cx = cos(x);
196 value = factor * (rc * (sx - cx) + y * rs * (sx + cx));
197 }
198
199 if (x < 0)
200 {
201 value *= -1; // odd function
202 }
203 return value;
204 }
205
206 }}} // namespaces
207
208 #endif // BOOST_MATH_BESSEL_J1_HPP
209
210