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_K1_HPP |
8 | #define BOOST_MATH_BESSEL_K1_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 | |
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 diagnostic ignored work :( |
27 | // |
28 | #pragma GCC system_header |
29 | #endif |
30 | |
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/2016/01/05/rational-approximations-for-the-modified-bessel-function-of-the-second-kind-k1-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_k1(const T& x); |
48 | |
49 | template <class T, class tag> |
50 | struct bessel_k1_initializer |
51 | { |
52 | struct init |
53 | { |
54 | init() |
55 | { |
56 | do_init(tag()); |
57 | } |
58 | static void do_init(const boost::integral_constant<int, 113>&) |
59 | { |
60 | bessel_k1(T(0.5)); |
61 | bessel_k1(T(2)); |
62 | bessel_k1(T(6)); |
63 | } |
64 | static void do_init(const boost::integral_constant<int, 64>&) |
65 | { |
66 | bessel_k1(T(0.5)); |
67 | bessel_k1(T(6)); |
68 | } |
69 | template <class U> |
70 | static void do_init(const U&) {} |
71 | void force_instantiate()const {} |
72 | }; |
73 | static const init initializer; |
74 | static void force_instantiate() |
75 | { |
76 | initializer.force_instantiate(); |
77 | } |
78 | }; |
79 | |
80 | template <class T, class tag> |
81 | const typename bessel_k1_initializer<T, tag>::init bessel_k1_initializer<T, tag>::initializer; |
82 | |
83 | |
84 | template <typename T, int N> |
85 | inline T bessel_k1_imp(const T& x, const boost::integral_constant<int, N>&) |
86 | { |
87 | BOOST_ASSERT(0); |
88 | return 0; |
89 | } |
90 | |
91 | template <typename T> |
92 | T bessel_k1_imp(const T& x, const boost::integral_constant<int, 24>&) |
93 | { |
94 | BOOST_MATH_STD_USING |
95 | if(x <= 1) |
96 | { |
97 | // Maximum Deviation Found: 3.090e-12 |
98 | // Expected Error Term : -3.053e-12 |
99 | // Maximum Relative Change in Control Points : 4.927e-02 |
100 | // Max Error found at float precision = Poly : 7.918347e-10 |
101 | static const T Y = 8.695471287e-02f; |
102 | static const T P[] = |
103 | { |
104 | -3.621379531e-03f, |
105 | 7.131781976e-03f, |
106 | -1.535278300e-05f |
107 | }; |
108 | static const T Q[] = |
109 | { |
110 | 1.000000000e+00f, |
111 | -5.173102701e-02f, |
112 | 9.203530671e-04f |
113 | }; |
114 | |
115 | T a = x * x / 4; |
116 | a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2; |
117 | |
118 | // Maximum Deviation Found: 3.556e-08 |
119 | // Expected Error Term : -3.541e-08 |
120 | // Maximum Relative Change in Control Points : 8.203e-02 |
121 | static const T P2[] = |
122 | { |
123 | -3.079657469e-01f, |
124 | -8.537108913e-02f, |
125 | -4.640275408e-03f, |
126 | -1.156442414e-04f |
127 | }; |
128 | |
129 | return tools::evaluate_polynomial(P2, T(x * x)) * x + 1 / x + log(x) * a; |
130 | } |
131 | else |
132 | { |
133 | // Maximum Deviation Found: 3.369e-08 |
134 | // Expected Error Term : -3.227e-08 |
135 | // Maximum Relative Change in Control Points : 9.917e-02 |
136 | // Max Error found at float precision = Poly : 6.084411e-08 |
137 | static const T Y = 1.450342178f; |
138 | static const T P[] = |
139 | { |
140 | -1.970280088e-01f, |
141 | 2.188747807e-02f, |
142 | 7.270394756e-01f, |
143 | 2.490678196e-01f |
144 | }; |
145 | static const T Q[] = |
146 | { |
147 | 1.000000000e+00f, |
148 | 2.274292882e+00f, |
149 | 9.904984851e-01f, |
150 | 4.585534549e-02f |
151 | }; |
152 | if(x < tools::log_max_value<T>()) |
153 | return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x)); |
154 | else |
155 | { |
156 | T ex = exp(-x / 2); |
157 | return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex; |
158 | } |
159 | } |
160 | } |
161 | |
162 | template <typename T> |
163 | T bessel_k1_imp(const T& x, const boost::integral_constant<int, 53>&) |
164 | { |
165 | BOOST_MATH_STD_USING |
166 | if(x <= 1) |
167 | { |
168 | // Maximum Deviation Found: 1.922e-17 |
169 | // Expected Error Term : 1.921e-17 |
170 | // Maximum Relative Change in Control Points : 5.287e-03 |
171 | // Max Error found at double precision = Poly : 2.004747e-17 |
172 | static const T Y = 8.69547128677368164e-02f; |
173 | static const T P[] = |
174 | { |
175 | -3.62137953440350228e-03, |
176 | 7.11842087490330300e-03, |
177 | 1.00302560256614306e-05, |
178 | 1.77231085381040811e-06 |
179 | }; |
180 | static const T Q[] = |
181 | { |
182 | 1.00000000000000000e+00, |
183 | -4.80414794429043831e-02, |
184 | 9.85972641934416525e-04, |
185 | -8.91196859397070326e-06 |
186 | }; |
187 | |
188 | T a = x * x / 4; |
189 | a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2; |
190 | |
191 | // Maximum Deviation Found: 4.053e-17 |
192 | // Expected Error Term : -4.053e-17 |
193 | // Maximum Relative Change in Control Points : 3.103e-04 |
194 | // Max Error found at double precision = Poly : 1.246698e-16 |
195 | |
196 | static const T P2[] = |
197 | { |
198 | -3.07965757829206184e-01, |
199 | -7.80929703673074907e-02, |
200 | -2.70619343754051620e-03, |
201 | -2.49549522229072008e-05 |
202 | }; |
203 | static const T Q2[] = |
204 | { |
205 | 1.00000000000000000e+00, |
206 | -2.36316836412163098e-02, |
207 | 2.64524577525962719e-04, |
208 | -1.49749618004162787e-06 |
209 | }; |
210 | |
211 | return tools::evaluate_rational(P2, Q2, T(x * x)) * x + 1 / x + log(x) * a; |
212 | } |
213 | else |
214 | { |
215 | // Maximum Deviation Found: 8.883e-17 |
216 | // Expected Error Term : -1.641e-17 |
217 | // Maximum Relative Change in Control Points : 2.786e-01 |
218 | // Max Error found at double precision = Poly : 1.258798e-16 |
219 | |
220 | static const T Y = 1.45034217834472656f; |
221 | static const T P[] = |
222 | { |
223 | -1.97028041029226295e-01, |
224 | -2.32408961548087617e+00, |
225 | -7.98269784507699938e+00, |
226 | -2.39968410774221632e+00, |
227 | 3.28314043780858713e+01, |
228 | 5.67713761158496058e+01, |
229 | 3.30907788466509823e+01, |
230 | 6.62582288933739787e+00, |
231 | 3.08851840645286691e-01 |
232 | }; |
233 | static const T Q[] = |
234 | { |
235 | 1.00000000000000000e+00, |
236 | 1.41811409298826118e+01, |
237 | 7.35979466317556420e+01, |
238 | 1.77821793937080859e+02, |
239 | 2.11014501598705982e+02, |
240 | 1.19425262951064454e+02, |
241 | 2.88448064302447607e+01, |
242 | 2.27912927104139732e+00, |
243 | 2.50358186953478678e-02 |
244 | }; |
245 | if(x < tools::log_max_value<T>()) |
246 | return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x)); |
247 | else |
248 | { |
249 | T ex = exp(-x / 2); |
250 | return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex; |
251 | } |
252 | } |
253 | } |
254 | |
255 | template <typename T> |
256 | T bessel_k1_imp(const T& x, const boost::integral_constant<int, 64>&) |
257 | { |
258 | BOOST_MATH_STD_USING |
259 | if(x <= 1) |
260 | { |
261 | // Maximum Deviation Found: 5.549e-23 |
262 | // Expected Error Term : -5.548e-23 |
263 | // Maximum Relative Change in Control Points : 2.002e-03 |
264 | // Max Error found at float80 precision = Poly : 9.352785e-22 |
265 | static const T Y = 8.695471286773681640625e-02f; |
266 | static const T P[] = |
267 | { |
268 | BOOST_MATH_BIG_CONSTANT(T, 64, -3.621379534403483072861e-03), |
269 | BOOST_MATH_BIG_CONSTANT(T, 64, 7.102135866103952705932e-03), |
270 | BOOST_MATH_BIG_CONSTANT(T, 64, 4.167545240236717601167e-05), |
271 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.537484002571894870830e-06), |
272 | BOOST_MATH_BIG_CONSTANT(T, 64, 6.603228256820000135990e-09) |
273 | }; |
274 | static const T Q[] = |
275 | { |
276 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00), |
277 | BOOST_MATH_BIG_CONSTANT(T, 64, -4.354457194045068370363e-02), |
278 | BOOST_MATH_BIG_CONSTANT(T, 64, 8.709137201220209072820e-04), |
279 | BOOST_MATH_BIG_CONSTANT(T, 64, -9.676151796359590545143e-06), |
280 | BOOST_MATH_BIG_CONSTANT(T, 64, 5.162715192766245311659e-08) |
281 | }; |
282 | |
283 | T a = x * x / 4; |
284 | a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2; |
285 | |
286 | // Maximum Deviation Found: 1.995e-23 |
287 | // Expected Error Term : 1.995e-23 |
288 | // Maximum Relative Change in Control Points : 8.174e-04 |
289 | // Max Error found at float80 precision = Poly : 4.137325e-20 |
290 | static const T P2[] = |
291 | { |
292 | BOOST_MATH_BIG_CONSTANT(T, 64, -3.079657578292062244054e-01), |
293 | BOOST_MATH_BIG_CONSTANT(T, 64, -7.963049154965966503231e-02), |
294 | BOOST_MATH_BIG_CONSTANT(T, 64, -3.103277523735639924895e-03), |
295 | BOOST_MATH_BIG_CONSTANT(T, 64, -4.023052834702215699504e-05), |
296 | BOOST_MATH_BIG_CONSTANT(T, 64, -1.719459155018493821839e-07) |
297 | }; |
298 | static const T Q2[] = |
299 | { |
300 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00), |
301 | BOOST_MATH_BIG_CONSTANT(T, 64, -1.863917670410152669768e-02), |
302 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.699367098849735298090e-04), |
303 | BOOST_MATH_BIG_CONSTANT(T, 64, -9.309358790546076298429e-07), |
304 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.708893480271612711933e-09) |
305 | }; |
306 | |
307 | return tools::evaluate_rational(P2, Q2, T(x * x)) * x + 1 / x + log(x) * a; |
308 | } |
309 | else |
310 | { |
311 | // Maximum Deviation Found: 9.785e-20 |
312 | // Expected Error Term : -3.302e-21 |
313 | // Maximum Relative Change in Control Points : 3.432e-01 |
314 | // Max Error found at float80 precision = Poly : 1.083755e-19 |
315 | static const T Y = 1.450342178344726562500e+00f; |
316 | static const T P[] = |
317 | { |
318 | BOOST_MATH_BIG_CONSTANT(T, 64, -1.970280410292263112917e-01), |
319 | BOOST_MATH_BIG_CONSTANT(T, 64, -4.058564803062959169322e+00), |
320 | BOOST_MATH_BIG_CONSTANT(T, 64, -3.036658174194917777473e+01), |
321 | BOOST_MATH_BIG_CONSTANT(T, 64, -9.576825392332820142173e+01), |
322 | BOOST_MATH_BIG_CONSTANT(T, 64, -6.706969489248020941949e+01), |
323 | BOOST_MATH_BIG_CONSTANT(T, 64, 3.264572499406168221382e+02), |
324 | BOOST_MATH_BIG_CONSTANT(T, 64, 8.584972047303151034100e+02), |
325 | BOOST_MATH_BIG_CONSTANT(T, 64, 8.422082733280017909550e+02), |
326 | BOOST_MATH_BIG_CONSTANT(T, 64, 3.738005441471368178383e+02), |
327 | BOOST_MATH_BIG_CONSTANT(T, 64, 7.016938390144121276609e+01), |
328 | BOOST_MATH_BIG_CONSTANT(T, 64, 4.319614662598089438939e+00), |
329 | BOOST_MATH_BIG_CONSTANT(T, 64, 3.710715864316521856193e-02) |
330 | }; |
331 | static const T Q[] = |
332 | { |
333 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00), |
334 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.298433045824439052398e+01), |
335 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.082047745067709230037e+02), |
336 | BOOST_MATH_BIG_CONSTANT(T, 64, 9.662367854250262046592e+02), |
337 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.504148628460454004686e+03), |
338 | BOOST_MATH_BIG_CONSTANT(T, 64, 3.712730364911389908905e+03), |
339 | BOOST_MATH_BIG_CONSTANT(T, 64, 3.108002081150068641112e+03), |
340 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.400149940532448553143e+03), |
341 | BOOST_MATH_BIG_CONSTANT(T, 64, 3.083303048095846226299e+02), |
342 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.748706060530351833346e+01), |
343 | BOOST_MATH_BIG_CONSTANT(T, 64, 6.321900849331506946977e-01), |
344 | }; |
345 | if(x < tools::log_max_value<T>()) |
346 | return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * exp(-x) / sqrt(x)); |
347 | else |
348 | { |
349 | T ex = exp(-x / 2); |
350 | return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex; |
351 | } |
352 | } |
353 | } |
354 | |
355 | template <typename T> |
356 | T bessel_k1_imp(const T& x, const boost::integral_constant<int, 113>&) |
357 | { |
358 | BOOST_MATH_STD_USING |
359 | if(x <= 1) |
360 | { |
361 | // Maximum Deviation Found: 7.120e-35 |
362 | // Expected Error Term : -7.119e-35 |
363 | // Maximum Relative Change in Control Points : 1.207e-03 |
364 | // Max Error found at float128 precision = Poly : 7.143688e-35 |
365 | static const T Y = 8.695471286773681640625000000000000000e-02f; |
366 | static const T P[] = |
367 | { |
368 | BOOST_MATH_BIG_CONSTANT(T, 113, -3.621379534403483072916666666666595475e-03), |
369 | BOOST_MATH_BIG_CONSTANT(T, 113, 7.074117676930975433219826471336547627e-03), |
370 | BOOST_MATH_BIG_CONSTANT(T, 113, 9.631337631362776369069668419033041661e-05), |
371 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.468935967870048731821071646104412775e-06), |
372 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.956705020559599861444492614737168261e-08), |
373 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.347140307321161346703214099534250263e-10), |
374 | BOOST_MATH_BIG_CONSTANT(T, 113, 5.569608494081482873946791086435679661e-13) |
375 | }; |
376 | static const T Q[] = |
377 | { |
378 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00), |
379 | BOOST_MATH_BIG_CONSTANT(T, 113, -3.580768910152105375615558920428350204e-02), |
380 | BOOST_MATH_BIG_CONSTANT(T, 113, 6.197467671701485365363068445534557369e-04), |
381 | BOOST_MATH_BIG_CONSTANT(T, 113, -6.707466533308630411966030561446666237e-06), |
382 | BOOST_MATH_BIG_CONSTANT(T, 113, 4.846687802282250112624373388491123527e-08), |
383 | BOOST_MATH_BIG_CONSTANT(T, 113, -2.248493131151981569517383040323900343e-10), |
384 | BOOST_MATH_BIG_CONSTANT(T, 113, 5.319279786372775264555728921709381080e-13) |
385 | }; |
386 | |
387 | T a = x * x / 4; |
388 | a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2; |
389 | |
390 | // Maximum Deviation Found: 4.473e-37 |
391 | // Expected Error Term : 4.473e-37 |
392 | // Maximum Relative Change in Control Points : 8.550e-04 |
393 | // Max Error found at float128 precision = Poly : 8.167701e-35 |
394 | static const T P2[] = |
395 | { |
396 | BOOST_MATH_BIG_CONSTANT(T, 113, -3.079657578292062244053600156878870690e-01), |
397 | BOOST_MATH_BIG_CONSTANT(T, 113, -8.133183745732467770755578848987414875e-02), |
398 | BOOST_MATH_BIG_CONSTANT(T, 113, -3.548968792764174773125420229299431951e-03), |
399 | BOOST_MATH_BIG_CONSTANT(T, 113, -5.886125468718182876076972186152445490e-05), |
400 | BOOST_MATH_BIG_CONSTANT(T, 113, -4.506712111733707245745396404449639865e-07), |
401 | BOOST_MATH_BIG_CONSTANT(T, 113, -1.632502325880313239698965376754406011e-09), |
402 | BOOST_MATH_BIG_CONSTANT(T, 113, -2.311973065898784812266544485665624227e-12) |
403 | }; |
404 | static const T Q2[] = |
405 | { |
406 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00), |
407 | BOOST_MATH_BIG_CONSTANT(T, 113, -1.311471216733781016657962995723287450e-02), |
408 | BOOST_MATH_BIG_CONSTANT(T, 113, 8.571876054797365417068164018709472969e-05), |
409 | BOOST_MATH_BIG_CONSTANT(T, 113, -3.630181215268238731442496851497901293e-07), |
410 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.070176111227805048604885986867484807e-09), |
411 | BOOST_MATH_BIG_CONSTANT(T, 113, -2.129046580769872602793220056461084761e-12), |
412 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.294906469421390890762001971790074432e-15) |
413 | }; |
414 | |
415 | return tools::evaluate_rational(P2, Q2, T(x * x)) * x + 1 / x + log(x) * a; |
416 | } |
417 | else if(x < 4) |
418 | { |
419 | // Max error in interpolated form: 5.307e-37 |
420 | // Max Error found at float128 precision = Poly: 7.087862e-35 |
421 | static const T Y = 1.5023040771484375f; |
422 | static const T P[] = |
423 | { |
424 | BOOST_MATH_BIG_CONSTANT(T, 113, -2.489899398329369710528254347931380044e-01), |
425 | BOOST_MATH_BIG_CONSTANT(T, 113, -6.819080211203854781858815596508456873e+00), |
426 | BOOST_MATH_BIG_CONSTANT(T, 113, -7.599915699069767382647695624952723034e+01), |
427 | BOOST_MATH_BIG_CONSTANT(T, 113, -4.450211910821295507926582231071300718e+02), |
428 | BOOST_MATH_BIG_CONSTANT(T, 113, -1.451374687870925175794150513723956533e+03), |
429 | BOOST_MATH_BIG_CONSTANT(T, 113, -2.405805746895098802803503988539098226e+03), |
430 | BOOST_MATH_BIG_CONSTANT(T, 113, -5.638808326778389656403861103277220518e+02), |
431 | BOOST_MATH_BIG_CONSTANT(T, 113, 5.513958744081268456191778822780865708e+03), |
432 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.121301640926540743072258116122834804e+04), |
433 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.080094900175649541266613109971296190e+04), |
434 | BOOST_MATH_BIG_CONSTANT(T, 113, 5.896531083639613332407534434915552429e+03), |
435 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.856602122319645694042555107114028437e+03), |
436 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.237121918853145421414003823957537419e+02), |
437 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.842072954561323076230238664623893504e+01), |
438 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.039705646510167437971862966128055524e+00), |
439 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.008418100718254816100425022904039530e-02) |
440 | }; |
441 | static const T Q[] = |
442 | { |
443 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00), |
444 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.927456835239137986889227412815459529e+01), |
445 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.598985593265577043711382994516531273e+02), |
446 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.449897377085510281395819892689690579e+03), |
447 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.025555887684561913263090023158085327e+04), |
448 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.774140447181062463181892531100679195e+04), |
449 | BOOST_MATH_BIG_CONSTANT(T, 113, 4.962055507843204417243602332246120418e+04), |
450 | BOOST_MATH_BIG_CONSTANT(T, 113, 5.908269326976180183216954452196772931e+04), |
451 | BOOST_MATH_BIG_CONSTANT(T, 113, 4.655160454422016855911700790722577942e+04), |
452 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.383586885019548163464418964577684608e+04), |
453 | BOOST_MATH_BIG_CONSTANT(T, 113, 7.679920375586960324298491662159976419e+03), |
454 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.478586421028842906987799049804565008e+03), |
455 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.565384974896746094224942654383537090e+02), |
456 | BOOST_MATH_BIG_CONSTANT(T, 113, 7.902617937084010911005732488607114511e+00), |
457 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.429293010387921526110949911029094926e-01), |
458 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.880342607911083143560111853491047663e-04) |
459 | }; |
460 | return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * exp(-x) / sqrt(x)); |
461 | } |
462 | else |
463 | { |
464 | // Maximum Deviation Found: 4.359e-37 |
465 | // Expected Error Term : -6.565e-40 |
466 | // Maximum Relative Change in Control Points : 1.880e-01 |
467 | // Max Error found at float128 precision = Poly : 2.943572e-35 |
468 | static const T Y = 1.308816909790039062500000000000000000f; |
469 | static const T P[] = |
470 | { |
471 | BOOST_MATH_BIG_CONSTANT(T, 113, -5.550277247453881129211735759447737350e-02), |
472 | BOOST_MATH_BIG_CONSTANT(T, 113, -3.485883080219574328217554864956175929e+00), |
473 | BOOST_MATH_BIG_CONSTANT(T, 113, -8.903760658131484239300875153154881958e+01), |
474 | BOOST_MATH_BIG_CONSTANT(T, 113, -1.144813672213626237418235110712293337e+03), |
475 | BOOST_MATH_BIG_CONSTANT(T, 113, -6.498400501156131446691826557494158173e+03), |
476 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.573531831870363502604119835922166116e+04), |
477 | BOOST_MATH_BIG_CONSTANT(T, 113, 5.417416550054632009958262596048841154e+05), |
478 | BOOST_MATH_BIG_CONSTANT(T, 113, 4.271266450613557412825896604269130661e+06), |
479 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.898386013314389952534433455681107783e+07), |
480 | BOOST_MATH_BIG_CONSTANT(T, 113, 5.353798784656436259250791761023512750e+07), |
481 | BOOST_MATH_BIG_CONSTANT(T, 113, 9.839619195427352438957774052763490067e+07), |
482 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.169246368651532232388152442538005637e+08), |
483 | BOOST_MATH_BIG_CONSTANT(T, 113, 8.696368884166831199967845883371116431e+07), |
484 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.810226630422736458064005843327500169e+07), |
485 | BOOST_MATH_BIG_CONSTANT(T, 113, 8.854996610560406127438950635716757614e+06), |
486 | BOOST_MATH_BIG_CONSTANT(T, 113, 8.981057433937398731355768088809437625e+05), |
487 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.519440069856232098711793483639792952e+04) |
488 | }; |
489 | static const T Q[] = |
490 | { |
491 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00), |
492 | BOOST_MATH_BIG_CONSTANT(T, 113, 7.127348248283623146544565916604103560e+01), |
493 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.205092684176906740104488180754982065e+03), |
494 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.911249195069050636298346469740075758e+04), |
495 | BOOST_MATH_BIG_CONSTANT(T, 113, 4.426103406579046249654548481377792614e+05), |
496 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.365861555422488771286500241966208541e+06), |
497 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.765377714160383676864913709252529840e+07), |
498 | BOOST_MATH_BIG_CONSTANT(T, 113, 6.453822726931857253365138260720815246e+07), |
499 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.643207885048369990391975749439783892e+08), |
500 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.882540678243694621895816336640877878e+08), |
501 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.410120808992380266174106812005338148e+08), |
502 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.628138016559335882019310900426773027e+08), |
503 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.250794693811010646965360198541047961e+08), |
504 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.378723408195485594610593014072950078e+07), |
505 | BOOST_MATH_BIG_CONSTANT(T, 113, 4.488253856312453816451380319061865560e+06), |
506 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.202167197882689873967723350537104582e+05), |
507 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.673233230356966539460728211412989843e+03) |
508 | }; |
509 | if(x < tools::log_max_value<T>()) |
510 | return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * exp(-x) / sqrt(x)); |
511 | else |
512 | { |
513 | T ex = exp(-x / 2); |
514 | return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex; |
515 | } |
516 | } |
517 | } |
518 | |
519 | template <typename T> |
520 | T bessel_k1_imp(const T& x, const boost::integral_constant<int, 0>&) |
521 | { |
522 | if(boost::math::tools::digits<T>() <= 24) |
523 | return bessel_k1_imp(x, boost::integral_constant<int, 24>()); |
524 | else if(boost::math::tools::digits<T>() <= 53) |
525 | return bessel_k1_imp(x, boost::integral_constant<int, 53>()); |
526 | else if(boost::math::tools::digits<T>() <= 64) |
527 | return bessel_k1_imp(x, boost::integral_constant<int, 64>()); |
528 | else if(boost::math::tools::digits<T>() <= 113) |
529 | return bessel_k1_imp(x, boost::integral_constant<int, 113>()); |
530 | BOOST_ASSERT(0); |
531 | return 0; |
532 | } |
533 | |
534 | template <typename T> |
535 | inline T bessel_k1(const T& x) |
536 | { |
537 | typedef boost::integral_constant<int, |
538 | ((std::numeric_limits<T>::digits == 0) || (std::numeric_limits<T>::radix != 2)) ? |
539 | 0 : |
540 | std::numeric_limits<T>::digits <= 24 ? |
541 | 24 : |
542 | std::numeric_limits<T>::digits <= 53 ? |
543 | 53 : |
544 | std::numeric_limits<T>::digits <= 64 ? |
545 | 64 : |
546 | std::numeric_limits<T>::digits <= 113 ? |
547 | 113 : -1 |
548 | > tag_type; |
549 | |
550 | bessel_k1_initializer<T, tag_type>::force_instantiate(); |
551 | return bessel_k1_imp(x, tag_type()); |
552 | } |
553 | |
554 | }}} // namespaces |
555 | |
556 | #ifdef _MSC_VER |
557 | #pragma warning(pop) |
558 | #endif |
559 | |
560 | #endif // BOOST_MATH_BESSEL_K1_HPP |
561 | |
562 | |