922 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::mpreal 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_MPREAL_BINDINGS_HPP12#define BOOST_MATH_MPREAL_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 <mpreal.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/config.hpp>39#ifndef BOOST_MATH_STANDALONE40#include <boost/lexical_cast.hpp>41#endif42 43namespace mpfr{44 45template <class T>46inline mpreal operator + (const mpreal& r, const T& t){ return r + mpreal(t); }47template <class T>48inline mpreal operator - (const mpreal& r, const T& t){ return r - mpreal(t); }49template <class T>50inline mpreal operator * (const mpreal& r, const T& t){ return r * mpreal(t); }51template <class T>52inline mpreal operator / (const mpreal& r, const T& t){ return r / mpreal(t); }53 54template <class T>55inline mpreal operator + (const T& t, const mpreal& r){ return mpreal(t) + r; }56template <class T>57inline mpreal operator - (const T& t, const mpreal& r){ return mpreal(t) - r; }58template <class T>59inline mpreal operator * (const T& t, const mpreal& r){ return mpreal(t) * r; }60template <class T>61inline mpreal operator / (const T& t, const mpreal& r){ return mpreal(t) / r; }62 63template <class T>64inline bool operator == (const mpreal& r, const T& t){ return r == mpreal(t); }65template <class T>66inline bool operator != (const mpreal& r, const T& t){ return r != mpreal(t); }67template <class T>68inline bool operator <= (const mpreal& r, const T& t){ return r <= mpreal(t); }69template <class T>70inline bool operator >= (const mpreal& r, const T& t){ return r >= mpreal(t); }71template <class T>72inline bool operator < (const mpreal& r, const T& t){ return r < mpreal(t); }73template <class T>74inline bool operator > (const mpreal& r, const T& t){ return r > mpreal(t); }75 76template <class T>77inline bool operator == (const T& t, const mpreal& r){ return mpreal(t) == r; }78template <class T>79inline bool operator != (const T& t, const mpreal& r){ return mpreal(t) != r; }80template <class T>81inline bool operator <= (const T& t, const mpreal& r){ return mpreal(t) <= r; }82template <class T>83inline bool operator >= (const T& t, const mpreal& r){ return mpreal(t) >= r; }84template <class T>85inline bool operator < (const T& t, const mpreal& r){ return mpreal(t) < r; }86template <class T>87inline bool operator > (const T& t, const mpreal& r){ return mpreal(t) > r; }88 89/*90inline mpfr::mpreal fabs(const mpfr::mpreal& v)91{92 return abs(v);93}94inline mpfr::mpreal pow(const mpfr::mpreal& b, const mpfr::mpreal e)95{96 mpfr::mpreal result;97 mpfr_pow(result.__get_mp(), b.__get_mp(), e.__get_mp(), GMP_RNDN);98 return result;99}100*/101inline mpfr::mpreal ldexp(const mpfr::mpreal& v, int e)102{103 return mpfr::ldexp(v, static_cast<mp_exp_t>(e));104}105 106inline mpfr::mpreal frexp(const mpfr::mpreal& v, int* expon)107{108 mp_exp_t e;109 mpfr::mpreal r = mpfr::frexp(v, &e);110 *expon = e;111 return r;112}113 114#if (MPFR_VERSION < MPFR_VERSION_NUM(2,4,0))115mpfr::mpreal fmod(const mpfr::mpreal& v1, const mpfr::mpreal& v2)116{117 mpfr::mpreal n;118 if(v1 < 0)119 n = ceil(v1 / v2);120 else121 n = floor(v1 / v2);122 return v1 - n * v2;123}124#endif125 126template <class Policy>127inline mpfr::mpreal modf(const mpfr::mpreal& v, long long* ipart, const Policy& pol)128{129 *ipart = lltrunc(v, pol);130 return v - boost::math::tools::real_cast<mpfr::mpreal>(*ipart);131}132template <class Policy>133inline int iround(mpfr::mpreal const& x, const Policy& pol)134{135 return boost::math::tools::real_cast<int>(boost::math::round(x, pol));136}137 138template <class Policy>139inline long lround(mpfr::mpreal const& x, const Policy& pol)140{141 return boost::math::tools::real_cast<long>(boost::math::round(x, pol));142}143 144template <class Policy>145inline long long llround(mpfr::mpreal const& x, const Policy& pol)146{147 return boost::math::tools::real_cast<long long>(boost::math::round(x, pol));148}149 150template <class Policy>151inline int itrunc(mpfr::mpreal const& x, const Policy& pol)152{153 return boost::math::tools::real_cast<int>(boost::math::trunc(x, pol));154}155 156template <class Policy>157inline long ltrunc(mpfr::mpreal const& x, const Policy& pol)158{159 return boost::math::tools::real_cast<long>(boost::math::trunc(x, pol));160}161 162template <class Policy>163inline long long lltrunc(mpfr::mpreal const& x, const Policy& pol)164{165 return boost::math::tools::real_cast<long long>(boost::math::trunc(x, pol));166}167 168}169 170namespace boost{ namespace math{171 172#if defined(__GNUC__) && (__GNUC__ < 4)173 using ::iround;174 using ::lround;175 using ::llround;176 using ::itrunc;177 using ::ltrunc;178 using ::lltrunc;179 using ::modf;180#endif181 182namespace lanczos{183 184struct mpreal_lanczos185{186 static mpfr::mpreal lanczos_sum(const mpfr::mpreal& z)187 {188 unsigned long p = z.get_default_prec();189 if(p <= 72)190 return lanczos13UDT::lanczos_sum(z);191 else if(p <= 120)192 return lanczos22UDT::lanczos_sum(z);193 else if(p <= 170)194 return lanczos31UDT::lanczos_sum(z);195 else //if(p <= 370) approx 100 digit precision:196 return lanczos61UDT::lanczos_sum(z);197 }198 static mpfr::mpreal lanczos_sum_expG_scaled(const mpfr::mpreal& z)199 {200 unsigned long p = z.get_default_prec();201 if(p <= 72)202 return lanczos13UDT::lanczos_sum_expG_scaled(z);203 else if(p <= 120)204 return lanczos22UDT::lanczos_sum_expG_scaled(z);205 else if(p <= 170)206 return lanczos31UDT::lanczos_sum_expG_scaled(z);207 else //if(p <= 370) approx 100 digit precision:208 return lanczos61UDT::lanczos_sum_expG_scaled(z);209 }210 static mpfr::mpreal lanczos_sum_near_1(const mpfr::mpreal& z)211 {212 unsigned long p = z.get_default_prec();213 if(p <= 72)214 return lanczos13UDT::lanczos_sum_near_1(z);215 else if(p <= 120)216 return lanczos22UDT::lanczos_sum_near_1(z);217 else if(p <= 170)218 return lanczos31UDT::lanczos_sum_near_1(z);219 else //if(p <= 370) approx 100 digit precision:220 return lanczos61UDT::lanczos_sum_near_1(z);221 }222 static mpfr::mpreal lanczos_sum_near_2(const mpfr::mpreal& z)223 {224 unsigned long p = z.get_default_prec();225 if(p <= 72)226 return lanczos13UDT::lanczos_sum_near_2(z);227 else if(p <= 120)228 return lanczos22UDT::lanczos_sum_near_2(z);229 else if(p <= 170)230 return lanczos31UDT::lanczos_sum_near_2(z);231 else //if(p <= 370) approx 100 digit precision:232 return lanczos61UDT::lanczos_sum_near_2(z);233 }234 static mpfr::mpreal g()235 {236 unsigned long p = mpfr::mpreal::get_default_prec();237 if(p <= 72)238 return lanczos13UDT::g();239 else if(p <= 120)240 return lanczos22UDT::g();241 else if(p <= 170)242 return lanczos31UDT::g();243 else //if(p <= 370) approx 100 digit precision:244 return lanczos61UDT::g();245 }246};247 248template<class Policy>249struct lanczos<mpfr::mpreal, Policy>250{251 typedef mpreal_lanczos type;252};253 254} // namespace lanczos255 256namespace tools257{258 259template<>260inline int digits<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))261{262 return mpfr::mpreal::get_default_prec();263}264 265namespace detail{266 267template<class Integer268void convert_to_long_result(mpfr::mpreal const& r, Integer& result)269{270 result = 0;271 I last_result(0);272 mpfr::mpreal t(r);273 double term;274 do275 {276 term = real_cast<double>(t);277 last_result = result;278 result += static_cast<I>(term);279 t -= term;280 }while(result != last_result);281}282 283}284 285template <>286inline mpfr::mpreal real_cast<mpfr::mpreal, long long>(long long t)287{288 mpfr::mpreal result;289 int expon = 0;290 int sign = 1;291 if(t < 0)292 {293 sign = -1;294 t = -t;295 }296 while(t)297 {298 result += ldexp(static_cast<double>(t & 0xffffL), expon);299 expon += 32;300 t >>= 32;301 }302 return result * sign;303}304/*305template <>306inline unsigned real_cast<unsigned, mpfr::mpreal>(mpfr::mpreal t)307{308 return t.get_ui();309}310template <>311inline int real_cast<int, mpfr::mpreal>(mpfr::mpreal t)312{313 return t.get_si();314}315template <>316inline double real_cast<double, mpfr::mpreal>(mpfr::mpreal t)317{318 return t.get_d();319}320template <>321inline float real_cast<float, mpfr::mpreal>(mpfr::mpreal t)322{323 return static_cast<float>(t.get_d());324}325template <>326inline long real_cast<long, mpfr::mpreal>(mpfr::mpreal t)327{328 long result;329 detail::convert_to_long_result(t, result);330 return result;331}332*/333template <>334inline long long real_cast<long long, mpfr::mpreal>(mpfr::mpreal t)335{336 long long result;337 detail::convert_to_long_result(t, result);338 return result;339}340 341template <>342inline mpfr::mpreal max_value<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))343{344 static bool has_init = false;345 static mpfr::mpreal val(0.5);346 if(!has_init)347 {348 val = ldexp(val, mpfr_get_emax());349 has_init = true;350 }351 return val;352}353 354template <>355inline mpfr::mpreal min_value<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))356{357 static bool has_init = false;358 static mpfr::mpreal val(0.5);359 if(!has_init)360 {361 val = ldexp(val, mpfr_get_emin());362 has_init = true;363 }364 return val;365}366 367template <>368inline mpfr::mpreal log_max_value<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))369{370 static bool has_init = false;371 static mpfr::mpreal val = max_value<mpfr::mpreal>();372 if(!has_init)373 {374 val = log(val);375 has_init = true;376 }377 return val;378}379 380template <>381inline mpfr::mpreal log_min_value<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))382{383 static bool has_init = false;384 static mpfr::mpreal val = max_value<mpfr::mpreal>();385 if(!has_init)386 {387 val = log(val);388 has_init = true;389 }390 return val;391}392 393template <>394inline mpfr::mpreal epsilon<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))395{396 return ldexp(mpfr::mpreal(1), 1-boost::math::policies::digits<mpfr::mpreal, boost::math::policies::policy<> >());397}398 399} // namespace tools400 401template <class Policy>402inline mpfr::mpreal skewness(const extreme_value_distribution<mpfr::mpreal, Policy>& /*dist*/)403{404 //405 // This is 12 * sqrt(6) * zeta(3) / pi^3:406 // See http://mathworld.wolfram.com/ExtremeValueDistribution.html407 //408 #ifdef BOOST_MATH_STANDALONE409 static_assert(sizeof(Policy) == 0, "mpreal skewness can not be calculated in standalone mode");410 #endif411 412 return boost::lexical_cast<mpfr::mpreal>("1.1395470994046486574927930193898461120875997958366");413}414 415template <class Policy>416inline mpfr::mpreal skewness(const rayleigh_distribution<mpfr::mpreal, Policy>& /*dist*/)417{418 // using namespace boost::math::constants;419 #ifdef BOOST_MATH_STANDALONE420 static_assert(sizeof(Policy) == 0, "mpreal skewness can not be calculated in standalone mode");421 #endif422 423 return boost::lexical_cast<mpfr::mpreal>("0.63111065781893713819189935154422777984404221106391");424 // Computed using NTL at 150 bit, about 50 decimal digits.425 // return 2 * root_pi<RealType>() * pi_minus_three<RealType>() / pow23_four_minus_pi<RealType>();426}427 428template <class Policy>429inline mpfr::mpreal kurtosis(const rayleigh_distribution<mpfr::mpreal, Policy>& /*dist*/)430{431 // using namespace boost::math::constants;432 #ifdef BOOST_MATH_STANDALONE433 static_assert(sizeof(Policy) == 0, "mpreal kurtosis can not be calculated in standalone mode");434 #endif435 436 return boost::lexical_cast<mpfr::mpreal>("3.2450893006876380628486604106197544154170667057995");437 // Computed using NTL at 150 bit, about 50 decimal digits.438 // return 3 - (6 * pi<RealType>() * pi<RealType>() - 24 * pi<RealType>() + 16) /439 // (four_minus_pi<RealType>() * four_minus_pi<RealType>());440}441 442template <class Policy>443inline mpfr::mpreal kurtosis_excess(const rayleigh_distribution<mpfr::mpreal, Policy>& /*dist*/)444{445 //using namespace boost::math::constants;446 // Computed using NTL at 150 bit, about 50 decimal digits.447 #ifdef BOOST_MATH_STANDALONE448 static_assert(sizeof(Policy) == 0, "mpreal excess kurtosis can not be calculated in standalone mode");449 #endif450 451 return boost::lexical_cast<mpfr::mpreal>("0.2450893006876380628486604106197544154170667057995");452 // return -(6 * pi<RealType>() * pi<RealType>() - 24 * pi<RealType>() + 16) /453 // (four_minus_pi<RealType>() * four_minus_pi<RealType>());454} // kurtosis455 456namespace detail{457 458//459// Version of Digamma accurate to ~100 decimal digits.460//461template <class Policy>462mpfr::mpreal digamma_imp(mpfr::mpreal x, const std::integral_constant<int, 0>* , const Policy& pol)463{464 //465 // This handles reflection of negative arguments, and all our466 // empfr_classor handling, then forwards to the T-specific approximation.467 //468 BOOST_MATH_STD_USING // ADL of std functions.469 470 mpfr::mpreal result = 0;471 //472 // Check for negative arguments and use reflection:473 //474 if(x < 0)475 {476 // Reflect:477 x = 1 - x;478 // Argument reduction for tan:479 mpfr::mpreal remainder = x - floor(x);480 // Shift to negative if > 0.5:481 if(remainder > 0.5)482 {483 remainder -= 1;484 }485 //486 // check for evaluation at a negative pole:487 //488 if(remainder == 0)489 {490 return policies::raise_pole_error<mpfr::mpreal>("boost::math::digamma<%1%>(%1%)", nullptr, (1-x), pol);491 }492 result = constants::pi<mpfr::mpreal>() / tan(constants::pi<mpfr::mpreal>() * remainder);493 }494 result += big_digamma(x);495 return result;496}497//498// Specialisations of this function provides the initial499// starting guess for Halley iteration:500//501template <class Policy>502mpfr::mpreal erf_inv_imp(const mpfr::mpreal& p, const mpfr::mpreal& q, const Policy&, const std::integral_constant<int, 64>*)503{504 BOOST_MATH_STD_USING // for ADL of std names.505 506 mpfr::mpreal result = 0;507 508 if(p <= 0.5)509 {510 //511 // Evaluate inverse erf using the rational approximation:512 //513 // x = p(p+10)(Y+R(p))514 //515 // Where Y is a constant, and R(p) is optimised for a low516 // absolute empfr_classor compared to |Y|.517 //518 // double: Max empfr_classor found: 2.001849e-18519 // long double: Max empfr_classor found: 1.017064e-20520 // Maximum Deviation Found (actual empfr_classor term at infinite precision) 8.030e-21521 //522 static const float Y = 0.0891314744949340820313f;523 static const mpfr::mpreal P[] = {524 -0.000508781949658280665617,525 -0.00836874819741736770379,526 0.0334806625409744615033,527 -0.0126926147662974029034,528 -0.0365637971411762664006,529 0.0219878681111168899165,530 0.00822687874676915743155,531 -0.00538772965071242932965532 };533 static const mpfr::mpreal Q[] = {534 1,535 -0.970005043303290640362,536 -1.56574558234175846809,537 1.56221558398423026363,538 0.662328840472002992063,539 -0.71228902341542847553,540 -0.0527396382340099713954,541 0.0795283687341571680018,542 -0.00233393759374190016776,543 0.000886216390456424707504544 };545 mpfr::mpreal g = p * (p + 10);546 mpfr::mpreal r = tools::evaluate_polynomial(P, p) / tools::evaluate_polynomial(Q, p);547 result = g * Y + g * r;548 }549 else if(q >= 0.25)550 {551 //552 // Rational approximation for 0.5 > q >= 0.25553 //554 // x = sqrt(-2*log(q)) / (Y + R(q))555 //556 // Where Y is a constant, and R(q) is optimised for a low557 // absolute empfr_classor compared to Y.558 //559 // double : Max empfr_classor found: 7.403372e-17560 // long double : Max empfr_classor found: 6.084616e-20561 // Maximum Deviation Found (empfr_classor term) 4.811e-20562 //563 static const float Y = 2.249481201171875f;564 static const mpfr::mpreal P[] = {565 -0.202433508355938759655,566 0.105264680699391713268,567 8.37050328343119927838,568 17.6447298408374015486,569 -18.8510648058714251895,570 -44.6382324441786960818,571 17.445385985570866523,572 21.1294655448340526258,573 -3.67192254707729348546574 };575 static const mpfr::mpreal Q[] = {576 1,577 6.24264124854247537712,578 3.9713437953343869095,579 -28.6608180499800029974,580 -20.1432634680485188801,581 48.5609213108739935468,582 10.8268667355460159008,583 -22.6436933413139721736,584 1.72114765761200282724585 };586 mpfr::mpreal g = sqrt(-2 * log(q));587 mpfr::mpreal xs = q - 0.25;588 mpfr::mpreal r = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);589 result = g / (Y + r);590 }591 else592 {593 //594 // For q < 0.25 we have a series of rational approximations all595 // of the general form:596 //597 // let: x = sqrt(-log(q))598 //599 // Then the result is given by:600 //601 // x(Y+R(x-B))602 //603 // where Y is a constant, B is the lowest value of x for which604 // the approximation is valid, and R(x-B) is optimised for a low605 // absolute empfr_classor compared to Y.606 //607 // Note that almost all code will really go through the first608 // or maybe second approximation. After than we're dealing with very609 // small input values indeed: 80 and 128 bit long double's go all the610 // way down to ~ 1e-5000 so the "tail" is rather long...611 //612 mpfr::mpreal x = sqrt(-log(q));613 if(x < 3)614 {615 // Max empfr_classor found: 1.089051e-20616 static const float Y = 0.807220458984375f;617 static const mpfr::mpreal P[] = {618 -0.131102781679951906451,619 -0.163794047193317060787,620 0.117030156341995252019,621 0.387079738972604337464,622 0.337785538912035898924,623 0.142869534408157156766,624 0.0290157910005329060432,625 0.00214558995388805277169,626 -0.679465575181126350155e-6,627 0.285225331782217055858e-7,628 -0.681149956853776992068e-9629 };630 static const mpfr::mpreal Q[] = {631 1,632 3.46625407242567245975,633 5.38168345707006855425,634 4.77846592945843778382,635 2.59301921623620271374,636 0.848854343457902036425,637 0.152264338295331783612,638 0.01105924229346489121639 };640 mpfr::mpreal xs = x - 1.125;641 mpfr::mpreal R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);642 result = Y * x + R * x;643 }644 else if(x < 6)645 {646 // Max empfr_classor found: 8.389174e-21647 static const float Y = 0.93995571136474609375f;648 static const mpfr::mpreal P[] = {649 -0.0350353787183177984712,650 -0.00222426529213447927281,651 0.0185573306514231072324,652 0.00950804701325919603619,653 0.00187123492819559223345,654 0.000157544617424960554631,655 0.460469890584317994083e-5,656 -0.230404776911882601748e-9,657 0.266339227425782031962e-11658 };659 static const mpfr::mpreal Q[] = {660 1,661 1.3653349817554063097,662 0.762059164553623404043,663 0.220091105764131249824,664 0.0341589143670947727934,665 0.00263861676657015992959,666 0.764675292302794483503e-4667 };668 mpfr::mpreal xs = x - 3;669 mpfr::mpreal R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);670 result = Y * x + R * x;671 }672 else if(x < 18)673 {674 // Max empfr_classor found: 1.481312e-19675 static const float Y = 0.98362827301025390625f;676 static const mpfr::mpreal P[] = {677 -0.0167431005076633737133,678 -0.00112951438745580278863,679 0.00105628862152492910091,680 0.000209386317487588078668,681 0.149624783758342370182e-4,682 0.449696789927706453732e-6,683 0.462596163522878599135e-8,684 -0.281128735628831791805e-13,685 0.99055709973310326855e-16686 };687 static const mpfr::mpreal Q[] = {688 1,689 0.591429344886417493481,690 0.138151865749083321638,691 0.0160746087093676504695,692 0.000964011807005165528527,693 0.275335474764726041141e-4,694 0.282243172016108031869e-6695 };696 mpfr::mpreal xs = x - 6;697 mpfr::mpreal R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);698 result = Y * x + R * x;699 }700 else if(x < 44)701 {702 // Max empfr_classor found: 5.697761e-20703 static const float Y = 0.99714565277099609375f;704 static const mpfr::mpreal P[] = {705 -0.0024978212791898131227,706 -0.779190719229053954292e-5,707 0.254723037413027451751e-4,708 0.162397777342510920873e-5,709 0.396341011304801168516e-7,710 0.411632831190944208473e-9,711 0.145596286718675035587e-11,712 -0.116765012397184275695e-17713 };714 static const mpfr::mpreal Q[] = {715 1,716 0.207123112214422517181,717 0.0169410838120975906478,718 0.000690538265622684595676,719 0.145007359818232637924e-4,720 0.144437756628144157666e-6,721 0.509761276599778486139e-9722 };723 mpfr::mpreal xs = x - 18;724 mpfr::mpreal R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);725 result = Y * x + R * x;726 }727 else728 {729 // Max empfr_classor found: 1.279746e-20730 static const float Y = 0.99941349029541015625f;731 static const mpfr::mpreal P[] = {732 -0.000539042911019078575891,733 -0.28398759004727721098e-6,734 0.899465114892291446442e-6,735 0.229345859265920864296e-7,736 0.225561444863500149219e-9,737 0.947846627503022684216e-12,738 0.135880130108924861008e-14,739 -0.348890393399948882918e-21740 };741 static const mpfr::mpreal Q[] = {742 1,743 0.0845746234001899436914,744 0.00282092984726264681981,745 0.468292921940894236786e-4,746 0.399968812193862100054e-6,747 0.161809290887904476097e-8,748 0.231558608310259605225e-11749 };750 mpfr::mpreal xs = x - 44;751 mpfr::mpreal R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);752 result = Y * x + R * x;753 }754 }755 return result;756}757 758inline mpfr::mpreal bessel_i0(mpfr::mpreal x)759{760 #ifdef BOOST_MATH_STANDALONE761 static_assert(sizeof(x) == 0, "mpreal bessel_i0 can not be calculated in standalone mode");762 #endif763 764 static const mpfr::mpreal P1[] = {765 boost::lexical_cast<mpfr::mpreal>("-2.2335582639474375249e+15"),766 boost::lexical_cast<mpfr::mpreal>("-5.5050369673018427753e+14"),767 boost::lexical_cast<mpfr::mpreal>("-3.2940087627407749166e+13"),768 boost::lexical_cast<mpfr::mpreal>("-8.4925101247114157499e+11"),769 boost::lexical_cast<mpfr::mpreal>("-1.1912746104985237192e+10"),770 boost::lexical_cast<mpfr::mpreal>("-1.0313066708737980747e+08"),771 boost::lexical_cast<mpfr::mpreal>("-5.9545626019847898221e+05"),772 boost::lexical_cast<mpfr::mpreal>("-2.4125195876041896775e+03"),773 boost::lexical_cast<mpfr::mpreal>("-7.0935347449210549190e+00"),774 boost::lexical_cast<mpfr::mpreal>("-1.5453977791786851041e-02"),775 boost::lexical_cast<mpfr::mpreal>("-2.5172644670688975051e-05"),776 boost::lexical_cast<mpfr::mpreal>("-3.0517226450451067446e-08"),777 boost::lexical_cast<mpfr::mpreal>("-2.6843448573468483278e-11"),778 boost::lexical_cast<mpfr::mpreal>("-1.5982226675653184646e-14"),779 boost::lexical_cast<mpfr::mpreal>("-5.2487866627945699800e-18"),780 };781 static const mpfr::mpreal Q1[] = {782 boost::lexical_cast<mpfr::mpreal>("-2.2335582639474375245e+15"),783 boost::lexical_cast<mpfr::mpreal>("7.8858692566751002988e+12"),784 boost::lexical_cast<mpfr::mpreal>("-1.2207067397808979846e+10"),785 boost::lexical_cast<mpfr::mpreal>("1.0377081058062166144e+07"),786 boost::lexical_cast<mpfr::mpreal>("-4.8527560179962773045e+03"),787 boost::lexical_cast<mpfr::mpreal>("1.0"),788 };789 static const mpfr::mpreal P2[] = {790 boost::lexical_cast<mpfr::mpreal>("-2.2210262233306573296e-04"),791 boost::lexical_cast<mpfr::mpreal>("1.3067392038106924055e-02"),792 boost::lexical_cast<mpfr::mpreal>("-4.4700805721174453923e-01"),793 boost::lexical_cast<mpfr::mpreal>("5.5674518371240761397e+00"),794 boost::lexical_cast<mpfr::mpreal>("-2.3517945679239481621e+01"),795 boost::lexical_cast<mpfr::mpreal>("3.1611322818701131207e+01"),796 boost::lexical_cast<mpfr::mpreal>("-9.6090021968656180000e+00"),797 };798 static const mpfr::mpreal Q2[] = {799 boost::lexical_cast<mpfr::mpreal>("-5.5194330231005480228e-04"),800 boost::lexical_cast<mpfr::mpreal>("3.2547697594819615062e-02"),801 boost::lexical_cast<mpfr::mpreal>("-1.1151759188741312645e+00"),802 boost::lexical_cast<mpfr::mpreal>("1.3982595353892851542e+01"),803 boost::lexical_cast<mpfr::mpreal>("-6.0228002066743340583e+01"),804 boost::lexical_cast<mpfr::mpreal>("8.5539563258012929600e+01"),805 boost::lexical_cast<mpfr::mpreal>("-3.1446690275135491500e+01"),806 boost::lexical_cast<mpfr::mpreal>("1.0"),807 };808 mpfr::mpreal value, factor, r;809 810 BOOST_MATH_STD_USING811 using namespace boost::math::tools;812 813 if (x < 0)814 {815 x = -x; // even function816 }817 if (x == 0)818 {819 return static_cast<mpfr::mpreal>(1);820 }821 if (x <= 15) // x in (0, 15]822 {823 mpfr::mpreal y = x * x;824 value = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);825 }826 else // x in (15, \infty)827 {828 mpfr::mpreal y = 1 / x - mpfr::mpreal(1) / 15;829 r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);830 factor = exp(x) / sqrt(x);831 value = factor * r;832 }833 834 return value;835}836 837inline mpfr::mpreal bessel_i1(mpfr::mpreal x)838{839 static const mpfr::mpreal P1[] = {840 static_cast<mpfr::mpreal>("-1.4577180278143463643e+15"),841 static_cast<mpfr::mpreal>("-1.7732037840791591320e+14"),842 static_cast<mpfr::mpreal>("-6.9876779648010090070e+12"),843 static_cast<mpfr::mpreal>("-1.3357437682275493024e+11"),844 static_cast<mpfr::mpreal>("-1.4828267606612366099e+09"),845 static_cast<mpfr::mpreal>("-1.0588550724769347106e+07"),846 static_cast<mpfr::mpreal>("-5.1894091982308017540e+04"),847 static_cast<mpfr::mpreal>("-1.8225946631657315931e+02"),848 static_cast<mpfr::mpreal>("-4.7207090827310162436e-01"),849 static_cast<mpfr::mpreal>("-9.1746443287817501309e-04"),850 static_cast<mpfr::mpreal>("-1.3466829827635152875e-06"),851 static_cast<mpfr::mpreal>("-1.4831904935994647675e-09"),852 static_cast<mpfr::mpreal>("-1.1928788903603238754e-12"),853 static_cast<mpfr::mpreal>("-6.5245515583151902910e-16"),854 static_cast<mpfr::mpreal>("-1.9705291802535139930e-19"),855 };856 static const mpfr::mpreal Q1[] = {857 static_cast<mpfr::mpreal>("-2.9154360556286927285e+15"),858 static_cast<mpfr::mpreal>("9.7887501377547640438e+12"),859 static_cast<mpfr::mpreal>("-1.4386907088588283434e+10"),860 static_cast<mpfr::mpreal>("1.1594225856856884006e+07"),861 static_cast<mpfr::mpreal>("-5.1326864679904189920e+03"),862 static_cast<mpfr::mpreal>("1.0"),863 };864 static const mpfr::mpreal P2[] = {865 static_cast<mpfr::mpreal>("1.4582087408985668208e-05"),866 static_cast<mpfr::mpreal>("-8.9359825138577646443e-04"),867 static_cast<mpfr::mpreal>("2.9204895411257790122e-02"),868 static_cast<mpfr::mpreal>("-3.4198728018058047439e-01"),869 static_cast<mpfr::mpreal>("1.3960118277609544334e+00"),870 static_cast<mpfr::mpreal>("-1.9746376087200685843e+00"),871 static_cast<mpfr::mpreal>("8.5591872901933459000e-01"),872 static_cast<mpfr::mpreal>("-6.0437159056137599999e-02"),873 };874 static const mpfr::mpreal Q2[] = {875 static_cast<mpfr::mpreal>("3.7510433111922824643e-05"),876 static_cast<mpfr::mpreal>("-2.2835624489492512649e-03"),877 static_cast<mpfr::mpreal>("7.4212010813186530069e-02"),878 static_cast<mpfr::mpreal>("-8.5017476463217924408e-01"),879 static_cast<mpfr::mpreal>("3.2593714889036996297e+00"),880 static_cast<mpfr::mpreal>("-3.8806586721556593450e+00"),881 static_cast<mpfr::mpreal>("1.0"),882 };883 mpfr::mpreal value, factor, r, w;884 885 BOOST_MATH_STD_USING886 using namespace boost::math::tools;887 888 w = abs(x);889 if (x == 0)890 {891 return static_cast<mpfr::mpreal>(0);892 }893 if (w <= 15) // w in (0, 15]894 {895 mpfr::mpreal y = x * x;896 r = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);897 factor = w;898 value = factor * r;899 }900 else // w in (15, \infty)901 {902 mpfr::mpreal y = 1 / w - mpfr::mpreal(1) / 15;903 r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);904 factor = exp(w) / sqrt(w);905 value = factor * r;906 }907 908 if (x < 0)909 {910 value *= -value; // odd function911 }912 return value;913}914 915} // namespace detail916} // namespace math917 918}919 920#endif // BOOST_MATH_MPLFR_BINDINGS_HPP921 922