ceph/src/boost/libs/math/include/boost/math/bindings/mpfr.hpp

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