]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | ||
92f5a8d4 TL |
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 dianostic ignored work :( | |
24 | // | |
25 | #pragma GCC system_header | |
26 | #endif | |
27 | ||
7c673cae FG |
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 | { | |
42 | init() | |
43 | { | |
44 | do_init(); | |
45 | } | |
46 | static void do_init() | |
47 | { | |
48 | bessel_j1(T(1)); | |
49 | } | |
50 | void force_instantiate()const{} | |
51 | }; | |
52 | static const init initializer; | |
53 | 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> | |
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 |