527 lines · plain
1// Copyright (c) 2006 Xiaogang Zhang2// Copyright (c) 2017 John Maddock 3// 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_K1_HPP8#define BOOST_MATH_BESSEL_K1_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/policies/error_handling.hpp>23#include <boost/math/tools/assert.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/2016/01/05/rational-approximations-for-the-modified-bessel-function-of-the-second-kind-k1-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 50 template <typename T, int N>51 BOOST_MATH_GPU_ENABLED inline T bessel_k1_imp(const T&, const boost::math::integral_constant<int, N>&)52 {53 BOOST_MATH_ASSERT(0);54 return 0;55 }56 57 template <typename T>58 BOOST_MATH_GPU_ENABLED T bessel_k1_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: 3.090e-1264 // Expected Error Term : -3.053e-1265 // Maximum Relative Change in Control Points : 4.927e-0266 // Max Error found at float precision = Poly : 7.918347e-1067 BOOST_MATH_STATIC const T Y = 8.695471287e-02f;68 BOOST_MATH_STATIC const T P[] =69 {70 -3.621379531e-03f,71 7.131781976e-03f,72 -1.535278300e-05f73 };74 BOOST_MATH_STATIC const T Q[] =75 {76 1.000000000e+00f,77 -5.173102701e-02f,78 9.203530671e-04f79 };80 81 T a = x * x / 4;82 a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2;83 84 // Maximum Deviation Found: 3.556e-0885 // Expected Error Term : -3.541e-0886 // Maximum Relative Change in Control Points : 8.203e-0287 BOOST_MATH_STATIC const T P2[] =88 {89 -3.079657469e-01f,90 -8.537108913e-02f,91 -4.640275408e-03f,92 -1.156442414e-04f93 };94 95 return tools::evaluate_polynomial(P2, T(x * x)) * x + 1 / x + log(x) * a;96 }97 else98 {99 // Maximum Deviation Found: 3.369e-08100 // Expected Error Term : -3.227e-08101 // Maximum Relative Change in Control Points : 9.917e-02102 // Max Error found at float precision = Poly : 6.084411e-08103 BOOST_MATH_STATIC const T Y = 1.450342178f;104 BOOST_MATH_STATIC const T P[] =105 {106 -1.970280088e-01f,107 2.188747807e-02f,108 7.270394756e-01f,109 2.490678196e-01f110 };111 BOOST_MATH_STATIC const T Q[] =112 {113 1.000000000e+00f,114 2.274292882e+00f,115 9.904984851e-01f,116 4.585534549e-02f117 };118 if(x < tools::log_max_value<T>())119 return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));120 else121 {122 T ex = exp(-x / 2);123 return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;124 }125 }126 }127 128 template <typename T>129 BOOST_MATH_GPU_ENABLED T bessel_k1_imp(const T& x, const boost::math::integral_constant<int, 53>&)130 {131 BOOST_MATH_STD_USING132 if(x <= 1)133 {134 // Maximum Deviation Found: 1.922e-17135 // Expected Error Term : 1.921e-17136 // Maximum Relative Change in Control Points : 5.287e-03137 // Max Error found at double precision = Poly : 2.004747e-17138 BOOST_MATH_STATIC const T Y = 8.69547128677368164e-02f;139 BOOST_MATH_STATIC const T P[] =140 {141 -3.62137953440350228e-03,142 7.11842087490330300e-03,143 1.00302560256614306e-05,144 1.77231085381040811e-06145 };146 BOOST_MATH_STATIC const T Q[] =147 {148 1.00000000000000000e+00,149 -4.80414794429043831e-02,150 9.85972641934416525e-04,151 -8.91196859397070326e-06152 };153 154 T a = x * x / 4;155 a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2;156 157 // Maximum Deviation Found: 4.053e-17158 // Expected Error Term : -4.053e-17159 // Maximum Relative Change in Control Points : 3.103e-04160 // Max Error found at double precision = Poly : 1.246698e-16161 162 BOOST_MATH_STATIC const T P2[] =163 {164 -3.07965757829206184e-01,165 -7.80929703673074907e-02,166 -2.70619343754051620e-03,167 -2.49549522229072008e-05168 };169 BOOST_MATH_STATIC const T Q2[] = 170 {171 1.00000000000000000e+00,172 -2.36316836412163098e-02,173 2.64524577525962719e-04,174 -1.49749618004162787e-06175 };176 177 return tools::evaluate_rational(P2, Q2, T(x * x)) * x + 1 / x + log(x) * a;178 }179 else180 {181 // Maximum Deviation Found: 8.883e-17182 // Expected Error Term : -1.641e-17183 // Maximum Relative Change in Control Points : 2.786e-01184 // Max Error found at double precision = Poly : 1.258798e-16185 186 BOOST_MATH_STATIC const T Y = 1.45034217834472656f;187 BOOST_MATH_STATIC const T P[] =188 {189 -1.97028041029226295e-01,190 -2.32408961548087617e+00,191 -7.98269784507699938e+00,192 -2.39968410774221632e+00,193 3.28314043780858713e+01,194 5.67713761158496058e+01,195 3.30907788466509823e+01,196 6.62582288933739787e+00,197 3.08851840645286691e-01198 };199 BOOST_MATH_STATIC const T Q[] =200 {201 1.00000000000000000e+00,202 1.41811409298826118e+01,203 7.35979466317556420e+01,204 1.77821793937080859e+02,205 2.11014501598705982e+02,206 1.19425262951064454e+02,207 2.88448064302447607e+01,208 2.27912927104139732e+00,209 2.50358186953478678e-02210 };211 if(x < tools::log_max_value<T>())212 return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));213 else214 {215 T ex = exp(-x / 2);216 return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;217 }218 }219 }220 221 template <typename T>222 BOOST_MATH_GPU_ENABLED T bessel_k1_imp(const T& x, const boost::math::integral_constant<int, 64>&)223 {224 BOOST_MATH_STD_USING225 if(x <= 1)226 {227 // Maximum Deviation Found: 5.549e-23228 // Expected Error Term : -5.548e-23229 // Maximum Relative Change in Control Points : 2.002e-03230 // Max Error found at float80 precision = Poly : 9.352785e-22231 BOOST_MATH_STATIC const T Y = 8.695471286773681640625e-02f;232 BOOST_MATH_STATIC const T P[] =233 {234 BOOST_MATH_BIG_CONSTANT(T, 64, -3.621379534403483072861e-03),235 BOOST_MATH_BIG_CONSTANT(T, 64, 7.102135866103952705932e-03),236 BOOST_MATH_BIG_CONSTANT(T, 64, 4.167545240236717601167e-05),237 BOOST_MATH_BIG_CONSTANT(T, 64, 2.537484002571894870830e-06),238 BOOST_MATH_BIG_CONSTANT(T, 64, 6.603228256820000135990e-09)239 };240 BOOST_MATH_STATIC const T Q[] =241 {242 BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),243 BOOST_MATH_BIG_CONSTANT(T, 64, -4.354457194045068370363e-02),244 BOOST_MATH_BIG_CONSTANT(T, 64, 8.709137201220209072820e-04),245 BOOST_MATH_BIG_CONSTANT(T, 64, -9.676151796359590545143e-06),246 BOOST_MATH_BIG_CONSTANT(T, 64, 5.162715192766245311659e-08)247 };248 249 T a = x * x / 4;250 a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2;251 252 // Maximum Deviation Found: 1.995e-23253 // Expected Error Term : 1.995e-23254 // Maximum Relative Change in Control Points : 8.174e-04255 // Max Error found at float80 precision = Poly : 4.137325e-20256 BOOST_MATH_STATIC const T P2[] =257 {258 BOOST_MATH_BIG_CONSTANT(T, 64, -3.079657578292062244054e-01),259 BOOST_MATH_BIG_CONSTANT(T, 64, -7.963049154965966503231e-02),260 BOOST_MATH_BIG_CONSTANT(T, 64, -3.103277523735639924895e-03),261 BOOST_MATH_BIG_CONSTANT(T, 64, -4.023052834702215699504e-05),262 BOOST_MATH_BIG_CONSTANT(T, 64, -1.719459155018493821839e-07)263 };264 BOOST_MATH_STATIC const T Q2[] = 265 {266 BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),267 BOOST_MATH_BIG_CONSTANT(T, 64, -1.863917670410152669768e-02),268 BOOST_MATH_BIG_CONSTANT(T, 64, 1.699367098849735298090e-04),269 BOOST_MATH_BIG_CONSTANT(T, 64, -9.309358790546076298429e-07),270 BOOST_MATH_BIG_CONSTANT(T, 64, 2.708893480271612711933e-09)271 };272 273 return tools::evaluate_rational(P2, Q2, T(x * x)) * x + 1 / x + log(x) * a;274 }275 else276 {277 // Maximum Deviation Found: 9.785e-20278 // Expected Error Term : -3.302e-21279 // Maximum Relative Change in Control Points : 3.432e-01280 // Max Error found at float80 precision = Poly : 1.083755e-19281 BOOST_MATH_STATIC const T Y = 1.450342178344726562500e+00f;282 BOOST_MATH_STATIC const T P[] =283 {284 BOOST_MATH_BIG_CONSTANT(T, 64, -1.970280410292263112917e-01),285 BOOST_MATH_BIG_CONSTANT(T, 64, -4.058564803062959169322e+00),286 BOOST_MATH_BIG_CONSTANT(T, 64, -3.036658174194917777473e+01),287 BOOST_MATH_BIG_CONSTANT(T, 64, -9.576825392332820142173e+01),288 BOOST_MATH_BIG_CONSTANT(T, 64, -6.706969489248020941949e+01),289 BOOST_MATH_BIG_CONSTANT(T, 64, 3.264572499406168221382e+02),290 BOOST_MATH_BIG_CONSTANT(T, 64, 8.584972047303151034100e+02),291 BOOST_MATH_BIG_CONSTANT(T, 64, 8.422082733280017909550e+02),292 BOOST_MATH_BIG_CONSTANT(T, 64, 3.738005441471368178383e+02),293 BOOST_MATH_BIG_CONSTANT(T, 64, 7.016938390144121276609e+01),294 BOOST_MATH_BIG_CONSTANT(T, 64, 4.319614662598089438939e+00),295 BOOST_MATH_BIG_CONSTANT(T, 64, 3.710715864316521856193e-02)296 };297 BOOST_MATH_STATIC const T Q[] =298 {299 BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),300 BOOST_MATH_BIG_CONSTANT(T, 64, 2.298433045824439052398e+01),301 BOOST_MATH_BIG_CONSTANT(T, 64, 2.082047745067709230037e+02),302 BOOST_MATH_BIG_CONSTANT(T, 64, 9.662367854250262046592e+02),303 BOOST_MATH_BIG_CONSTANT(T, 64, 2.504148628460454004686e+03),304 BOOST_MATH_BIG_CONSTANT(T, 64, 3.712730364911389908905e+03),305 BOOST_MATH_BIG_CONSTANT(T, 64, 3.108002081150068641112e+03),306 BOOST_MATH_BIG_CONSTANT(T, 64, 1.400149940532448553143e+03),307 BOOST_MATH_BIG_CONSTANT(T, 64, 3.083303048095846226299e+02),308 BOOST_MATH_BIG_CONSTANT(T, 64, 2.748706060530351833346e+01),309 BOOST_MATH_BIG_CONSTANT(T, 64, 6.321900849331506946977e-01),310 };311 if(x < tools::log_max_value<T>())312 return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));313 else314 {315 T ex = exp(-x / 2);316 return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;317 }318 }319 }320 321 template <typename T>322 BOOST_MATH_GPU_ENABLED T bessel_k1_imp(const T& x, const boost::math::integral_constant<int, 113>&)323 {324 BOOST_MATH_STD_USING325 if(x <= 1)326 {327 // Maximum Deviation Found: 7.120e-35328 // Expected Error Term : -7.119e-35329 // Maximum Relative Change in Control Points : 1.207e-03330 // Max Error found at float128 precision = Poly : 7.143688e-35331 BOOST_MATH_STATIC const T Y = 8.695471286773681640625000000000000000e-02f;332 BOOST_MATH_STATIC const T P[] =333 {334 BOOST_MATH_BIG_CONSTANT(T, 113, -3.621379534403483072916666666666595475e-03),335 BOOST_MATH_BIG_CONSTANT(T, 113, 7.074117676930975433219826471336547627e-03),336 BOOST_MATH_BIG_CONSTANT(T, 113, 9.631337631362776369069668419033041661e-05),337 BOOST_MATH_BIG_CONSTANT(T, 113, 3.468935967870048731821071646104412775e-06),338 BOOST_MATH_BIG_CONSTANT(T, 113, 2.956705020559599861444492614737168261e-08),339 BOOST_MATH_BIG_CONSTANT(T, 113, 2.347140307321161346703214099534250263e-10),340 BOOST_MATH_BIG_CONSTANT(T, 113, 5.569608494081482873946791086435679661e-13)341 };342 BOOST_MATH_STATIC const T Q[] =343 {344 BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),345 BOOST_MATH_BIG_CONSTANT(T, 113, -3.580768910152105375615558920428350204e-02),346 BOOST_MATH_BIG_CONSTANT(T, 113, 6.197467671701485365363068445534557369e-04),347 BOOST_MATH_BIG_CONSTANT(T, 113, -6.707466533308630411966030561446666237e-06),348 BOOST_MATH_BIG_CONSTANT(T, 113, 4.846687802282250112624373388491123527e-08),349 BOOST_MATH_BIG_CONSTANT(T, 113, -2.248493131151981569517383040323900343e-10),350 BOOST_MATH_BIG_CONSTANT(T, 113, 5.319279786372775264555728921709381080e-13)351 };352 353 T a = x * x / 4;354 a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2;355 356 // Maximum Deviation Found: 4.473e-37357 // Expected Error Term : 4.473e-37358 // Maximum Relative Change in Control Points : 8.550e-04359 // Max Error found at float128 precision = Poly : 8.167701e-35360 BOOST_MATH_STATIC const T P2[] =361 {362 BOOST_MATH_BIG_CONSTANT(T, 113, -3.079657578292062244053600156878870690e-01),363 BOOST_MATH_BIG_CONSTANT(T, 113, -8.133183745732467770755578848987414875e-02),364 BOOST_MATH_BIG_CONSTANT(T, 113, -3.548968792764174773125420229299431951e-03),365 BOOST_MATH_BIG_CONSTANT(T, 113, -5.886125468718182876076972186152445490e-05),366 BOOST_MATH_BIG_CONSTANT(T, 113, -4.506712111733707245745396404449639865e-07),367 BOOST_MATH_BIG_CONSTANT(T, 113, -1.632502325880313239698965376754406011e-09),368 BOOST_MATH_BIG_CONSTANT(T, 113, -2.311973065898784812266544485665624227e-12)369 };370 BOOST_MATH_STATIC const T Q2[] = 371 {372 BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),373 BOOST_MATH_BIG_CONSTANT(T, 113, -1.311471216733781016657962995723287450e-02),374 BOOST_MATH_BIG_CONSTANT(T, 113, 8.571876054797365417068164018709472969e-05),375 BOOST_MATH_BIG_CONSTANT(T, 113, -3.630181215268238731442496851497901293e-07),376 BOOST_MATH_BIG_CONSTANT(T, 113, 1.070176111227805048604885986867484807e-09),377 BOOST_MATH_BIG_CONSTANT(T, 113, -2.129046580769872602793220056461084761e-12),378 BOOST_MATH_BIG_CONSTANT(T, 113, 2.294906469421390890762001971790074432e-15)379 };380 381 return tools::evaluate_rational(P2, Q2, T(x * x)) * x + 1 / x + log(x) * a;382 }383 else if(x < 4)384 {385 // Max error in interpolated form: 5.307e-37386 // Max Error found at float128 precision = Poly: 7.087862e-35387 BOOST_MATH_STATIC const T Y = 1.5023040771484375f;388 BOOST_MATH_STATIC const T P[] =389 {390 BOOST_MATH_BIG_CONSTANT(T, 113, -2.489899398329369710528254347931380044e-01),391 BOOST_MATH_BIG_CONSTANT(T, 113, -6.819080211203854781858815596508456873e+00),392 BOOST_MATH_BIG_CONSTANT(T, 113, -7.599915699069767382647695624952723034e+01),393 BOOST_MATH_BIG_CONSTANT(T, 113, -4.450211910821295507926582231071300718e+02),394 BOOST_MATH_BIG_CONSTANT(T, 113, -1.451374687870925175794150513723956533e+03),395 BOOST_MATH_BIG_CONSTANT(T, 113, -2.405805746895098802803503988539098226e+03),396 BOOST_MATH_BIG_CONSTANT(T, 113, -5.638808326778389656403861103277220518e+02),397 BOOST_MATH_BIG_CONSTANT(T, 113, 5.513958744081268456191778822780865708e+03),398 BOOST_MATH_BIG_CONSTANT(T, 113, 1.121301640926540743072258116122834804e+04),399 BOOST_MATH_BIG_CONSTANT(T, 113, 1.080094900175649541266613109971296190e+04),400 BOOST_MATH_BIG_CONSTANT(T, 113, 5.896531083639613332407534434915552429e+03),401 BOOST_MATH_BIG_CONSTANT(T, 113, 1.856602122319645694042555107114028437e+03),402 BOOST_MATH_BIG_CONSTANT(T, 113, 3.237121918853145421414003823957537419e+02),403 BOOST_MATH_BIG_CONSTANT(T, 113, 2.842072954561323076230238664623893504e+01),404 BOOST_MATH_BIG_CONSTANT(T, 113, 1.039705646510167437971862966128055524e+00),405 BOOST_MATH_BIG_CONSTANT(T, 113, 1.008418100718254816100425022904039530e-02)406 };407 BOOST_MATH_STATIC const T Q[] =408 {409 BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),410 BOOST_MATH_BIG_CONSTANT(T, 113, 2.927456835239137986889227412815459529e+01),411 BOOST_MATH_BIG_CONSTANT(T, 113, 3.598985593265577043711382994516531273e+02),412 BOOST_MATH_BIG_CONSTANT(T, 113, 2.449897377085510281395819892689690579e+03),413 BOOST_MATH_BIG_CONSTANT(T, 113, 1.025555887684561913263090023158085327e+04),414 BOOST_MATH_BIG_CONSTANT(T, 113, 2.774140447181062463181892531100679195e+04),415 BOOST_MATH_BIG_CONSTANT(T, 113, 4.962055507843204417243602332246120418e+04),416 BOOST_MATH_BIG_CONSTANT(T, 113, 5.908269326976180183216954452196772931e+04),417 BOOST_MATH_BIG_CONSTANT(T, 113, 4.655160454422016855911700790722577942e+04),418 BOOST_MATH_BIG_CONSTANT(T, 113, 2.383586885019548163464418964577684608e+04),419 BOOST_MATH_BIG_CONSTANT(T, 113, 7.679920375586960324298491662159976419e+03),420 BOOST_MATH_BIG_CONSTANT(T, 113, 1.478586421028842906987799049804565008e+03),421 BOOST_MATH_BIG_CONSTANT(T, 113, 1.565384974896746094224942654383537090e+02),422 BOOST_MATH_BIG_CONSTANT(T, 113, 7.902617937084010911005732488607114511e+00),423 BOOST_MATH_BIG_CONSTANT(T, 113, 1.429293010387921526110949911029094926e-01),424 BOOST_MATH_BIG_CONSTANT(T, 113, 3.880342607911083143560111853491047663e-04)425 };426 return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));427 }428 else429 {430 // Maximum Deviation Found: 4.359e-37431 // Expected Error Term : -6.565e-40432 // Maximum Relative Change in Control Points : 1.880e-01433 // Max Error found at float128 precision = Poly : 2.943572e-35434 BOOST_MATH_STATIC const T Y = 1.308816909790039062500000000000000000f;435 BOOST_MATH_STATIC const T P[] =436 {437 BOOST_MATH_BIG_CONSTANT(T, 113, -5.550277247453881129211735759447737350e-02),438 BOOST_MATH_BIG_CONSTANT(T, 113, -3.485883080219574328217554864956175929e+00),439 BOOST_MATH_BIG_CONSTANT(T, 113, -8.903760658131484239300875153154881958e+01),440 BOOST_MATH_BIG_CONSTANT(T, 113, -1.144813672213626237418235110712293337e+03),441 BOOST_MATH_BIG_CONSTANT(T, 113, -6.498400501156131446691826557494158173e+03),442 BOOST_MATH_BIG_CONSTANT(T, 113, 1.573531831870363502604119835922166116e+04),443 BOOST_MATH_BIG_CONSTANT(T, 113, 5.417416550054632009958262596048841154e+05),444 BOOST_MATH_BIG_CONSTANT(T, 113, 4.271266450613557412825896604269130661e+06),445 BOOST_MATH_BIG_CONSTANT(T, 113, 1.898386013314389952534433455681107783e+07),446 BOOST_MATH_BIG_CONSTANT(T, 113, 5.353798784656436259250791761023512750e+07),447 BOOST_MATH_BIG_CONSTANT(T, 113, 9.839619195427352438957774052763490067e+07),448 BOOST_MATH_BIG_CONSTANT(T, 113, 1.169246368651532232388152442538005637e+08),449 BOOST_MATH_BIG_CONSTANT(T, 113, 8.696368884166831199967845883371116431e+07),450 BOOST_MATH_BIG_CONSTANT(T, 113, 3.810226630422736458064005843327500169e+07),451 BOOST_MATH_BIG_CONSTANT(T, 113, 8.854996610560406127438950635716757614e+06),452 BOOST_MATH_BIG_CONSTANT(T, 113, 8.981057433937398731355768088809437625e+05),453 BOOST_MATH_BIG_CONSTANT(T, 113, 2.519440069856232098711793483639792952e+04)454 };455 BOOST_MATH_STATIC const T Q[] =456 {457 BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),458 BOOST_MATH_BIG_CONSTANT(T, 113, 7.127348248283623146544565916604103560e+01),459 BOOST_MATH_BIG_CONSTANT(T, 113, 2.205092684176906740104488180754982065e+03),460 BOOST_MATH_BIG_CONSTANT(T, 113, 3.911249195069050636298346469740075758e+04),461 BOOST_MATH_BIG_CONSTANT(T, 113, 4.426103406579046249654548481377792614e+05),462 BOOST_MATH_BIG_CONSTANT(T, 113, 3.365861555422488771286500241966208541e+06),463 BOOST_MATH_BIG_CONSTANT(T, 113, 1.765377714160383676864913709252529840e+07),464 BOOST_MATH_BIG_CONSTANT(T, 113, 6.453822726931857253365138260720815246e+07),465 BOOST_MATH_BIG_CONSTANT(T, 113, 1.643207885048369990391975749439783892e+08),466 BOOST_MATH_BIG_CONSTANT(T, 113, 2.882540678243694621895816336640877878e+08),467 BOOST_MATH_BIG_CONSTANT(T, 113, 3.410120808992380266174106812005338148e+08),468 BOOST_MATH_BIG_CONSTANT(T, 113, 2.628138016559335882019310900426773027e+08),469 BOOST_MATH_BIG_CONSTANT(T, 113, 1.250794693811010646965360198541047961e+08),470 BOOST_MATH_BIG_CONSTANT(T, 113, 3.378723408195485594610593014072950078e+07),471 BOOST_MATH_BIG_CONSTANT(T, 113, 4.488253856312453816451380319061865560e+06),472 BOOST_MATH_BIG_CONSTANT(T, 113, 2.202167197882689873967723350537104582e+05),473 BOOST_MATH_BIG_CONSTANT(T, 113, 1.673233230356966539460728211412989843e+03)474 };475 if(x < tools::log_max_value<T>())476 return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));477 else478 {479 T ex = exp(-x / 2);480 return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;481 }482 }483 }484 485 template <typename T>486 BOOST_MATH_GPU_ENABLED T bessel_k1_imp(const T& x, const boost::math::integral_constant<int, 0>&)487 {488 if(boost::math::tools::digits<T>() <= 24)489 return bessel_k1_imp(x, boost::math::integral_constant<int, 24>());490 else if(boost::math::tools::digits<T>() <= 53)491 return bessel_k1_imp(x, boost::math::integral_constant<int, 53>());492 else if(boost::math::tools::digits<T>() <= 64)493 return bessel_k1_imp(x, boost::math::integral_constant<int, 64>());494 else if(boost::math::tools::digits<T>() <= 113)495 return bessel_k1_imp(x, boost::math::integral_constant<int, 113>());496 BOOST_MATH_ASSERT(0);497 return 0;498 }499 500 template <typename T>501 BOOST_MATH_GPU_ENABLED inline T bessel_k1(const T& x)502 {503 typedef boost::math::integral_constant<int,504 ((boost::math::numeric_limits<T>::digits == 0) || (boost::math::numeric_limits<T>::radix != 2)) ?505 0 :506 boost::math::numeric_limits<T>::digits <= 24 ?507 24 :508 boost::math::numeric_limits<T>::digits <= 53 ?509 53 :510 boost::math::numeric_limits<T>::digits <= 64 ?511 64 :512 boost::math::numeric_limits<T>::digits <= 113 ?513 113 : -1514 > tag_type;515 516 return bessel_k1_imp(x, tag_type());517 }518 519}}} // namespaces520 521#ifdef _MSC_VER522#pragma warning(pop)523#endif524 525#endif // BOOST_MATH_BESSEL_K1_HPP526 527