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
44namespace 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

source code of include/boost/math/special_functions/detail/bessel_k1.hpp