math/bindings/mpreal.hpp

0001 //  Copyright John Maddock 2008.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005 //
0006 // Wrapper that works with mpfr::mpreal defined in gmpfrxx.h
0007 // See http://math.berkeley.edu/~wilken/code/gmpfrxx/
0008 // Also requires the gmp and mpfr libraries.
0009 //
0010
0011 #ifndef BOOST_MATH_MPREAL_BINDINGS_HPP
0012 #define BOOST_MATH_MPREAL_BINDINGS_HPP
0013
0014 #include <type_traits>
0015
0016 #ifdef _MSC_VER
0017 //
0018 // We get a lot of warnings from the gmp, mpfr and gmpfrxx headers,
0019 // disable them here, so we only see warnings from *our* code:
0020 //
0021 #pragma warning(push)
0022 #pragma warning(disable: 4127 4800 4512)
0023 #endif
0024
0025 #include <mpreal.h>
0026
0027 #ifdef _MSC_VER
0028 #pragma warning(pop)
0029 #endif
0030
0031 #include <boost/math/tools/precision.hpp>
0032 #include <boost/math/tools/real_cast.hpp>
0033 #include <boost/math/policies/policy.hpp>
0034 #include <boost/math/distributions/fwd.hpp>
0035 #include <boost/math/special_functions/math_fwd.hpp>
0036 #include <boost/math/bindings/detail/big_digamma.hpp>
0037 #include <boost/math/bindings/detail/big_lanczos.hpp>
0038 #include <boost/math/tools/config.hpp>
0039 #ifndef BOOST_MATH_STANDALONE
0040 #include <boost/lexical_cast.hpp>
0041 #endif
0042
0043 namespace mpfr{
0044
0045 template <class T>
0046 inline mpreal operator + (const mpreal& r, const T& t){ return r + mpreal(t); }
0047 template <class T>
0048 inline mpreal operator - (const mpreal& r, const T& t){ return r - mpreal(t); }
0049 template <class T>
0050 inline mpreal operator * (const mpreal& r, const T& t){ return r * mpreal(t); }
0051 template <class T>
0052 inline mpreal operator / (const mpreal& r, const T& t){ return r / mpreal(t); }
0053
0054 template <class T>
0055 inline mpreal operator + (const T& t, const mpreal& r){ return mpreal(t) + r; }
0056 template <class T>
0057 inline mpreal operator - (const T& t, const mpreal& r){ return mpreal(t) - r; }
0058 template <class T>
0059 inline mpreal operator * (const T& t, const mpreal& r){ return mpreal(t) * r; }
0060 template <class T>
0061 inline mpreal operator / (const T& t, const mpreal& r){ return mpreal(t) / r; }
0062
0063 template <class T>
0064 inline bool operator == (const mpreal& r, const T& t){ return r == mpreal(t); }
0065 template <class T>
0066 inline bool operator != (const mpreal& r, const T& t){ return r != mpreal(t); }
0067 template <class T>
0068 inline bool operator <= (const mpreal& r, const T& t){ return r <= mpreal(t); }
0069 template <class T>
0070 inline bool operator >= (const mpreal& r, const T& t){ return r >= mpreal(t); }
0071 template <class T>
0072 inline bool operator < (const mpreal& r, const T& t){ return r < mpreal(t); }
0073 template <class T>
0074 inline bool operator > (const mpreal& r, const T& t){ return r > mpreal(t); }
0075
0076 template <class T>
0077 inline bool operator == (const T& t, const mpreal& r){ return mpreal(t) == r; }
0078 template <class T>
0079 inline bool operator != (const T& t, const mpreal& r){ return mpreal(t) != r; }
0080 template <class T>
0081 inline bool operator <= (const T& t, const mpreal& r){ return mpreal(t) <= r; }
0082 template <class T>
0083 inline bool operator >= (const T& t, const mpreal& r){ return mpreal(t) >= r; }
0084 template <class T>
0085 inline bool operator < (const T& t, const mpreal& r){ return mpreal(t) < r; }
0086 template <class T>
0087 inline bool operator > (const T& t, const mpreal& r){ return mpreal(t) > r; }
0088
0089 /*
0090 inline mpfr::mpreal fabs(const mpfr::mpreal& v)
0091 {
0092    return abs(v);
0093 }
0094 inline mpfr::mpreal pow(const mpfr::mpreal& b, const mpfr::mpreal e)
0095 {
0096    mpfr::mpreal result;
0097    mpfr_pow(result.__get_mp(), b.__get_mp(), e.__get_mp(), GMP_RNDN);
0098    return result;
0099 }
0100 */
0101 inline mpfr::mpreal ldexp(const mpfr::mpreal& v, int e)
0102 {
0103    return mpfr::ldexp(v, static_cast<mp_exp_t>(e));
0104 }
0105
0106 inline mpfr::mpreal frexp(const mpfr::mpreal& v, int* expon)
0107 {
0108    mp_exp_t e;
0109    mpfr::mpreal r = mpfr::frexp(v, &e);
0110    *expon = e;
0111    return r;
0112 }
0113
0114 #if (MPFR_VERSION < MPFR_VERSION_NUM(2,4,0))
0115 mpfr::mpreal fmod(const mpfr::mpreal& v1, const mpfr::mpreal& v2)
0116 {
0117    mpfr::mpreal n;
0118    if(v1 < 0)
0119       n = ceil(v1 / v2);
0120    else
0121       n = floor(v1 / v2);
0122    return v1 - n * v2;
0123 }
0124 #endif
0125
0126 template <class Policy>
0127 inline mpfr::mpreal modf(const mpfr::mpreal& v, long long* ipart, const Policy& pol)
0128 {
0129    *ipart = lltrunc(v, pol);
0130    return v - boost::math::tools::real_cast<mpfr::mpreal>(*ipart);
0131 }
0132 template <class Policy>
0133 inline int iround(mpfr::mpreal const& x, const Policy& pol)
0134 {
0135    return boost::math::tools::real_cast<int>(boost::math::round(x, pol));
0136 }
0137
0138 template <class Policy>
0139 inline long lround(mpfr::mpreal const& x, const Policy& pol)
0140 {
0141    return boost::math::tools::real_cast<long>(boost::math::round(x, pol));
0142 }
0143
0144 template <class Policy>
0145 inline long long llround(mpfr::mpreal const& x, const Policy& pol)
0146 {
0147    return boost::math::tools::real_cast<long long>(boost::math::round(x, pol));
0148 }
0149
0150 template <class Policy>
0151 inline int itrunc(mpfr::mpreal const& x, const Policy& pol)
0152 {
0153    return boost::math::tools::real_cast<int>(boost::math::trunc(x, pol));
0154 }
0155
0156 template <class Policy>
0157 inline long ltrunc(mpfr::mpreal const& x, const Policy& pol)
0158 {
0159    return boost::math::tools::real_cast<long>(boost::math::trunc(x, pol));
0160 }
0161
0162 template <class Policy>
0163 inline long long lltrunc(mpfr::mpreal const& x, const Policy& pol)
0164 {
0165    return boost::math::tools::real_cast<long long>(boost::math::trunc(x, pol));
0166 }
0167
0168 }
0169
0170 namespace boost{ namespace math{
0171
0172 #if defined(__GNUC__) && (__GNUC__ < 4)
0173    using ::iround;
0174    using ::lround;
0175    using ::llround;
0176    using ::itrunc;
0177    using ::ltrunc;
0178    using ::lltrunc;
0179    using ::modf;
0180 #endif
0181
0182 namespace lanczos{
0183
0184 struct mpreal_lanczos
0185 {
0186    static mpfr::mpreal lanczos_sum(const mpfr::mpreal& z)
0187    {
0188       unsigned long p = z.get_default_prec();
0189       if(p <= 72)
0190          return lanczos13UDT::lanczos_sum(z);
0191       else if(p <= 120)
0192          return lanczos22UDT::lanczos_sum(z);
0193       else if(p <= 170)
0194          return lanczos31UDT::lanczos_sum(z);
0195       else //if(p <= 370) approx 100 digit precision:
0196          return lanczos61UDT::lanczos_sum(z);
0197    }
0198    static mpfr::mpreal lanczos_sum_expG_scaled(const mpfr::mpreal& z)
0199    {
0200       unsigned long p = z.get_default_prec();
0201       if(p <= 72)
0202          return lanczos13UDT::lanczos_sum_expG_scaled(z);
0203       else if(p <= 120)
0204          return lanczos22UDT::lanczos_sum_expG_scaled(z);
0205       else if(p <= 170)
0206          return lanczos31UDT::lanczos_sum_expG_scaled(z);
0207       else //if(p <= 370) approx 100 digit precision:
0208          return lanczos61UDT::lanczos_sum_expG_scaled(z);
0209    }
0210    static mpfr::mpreal lanczos_sum_near_1(const mpfr::mpreal& z)
0211    {
0212       unsigned long p = z.get_default_prec();
0213       if(p <= 72)
0214          return lanczos13UDT::lanczos_sum_near_1(z);
0215       else if(p <= 120)
0216          return lanczos22UDT::lanczos_sum_near_1(z);
0217       else if(p <= 170)
0218          return lanczos31UDT::lanczos_sum_near_1(z);
0219       else //if(p <= 370) approx 100 digit precision:
0220          return lanczos61UDT::lanczos_sum_near_1(z);
0221    }
0222    static mpfr::mpreal lanczos_sum_near_2(const mpfr::mpreal& z)
0223    {
0224       unsigned long p = z.get_default_prec();
0225       if(p <= 72)
0226          return lanczos13UDT::lanczos_sum_near_2(z);
0227       else if(p <= 120)
0228          return lanczos22UDT::lanczos_sum_near_2(z);
0229       else if(p <= 170)
0230          return lanczos31UDT::lanczos_sum_near_2(z);
0231       else //if(p <= 370) approx 100 digit precision:
0232          return lanczos61UDT::lanczos_sum_near_2(z);
0233    }
0234    static mpfr::mpreal g()
0235    {
0236       unsigned long p = mpfr::mpreal::get_default_prec();
0237       if(p <= 72)
0238          return lanczos13UDT::g();
0239       else if(p <= 120)
0240          return lanczos22UDT::g();
0241       else if(p <= 170)
0242          return lanczos31UDT::g();
0243       else //if(p <= 370) approx 100 digit precision:
0244          return lanczos61UDT::g();
0245    }
0246 };
0247
0248 template<class Policy>
0249 struct lanczos<mpfr::mpreal, Policy>
0250 {
0251    typedef mpreal_lanczos type;
0252 };
0253
0254 } // namespace lanczos
0255
0256 namespace tools
0257 {
0258
0259 template<>
0260 inline int digits<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))
0261 {
0262    return mpfr::mpreal::get_default_prec();
0263 }
0264
0265 namespace detail{
0266
0267 template<class I>
0268 void convert_to_long_result(mpfr::mpreal const& r, I& result)
0269 {
0270    result = 0;
0271    I last_result(0);
0272    mpfr::mpreal t(r);
0273    double term;
0274    do
0275    {
0276       term = real_cast<double>(t);
0277       last_result = result;
0278       result += static_cast<I>(term);
0279       t -= term;
0280    }while(result != last_result);
0281 }
0282
0283 }
0284
0285 template <>
0286 inline mpfr::mpreal real_cast<mpfr::mpreal, long long>(long long t)
0287 {
0288    mpfr::mpreal result;
0289    int expon = 0;
0290    int sign = 1;
0291    if(t < 0)
0292    {
0293       sign = -1;
0294       t = -t;
0295    }
0296    while(t)
0297    {
0298       result += ldexp(static_cast<double>(t & 0xffffL), expon);
0299       expon += 32;
0300       t >>= 32;
0301    }
0302    return result * sign;
0303 }
0304 /*
0305 template <>
0306 inline unsigned real_cast<unsigned, mpfr::mpreal>(mpfr::mpreal t)
0307 {
0308    return t.get_ui();
0309 }
0310 template <>
0311 inline int real_cast<int, mpfr::mpreal>(mpfr::mpreal t)
0312 {
0313    return t.get_si();
0314 }
0315 template <>
0316 inline double real_cast<double, mpfr::mpreal>(mpfr::mpreal t)
0317 {
0318    return t.get_d();
0319 }
0320 template <>
0321 inline float real_cast<float, mpfr::mpreal>(mpfr::mpreal t)
0322 {
0323    return static_cast<float>(t.get_d());
0324 }
0325 template <>
0326 inline long real_cast<long, mpfr::mpreal>(mpfr::mpreal t)
0327 {
0328    long result;
0329    detail::convert_to_long_result(t, result);
0330    return result;
0331 }
0332 */
0333 template <>
0334 inline long long real_cast<long long, mpfr::mpreal>(mpfr::mpreal t)
0335 {
0336    long long result;
0337    detail::convert_to_long_result(t, result);
0338    return result;
0339 }
0340
0341 template <>
0342 inline mpfr::mpreal max_value<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))
0343 {
0344    static bool has_init = false;
0345    static mpfr::mpreal val(0.5);
0346    if(!has_init)
0347    {
0348       val = ldexp(val, mpfr_get_emax());
0349       has_init = true;
0350    }
0351    return val;
0352 }
0353
0354 template <>
0355 inline mpfr::mpreal min_value<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))
0356 {
0357    static bool has_init = false;
0358    static mpfr::mpreal val(0.5);
0359    if(!has_init)
0360    {
0361       val = ldexp(val, mpfr_get_emin());
0362       has_init = true;
0363    }
0364    return val;
0365 }
0366
0367 template <>
0368 inline mpfr::mpreal log_max_value<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))
0369 {
0370    static bool has_init = false;
0371    static mpfr::mpreal val = max_value<mpfr::mpreal>();
0372    if(!has_init)
0373    {
0374       val = log(val);
0375       has_init = true;
0376    }
0377    return val;
0378 }
0379
0380 template <>
0381 inline mpfr::mpreal log_min_value<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))
0382 {
0383    static bool has_init = false;
0384    static mpfr::mpreal val = max_value<mpfr::mpreal>();
0385    if(!has_init)
0386    {
0387       val = log(val);
0388       has_init = true;
0389    }
0390    return val;
0391 }
0392
0393 template <>
0394 inline mpfr::mpreal epsilon<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))
0395 {
0396    return ldexp(mpfr::mpreal(1), 1-boost::math::policies::digits<mpfr::mpreal, boost::math::policies::policy<> >());
0397 }
0398
0399 } // namespace tools
0400
0401 template <class Policy>
0402 inline mpfr::mpreal skewness(const extreme_value_distribution<mpfr::mpreal, Policy>& /*dist*/)
0403 {
0404    //
0405    // This is 12 * sqrt(6) * zeta(3) / pi^3:
0406    // See http://mathworld.wolfram.com/ExtremeValueDistribution.html
0407    //
0408    #ifdef BOOST_MATH_STANDALONE
0409    static_assert(sizeof(Policy) == 0, "mpreal skewness can not be calculated in standalone mode");
0410    #endif
0411
0412    return boost::lexical_cast<mpfr::mpreal>("1.1395470994046486574927930193898461120875997958366");
0413 }
0414
0415 template <class Policy>
0416 inline mpfr::mpreal skewness(const rayleigh_distribution<mpfr::mpreal, Policy>& /*dist*/)
0417 {
0418   // using namespace boost::math::constants;
0419   #ifdef BOOST_MATH_STANDALONE
0420   static_assert(sizeof(Policy) == 0, "mpreal skewness can not be calculated in standalone mode");
0421   #endif
0422
0423   return boost::lexical_cast<mpfr::mpreal>("0.63111065781893713819189935154422777984404221106391");
0424   // Computed using NTL at 150 bit, about 50 decimal digits.
0425   // return 2 * root_pi<RealType>() * pi_minus_three<RealType>() / pow23_four_minus_pi<RealType>();
0426 }
0427
0428 template <class Policy>
0429 inline mpfr::mpreal kurtosis(const rayleigh_distribution<mpfr::mpreal, Policy>& /*dist*/)
0430 {
0431   // using namespace boost::math::constants;
0432   #ifdef BOOST_MATH_STANDALONE
0433   static_assert(sizeof(Policy) == 0, "mpreal kurtosis can not be calculated in standalone mode");
0434   #endif
0435
0436   return boost::lexical_cast<mpfr::mpreal>("3.2450893006876380628486604106197544154170667057995");
0437   // Computed using NTL at 150 bit, about 50 decimal digits.
0438   // return 3 - (6 * pi<RealType>() * pi<RealType>() - 24 * pi<RealType>() + 16) /
0439   // (four_minus_pi<RealType>() * four_minus_pi<RealType>());
0440 }
0441
0442 template <class Policy>
0443 inline mpfr::mpreal kurtosis_excess(const rayleigh_distribution<mpfr::mpreal, Policy>& /*dist*/)
0444 {
0445   //using namespace boost::math::constants;
0446   // Computed using NTL at 150 bit, about 50 decimal digits.
0447   #ifdef BOOST_MATH_STANDALONE
0448   static_assert(sizeof(Policy) == 0, "mpreal excess kurtosis can not be calculated in standalone mode");
0449   #endif
0450
0451   return boost::lexical_cast<mpfr::mpreal>("0.2450893006876380628486604106197544154170667057995");
0452   // return -(6 * pi<RealType>() * pi<RealType>() - 24 * pi<RealType>() + 16) /
0453   //   (four_minus_pi<RealType>() * four_minus_pi<RealType>());
0454 } // kurtosis
0455
0456 namespace detail{
0457
0458 //
0459 // Version of Digamma accurate to ~100 decimal digits.
0460 //
0461 template <class Policy>
0462 mpfr::mpreal digamma_imp(mpfr::mpreal x, const std::integral_constant<int, 0>* , const Policy& pol)
0463 {
0464    //
0465    // This handles reflection of negative arguments, and all our
0466    // empfr_classor handling, then forwards to the T-specific approximation.
0467    //
0468    BOOST_MATH_STD_USING // ADL of std functions.
0469
0470    mpfr::mpreal result = 0;
0471    //
0472    // Check for negative arguments and use reflection:
0473    //
0474    if(x < 0)
0475    {
0476       // Reflect:
0477       x = 1 - x;
0478       // Argument reduction for tan:
0479       mpfr::mpreal remainder = x - floor(x);
0480       // Shift to negative if > 0.5:
0481       if(remainder > 0.5)
0482       {
0483          remainder -= 1;
0484       }
0485       //
0486       // check for evaluation at a negative pole:
0487       //
0488       if(remainder == 0)
0489       {
0490          return policies::raise_pole_error<mpfr::mpreal>("boost::math::digamma<%1%>(%1%)", nullptr, (1-x), pol);
0491       }
0492       result = constants::pi<mpfr::mpreal>() / tan(constants::pi<mpfr::mpreal>() * remainder);
0493    }
0494    result += big_digamma(x);
0495    return result;
0496 }
0497 //
0498 // Specialisations of this function provides the initial
0499 // starting guess for Halley iteration:
0500 //
0501 template <class Policy>
0502 mpfr::mpreal erf_inv_imp(const mpfr::mpreal& p, const mpfr::mpreal& q, const Policy&, const std::integral_constant<int, 64>*)
0503 {
0504    BOOST_MATH_STD_USING // for ADL of std names.
0505
0506    mpfr::mpreal result = 0;
0507
0508    if(p <= 0.5)
0509    {
0510       //
0511       // Evaluate inverse erf using the rational approximation:
0512       //
0513       // x = p(p+10)(Y+R(p))
0514       //
0515       // Where Y is a constant, and R(p) is optimised for a low
0516       // absolute empfr_classor compared to |Y|.
0517       //
0518       // double: Max empfr_classor found: 2.001849e-18
0519       // long double: Max empfr_classor found: 1.017064e-20
0520       // Maximum Deviation Found (actual empfr_classor term at infinite precision) 8.030e-21
0521       //
0522       static const float Y = 0.0891314744949340820313f;
0523       static const mpfr::mpreal P[] = {
0524          -0.000508781949658280665617,
0525          -0.00836874819741736770379,
0526          0.0334806625409744615033,
0527          -0.0126926147662974029034,
0528          -0.0365637971411762664006,
0529          0.0219878681111168899165,
0530          0.00822687874676915743155,
0531          -0.00538772965071242932965
0532       };
0533       static const mpfr::mpreal Q[] = {
0534          1,
0535          -0.970005043303290640362,
0536          -1.56574558234175846809,
0537          1.56221558398423026363,
0538          0.662328840472002992063,
0539          -0.71228902341542847553,
0540          -0.0527396382340099713954,
0541          0.0795283687341571680018,
0542          -0.00233393759374190016776,
0543          0.000886216390456424707504
0544       };
0545       mpfr::mpreal g = p * (p + 10);
0546       mpfr::mpreal r = tools::evaluate_polynomial(P, p) / tools::evaluate_polynomial(Q, p);
0547       result = g * Y + g * r;
0548    }
0549    else if(q >= 0.25)
0550    {
0551       //
0552       // Rational approximation for 0.5 > q >= 0.25
0553       //
0554       // x = sqrt(-2*log(q)) / (Y + R(q))
0555       //
0556       // Where Y is a constant, and R(q) is optimised for a low
0557       // absolute empfr_classor compared to Y.
0558       //
0559       // double : Max empfr_classor found: 7.403372e-17
0560       // long double : Max empfr_classor found: 6.084616e-20
0561       // Maximum Deviation Found (empfr_classor term) 4.811e-20
0562       //
0563       static const float Y = 2.249481201171875f;
0564       static const mpfr::mpreal P[] = {
0565          -0.202433508355938759655,
0566          0.105264680699391713268,
0567          8.37050328343119927838,
0568          17.6447298408374015486,
0569          -18.8510648058714251895,
0570          -44.6382324441786960818,
0571          17.445385985570866523,
0572          21.1294655448340526258,
0573          -3.67192254707729348546
0574       };
0575       static const mpfr::mpreal Q[] = {
0576          1,
0577          6.24264124854247537712,
0578          3.9713437953343869095,
0579          -28.6608180499800029974,
0580          -20.1432634680485188801,
0581          48.5609213108739935468,
0582          10.8268667355460159008,
0583          -22.6436933413139721736,
0584          1.72114765761200282724
0585       };
0586       mpfr::mpreal g = sqrt(-2 * log(q));
0587       mpfr::mpreal xs = q - 0.25;
0588       mpfr::mpreal r = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0589       result = g / (Y + r);
0590    }
0591    else
0592    {
0593       //
0594       // For q < 0.25 we have a series of rational approximations all
0595       // of the general form:
0596       //
0597       // let: x = sqrt(-log(q))
0598       //
0599       // Then the result is given by:
0600       //
0601       // x(Y+R(x-B))
0602       //
0603       // where Y is a constant, B is the lowest value of x for which
0604       // the approximation is valid, and R(x-B) is optimised for a low
0605       // absolute empfr_classor compared to Y.
0606       //
0607       // Note that almost all code will really go through the first
0608       // or maybe second approximation.  After than we're dealing with very
0609       // small input values indeed: 80 and 128 bit long double's go all the
0610       // way down to ~ 1e-5000 so the "tail" is rather long...
0611       //
0612       mpfr::mpreal x = sqrt(-log(q));
0613       if(x < 3)
0614       {
0615          // Max empfr_classor found: 1.089051e-20
0616          static const float Y = 0.807220458984375f;
0617          static const mpfr::mpreal P[] = {
0618             -0.131102781679951906451,
0619             -0.163794047193317060787,
0620             0.117030156341995252019,
0621             0.387079738972604337464,
0622             0.337785538912035898924,
0623             0.142869534408157156766,
0624             0.0290157910005329060432,
0625             0.00214558995388805277169,
0626             -0.679465575181126350155e-6,
0627             0.285225331782217055858e-7,
0628             -0.681149956853776992068e-9
0629          };
0630          static const mpfr::mpreal Q[] = {
0631             1,
0632             3.46625407242567245975,
0633             5.38168345707006855425,
0634             4.77846592945843778382,
0635             2.59301921623620271374,
0636             0.848854343457902036425,
0637             0.152264338295331783612,
0638             0.01105924229346489121
0639          };
0640          mpfr::mpreal xs = x - 1.125;
0641          mpfr::mpreal R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0642          result = Y * x + R * x;
0643       }
0644       else if(x < 6)
0645       {
0646          // Max empfr_classor found: 8.389174e-21
0647          static const float Y = 0.93995571136474609375f;
0648          static const mpfr::mpreal P[] = {
0649             -0.0350353787183177984712,
0650             -0.00222426529213447927281,
0651             0.0185573306514231072324,
0652             0.00950804701325919603619,
0653             0.00187123492819559223345,
0654             0.000157544617424960554631,
0655             0.460469890584317994083e-5,
0656             -0.230404776911882601748e-9,
0657             0.266339227425782031962e-11
0658          };
0659          static const mpfr::mpreal Q[] = {
0660             1,
0661             1.3653349817554063097,
0662             0.762059164553623404043,
0663             0.220091105764131249824,
0664             0.0341589143670947727934,
0665             0.00263861676657015992959,
0666             0.764675292302794483503e-4
0667          };
0668          mpfr::mpreal xs = x - 3;
0669          mpfr::mpreal R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0670          result = Y * x + R * x;
0671       }
0672       else if(x < 18)
0673       {
0674          // Max empfr_classor found: 1.481312e-19
0675          static const float Y = 0.98362827301025390625f;
0676          static const mpfr::mpreal P[] = {
0677             -0.0167431005076633737133,
0678             -0.00112951438745580278863,
0679             0.00105628862152492910091,
0680             0.000209386317487588078668,
0681             0.149624783758342370182e-4,
0682             0.449696789927706453732e-6,
0683             0.462596163522878599135e-8,
0684             -0.281128735628831791805e-13,
0685             0.99055709973310326855e-16
0686          };
0687          static const mpfr::mpreal Q[] = {
0688             1,
0689             0.591429344886417493481,
0690             0.138151865749083321638,
0691             0.0160746087093676504695,
0692             0.000964011807005165528527,
0693             0.275335474764726041141e-4,
0694             0.282243172016108031869e-6
0695          };
0696          mpfr::mpreal xs = x - 6;
0697          mpfr::mpreal R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0698          result = Y * x + R * x;
0699       }
0700       else if(x < 44)
0701       {
0702          // Max empfr_classor found: 5.697761e-20
0703          static const float Y = 0.99714565277099609375f;
0704          static const mpfr::mpreal P[] = {
0705             -0.0024978212791898131227,
0706             -0.779190719229053954292e-5,
0707             0.254723037413027451751e-4,
0708             0.162397777342510920873e-5,
0709             0.396341011304801168516e-7,
0710             0.411632831190944208473e-9,
0711             0.145596286718675035587e-11,
0712             -0.116765012397184275695e-17
0713          };
0714          static const mpfr::mpreal Q[] = {
0715             1,
0716             0.207123112214422517181,
0717             0.0169410838120975906478,
0718             0.000690538265622684595676,
0719             0.145007359818232637924e-4,
0720             0.144437756628144157666e-6,
0721             0.509761276599778486139e-9
0722          };
0723          mpfr::mpreal xs = x - 18;
0724          mpfr::mpreal R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0725          result = Y * x + R * x;
0726       }
0727       else
0728       {
0729          // Max empfr_classor found: 1.279746e-20
0730          static const float Y = 0.99941349029541015625f;
0731          static const mpfr::mpreal P[] = {
0732             -0.000539042911019078575891,
0733             -0.28398759004727721098e-6,
0734             0.899465114892291446442e-6,
0735             0.229345859265920864296e-7,
0736             0.225561444863500149219e-9,
0737             0.947846627503022684216e-12,
0738             0.135880130108924861008e-14,
0739             -0.348890393399948882918e-21
0740          };
0741          static const mpfr::mpreal Q[] = {
0742             1,
0743             0.0845746234001899436914,
0744             0.00282092984726264681981,
0745             0.468292921940894236786e-4,
0746             0.399968812193862100054e-6,
0747             0.161809290887904476097e-8,
0748             0.231558608310259605225e-11
0749          };
0750          mpfr::mpreal xs = x - 44;
0751          mpfr::mpreal R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0752          result = Y * x + R * x;
0753       }
0754    }
0755    return result;
0756 }
0757
0758 inline mpfr::mpreal bessel_i0(mpfr::mpreal x)
0759 {
0760    #ifdef BOOST_MATH_STANDALONE
0761    static_assert(sizeof(x) == 0, "mpreal bessel_i0 can not be calculated in standalone mode");
0762    #endif
0763
0764     static const mpfr::mpreal P1[] = {
0765         boost::lexical_cast<mpfr::mpreal>("-2.2335582639474375249e+15"),
0766         boost::lexical_cast<mpfr::mpreal>("-5.5050369673018427753e+14"),
0767         boost::lexical_cast<mpfr::mpreal>("-3.2940087627407749166e+13"),
0768         boost::lexical_cast<mpfr::mpreal>("-8.4925101247114157499e+11"),
0769         boost::lexical_cast<mpfr::mpreal>("-1.1912746104985237192e+10"),
0770         boost::lexical_cast<mpfr::mpreal>("-1.0313066708737980747e+08"),
0771         boost::lexical_cast<mpfr::mpreal>("-5.9545626019847898221e+05"),
0772         boost::lexical_cast<mpfr::mpreal>("-2.4125195876041896775e+03"),
0773         boost::lexical_cast<mpfr::mpreal>("-7.0935347449210549190e+00"),
0774         boost::lexical_cast<mpfr::mpreal>("-1.5453977791786851041e-02"),
0775         boost::lexical_cast<mpfr::mpreal>("-2.5172644670688975051e-05"),
0776         boost::lexical_cast<mpfr::mpreal>("-3.0517226450451067446e-08"),
0777         boost::lexical_cast<mpfr::mpreal>("-2.6843448573468483278e-11"),
0778         boost::lexical_cast<mpfr::mpreal>("-1.5982226675653184646e-14"),
0779         boost::lexical_cast<mpfr::mpreal>("-5.2487866627945699800e-18"),
0780     };
0781     static const mpfr::mpreal Q1[] = {
0782         boost::lexical_cast<mpfr::mpreal>("-2.2335582639474375245e+15"),
0783         boost::lexical_cast<mpfr::mpreal>("7.8858692566751002988e+12"),
0784         boost::lexical_cast<mpfr::mpreal>("-1.2207067397808979846e+10"),
0785         boost::lexical_cast<mpfr::mpreal>("1.0377081058062166144e+07"),
0786         boost::lexical_cast<mpfr::mpreal>("-4.8527560179962773045e+03"),
0787         boost::lexical_cast<mpfr::mpreal>("1.0"),
0788     };
0789     static const mpfr::mpreal P2[] = {
0790         boost::lexical_cast<mpfr::mpreal>("-2.2210262233306573296e-04"),
0791         boost::lexical_cast<mpfr::mpreal>("1.3067392038106924055e-02"),
0792         boost::lexical_cast<mpfr::mpreal>("-4.4700805721174453923e-01"),
0793         boost::lexical_cast<mpfr::mpreal>("5.5674518371240761397e+00"),
0794         boost::lexical_cast<mpfr::mpreal>("-2.3517945679239481621e+01"),
0795         boost::lexical_cast<mpfr::mpreal>("3.1611322818701131207e+01"),
0796         boost::lexical_cast<mpfr::mpreal>("-9.6090021968656180000e+00"),
0797     };
0798     static const mpfr::mpreal Q2[] = {
0799         boost::lexical_cast<mpfr::mpreal>("-5.5194330231005480228e-04"),
0800         boost::lexical_cast<mpfr::mpreal>("3.2547697594819615062e-02"),
0801         boost::lexical_cast<mpfr::mpreal>("-1.1151759188741312645e+00"),
0802         boost::lexical_cast<mpfr::mpreal>("1.3982595353892851542e+01"),
0803         boost::lexical_cast<mpfr::mpreal>("-6.0228002066743340583e+01"),
0804         boost::lexical_cast<mpfr::mpreal>("8.5539563258012929600e+01"),
0805         boost::lexical_cast<mpfr::mpreal>("-3.1446690275135491500e+01"),
0806         boost::lexical_cast<mpfr::mpreal>("1.0"),
0807     };
0808     mpfr::mpreal value, factor, r;
0809
0810     BOOST_MATH_STD_USING
0811     using namespace boost::math::tools;
0812
0813     if (x < 0)
0814     {
0815         x = -x;                         // even function
0816     }
0817     if (x == 0)
0818     {
0819         return static_cast<mpfr::mpreal>(1);
0820     }
0821     if (x <= 15)                        // x in (0, 15]
0822     {
0823         mpfr::mpreal y = x * x;
0824         value = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);
0825     }
0826     else                                // x in (15, \infty)
0827     {
0828         mpfr::mpreal y = 1 / x - mpfr::mpreal(1) / 15;
0829         r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);
0830         factor = exp(x) / sqrt(x);
0831         value = factor * r;
0832     }
0833
0834     return value;
0835 }
0836
0837 inline mpfr::mpreal bessel_i1(mpfr::mpreal x)
0838 {
0839     static const mpfr::mpreal P1[] = {
0840         static_cast<mpfr::mpreal>("-1.4577180278143463643e+15"),
0841         static_cast<mpfr::mpreal>("-1.7732037840791591320e+14"),
0842         static_cast<mpfr::mpreal>("-6.9876779648010090070e+12"),
0843         static_cast<mpfr::mpreal>("-1.3357437682275493024e+11"),
0844         static_cast<mpfr::mpreal>("-1.4828267606612366099e+09"),
0845         static_cast<mpfr::mpreal>("-1.0588550724769347106e+07"),
0846         static_cast<mpfr::mpreal>("-5.1894091982308017540e+04"),
0847         static_cast<mpfr::mpreal>("-1.8225946631657315931e+02"),
0848         static_cast<mpfr::mpreal>("-4.7207090827310162436e-01"),
0849         static_cast<mpfr::mpreal>("-9.1746443287817501309e-04"),
0850         static_cast<mpfr::mpreal>("-1.3466829827635152875e-06"),
0851         static_cast<mpfr::mpreal>("-1.4831904935994647675e-09"),
0852         static_cast<mpfr::mpreal>("-1.1928788903603238754e-12"),
0853         static_cast<mpfr::mpreal>("-6.5245515583151902910e-16"),
0854         static_cast<mpfr::mpreal>("-1.9705291802535139930e-19"),
0855     };
0856     static const mpfr::mpreal Q1[] = {
0857         static_cast<mpfr::mpreal>("-2.9154360556286927285e+15"),
0858         static_cast<mpfr::mpreal>("9.7887501377547640438e+12"),
0859         static_cast<mpfr::mpreal>("-1.4386907088588283434e+10"),
0860         static_cast<mpfr::mpreal>("1.1594225856856884006e+07"),
0861         static_cast<mpfr::mpreal>("-5.1326864679904189920e+03"),
0862         static_cast<mpfr::mpreal>("1.0"),
0863     };
0864     static const mpfr::mpreal P2[] = {
0865         static_cast<mpfr::mpreal>("1.4582087408985668208e-05"),
0866         static_cast<mpfr::mpreal>("-8.9359825138577646443e-04"),
0867         static_cast<mpfr::mpreal>("2.9204895411257790122e-02"),
0868         static_cast<mpfr::mpreal>("-3.4198728018058047439e-01"),
0869         static_cast<mpfr::mpreal>("1.3960118277609544334e+00"),
0870         static_cast<mpfr::mpreal>("-1.9746376087200685843e+00"),
0871         static_cast<mpfr::mpreal>("8.5591872901933459000e-01"),
0872         static_cast<mpfr::mpreal>("-6.0437159056137599999e-02"),
0873     };
0874     static const mpfr::mpreal Q2[] = {
0875         static_cast<mpfr::mpreal>("3.7510433111922824643e-05"),
0876         static_cast<mpfr::mpreal>("-2.2835624489492512649e-03"),
0877         static_cast<mpfr::mpreal>("7.4212010813186530069e-02"),
0878         static_cast<mpfr::mpreal>("-8.5017476463217924408e-01"),
0879         static_cast<mpfr::mpreal>("3.2593714889036996297e+00"),
0880         static_cast<mpfr::mpreal>("-3.8806586721556593450e+00"),
0881         static_cast<mpfr::mpreal>("1.0"),
0882     };
0883     mpfr::mpreal value, factor, r, w;
0884
0885     BOOST_MATH_STD_USING
0886     using namespace boost::math::tools;
0887
0888     w = abs(x);
0889     if (x == 0)
0890     {
0891         return static_cast<mpfr::mpreal>(0);
0892     }
0893     if (w <= 15)                        // w in (0, 15]
0894     {
0895         mpfr::mpreal y = x * x;
0896         r = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);
0897         factor = w;
0898         value = factor * r;
0899     }
0900     else                                // w in (15, \infty)
0901     {
0902         mpfr::mpreal y = 1 / w - mpfr::mpreal(1) / 15;
0903         r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);
0904         factor = exp(w) / sqrt(w);
0905         value = factor * r;
0906     }
0907
0908     if (x < 0)
0909     {
0910         value *= -value;                 // odd function
0911     }
0912     return value;
0913 }
0914
0915 } // namespace detail
0916 } // namespace math
0917
0918 }
0919
0920 #endif // BOOST_MATH_MPLFR_BINDINGS_HPP
0921