]>
Commit | Line | Data |
---|---|---|
b32b8144 FG |
1 | // Copyright (c) 2006 Xiaogang Zhang |
2 | // Copyright (c) 2017 John Maddock | |
3 | // Use, modification and distribution are subject to the | |
4 | // Boost Software License, Version 1.0. (See accompanying file | |
5 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
6 | ||
7 | #ifndef BOOST_MATH_BESSEL_K0_HPP | |
8 | #define BOOST_MATH_BESSEL_K0_HPP | |
9 | ||
10 | #ifdef _MSC_VER | |
11 | #pragma once | |
12 | #pragma warning(push) | |
13 | #pragma warning(disable:4702) // Unreachable code (release mode only warning) | |
14 | #endif | |
15 | ||
16 | #include <boost/math/tools/rational.hpp> | |
17 | #include <boost/math/tools/big_constant.hpp> | |
18 | #include <boost/math/policies/error_handling.hpp> | |
19 | #include <boost/assert.hpp> | |
20 | ||
92f5a8d4 TL |
21 | #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128) |
22 | // | |
23 | // This is the only way we can avoid | |
24 | // warning: non-standard suffix on floating constant [-Wpedantic] | |
25 | // when building with -Wall -pedantic. Neither __extension__ | |
26 | // nor #pragma dianostic ignored work :( | |
27 | // | |
28 | #pragma GCC system_header | |
29 | #endif | |
30 | ||
b32b8144 FG |
31 | // Modified Bessel function of the second kind of order zero |
32 | // minimax rational approximations on intervals, see | |
33 | // Russon and Blair, Chalk River Report AECL-3461, 1969, | |
34 | // as revised by Pavel Holoborodko in "Rational Approximations | |
35 | // for the Modified Bessel Function of the Second Kind - K0(x) | |
36 | // for Computations with Double Precision", see | |
37 | // http://www.advanpix.com/2015/11/25/rational-approximations-for-the-modified-bessel-function-of-the-second-kind-k0-for-computations-with-double-precision/ | |
38 | // | |
39 | // The actual coefficients used are our own derivation (by JM) | |
40 | // since we extend to both greater and lesser precision than the | |
41 | // references above. We can also improve performance WRT to | |
42 | // Holoborodko without loss of precision. | |
43 | ||
44 | namespace boost { namespace math { namespace detail{ | |
45 | ||
46 | template <typename T> | |
47 | T bessel_k0(const T& x); | |
48 | ||
49 | template <class T, class tag> | |
50 | struct bessel_k0_initializer | |
51 | { | |
52 | struct init | |
53 | { | |
54 | init() | |
55 | { | |
56 | do_init(tag()); | |
57 | } | |
58 | static void do_init(const mpl::int_<113>&) | |
59 | { | |
60 | bessel_k0(T(0.5)); | |
61 | bessel_k0(T(1.5)); | |
62 | } | |
63 | static void do_init(const mpl::int_<64>&) | |
64 | { | |
65 | bessel_k0(T(0.5)); | |
66 | bessel_k0(T(1.5)); | |
67 | } | |
68 | template <class U> | |
69 | static void do_init(const U&){} | |
70 | void force_instantiate()const{} | |
71 | }; | |
72 | static const init initializer; | |
73 | static void force_instantiate() | |
74 | { | |
75 | initializer.force_instantiate(); | |
76 | } | |
77 | }; | |
78 | ||
79 | template <class T, class tag> | |
80 | const typename bessel_k0_initializer<T, tag>::init bessel_k0_initializer<T, tag>::initializer; | |
81 | ||
82 | ||
83 | template <typename T, int N> | |
84 | T bessel_k0_imp(const T& x, const mpl::int_<N>&) | |
85 | { | |
86 | BOOST_ASSERT(0); | |
87 | return 0; | |
88 | } | |
89 | ||
90 | template <typename T> | |
91 | T bessel_k0_imp(const T& x, const mpl::int_<24>&) | |
92 | { | |
93 | BOOST_MATH_STD_USING | |
94 | if(x <= 1) | |
95 | { | |
96 | // Maximum Deviation Found : 2.358e-09 | |
97 | // Expected Error Term : -2.358e-09 | |
98 | // Maximum Relative Change in Control Points : 9.552e-02 | |
99 | // Max Error found at float precision = Poly : 4.448220e-08 | |
100 | static const T Y = 1.137250900268554688f; | |
101 | static const T P[] = | |
102 | { | |
103 | -1.372508979104259711e-01f, | |
104 | 2.622545986273687617e-01f, | |
105 | 5.047103728247919836e-03f | |
106 | }; | |
107 | static const T Q[] = | |
108 | { | |
109 | 1.000000000000000000e+00f, | |
110 | -8.928694018000029415e-02f, | |
111 | 2.985980684180969241e-03f | |
112 | }; | |
113 | T a = x * x / 4; | |
114 | a = (tools::evaluate_rational(P, Q, a) + Y) * a + 1; | |
115 | ||
116 | // Maximum Deviation Found: 1.346e-09 | |
117 | // Expected Error Term : -1.343e-09 | |
118 | // Maximum Relative Change in Control Points : 2.405e-02 | |
119 | // Max Error found at float precision = Poly : 1.354814e-07 | |
120 | static const T P2[] = { | |
121 | 1.159315158e-01f, | |
122 | 2.789828686e-01f, | |
123 | 2.524902861e-02f, | |
124 | 8.457241514e-04f, | |
125 | 1.530051997e-05f | |
126 | }; | |
127 | return tools::evaluate_polynomial(P2, T(x * x)) - log(x) * a; | |
128 | } | |
129 | else | |
130 | { | |
131 | // Maximum Deviation Found: 1.587e-08 | |
132 | // Expected Error Term : 1.531e-08 | |
133 | // Maximum Relative Change in Control Points : 9.064e-02 | |
134 | // Max Error found at float precision = Poly : 5.065020e-08 | |
135 | ||
136 | static const T P[] = | |
137 | { | |
92f5a8d4 TL |
138 | 2.533141220e-01f, |
139 | 5.221502603e-01f, | |
140 | 6.380180669e-02f, | |
141 | -5.934976547e-02f | |
b32b8144 FG |
142 | }; |
143 | static const T Q[] = | |
144 | { | |
92f5a8d4 TL |
145 | 1.000000000e+00f, |
146 | 2.679722431e+00f, | |
147 | 1.561635813e+00f, | |
148 | 1.573660661e-01f | |
b32b8144 FG |
149 | }; |
150 | if(x < tools::log_max_value<T>()) | |
151 | return ((tools::evaluate_rational(P, Q, T(1 / x)) + 1) * exp(-x) / sqrt(x)); | |
152 | else | |
153 | { | |
154 | T ex = exp(-x / 2); | |
155 | return ((tools::evaluate_rational(P, Q, T(1 / x)) + 1) * ex / sqrt(x)) * ex; | |
156 | } | |
157 | } | |
158 | } | |
159 | ||
160 | template <typename T> | |
161 | T bessel_k0_imp(const T& x, const mpl::int_<53>&) | |
162 | { | |
163 | BOOST_MATH_STD_USING | |
164 | if(x <= 1) | |
165 | { | |
166 | // Maximum Deviation Found: 6.077e-17 | |
167 | // Expected Error Term : -6.077e-17 | |
168 | // Maximum Relative Change in Control Points : 7.797e-02 | |
169 | // Max Error found at double precision = Poly : 1.003156e-16 | |
170 | static const T Y = 1.137250900268554688; | |
171 | static const T P[] = | |
172 | { | |
173 | -1.372509002685546267e-01, | |
174 | 2.574916117833312855e-01, | |
175 | 1.395474602146869316e-02, | |
176 | 5.445476986653926759e-04, | |
177 | 7.125159422136622118e-06 | |
178 | }; | |
179 | static const T Q[] = | |
180 | { | |
181 | 1.000000000000000000e+00, | |
182 | -5.458333438017788530e-02, | |
183 | 1.291052816975251298e-03, | |
184 | -1.367653946978586591e-05 | |
185 | }; | |
186 | ||
187 | T a = x * x / 4; | |
188 | a = (tools::evaluate_polynomial(P, a) / tools::evaluate_polynomial(Q, a) + Y) * a + 1; | |
189 | ||
190 | // Maximum Deviation Found: 3.429e-18 | |
191 | // Expected Error Term : 3.392e-18 | |
192 | // Maximum Relative Change in Control Points : 2.041e-02 | |
193 | // Max Error found at double precision = Poly : 2.513112e-16 | |
194 | static const T P2[] = | |
195 | { | |
196 | 1.159315156584124484e-01, | |
197 | 2.789828789146031732e-01, | |
198 | 2.524892993216121934e-02, | |
199 | 8.460350907213637784e-04, | |
200 | 1.491471924309617534e-05, | |
201 | 1.627106892422088488e-07, | |
202 | 1.208266102392756055e-09, | |
203 | 6.611686391749704310e-12 | |
204 | }; | |
205 | ||
206 | return tools::evaluate_polynomial(P2, T(x * x)) - log(x) * a; | |
207 | } | |
208 | else | |
209 | { | |
210 | // Maximum Deviation Found: 4.316e-17 | |
211 | // Expected Error Term : 9.570e-18 | |
212 | // Maximum Relative Change in Control Points : 2.757e-01 | |
213 | // Max Error found at double precision = Poly : 1.001560e-16 | |
214 | ||
215 | static const T Y = 1; | |
216 | static const T P[] = | |
217 | { | |
218 | 2.533141373155002416e-01, | |
219 | 3.628342133984595192e+00, | |
220 | 1.868441889406606057e+01, | |
221 | 4.306243981063412784e+01, | |
222 | 4.424116209627428189e+01, | |
223 | 1.562095339356220468e+01, | |
224 | -1.810138978229410898e+00, | |
225 | -1.414237994269995877e+00, | |
226 | -9.369168119754924625e-02 | |
227 | }; | |
228 | static const T Q[] = | |
229 | { | |
230 | 1.000000000000000000e+00, | |
231 | 1.494194694879908328e+01, | |
232 | 8.265296455388554217e+01, | |
233 | 2.162779506621866970e+02, | |
234 | 2.845145155184222157e+02, | |
235 | 1.851714491916334995e+02, | |
236 | 5.486540717439723515e+01, | |
237 | 6.118075837628957015e+00, | |
238 | 1.586261269326235053e-01 | |
239 | }; | |
240 | if(x < tools::log_max_value<T>()) | |
241 | return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x)); | |
242 | else | |
243 | { | |
244 | T ex = exp(-x / 2); | |
245 | return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex; | |
246 | } | |
247 | } | |
248 | } | |
249 | ||
250 | template <typename T> | |
251 | T bessel_k0_imp(const T& x, const mpl::int_<64>&) | |
252 | { | |
253 | BOOST_MATH_STD_USING | |
254 | if(x <= 1) | |
255 | { | |
256 | // Maximum Deviation Found: 2.180e-22 | |
257 | // Expected Error Term : 2.180e-22 | |
258 | // Maximum Relative Change in Control Points : 2.943e-01 | |
259 | // Max Error found at float80 precision = Poly : 3.923207e-20 | |
260 | static const T Y = 1.137250900268554687500e+00; | |
261 | static const T P[] = | |
262 | { | |
263 | BOOST_MATH_BIG_CONSTANT(T, 64, -1.372509002685546875002e-01), | |
264 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.566481981037407600436e-01), | |
265 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.551881122448948854873e-02), | |
266 | BOOST_MATH_BIG_CONSTANT(T, 64, 6.646112454323276529650e-04), | |
267 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.213747930378196492543e-05), | |
268 | BOOST_MATH_BIG_CONSTANT(T, 64, 9.423709328020389560844e-08) | |
269 | }; | |
270 | static const T Q[] = | |
271 | { | |
272 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00), | |
273 | BOOST_MATH_BIG_CONSTANT(T, 64, -4.843828412587773008342e-02), | |
274 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.088484822515098936140e-03), | |
275 | BOOST_MATH_BIG_CONSTANT(T, 64, -1.374724008530702784829e-05), | |
276 | BOOST_MATH_BIG_CONSTANT(T, 64, 8.452665455952581680339e-08) | |
277 | }; | |
278 | ||
279 | ||
280 | T a = x * x / 4; | |
281 | a = (tools::evaluate_polynomial(P, a) / tools::evaluate_polynomial(Q, a) + Y) * a + 1; | |
282 | ||
283 | // Maximum Deviation Found: 2.440e-21 | |
284 | // Expected Error Term : -2.434e-21 | |
285 | // Maximum Relative Change in Control Points : 2.459e-02 | |
286 | // Max Error found at float80 precision = Poly : 1.482487e-19 | |
287 | static const T P2[] = | |
288 | { | |
289 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.159315156584124488110e-01), | |
290 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.764832791416047889734e-01), | |
291 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.926062887220923354112e-02), | |
292 | BOOST_MATH_BIG_CONSTANT(T, 64, 3.660777862036966089410e-04), | |
293 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.094942446930673386849e-06) | |
294 | }; | |
295 | static const T Q2[] = | |
296 | { | |
297 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00), | |
298 | BOOST_MATH_BIG_CONSTANT(T, 64, -2.156100313881251616320e-02), | |
299 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.315993873344905957033e-04), | |
300 | BOOST_MATH_BIG_CONSTANT(T, 64, -1.529444499350703363451e-06), | |
301 | BOOST_MATH_BIG_CONSTANT(T, 64, 5.524988589917857531177e-09) | |
302 | }; | |
303 | return tools::evaluate_rational(P2, Q2, T(x * x)) - log(x) * a; | |
304 | } | |
305 | else | |
306 | { | |
307 | // Maximum Deviation Found: 4.291e-20 | |
308 | // Expected Error Term : 2.236e-21 | |
309 | // Maximum Relative Change in Control Points : 3.021e-01 | |
310 | //Max Error found at float80 precision = Poly : 8.727378e-20 | |
311 | static const T Y = 1; | |
312 | static const T P[] = | |
313 | { | |
314 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.533141373155002512056e-01), | |
315 | BOOST_MATH_BIG_CONSTANT(T, 64, 5.417942070721928652715e+00), | |
316 | BOOST_MATH_BIG_CONSTANT(T, 64, 4.477464607463971754433e+01), | |
317 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.838745728725943889876e+02), | |
318 | BOOST_MATH_BIG_CONSTANT(T, 64, 4.009736314927811202517e+02), | |
319 | BOOST_MATH_BIG_CONSTANT(T, 64, 4.557411293123609803452e+02), | |
320 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.360222564015361268955e+02), | |
321 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.385435333168505701022e+01), | |
322 | BOOST_MATH_BIG_CONSTANT(T, 64, -1.750195760942181592050e+01), | |
323 | BOOST_MATH_BIG_CONSTANT(T, 64, -4.059789241612946683713e+00), | |
324 | BOOST_MATH_BIG_CONSTANT(T, 64, -1.612783121537333908889e-01) | |
325 | }; | |
326 | static const T Q[] = | |
327 | { | |
328 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00), | |
329 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.200669254769325861404e+01), | |
330 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.900177593527144126549e+02), | |
331 | BOOST_MATH_BIG_CONSTANT(T, 64, 8.361003989965786932682e+02), | |
332 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.041319870804843395893e+03), | |
333 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.828491555113790345068e+03), | |
334 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.190342229261529076624e+03), | |
335 | BOOST_MATH_BIG_CONSTANT(T, 64, 9.003330795963812219852e+02), | |
336 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.773371397243777891569e+02), | |
337 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.368634935531158398439e+01), | |
338 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.543310879400359967327e-01) | |
339 | }; | |
340 | if(x < tools::log_max_value<T>()) | |
341 | return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x)); | |
342 | else | |
343 | { | |
344 | T ex = exp(-x / 2); | |
345 | return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex; | |
346 | } | |
347 | } | |
348 | } | |
349 | ||
350 | template <typename T> | |
351 | T bessel_k0_imp(const T& x, const mpl::int_<113>&) | |
352 | { | |
353 | BOOST_MATH_STD_USING | |
354 | if(x <= 1) | |
355 | { | |
356 | // Maximum Deviation Found: 5.682e-37 | |
357 | // Expected Error Term : 5.682e-37 | |
358 | // Maximum Relative Change in Control Points : 6.094e-04 | |
359 | // Max Error found at float128 precision = Poly : 5.338213e-35 | |
360 | static const T Y = 1.137250900268554687500000000000000000e+00f; | |
361 | static const T P[] = | |
362 | { | |
363 | BOOST_MATH_BIG_CONSTANT(T, 113, -1.372509002685546875000000000000000006e-01), | |
364 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.556212905071072782462974351698081303e-01), | |
365 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.742459135264203478530904179889103929e-02), | |
366 | BOOST_MATH_BIG_CONSTANT(T, 113, 8.077860530453688571555479526961318918e-04), | |
367 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.868173911669241091399374307788635148e-05), | |
368 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.496405768838992243478709145123306602e-07), | |
369 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.752489221949580551692915881999762125e-09), | |
370 | BOOST_MATH_BIG_CONSTANT(T, 113, 5.243010555737173524710512824955368526e-12) | |
371 | }; | |
372 | static const T Q[] = | |
373 | { | |
374 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00), | |
375 | BOOST_MATH_BIG_CONSTANT(T, 113, -4.095631064064621099785696980653193721e-02), | |
376 | BOOST_MATH_BIG_CONSTANT(T, 113, 8.313880983725212151967078809725835532e-04), | |
377 | BOOST_MATH_BIG_CONSTANT(T, 113, -1.095229912293480063501285562382835142e-05), | |
378 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.022828799511943141130509410251996277e-07), | |
379 | BOOST_MATH_BIG_CONSTANT(T, 113, -6.860874007419812445494782795829046836e-10), | |
380 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.107297802344970725756092082686799037e-12), | |
381 | BOOST_MATH_BIG_CONSTANT(T, 113, -7.460529579244623559164763757787600944e-15) | |
382 | }; | |
383 | T a = x * x / 4; | |
384 | a = (tools::evaluate_rational(P, Q, a) + Y) * a + 1; | |
385 | ||
386 | // Maximum Deviation Found: 5.173e-38 | |
387 | // Expected Error Term : 5.105e-38 | |
388 | // Maximum Relative Change in Control Points : 9.734e-03 | |
389 | // Max Error found at float128 precision = Poly : 1.688806e-34 | |
390 | static const T P2[] = | |
391 | { | |
392 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.159315156584124488107200313757741370e-01), | |
393 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.789828789146031122026800078439435369e-01), | |
394 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.524892993216269451266750049024628432e-02), | |
395 | BOOST_MATH_BIG_CONSTANT(T, 113, 8.460350907082229957222453839935101823e-04), | |
396 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.491471929926042875260452849503857976e-05), | |
397 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.627105610481598430816014719558896866e-07), | |
398 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.208426165007797264194914898538250281e-09), | |
399 | BOOST_MATH_BIG_CONSTANT(T, 113, 6.508697838747354949164182457073784117e-12), | |
400 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.659784680639805301101014383907273109e-14), | |
401 | BOOST_MATH_BIG_CONSTANT(T, 113, 8.531090131964391104248859415958109654e-17), | |
402 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.205195117066478034260323124669936314e-19), | |
403 | BOOST_MATH_BIG_CONSTANT(T, 113, 4.692219280289030165761119775783115426e-22), | |
404 | BOOST_MATH_BIG_CONSTANT(T, 113, 8.362350161092532344171965861545860747e-25), | |
405 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.277990623924628999539014980773738258e-27) | |
406 | }; | |
407 | ||
408 | return tools::evaluate_polynomial(P2, T(x * x)) - log(x) * a; | |
409 | } | |
410 | else | |
411 | { | |
412 | // Maximum Deviation Found: 1.462e-34 | |
413 | // Expected Error Term : 4.917e-40 | |
414 | // Maximum Relative Change in Control Points : 3.385e-01 | |
415 | // Max Error found at float128 precision = Poly : 1.567573e-34 | |
416 | static const T Y = 1; | |
417 | static const T P[] = | |
418 | { | |
419 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.533141373155002512078826424055226265e-01), | |
420 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.001949740768235770078339977110749204e+01), | |
421 | BOOST_MATH_BIG_CONSTANT(T, 113, 6.991516715983883248363351472378349986e+02), | |
422 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.429587951594593159075690819360687720e+04), | |
423 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.911933815201948768044660065771258450e+05), | |
424 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.769943016204926614862175317962439875e+06), | |
425 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.170866154649560750500954150401105606e+07), | |
426 | BOOST_MATH_BIG_CONSTANT(T, 113, 5.634687099724383996792011977705727661e+07), | |
427 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.989524036456492581597607246664394014e+08), | |
428 | BOOST_MATH_BIG_CONSTANT(T, 113, 5.160394785715328062088529400178080360e+08), | |
429 | BOOST_MATH_BIG_CONSTANT(T, 113, 9.778173054417826368076483100902201433e+08), | |
430 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.335667778588806892764139643950439733e+09), | |
431 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.283635100080306980206494425043706838e+09), | |
432 | BOOST_MATH_BIG_CONSTANT(T, 113, 8.300616188213640626577036321085025855e+08), | |
433 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.277591957076162984986406540894621482e+08), | |
434 | BOOST_MATH_BIG_CONSTANT(T, 113, 5.564360536834214058158565361486115932e+07), | |
435 | BOOST_MATH_BIG_CONSTANT(T, 113, -1.043505161612403359098596828115690596e+07), | |
436 | BOOST_MATH_BIG_CONSTANT(T, 113, -7.217035248223503605127967970903027314e+06), | |
437 | BOOST_MATH_BIG_CONSTANT(T, 113, -1.422938158797326748375799596769964430e+06), | |
438 | BOOST_MATH_BIG_CONSTANT(T, 113, -1.229125746200586805278634786674745210e+05), | |
439 | BOOST_MATH_BIG_CONSTANT(T, 113, -4.201632288615609937883545928660649813e+03), | |
440 | BOOST_MATH_BIG_CONSTANT(T, 113, -3.690820607338480548346746717311811406e+01) | |
441 | }; | |
442 | static const T Q[] = | |
443 | { | |
444 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00), | |
445 | BOOST_MATH_BIG_CONSTANT(T, 113, 7.964877874035741452203497983642653107e+01), | |
446 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.808929943826193766839360018583294769e+03), | |
447 | BOOST_MATH_BIG_CONSTANT(T, 113, 5.814524004679994110944366890912384139e+04), | |
448 | BOOST_MATH_BIG_CONSTANT(T, 113, 7.897794522506725610540209610337355118e+05), | |
449 | BOOST_MATH_BIG_CONSTANT(T, 113, 7.456339470955813675629523617440433672e+06), | |
450 | BOOST_MATH_BIG_CONSTANT(T, 113, 5.057818717813969772198911392875127212e+07), | |
451 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.513821619536852436424913886081133209e+08), | |
452 | BOOST_MATH_BIG_CONSTANT(T, 113, 9.255938846873380596038513316919990776e+08), | |
453 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.537077551699028079347581816919572141e+09), | |
454 | BOOST_MATH_BIG_CONSTANT(T, 113, 5.176769339768120752974843214652367321e+09), | |
455 | BOOST_MATH_BIG_CONSTANT(T, 113, 7.828722317390455845253191337207432060e+09), | |
456 | BOOST_MATH_BIG_CONSTANT(T, 113, 8.698864296569996402006511705803675890e+09), | |
457 | BOOST_MATH_BIG_CONSTANT(T, 113, 7.007803261356636409943826918468544629e+09), | |
458 | BOOST_MATH_BIG_CONSTANT(T, 113, 4.016564631288740308993071395104715469e+09), | |
459 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.595893010619754750655947035567624730e+09), | |
460 | BOOST_MATH_BIG_CONSTANT(T, 113, 4.241241839120481076862742189989406856e+08), | |
461 | BOOST_MATH_BIG_CONSTANT(T, 113, 7.168778094393076220871007550235840858e+07), | |
462 | BOOST_MATH_BIG_CONSTANT(T, 113, 7.156200301360388147635052029404211109e+06), | |
463 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.752130382550379886741949463587008794e+05), | |
464 | BOOST_MATH_BIG_CONSTANT(T, 113, 8.370574966987293592457152146806662562e+03), | |
465 | BOOST_MATH_BIG_CONSTANT(T, 113, 4.871254714311063594080644835895740323e+01) | |
466 | }; | |
467 | if(x < tools::log_max_value<T>()) | |
468 | return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x)); | |
469 | else | |
470 | { | |
471 | T ex = exp(-x / 2); | |
472 | return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex; | |
473 | } | |
474 | } | |
475 | } | |
476 | ||
477 | template <typename T> | |
478 | T bessel_k0_imp(const T& x, const mpl::int_<0>&) | |
479 | { | |
480 | if(boost::math::tools::digits<T>() <= 24) | |
481 | return bessel_k0_imp(x, mpl::int_<24>()); | |
482 | else if(boost::math::tools::digits<T>() <= 53) | |
483 | return bessel_k0_imp(x, mpl::int_<53>()); | |
484 | else if(boost::math::tools::digits<T>() <= 64) | |
485 | return bessel_k0_imp(x, mpl::int_<64>()); | |
486 | else if(boost::math::tools::digits<T>() <= 113) | |
487 | return bessel_k0_imp(x, mpl::int_<113>()); | |
488 | BOOST_ASSERT(0); | |
489 | return 0; | |
490 | } | |
491 | ||
492 | template <typename T> | |
493 | inline T bessel_k0(const T& x) | |
494 | { | |
495 | typedef mpl::int_< | |
496 | ((std::numeric_limits<T>::digits == 0) || (std::numeric_limits<T>::radix != 2)) ? | |
497 | 0 : | |
498 | std::numeric_limits<T>::digits <= 24 ? | |
499 | 24 : | |
500 | std::numeric_limits<T>::digits <= 53 ? | |
501 | 53 : | |
502 | std::numeric_limits<T>::digits <= 64 ? | |
503 | 64 : | |
504 | std::numeric_limits<T>::digits <= 113 ? | |
505 | 113 : -1 | |
506 | > tag_type; | |
507 | ||
508 | bessel_k0_initializer<T, tag_type>::force_instantiate(); | |
509 | return bessel_k0_imp(x, tag_type()); | |
510 | } | |
511 | ||
512 | }}} // namespaces | |
513 | ||
514 | #ifdef _MSC_VER | |
515 | #pragma warning(pop) | |
516 | #endif | |
517 | ||
518 | #endif // BOOST_MATH_BESSEL_K0_HPP | |
519 |