486 lines · plain
1// Copyright (c) 2006 Xiaogang Zhang2// Copyright (c) 2017 John Maddock3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_BESSEL_K0_HPP8#define BOOST_MATH_BESSEL_K0_HPP9 10#ifdef _MSC_VER11#pragma once12#pragma warning(push)13#pragma warning(disable:4702) // Unreachable code (release mode only warning)14#endif15 16#include <boost/math/tools/config.hpp>17#include <boost/math/tools/type_traits.hpp>18#include <boost/math/tools/numeric_limits.hpp>19#include <boost/math/tools/precision.hpp>20#include <boost/math/tools/rational.hpp>21#include <boost/math/tools/big_constant.hpp>22#include <boost/math/tools/assert.hpp>23#include <boost/math/policies/error_handling.hpp>24 25#if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)26//27// This is the only way we can avoid28// warning: non-standard suffix on floating constant [-Wpedantic]29// when building with -Wall -pedantic. Neither __extension__30// nor #pragma diagnostic ignored work :(31//32#pragma GCC system_header33#endif34 35// Modified Bessel function of the second kind of order zero36// minimax rational approximations on intervals, see37// Russon and Blair, Chalk River Report AECL-3461, 1969,38// as revised by Pavel Holoborodko in "Rational Approximations 39// for the Modified Bessel Function of the Second Kind - K0(x) 40// for Computations with Double Precision", see 41// http://www.advanpix.com/2015/11/25/rational-approximations-for-the-modified-bessel-function-of-the-second-kind-k0-for-computations-with-double-precision/42//43// The actual coefficients used are our own derivation (by JM)44// since we extend to both greater and lesser precision than the45// references above. We can also improve performance WRT to46// Holoborodko without loss of precision.47 48namespace boost { namespace math { namespace detail{49 50template <typename T, int N>51BOOST_MATH_GPU_ENABLED T bessel_k0_imp(const T&, const boost::math::integral_constant<int, N>&)52{53 BOOST_MATH_ASSERT(0);54 return 0;55}56 57template <typename T>58BOOST_MATH_GPU_ENABLED T bessel_k0_imp(const T& x, const boost::math::integral_constant<int, 24>&)59{60 BOOST_MATH_STD_USING61 if(x <= 1)62 {63 // Maximum Deviation Found : 2.358e-0964 // Expected Error Term : -2.358e-0965 // Maximum Relative Change in Control Points : 9.552e-0266 // Max Error found at float precision = Poly : 4.448220e-0867 BOOST_MATH_STATIC const T Y = 1.137250900268554688f;68 BOOST_MATH_STATIC const T P[] = 69 {70 -1.372508979104259711e-01f,71 2.622545986273687617e-01f,72 5.047103728247919836e-03f73 };74 BOOST_MATH_STATIC const T Q[] = 75 {76 1.000000000000000000e+00f,77 -8.928694018000029415e-02f,78 2.985980684180969241e-03f79 };80 T a = x * x / 4;81 a = (tools::evaluate_rational(P, Q, a) + Y) * a + 1;82 83 // Maximum Deviation Found: 1.346e-0984 // Expected Error Term : -1.343e-0985 // Maximum Relative Change in Control Points : 2.405e-0286 // Max Error found at float precision = Poly : 1.354814e-0787 BOOST_MATH_STATIC const T P2[] = {88 1.159315158e-01f,89 2.789828686e-01f,90 2.524902861e-02f,91 8.457241514e-04f,92 1.530051997e-05f93 };94 return tools::evaluate_polynomial(P2, T(x * x)) - log(x) * a;95 }96 else97 {98 // Maximum Deviation Found: 1.587e-0899 // Expected Error Term : 1.531e-08100 // Maximum Relative Change in Control Points : 9.064e-02101 // Max Error found at float precision = Poly : 5.065020e-08102 103 BOOST_MATH_STATIC const T P[] =104 {105 2.533141220e-01f,106 5.221502603e-01f,107 6.380180669e-02f,108 -5.934976547e-02f109 };110 BOOST_MATH_STATIC const T Q[] =111 {112 1.000000000e+00f,113 2.679722431e+00f,114 1.561635813e+00f,115 1.573660661e-01f116 };117 if(x < tools::log_max_value<T>())118 return ((tools::evaluate_rational(P, Q, T(1 / x)) + 1) * exp(-x) / sqrt(x));119 else120 {121 T ex = exp(-x / 2);122 return ((tools::evaluate_rational(P, Q, T(1 / x)) + 1) * ex / sqrt(x)) * ex;123 }124 }125}126 127template <typename T>128BOOST_MATH_GPU_ENABLED T bessel_k0_imp(const T& x, const boost::math::integral_constant<int, 53>&)129{130 BOOST_MATH_STD_USING131 if(x <= 1)132 {133 // Maximum Deviation Found: 6.077e-17134 // Expected Error Term : -6.077e-17135 // Maximum Relative Change in Control Points : 7.797e-02136 // Max Error found at double precision = Poly : 1.003156e-16137 BOOST_MATH_STATIC const T Y = 1.137250900268554688;138 BOOST_MATH_STATIC const T P[] =139 {140 -1.372509002685546267e-01,141 2.574916117833312855e-01,142 1.395474602146869316e-02,143 5.445476986653926759e-04,144 7.125159422136622118e-06145 };146 BOOST_MATH_STATIC const T Q[] =147 {148 1.000000000000000000e+00,149 -5.458333438017788530e-02,150 1.291052816975251298e-03,151 -1.367653946978586591e-05152 };153 154 T a = x * x / 4;155 a = (tools::evaluate_polynomial(P, a) / tools::evaluate_polynomial(Q, a) + Y) * a + 1;156 157 // Maximum Deviation Found: 3.429e-18158 // Expected Error Term : 3.392e-18159 // Maximum Relative Change in Control Points : 2.041e-02160 // Max Error found at double precision = Poly : 2.513112e-16161 BOOST_MATH_STATIC const T P2[] =162 {163 1.159315156584124484e-01,164 2.789828789146031732e-01,165 2.524892993216121934e-02,166 8.460350907213637784e-04,167 1.491471924309617534e-05,168 1.627106892422088488e-07,169 1.208266102392756055e-09,170 6.611686391749704310e-12171 };172 173 return tools::evaluate_polynomial(P2, T(x * x)) - log(x) * a;174 }175 else176 {177 // Maximum Deviation Found: 4.316e-17178 // Expected Error Term : 9.570e-18179 // Maximum Relative Change in Control Points : 2.757e-01180 // Max Error found at double precision = Poly : 1.001560e-16181 182 BOOST_MATH_STATIC const T Y = 1;183 BOOST_MATH_STATIC const T P[] =184 {185 2.533141373155002416e-01,186 3.628342133984595192e+00,187 1.868441889406606057e+01,188 4.306243981063412784e+01,189 4.424116209627428189e+01,190 1.562095339356220468e+01,191 -1.810138978229410898e+00,192 -1.414237994269995877e+00,193 -9.369168119754924625e-02194 };195 BOOST_MATH_STATIC const T Q[] =196 {197 1.000000000000000000e+00,198 1.494194694879908328e+01,199 8.265296455388554217e+01,200 2.162779506621866970e+02,201 2.845145155184222157e+02,202 1.851714491916334995e+02,203 5.486540717439723515e+01,204 6.118075837628957015e+00,205 1.586261269326235053e-01206 };207 if(x < tools::log_max_value<T>())208 return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));209 else210 {211 T ex = exp(-x / 2);212 return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;213 }214 }215}216 217template <typename T>218BOOST_MATH_GPU_ENABLED T bessel_k0_imp(const T& x, const boost::math::integral_constant<int, 64>&)219{220 BOOST_MATH_STD_USING221 if(x <= 1)222 {223 // Maximum Deviation Found: 2.180e-22224 // Expected Error Term : 2.180e-22225 // Maximum Relative Change in Control Points : 2.943e-01226 // Max Error found at float80 precision = Poly : 3.923207e-20227 BOOST_MATH_STATIC const T Y = 1.137250900268554687500e+00;228 BOOST_MATH_STATIC const T P[] =229 {230 BOOST_MATH_BIG_CONSTANT(T, 64, -1.372509002685546875002e-01),231 BOOST_MATH_BIG_CONSTANT(T, 64, 2.566481981037407600436e-01),232 BOOST_MATH_BIG_CONSTANT(T, 64, 1.551881122448948854873e-02),233 BOOST_MATH_BIG_CONSTANT(T, 64, 6.646112454323276529650e-04),234 BOOST_MATH_BIG_CONSTANT(T, 64, 1.213747930378196492543e-05),235 BOOST_MATH_BIG_CONSTANT(T, 64, 9.423709328020389560844e-08)236 };237 BOOST_MATH_STATIC const T Q[] =238 {239 BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),240 BOOST_MATH_BIG_CONSTANT(T, 64, -4.843828412587773008342e-02),241 BOOST_MATH_BIG_CONSTANT(T, 64, 1.088484822515098936140e-03),242 BOOST_MATH_BIG_CONSTANT(T, 64, -1.374724008530702784829e-05),243 BOOST_MATH_BIG_CONSTANT(T, 64, 8.452665455952581680339e-08)244 };245 246 247 T a = x * x / 4;248 a = (tools::evaluate_polynomial(P, a) / tools::evaluate_polynomial(Q, a) + Y) * a + 1;249 250 // Maximum Deviation Found: 2.440e-21251 // Expected Error Term : -2.434e-21252 // Maximum Relative Change in Control Points : 2.459e-02253 // Max Error found at float80 precision = Poly : 1.482487e-19254 BOOST_MATH_STATIC const T P2[] =255 {256 BOOST_MATH_BIG_CONSTANT(T, 64, 1.159315156584124488110e-01),257 BOOST_MATH_BIG_CONSTANT(T, 64, 2.764832791416047889734e-01),258 BOOST_MATH_BIG_CONSTANT(T, 64, 1.926062887220923354112e-02),259 BOOST_MATH_BIG_CONSTANT(T, 64, 3.660777862036966089410e-04),260 BOOST_MATH_BIG_CONSTANT(T, 64, 2.094942446930673386849e-06)261 };262 BOOST_MATH_STATIC const T Q2[] =263 {264 BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),265 BOOST_MATH_BIG_CONSTANT(T, 64, -2.156100313881251616320e-02),266 BOOST_MATH_BIG_CONSTANT(T, 64, 2.315993873344905957033e-04),267 BOOST_MATH_BIG_CONSTANT(T, 64, -1.529444499350703363451e-06),268 BOOST_MATH_BIG_CONSTANT(T, 64, 5.524988589917857531177e-09)269 };270 return tools::evaluate_rational(P2, Q2, T(x * x)) - log(x) * a;271 }272 else273 {274 // Maximum Deviation Found: 4.291e-20275 // Expected Error Term : 2.236e-21276 // Maximum Relative Change in Control Points : 3.021e-01277 //Max Error found at float80 precision = Poly : 8.727378e-20278 BOOST_MATH_STATIC const T Y = 1;279 BOOST_MATH_STATIC const T P[] =280 {281 BOOST_MATH_BIG_CONSTANT(T, 64, 2.533141373155002512056e-01),282 BOOST_MATH_BIG_CONSTANT(T, 64, 5.417942070721928652715e+00),283 BOOST_MATH_BIG_CONSTANT(T, 64, 4.477464607463971754433e+01),284 BOOST_MATH_BIG_CONSTANT(T, 64, 1.838745728725943889876e+02),285 BOOST_MATH_BIG_CONSTANT(T, 64, 4.009736314927811202517e+02),286 BOOST_MATH_BIG_CONSTANT(T, 64, 4.557411293123609803452e+02),287 BOOST_MATH_BIG_CONSTANT(T, 64, 2.360222564015361268955e+02),288 BOOST_MATH_BIG_CONSTANT(T, 64, 2.385435333168505701022e+01),289 BOOST_MATH_BIG_CONSTANT(T, 64, -1.750195760942181592050e+01),290 BOOST_MATH_BIG_CONSTANT(T, 64, -4.059789241612946683713e+00),291 BOOST_MATH_BIG_CONSTANT(T, 64, -1.612783121537333908889e-01)292 };293 BOOST_MATH_STATIC const T Q[] =294 {295 BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),296 BOOST_MATH_BIG_CONSTANT(T, 64, 2.200669254769325861404e+01),297 BOOST_MATH_BIG_CONSTANT(T, 64, 1.900177593527144126549e+02),298 BOOST_MATH_BIG_CONSTANT(T, 64, 8.361003989965786932682e+02),299 BOOST_MATH_BIG_CONSTANT(T, 64, 2.041319870804843395893e+03),300 BOOST_MATH_BIG_CONSTANT(T, 64, 2.828491555113790345068e+03),301 BOOST_MATH_BIG_CONSTANT(T, 64, 2.190342229261529076624e+03),302 BOOST_MATH_BIG_CONSTANT(T, 64, 9.003330795963812219852e+02),303 BOOST_MATH_BIG_CONSTANT(T, 64, 1.773371397243777891569e+02),304 BOOST_MATH_BIG_CONSTANT(T, 64, 1.368634935531158398439e+01),305 BOOST_MATH_BIG_CONSTANT(T, 64, 2.543310879400359967327e-01)306 };307 if(x < tools::log_max_value<T>())308 return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));309 else310 {311 T ex = exp(-x / 2);312 return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;313 }314 }315}316 317template <typename T>318BOOST_MATH_GPU_ENABLED T bessel_k0_imp(const T& x, const boost::math::integral_constant<int, 113>&)319{320 BOOST_MATH_STD_USING321 if(x <= 1)322 {323 // Maximum Deviation Found: 5.682e-37324 // Expected Error Term : 5.682e-37325 // Maximum Relative Change in Control Points : 6.094e-04326 // Max Error found at float128 precision = Poly : 5.338213e-35327 BOOST_MATH_STATIC const T Y = 1.137250900268554687500000000000000000e+00f;328 BOOST_MATH_STATIC const T P[] =329 {330 BOOST_MATH_BIG_CONSTANT(T, 113, -1.372509002685546875000000000000000006e-01),331 BOOST_MATH_BIG_CONSTANT(T, 113, 2.556212905071072782462974351698081303e-01),332 BOOST_MATH_BIG_CONSTANT(T, 113, 1.742459135264203478530904179889103929e-02),333 BOOST_MATH_BIG_CONSTANT(T, 113, 8.077860530453688571555479526961318918e-04),334 BOOST_MATH_BIG_CONSTANT(T, 113, 1.868173911669241091399374307788635148e-05),335 BOOST_MATH_BIG_CONSTANT(T, 113, 2.496405768838992243478709145123306602e-07),336 BOOST_MATH_BIG_CONSTANT(T, 113, 1.752489221949580551692915881999762125e-09),337 BOOST_MATH_BIG_CONSTANT(T, 113, 5.243010555737173524710512824955368526e-12)338 };339 BOOST_MATH_STATIC const T Q[] =340 {341 BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),342 BOOST_MATH_BIG_CONSTANT(T, 113, -4.095631064064621099785696980653193721e-02),343 BOOST_MATH_BIG_CONSTANT(T, 113, 8.313880983725212151967078809725835532e-04),344 BOOST_MATH_BIG_CONSTANT(T, 113, -1.095229912293480063501285562382835142e-05),345 BOOST_MATH_BIG_CONSTANT(T, 113, 1.022828799511943141130509410251996277e-07),346 BOOST_MATH_BIG_CONSTANT(T, 113, -6.860874007419812445494782795829046836e-10),347 BOOST_MATH_BIG_CONSTANT(T, 113, 3.107297802344970725756092082686799037e-12),348 BOOST_MATH_BIG_CONSTANT(T, 113, -7.460529579244623559164763757787600944e-15)349 };350 T a = x * x / 4;351 a = (tools::evaluate_rational(P, Q, a) + Y) * a + 1;352 353 // Maximum Deviation Found: 5.173e-38354 // Expected Error Term : 5.105e-38355 // Maximum Relative Change in Control Points : 9.734e-03356 // Max Error found at float128 precision = Poly : 1.688806e-34357 BOOST_MATH_STATIC const T P2[] =358 {359 BOOST_MATH_BIG_CONSTANT(T, 113, 1.159315156584124488107200313757741370e-01),360 BOOST_MATH_BIG_CONSTANT(T, 113, 2.789828789146031122026800078439435369e-01),361 BOOST_MATH_BIG_CONSTANT(T, 113, 2.524892993216269451266750049024628432e-02),362 BOOST_MATH_BIG_CONSTANT(T, 113, 8.460350907082229957222453839935101823e-04),363 BOOST_MATH_BIG_CONSTANT(T, 113, 1.491471929926042875260452849503857976e-05),364 BOOST_MATH_BIG_CONSTANT(T, 113, 1.627105610481598430816014719558896866e-07),365 BOOST_MATH_BIG_CONSTANT(T, 113, 1.208426165007797264194914898538250281e-09),366 BOOST_MATH_BIG_CONSTANT(T, 113, 6.508697838747354949164182457073784117e-12),367 BOOST_MATH_BIG_CONSTANT(T, 113, 2.659784680639805301101014383907273109e-14),368 BOOST_MATH_BIG_CONSTANT(T, 113, 8.531090131964391104248859415958109654e-17),369 BOOST_MATH_BIG_CONSTANT(T, 113, 2.205195117066478034260323124669936314e-19),370 BOOST_MATH_BIG_CONSTANT(T, 113, 4.692219280289030165761119775783115426e-22),371 BOOST_MATH_BIG_CONSTANT(T, 113, 8.362350161092532344171965861545860747e-25),372 BOOST_MATH_BIG_CONSTANT(T, 113, 1.277990623924628999539014980773738258e-27)373 };374 375 return tools::evaluate_polynomial(P2, T(x * x)) - log(x) * a;376 }377 else378 {379 // Maximum Deviation Found: 1.462e-34380 // Expected Error Term : 4.917e-40381 // Maximum Relative Change in Control Points : 3.385e-01382 // Max Error found at float128 precision = Poly : 1.567573e-34383 BOOST_MATH_STATIC const T Y = 1;384 BOOST_MATH_STATIC const T P[] =385 {386 BOOST_MATH_BIG_CONSTANT(T, 113, 2.533141373155002512078826424055226265e-01),387 BOOST_MATH_BIG_CONSTANT(T, 113, 2.001949740768235770078339977110749204e+01),388 BOOST_MATH_BIG_CONSTANT(T, 113, 6.991516715983883248363351472378349986e+02),389 BOOST_MATH_BIG_CONSTANT(T, 113, 1.429587951594593159075690819360687720e+04),390 BOOST_MATH_BIG_CONSTANT(T, 113, 1.911933815201948768044660065771258450e+05),391 BOOST_MATH_BIG_CONSTANT(T, 113, 1.769943016204926614862175317962439875e+06),392 BOOST_MATH_BIG_CONSTANT(T, 113, 1.170866154649560750500954150401105606e+07),393 BOOST_MATH_BIG_CONSTANT(T, 113, 5.634687099724383996792011977705727661e+07),394 BOOST_MATH_BIG_CONSTANT(T, 113, 1.989524036456492581597607246664394014e+08),395 BOOST_MATH_BIG_CONSTANT(T, 113, 5.160394785715328062088529400178080360e+08),396 BOOST_MATH_BIG_CONSTANT(T, 113, 9.778173054417826368076483100902201433e+08),397 BOOST_MATH_BIG_CONSTANT(T, 113, 1.335667778588806892764139643950439733e+09),398 BOOST_MATH_BIG_CONSTANT(T, 113, 1.283635100080306980206494425043706838e+09),399 BOOST_MATH_BIG_CONSTANT(T, 113, 8.300616188213640626577036321085025855e+08),400 BOOST_MATH_BIG_CONSTANT(T, 113, 3.277591957076162984986406540894621482e+08),401 BOOST_MATH_BIG_CONSTANT(T, 113, 5.564360536834214058158565361486115932e+07),402 BOOST_MATH_BIG_CONSTANT(T, 113, -1.043505161612403359098596828115690596e+07),403 BOOST_MATH_BIG_CONSTANT(T, 113, -7.217035248223503605127967970903027314e+06),404 BOOST_MATH_BIG_CONSTANT(T, 113, -1.422938158797326748375799596769964430e+06),405 BOOST_MATH_BIG_CONSTANT(T, 113, -1.229125746200586805278634786674745210e+05),406 BOOST_MATH_BIG_CONSTANT(T, 113, -4.201632288615609937883545928660649813e+03),407 BOOST_MATH_BIG_CONSTANT(T, 113, -3.690820607338480548346746717311811406e+01)408 };409 BOOST_MATH_STATIC const T Q[] =410 {411 BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),412 BOOST_MATH_BIG_CONSTANT(T, 113, 7.964877874035741452203497983642653107e+01),413 BOOST_MATH_BIG_CONSTANT(T, 113, 2.808929943826193766839360018583294769e+03),414 BOOST_MATH_BIG_CONSTANT(T, 113, 5.814524004679994110944366890912384139e+04),415 BOOST_MATH_BIG_CONSTANT(T, 113, 7.897794522506725610540209610337355118e+05),416 BOOST_MATH_BIG_CONSTANT(T, 113, 7.456339470955813675629523617440433672e+06),417 BOOST_MATH_BIG_CONSTANT(T, 113, 5.057818717813969772198911392875127212e+07),418 BOOST_MATH_BIG_CONSTANT(T, 113, 2.513821619536852436424913886081133209e+08),419 BOOST_MATH_BIG_CONSTANT(T, 113, 9.255938846873380596038513316919990776e+08),420 BOOST_MATH_BIG_CONSTANT(T, 113, 2.537077551699028079347581816919572141e+09),421 BOOST_MATH_BIG_CONSTANT(T, 113, 5.176769339768120752974843214652367321e+09),422 BOOST_MATH_BIG_CONSTANT(T, 113, 7.828722317390455845253191337207432060e+09),423 BOOST_MATH_BIG_CONSTANT(T, 113, 8.698864296569996402006511705803675890e+09),424 BOOST_MATH_BIG_CONSTANT(T, 113, 7.007803261356636409943826918468544629e+09),425 BOOST_MATH_BIG_CONSTANT(T, 113, 4.016564631288740308993071395104715469e+09),426 BOOST_MATH_BIG_CONSTANT(T, 113, 1.595893010619754750655947035567624730e+09),427 BOOST_MATH_BIG_CONSTANT(T, 113, 4.241241839120481076862742189989406856e+08),428 BOOST_MATH_BIG_CONSTANT(T, 113, 7.168778094393076220871007550235840858e+07),429 BOOST_MATH_BIG_CONSTANT(T, 113, 7.156200301360388147635052029404211109e+06),430 BOOST_MATH_BIG_CONSTANT(T, 113, 3.752130382550379886741949463587008794e+05),431 BOOST_MATH_BIG_CONSTANT(T, 113, 8.370574966987293592457152146806662562e+03),432 BOOST_MATH_BIG_CONSTANT(T, 113, 4.871254714311063594080644835895740323e+01)433 };434 if(-x > tools::log_min_value<T>())435 return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));436 else437 {438 T ex = exp(-x / 2);439 return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;440 }441 }442}443 444template <typename T>445BOOST_MATH_GPU_ENABLED T bessel_k0_imp(const T& x, const boost::math::integral_constant<int, 0>&)446{447 if(boost::math::tools::digits<T>() <= 24)448 return bessel_k0_imp(x, boost::math::integral_constant<int, 24>());449 else if(boost::math::tools::digits<T>() <= 53)450 return bessel_k0_imp(x, boost::math::integral_constant<int, 53>());451 else if(boost::math::tools::digits<T>() <= 64)452 return bessel_k0_imp(x, boost::math::integral_constant<int, 64>());453 else if(boost::math::tools::digits<T>() <= 113)454 return bessel_k0_imp(x, boost::math::integral_constant<int, 113>());455 BOOST_MATH_ASSERT(0);456 return 0;457}458 459template <typename T>460BOOST_MATH_GPU_ENABLED inline T bessel_k0(const T& x)461{462 typedef boost::math::integral_constant<int,463 ((boost::math::numeric_limits<T>::digits == 0) || (boost::math::numeric_limits<T>::radix != 2)) ?464 0 :465 boost::math::numeric_limits<T>::digits <= 24 ?466 24 :467 boost::math::numeric_limits<T>::digits <= 53 ?468 53 :469 boost::math::numeric_limits<T>::digits <= 64 ?470 64 :471 boost::math::numeric_limits<T>::digits <= 113 ?472 113 : -1473 > tag_type;474 475 return bessel_k0_imp(x, tag_type());476}477 478}}} // namespaces479 480#ifdef _MSC_VER481#pragma warning(pop)482#endif483 484#endif // BOOST_MATH_BESSEL_K0_HPP485 486