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
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/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
44namespace boost { namespace math { namespace detail{
45
46template <typename T>
47T bessel_k0(const T& x);
48
49template <class T, class tag>
50struct bessel_k0_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_k0(T(0.5));
61 bessel_k0(T(1.5));
62 }
63 static void do_init(const boost::integral_constant<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
79template <class T, class tag>
80const typename bessel_k0_initializer<T, tag>::init bessel_k0_initializer<T, tag>::initializer;
81
82
83template <typename T, int N>
84T bessel_k0_imp(const T& x, const boost::integral_constant<int, N>&)
85{
86 BOOST_ASSERT(0);
87 return 0;
88}
89
90template <typename T>
91T bessel_k0_imp(const T& x, const boost::integral_constant<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 {
138 2.533141220e-01f,
139 5.221502603e-01f,
140 6.380180669e-02f,
141 -5.934976547e-02f
142 };
143 static const T Q[] =
144 {
145 1.000000000e+00f,
146 2.679722431e+00f,
147 1.561635813e+00f,
148 1.573660661e-01f
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
160template <typename T>
161T bessel_k0_imp(const T& x, const boost::integral_constant<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
250template <typename T>
251T bessel_k0_imp(const T& x, const boost::integral_constant<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
350template <typename T>
351T bessel_k0_imp(const T& x, const boost::integral_constant<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
477template <typename T>
478T bessel_k0_imp(const T& x, const boost::integral_constant<int, 0>&)
479{
480 if(boost::math::tools::digits<T>() <= 24)
481 return bessel_k0_imp(x, boost::integral_constant<int, 24>());
482 else if(boost::math::tools::digits<T>() <= 53)
483 return bessel_k0_imp(x, boost::integral_constant<int, 53>());
484 else if(boost::math::tools::digits<T>() <= 64)
485 return bessel_k0_imp(x, boost::integral_constant<int, 64>());
486 else if(boost::math::tools::digits<T>() <= 113)
487 return bessel_k0_imp(x, boost::integral_constant<int, 113>());
488 BOOST_ASSERT(0);
489 return 0;
490}
491
492template <typename T>
493inline T bessel_k0(const T& x)
494{
495 typedef boost::integral_constant<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
520

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