975 lines · plain
1// Copyright John Maddock 2008.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5//6// Wrapper that works with mpfr_class defined in gmpfrxx.h7// See http://math.berkeley.edu/~wilken/code/gmpfrxx/8// Also requires the gmp and mpfr libraries.9//10 11#ifndef BOOST_MATH_MPLFR_BINDINGS_HPP12#define BOOST_MATH_MPLFR_BINDINGS_HPP13 14#include <type_traits>15 16#ifdef _MSC_VER17//18// We get a lot of warnings from the gmp, mpfr and gmpfrxx headers,19// disable them here, so we only see warnings from *our* code:20//21#pragma warning(push)22#pragma warning(disable: 4127 4800 4512)23#endif24 25#include <gmpfrxx.h>26 27#ifdef _MSC_VER28#pragma warning(pop)29#endif30 31#include <boost/math/tools/precision.hpp>32#include <boost/math/tools/real_cast.hpp>33#include <boost/math/policies/policy.hpp>34#include <boost/math/distributions/fwd.hpp>35#include <boost/math/special_functions/math_fwd.hpp>36#include <boost/math/bindings/detail/big_digamma.hpp>37#include <boost/math/bindings/detail/big_lanczos.hpp>38#include <boost/math/tools/big_constant.hpp>39#include <boost/math/tools/config.hpp>40 41inline mpfr_class fabs(const mpfr_class& v)42{43 return abs(v);44}45template <class T, class U>46inline mpfr_class fabs(const __gmp_expr<T,U>& v)47{48 return abs(static_cast<mpfr_class>(v));49}50 51inline mpfr_class pow(const mpfr_class& b, const mpfr_class& e)52{53 mpfr_class result;54 mpfr_pow(result.__get_mp(), b.__get_mp(), e.__get_mp(), GMP_RNDN);55 return result;56}57/*58template <class T, class U, class V, class W>59inline mpfr_class pow(const __gmp_expr<T,U>& b, const __gmp_expr<V,W>& e)60{61 return pow(static_cast<mpfr_class>(b), static_cast<mpfr_class>(e));62}63*/64inline mpfr_class ldexp(const mpfr_class& v, int e)65{66 //int e = mpfr_get_exp(*v.__get_mp());67 mpfr_class result(v);68 mpfr_set_exp(result.__get_mp(), e);69 return result;70}71template <class T, class U>72inline mpfr_class ldexp(const __gmp_expr<T,U>& v, int e)73{74 return ldexp(static_cast<mpfr_class>(v), e);75}76 77inline mpfr_class frexp(const mpfr_class& v, int* expon)78{79 int e = mpfr_get_exp(v.__get_mp());80 mpfr_class result(v);81 mpfr_set_exp(result.__get_mp(), 0);82 *expon = e;83 return result;84}85template <class T, class U>86inline mpfr_class frexp(const __gmp_expr<T,U>& v, int* expon)87{88 return frexp(static_cast<mpfr_class>(v), expon);89}90 91inline mpfr_class fmod(const mpfr_class& v1, const mpfr_class& v2)92{93 mpfr_class n;94 if(v1 < 0)95 n = ceil(v1 / v2);96 else97 n = floor(v1 / v2);98 return v1 - n * v2;99}100template <class T, class U, class V, class W>101inline mpfr_class fmod(const __gmp_expr<T,U>& v1, const __gmp_expr<V,W>& v2)102{103 return fmod(static_cast<mpfr_class>(v1), static_cast<mpfr_class>(v2));104}105 106template <class Policy>107inline mpfr_class modf(const mpfr_class& v, long long* ipart, const Policy& pol)108{109 *ipart = lltrunc(v, pol);110 return v - boost::math::tools::real_cast<mpfr_class>(*ipart);111}112template <class T, class U, class Policy>113inline mpfr_class modf(const __gmp_expr<T,U>& v, long long* ipart, const Policy& pol)114{115 return modf(static_cast<mpfr_class>(v), ipart, pol);116}117 118template <class Policy>119inline int iround(mpfr_class const& x, const Policy&)120{121 return boost::math::tools::real_cast<int>(boost::math::round(x, typename boost::math::policies::normalise<Policy, boost::math::policies::rounding_error< boost::math::policies::throw_on_error> >::type()));122}123template <class T, class U, class Policy>124inline int iround(__gmp_expr<T,U> const& x, const Policy& pol)125{126 return iround(static_cast<mpfr_class>(x), pol);127}128 129template <class Policy>130inline long lround(mpfr_class const& x, const Policy&)131{132 return boost::math::tools::real_cast<long>(boost::math::round(x, typename boost::math::policies::normalise<Policy, boost::math::policies::rounding_error< boost::math::policies::throw_on_error> >::type()));133}134template <class T, class U, class Policy>135inline long lround(__gmp_expr<T,U> const& x, const Policy& pol)136{137 return lround(static_cast<mpfr_class>(x), pol);138}139 140template <class Policy>141inline long long llround(mpfr_class const& x, const Policy&)142{143 return boost::math::tools::real_cast<long long>(boost::math::round(x, typename boost::math::policies::normalise<Policy, boost::math::policies::rounding_error< boost::math::policies::throw_on_error> >::type()));144}145template <class T, class U, class Policy>146inline long long llround(__gmp_expr<T,U> const& x, const Policy& pol)147{148 return llround(static_cast<mpfr_class>(x), pol);149}150 151template <class Policy>152inline int itrunc(mpfr_class const& x, const Policy&)153{154 return boost::math::tools::real_cast<int>(boost::math::trunc(x, typename boost::math::policies::normalise<Policy, boost::math::policies::rounding_error< boost::math::policies::throw_on_error> >::type()));155}156template <class T, class U, class Policy>157inline int itrunc(__gmp_expr<T,U> const& x, const Policy& pol)158{159 return itrunc(static_cast<mpfr_class>(x), pol);160}161 162template <class Policy>163inline long ltrunc(mpfr_class const& x, const Policy&)164{165 return boost::math::tools::real_cast<long>(boost::math::trunc(x, typename boost::math::policies::normalise<Policy, boost::math::policies::rounding_error< boost::math::policies::throw_on_error> >::type()));166}167template <class T, class U, class Policy>168inline long ltrunc(__gmp_expr<T,U> const& x, const Policy& pol)169{170 return ltrunc(static_cast<mpfr_class>(x), pol);171}172 173template <class Policy>174inline long long lltrunc(mpfr_class const& x, const Policy&)175{176 return boost::math::tools::real_cast<long long>(boost::math::trunc(x, typename boost::math::policies::normalise<Policy, boost::math::policies::rounding_error< boost::math::policies::throw_on_error> >::type()));177}178template <class T, class U, class Policy>179inline long long lltrunc(__gmp_expr<T,U> const& x, const Policy& pol)180{181 return lltrunc(static_cast<mpfr_class>(x), pol);182}183 184namespace boost{185 186#ifdef BOOST_MATH_USE_FLOAT128187 template<> struct std::is_convertible<BOOST_MATH_FLOAT128_TYPE, mpfr_class> : public std::integral_constant<bool, false>{};188#endif189 template<> struct std::is_convertible<long long, mpfr_class> : public std::integral_constant<bool, false>{};190 191namespace math{192 193#if defined(__GNUC__) && (__GNUC__ < 4)194 using ::iround;195 using ::lround;196 using ::llround;197 using ::itrunc;198 using ::ltrunc;199 using ::lltrunc;200 using ::modf;201#endif202 203namespace lanczos{204 205struct mpfr_lanczos206{207 static mpfr_class lanczos_sum(const mpfr_class& z)208 {209 unsigned long p = z.get_dprec();210 if(p <= 72)211 return lanczos13UDT::lanczos_sum(z);212 else if(p <= 120)213 return lanczos22UDT::lanczos_sum(z);214 else if(p <= 170)215 return lanczos31UDT::lanczos_sum(z);216 else //if(p <= 370) approx 100 digit precision:217 return lanczos61UDT::lanczos_sum(z);218 }219 static mpfr_class lanczos_sum_expG_scaled(const mpfr_class& z)220 {221 unsigned long p = z.get_dprec();222 if(p <= 72)223 return lanczos13UDT::lanczos_sum_expG_scaled(z);224 else if(p <= 120)225 return lanczos22UDT::lanczos_sum_expG_scaled(z);226 else if(p <= 170)227 return lanczos31UDT::lanczos_sum_expG_scaled(z);228 else //if(p <= 370) approx 100 digit precision:229 return lanczos61UDT::lanczos_sum_expG_scaled(z);230 }231 static mpfr_class lanczos_sum_near_1(const mpfr_class& z)232 {233 unsigned long p = z.get_dprec();234 if(p <= 72)235 return lanczos13UDT::lanczos_sum_near_1(z);236 else if(p <= 120)237 return lanczos22UDT::lanczos_sum_near_1(z);238 else if(p <= 170)239 return lanczos31UDT::lanczos_sum_near_1(z);240 else //if(p <= 370) approx 100 digit precision:241 return lanczos61UDT::lanczos_sum_near_1(z);242 }243 static mpfr_class lanczos_sum_near_2(const mpfr_class& z)244 {245 unsigned long p = z.get_dprec();246 if(p <= 72)247 return lanczos13UDT::lanczos_sum_near_2(z);248 else if(p <= 120)249 return lanczos22UDT::lanczos_sum_near_2(z);250 else if(p <= 170)251 return lanczos31UDT::lanczos_sum_near_2(z);252 else //if(p <= 370) approx 100 digit precision:253 return lanczos61UDT::lanczos_sum_near_2(z);254 }255 static mpfr_class g()256 {257 unsigned long p = mpfr_class::get_dprec();258 if(p <= 72)259 return lanczos13UDT::g();260 else if(p <= 120)261 return lanczos22UDT::g();262 else if(p <= 170)263 return lanczos31UDT::g();264 else //if(p <= 370) approx 100 digit precision:265 return lanczos61UDT::g();266 }267};268 269template<class Policy>270struct lanczos<mpfr_class, Policy>271{272 typedef mpfr_lanczos type;273};274 275} // namespace lanczos276 277namespace constants{278 279template <class Real, class Policy>280struct construction_traits;281 282template <class Policy>283struct construction_traits<mpfr_class, Policy>284{285 typedef std::integral_constant<int, 0> type;286};287 288}289 290namespace tools291{292 293template <class T, class U>294struct promote_arg<__gmp_expr<T,U> >295{ // If T is integral type, then promote to double.296 typedef mpfr_class type;297};298 299template<>300inline int digits<mpfr_class>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr_class)) noexcept301{302 return mpfr_class::get_dprec();303}304 305namespace detail{306 307template<class Integer>308void convert_to_long_result(mpfr_class const& r, Integer& result)309{310 result = 0;311 I last_result(0);312 mpfr_class t(r);313 double term;314 do315 {316 term = real_cast<double>(t);317 last_result = result;318 result += static_cast<I>(term);319 t -= term;320 }while(result != last_result);321}322 323}324 325template <>326inline mpfr_class real_cast<mpfr_class, long long>(long long t)327{328 mpfr_class result;329 int expon = 0;330 int sign = 1;331 if(t < 0)332 {333 sign = -1;334 t = -t;335 }336 while(t)337 {338 result += ldexp(static_cast<double>(t & 0xffffL), expon);339 expon += 32;340 t >>= 32;341 }342 return result * sign;343}344template <>345inline unsigned real_cast<unsigned, mpfr_class>(mpfr_class t)346{347 return t.get_ui();348}349template <>350inline int real_cast<int, mpfr_class>(mpfr_class t)351{352 return t.get_si();353}354template <>355inline double real_cast<double, mpfr_class>(mpfr_class t)356{357 return t.get_d();358}359template <>360inline float real_cast<float, mpfr_class>(mpfr_class t)361{362 return static_cast<float>(t.get_d());363}364template <>365inline long real_cast<long, mpfr_class>(mpfr_class t)366{367 long result;368 detail::convert_to_long_result(t, result);369 return result;370}371template <>372inline long long real_cast<long long, mpfr_class>(mpfr_class t)373{374 long long result;375 detail::convert_to_long_result(t, result);376 return result;377}378 379template <>380inline mpfr_class max_value<mpfr_class>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr_class))381{382 static bool has_init = false;383 static mpfr_class val;384 if(!has_init)385 {386 val = 0.5;387 mpfr_set_exp(val.__get_mp(), mpfr_get_emax());388 has_init = true;389 }390 return val;391}392 393template <>394inline mpfr_class min_value<mpfr_class>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr_class))395{396 static bool has_init = false;397 static mpfr_class val;398 if(!has_init)399 {400 val = 0.5;401 mpfr_set_exp(val.__get_mp(), mpfr_get_emin());402 has_init = true;403 }404 return val;405}406 407template <>408inline mpfr_class log_max_value<mpfr_class>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr_class))409{410 static bool has_init = false;411 static mpfr_class val = max_value<mpfr_class>();412 if(!has_init)413 {414 val = log(val);415 has_init = true;416 }417 return val;418}419 420template <>421inline mpfr_class log_min_value<mpfr_class>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr_class))422{423 static bool has_init = false;424 static mpfr_class val = max_value<mpfr_class>();425 if(!has_init)426 {427 val = log(val);428 has_init = true;429 }430 return val;431}432 433template <>434inline mpfr_class epsilon<mpfr_class>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr_class))435{436 return ldexp(mpfr_class(1), 1-boost::math::policies::digits<mpfr_class, boost::math::policies::policy<> >());437}438 439} // namespace tools440 441namespace policies{442 443template <class T, class U, class Policy>444struct evaluation<__gmp_expr<T, U>, Policy>445{446 typedef mpfr_class type;447};448 449}450 451template <class Policy>452inline mpfr_class skewness(const extreme_value_distribution<mpfr_class, Policy>& /*dist*/)453{454 //455 // This is 12 * sqrt(6) * zeta(3) / pi^3:456 // See http://mathworld.wolfram.com/ExtremeValueDistribution.html457 //458 #ifdef BOOST_MATH_STANDALONE459 static_assert(sizeof(Policy) == 0, "mpfr skewness can not be calculated in standalone mode");460 #endif461 462 return static_cast<mpfr_class>("1.1395470994046486574927930193898461120875997958366");463}464 465template <class Policy>466inline mpfr_class skewness(const rayleigh_distribution<mpfr_class, Policy>& /*dist*/)467{468 // using namespace boost::math::constants;469 #ifdef BOOST_MATH_STANDALONE470 static_assert(sizeof(Policy) == 0, "mpfr skewness can not be calculated in standalone mode");471 #endif472 473 return static_cast<mpfr_class>("0.63111065781893713819189935154422777984404221106391");474 // Computed using NTL at 150 bit, about 50 decimal digits.475 // return 2 * root_pi<RealType>() * pi_minus_three<RealType>() / pow23_four_minus_pi<RealType>();476}477 478template <class Policy>479inline mpfr_class kurtosis(const rayleigh_distribution<mpfr_class, Policy>& /*dist*/)480{481 // using namespace boost::math::constants;482 #ifdef BOOST_MATH_STANDALONE483 static_assert(sizeof(Policy) == 0, "mpfr kurtosis can not be calculated in standalone mode");484 #endif485 486 return static_cast<mpfr_class>("3.2450893006876380628486604106197544154170667057995");487 // Computed using NTL at 150 bit, about 50 decimal digits.488 // return 3 - (6 * pi<RealType>() * pi<RealType>() - 24 * pi<RealType>() + 16) /489 // (four_minus_pi<RealType>() * four_minus_pi<RealType>());490}491 492template <class Policy>493inline mpfr_class kurtosis_excess(const rayleigh_distribution<mpfr_class, Policy>& /*dist*/)494{495 //using namespace boost::math::constants;496 // Computed using NTL at 150 bit, about 50 decimal digits.497 #ifdef BOOST_MATH_STANDALONE498 static_assert(sizeof(Policy) == 0, "mpfr excess kurtosis can not be calculated in standalone mode");499 #endif500 501 return static_cast<mpfr_class>("0.2450893006876380628486604106197544154170667057995");502 // return -(6 * pi<RealType>() * pi<RealType>() - 24 * pi<RealType>() + 16) /503 // (four_minus_pi<RealType>() * four_minus_pi<RealType>());504} // kurtosis505 506namespace detail{507 508//509// Version of Digamma accurate to ~100 decimal digits.510//511template <class Policy>512mpfr_class digamma_imp(mpfr_class x, const std::integral_constant<int, 0>* , const Policy& pol)513{514 //515 // This handles reflection of negative arguments, and all our516 // empfr_classor handling, then forwards to the T-specific approximation.517 //518 BOOST_MATH_STD_USING // ADL of std functions.519 520 mpfr_class result = 0;521 //522 // Check for negative arguments and use reflection:523 //524 if(x < 0)525 {526 // Reflect:527 x = 1 - x;528 // Argument reduction for tan:529 mpfr_class remainder = x - floor(x);530 // Shift to negative if > 0.5:531 if(remainder > 0.5)532 {533 remainder -= 1;534 }535 //536 // check for evaluation at a negative pole:537 //538 if(remainder == 0)539 {540 return policies::raise_pole_error<mpfr_class>("boost::math::digamma<%1%>(%1%)", nullptr, (1-x), pol);541 }542 result = constants::pi<mpfr_class>() / tan(constants::pi<mpfr_class>() * remainder);543 }544 result += big_digamma(x);545 return result;546}547//548// Specialisations of this function provides the initial549// starting guess for Halley iteration:550//551template <class Policy>552inline mpfr_class erf_inv_imp(const mpfr_class& p, const mpfr_class& q, const Policy&, const std::integral_constant<int, 64>*)553{554 BOOST_MATH_STD_USING // for ADL of std names.555 556 mpfr_class result = 0;557 558 if(p <= 0.5)559 {560 //561 // Evaluate inverse erf using the rational approximation:562 //563 // x = p(p+10)(Y+R(p))564 //565 // Where Y is a constant, and R(p) is optimised for a low566 // absolute empfr_classor compared to |Y|.567 //568 // double: Max empfr_classor found: 2.001849e-18569 // long double: Max empfr_classor found: 1.017064e-20570 // Maximum Deviation Found (actual empfr_classor term at infinite precision) 8.030e-21571 //572 static const float Y = 0.0891314744949340820313f;573 static const mpfr_class P[] = {574 -0.000508781949658280665617,575 -0.00836874819741736770379,576 0.0334806625409744615033,577 -0.0126926147662974029034,578 -0.0365637971411762664006,579 0.0219878681111168899165,580 0.00822687874676915743155,581 -0.00538772965071242932965582 };583 static const mpfr_class Q[] = {584 1,585 -0.970005043303290640362,586 -1.56574558234175846809,587 1.56221558398423026363,588 0.662328840472002992063,589 -0.71228902341542847553,590 -0.0527396382340099713954,591 0.0795283687341571680018,592 -0.00233393759374190016776,593 0.000886216390456424707504594 };595 mpfr_class g = p * (p + 10);596 mpfr_class r = tools::evaluate_polynomial(P, p) / tools::evaluate_polynomial(Q, p);597 result = g * Y + g * r;598 }599 else if(q >= 0.25)600 {601 //602 // Rational approximation for 0.5 > q >= 0.25603 //604 // x = sqrt(-2*log(q)) / (Y + R(q))605 //606 // Where Y is a constant, and R(q) is optimised for a low607 // absolute empfr_classor compared to Y.608 //609 // double : Max empfr_classor found: 7.403372e-17610 // long double : Max empfr_classor found: 6.084616e-20611 // Maximum Deviation Found (empfr_classor term) 4.811e-20612 //613 static const float Y = 2.249481201171875f;614 static const mpfr_class P[] = {615 -0.202433508355938759655,616 0.105264680699391713268,617 8.37050328343119927838,618 17.6447298408374015486,619 -18.8510648058714251895,620 -44.6382324441786960818,621 17.445385985570866523,622 21.1294655448340526258,623 -3.67192254707729348546624 };625 static const mpfr_class Q[] = {626 1,627 6.24264124854247537712,628 3.9713437953343869095,629 -28.6608180499800029974,630 -20.1432634680485188801,631 48.5609213108739935468,632 10.8268667355460159008,633 -22.6436933413139721736,634 1.72114765761200282724635 };636 mpfr_class g = sqrt(-2 * log(q));637 mpfr_class xs = q - 0.25;638 mpfr_class r = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);639 result = g / (Y + r);640 }641 else642 {643 //644 // For q < 0.25 we have a series of rational approximations all645 // of the general form:646 //647 // let: x = sqrt(-log(q))648 //649 // Then the result is given by:650 //651 // x(Y+R(x-B))652 //653 // where Y is a constant, B is the lowest value of x for which654 // the approximation is valid, and R(x-B) is optimised for a low655 // absolute empfr_classor compared to Y.656 //657 // Note that almost all code will really go through the first658 // or maybe second approximation. After than we're dealing with very659 // small input values indeed: 80 and 128 bit long double's go all the660 // way down to ~ 1e-5000 so the "tail" is rather long...661 //662 mpfr_class x = sqrt(-log(q));663 if(x < 3)664 {665 // Max empfr_classor found: 1.089051e-20666 static const float Y = 0.807220458984375f;667 static const mpfr_class P[] = {668 -0.131102781679951906451,669 -0.163794047193317060787,670 0.117030156341995252019,671 0.387079738972604337464,672 0.337785538912035898924,673 0.142869534408157156766,674 0.0290157910005329060432,675 0.00214558995388805277169,676 -0.679465575181126350155e-6,677 0.285225331782217055858e-7,678 -0.681149956853776992068e-9679 };680 static const mpfr_class Q[] = {681 1,682 3.46625407242567245975,683 5.38168345707006855425,684 4.77846592945843778382,685 2.59301921623620271374,686 0.848854343457902036425,687 0.152264338295331783612,688 0.01105924229346489121689 };690 mpfr_class xs = x - 1.125;691 mpfr_class R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);692 result = Y * x + R * x;693 }694 else if(x < 6)695 {696 // Max empfr_classor found: 8.389174e-21697 static const float Y = 0.93995571136474609375f;698 static const mpfr_class P[] = {699 -0.0350353787183177984712,700 -0.00222426529213447927281,701 0.0185573306514231072324,702 0.00950804701325919603619,703 0.00187123492819559223345,704 0.000157544617424960554631,705 0.460469890584317994083e-5,706 -0.230404776911882601748e-9,707 0.266339227425782031962e-11708 };709 static const mpfr_class Q[] = {710 1,711 1.3653349817554063097,712 0.762059164553623404043,713 0.220091105764131249824,714 0.0341589143670947727934,715 0.00263861676657015992959,716 0.764675292302794483503e-4717 };718 mpfr_class xs = x - 3;719 mpfr_class R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);720 result = Y * x + R * x;721 }722 else if(x < 18)723 {724 // Max empfr_classor found: 1.481312e-19725 static const float Y = 0.98362827301025390625f;726 static const mpfr_class P[] = {727 -0.0167431005076633737133,728 -0.00112951438745580278863,729 0.00105628862152492910091,730 0.000209386317487588078668,731 0.149624783758342370182e-4,732 0.449696789927706453732e-6,733 0.462596163522878599135e-8,734 -0.281128735628831791805e-13,735 0.99055709973310326855e-16736 };737 static const mpfr_class Q[] = {738 1,739 0.591429344886417493481,740 0.138151865749083321638,741 0.0160746087093676504695,742 0.000964011807005165528527,743 0.275335474764726041141e-4,744 0.282243172016108031869e-6745 };746 mpfr_class xs = x - 6;747 mpfr_class R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);748 result = Y * x + R * x;749 }750 else if(x < 44)751 {752 // Max empfr_classor found: 5.697761e-20753 static const float Y = 0.99714565277099609375f;754 static const mpfr_class P[] = {755 -0.0024978212791898131227,756 -0.779190719229053954292e-5,757 0.254723037413027451751e-4,758 0.162397777342510920873e-5,759 0.396341011304801168516e-7,760 0.411632831190944208473e-9,761 0.145596286718675035587e-11,762 -0.116765012397184275695e-17763 };764 static const mpfr_class Q[] = {765 1,766 0.207123112214422517181,767 0.0169410838120975906478,768 0.000690538265622684595676,769 0.145007359818232637924e-4,770 0.144437756628144157666e-6,771 0.509761276599778486139e-9772 };773 mpfr_class xs = x - 18;774 mpfr_class R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);775 result = Y * x + R * x;776 }777 else778 {779 // Max empfr_classor found: 1.279746e-20780 static const float Y = 0.99941349029541015625f;781 static const mpfr_class P[] = {782 -0.000539042911019078575891,783 -0.28398759004727721098e-6,784 0.899465114892291446442e-6,785 0.229345859265920864296e-7,786 0.225561444863500149219e-9,787 0.947846627503022684216e-12,788 0.135880130108924861008e-14,789 -0.348890393399948882918e-21790 };791 static const mpfr_class Q[] = {792 1,793 0.0845746234001899436914,794 0.00282092984726264681981,795 0.468292921940894236786e-4,796 0.399968812193862100054e-6,797 0.161809290887904476097e-8,798 0.231558608310259605225e-11799 };800 mpfr_class xs = x - 44;801 mpfr_class R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);802 result = Y * x + R * x;803 }804 }805 return result;806}807 808inline mpfr_class bessel_i0(mpfr_class x)809{810 #ifdef BOOST_MATH_STANDALONE811 static_assert(sizeof(x) == 0, "mpfr bessel_i0 can not be calculated in standalone mode");812 #endif813 814 static const mpfr_class P1[] = {815 static_cast<mpfr_class>("-2.2335582639474375249e+15"),816 static_cast<mpfr_class>("-5.5050369673018427753e+14"),817 static_cast<mpfr_class>("-3.2940087627407749166e+13"),818 static_cast<mpfr_class>("-8.4925101247114157499e+11"),819 static_cast<mpfr_class>("-1.1912746104985237192e+10"),820 static_cast<mpfr_class>("-1.0313066708737980747e+08"),821 static_cast<mpfr_class>("-5.9545626019847898221e+05"),822 static_cast<mpfr_class>("-2.4125195876041896775e+03"),823 static_cast<mpfr_class>("-7.0935347449210549190e+00"),824 static_cast<mpfr_class>("-1.5453977791786851041e-02"),825 static_cast<mpfr_class>("-2.5172644670688975051e-05"),826 static_cast<mpfr_class>("-3.0517226450451067446e-08"),827 static_cast<mpfr_class>("-2.6843448573468483278e-11"),828 static_cast<mpfr_class>("-1.5982226675653184646e-14"),829 static_cast<mpfr_class>("-5.2487866627945699800e-18"),830 };831 static const mpfr_class Q1[] = {832 static_cast<mpfr_class>("-2.2335582639474375245e+15"),833 static_cast<mpfr_class>("7.8858692566751002988e+12"),834 static_cast<mpfr_class>("-1.2207067397808979846e+10"),835 static_cast<mpfr_class>("1.0377081058062166144e+07"),836 static_cast<mpfr_class>("-4.8527560179962773045e+03"),837 static_cast<mpfr_class>("1.0"),838 };839 static const mpfr_class P2[] = {840 static_cast<mpfr_class>("-2.2210262233306573296e-04"),841 static_cast<mpfr_class>("1.3067392038106924055e-02"),842 static_cast<mpfr_class>("-4.4700805721174453923e-01"),843 static_cast<mpfr_class>("5.5674518371240761397e+00"),844 static_cast<mpfr_class>("-2.3517945679239481621e+01"),845 static_cast<mpfr_class>("3.1611322818701131207e+01"),846 static_cast<mpfr_class>("-9.6090021968656180000e+00"),847 };848 static const mpfr_class Q2[] = {849 static_cast<mpfr_class>("-5.5194330231005480228e-04"),850 static_cast<mpfr_class>("3.2547697594819615062e-02"),851 static_cast<mpfr_class>("-1.1151759188741312645e+00"),852 static_cast<mpfr_class>("1.3982595353892851542e+01"),853 static_cast<mpfr_class>("-6.0228002066743340583e+01"),854 static_cast<mpfr_class>("8.5539563258012929600e+01"),855 static_cast<mpfr_class>("-3.1446690275135491500e+01"),856 static_cast<mpfr_class>("1.0"),857 };858 mpfr_class value, factor, r;859 860 BOOST_MATH_STD_USING861 using namespace boost::math::tools;862 863 if (x < 0)864 {865 x = -x; // even function866 }867 if (x == 0)868 {869 return static_cast<mpfr_class>(1);870 }871 if (x <= 15) // x in (0, 15]872 {873 mpfr_class y = x * x;874 value = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);875 }876 else // x in (15, \infty)877 {878 mpfr_class y = 1 / x - mpfr_class(1) / 15;879 r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);880 factor = exp(x) / sqrt(x);881 value = factor * r;882 }883 884 return value;885}886 887inline mpfr_class bessel_i1(mpfr_class x)888{889 static const mpfr_class P1[] = {890 static_cast<mpfr_class>("-1.4577180278143463643e+15"),891 static_cast<mpfr_class>("-1.7732037840791591320e+14"),892 static_cast<mpfr_class>("-6.9876779648010090070e+12"),893 static_cast<mpfr_class>("-1.3357437682275493024e+11"),894 static_cast<mpfr_class>("-1.4828267606612366099e+09"),895 static_cast<mpfr_class>("-1.0588550724769347106e+07"),896 static_cast<mpfr_class>("-5.1894091982308017540e+04"),897 static_cast<mpfr_class>("-1.8225946631657315931e+02"),898 static_cast<mpfr_class>("-4.7207090827310162436e-01"),899 static_cast<mpfr_class>("-9.1746443287817501309e-04"),900 static_cast<mpfr_class>("-1.3466829827635152875e-06"),901 static_cast<mpfr_class>("-1.4831904935994647675e-09"),902 static_cast<mpfr_class>("-1.1928788903603238754e-12"),903 static_cast<mpfr_class>("-6.5245515583151902910e-16"),904 static_cast<mpfr_class>("-1.9705291802535139930e-19"),905 };906 static const mpfr_class Q1[] = {907 static_cast<mpfr_class>("-2.9154360556286927285e+15"),908 static_cast<mpfr_class>("9.7887501377547640438e+12"),909 static_cast<mpfr_class>("-1.4386907088588283434e+10"),910 static_cast<mpfr_class>("1.1594225856856884006e+07"),911 static_cast<mpfr_class>("-5.1326864679904189920e+03"),912 static_cast<mpfr_class>("1.0"),913 };914 static const mpfr_class P2[] = {915 static_cast<mpfr_class>("1.4582087408985668208e-05"),916 static_cast<mpfr_class>("-8.9359825138577646443e-04"),917 static_cast<mpfr_class>("2.9204895411257790122e-02"),918 static_cast<mpfr_class>("-3.4198728018058047439e-01"),919 static_cast<mpfr_class>("1.3960118277609544334e+00"),920 static_cast<mpfr_class>("-1.9746376087200685843e+00"),921 static_cast<mpfr_class>("8.5591872901933459000e-01"),922 static_cast<mpfr_class>("-6.0437159056137599999e-02"),923 };924 static const mpfr_class Q2[] = {925 static_cast<mpfr_class>("3.7510433111922824643e-05"),926 static_cast<mpfr_class>("-2.2835624489492512649e-03"),927 static_cast<mpfr_class>("7.4212010813186530069e-02"),928 static_cast<mpfr_class>("-8.5017476463217924408e-01"),929 static_cast<mpfr_class>("3.2593714889036996297e+00"),930 static_cast<mpfr_class>("-3.8806586721556593450e+00"),931 static_cast<mpfr_class>("1.0"),932 };933 mpfr_class value, factor, r, w;934 935 BOOST_MATH_STD_USING936 using namespace boost::math::tools;937 938 w = abs(x);939 if (x == 0)940 {941 return static_cast<mpfr_class>(0);942 }943 if (w <= 15) // w in (0, 15]944 {945 mpfr_class y = x * x;946 r = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);947 factor = w;948 value = factor * r;949 }950 else // w in (15, \infty)951 {952 mpfr_class y = 1 / w - mpfr_class(1) / 15;953 r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);954 factor = exp(w) / sqrt(w);955 value = factor * r;956 }957 958 if (x < 0)959 {960 value *= -value; // odd function961 }962 return value;963}964 965} // namespace detail966 967}968 969template<> struct std::is_convertible<long double, mpfr_class> : public std::false_type{};970 971}972 973#endif // BOOST_MATH_MPLFR_BINDINGS_HPP974 975