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