math/constants/calculate_constants.hpp

0001 //  Copyright John Maddock 2010, 2012.
0002 //  Copyright Paul A. Bristow 2011, 2012.
0003
0004 //  Use, modification and distribution are subject to the
0005 //  Boost Software License, Version 1.0. (See accompanying file
0006 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0007
0008 #ifndef BOOST_MATH_CALCULATE_CONSTANTS_CONSTANTS_INCLUDED
0009 #define BOOST_MATH_CALCULATE_CONSTANTS_CONSTANTS_INCLUDED
0010 #include <type_traits>
0011
0012 namespace boost{ namespace math{ namespace constants{ namespace detail{
0013
0014 template <class T>
0015 template<int N>
0016 inline T constant_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0017 {
0018    BOOST_MATH_STD_USING
0019
0020    return ldexp(acos(T(0)), 1);
0021
0022    /*
0023    // Although this code works well, it's usually more accurate to just call acos
0024    // and access the number types own representation of PI which is usually calculated
0025    // at slightly higher precision...
0026
0027    T result;
0028    T a = 1;
0029    T b;
0030    T A(a);
0031    T B = 0.5f;
0032    T D = 0.25f;
0033
0034    T lim;
0035    lim = boost::math::tools::epsilon<T>();
0036
0037    unsigned k = 1;
0038
0039    do
0040    {
0041       result = A + B;
0042       result = ldexp(result, -2);
0043       b = sqrt(B);
0044       a += b;
0045       a = ldexp(a, -1);
0046       A = a * a;
0047       B = A - result;
0048       B = ldexp(B, 1);
0049       result = A - B;
0050       bool neg = boost::math::sign(result) < 0;
0051       if(neg)
0052          result = -result;
0053       if(result <= lim)
0054          break;
0055       if(neg)
0056          result = -result;
0057       result = ldexp(result, k - 1);
0058       D -= result;
0059       ++k;
0060       lim = ldexp(lim, 1);
0061    }
0062    while(true);
0063
0064    result = B / D;
0065    return result;
0066    */
0067 }
0068
0069 template <class T>
0070 template<int N>
0071 inline T constant_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0072 {
0073    return 2 * pi<T, policies::policy<policies::digits2<N> > >();
0074 }
0075
0076 template <class T> // 2 / pi
0077 template<int N>
0078 inline T constant_two_div_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0079 {
0080    return 2 / pi<T, policies::policy<policies::digits2<N> > >();
0081 }
0082
0083 template <class T>  // sqrt(2/pi)
0084 template <int N>
0085 inline T constant_root_two_div_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0086 {
0087    BOOST_MATH_STD_USING
0088    return sqrt((2 / pi<T, policies::policy<policies::digits2<N> > >()));
0089 }
0090
0091 template <class T>
0092 template<int N>
0093 inline T constant_one_div_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0094 {
0095    return 1 / two_pi<T, policies::policy<policies::digits2<N> > >();
0096 }
0097
0098 template <class T>
0099 template<int N>
0100 inline T constant_root_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0101 {
0102    BOOST_MATH_STD_USING
0103    return sqrt(pi<T, policies::policy<policies::digits2<N> > >());
0104 }
0105
0106 template <class T>
0107 template<int N>
0108 inline T constant_root_half_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0109 {
0110    BOOST_MATH_STD_USING
0111    return sqrt(pi<T, policies::policy<policies::digits2<N> > >() / 2);
0112 }
0113
0114 template <class T>
0115 template<int N>
0116 inline T constant_root_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0117 {
0118    BOOST_MATH_STD_USING
0119    return sqrt(two_pi<T, policies::policy<policies::digits2<N> > >());
0120 }
0121
0122 template <class T>
0123 template<int N>
0124 inline T constant_log_root_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0125 {
0126    BOOST_MATH_STD_USING
0127    return log(root_two_pi<T, policies::policy<policies::digits2<N> > >());
0128 }
0129
0130 template <class T>
0131 template<int N>
0132 inline T constant_root_ln_four<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0133 {
0134    BOOST_MATH_STD_USING
0135    return sqrt(log(static_cast<T>(4)));
0136 }
0137
0138 template <class T>
0139 template<int N>
0140 inline T constant_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0141 {
0142    //
0143    // Although we can clearly calculate this from first principles, this hooks into
0144    // T's own notion of e, which hopefully will more accurate than one calculated to
0145    // a few epsilon:
0146    //
0147    BOOST_MATH_STD_USING
0148    return exp(static_cast<T>(1));
0149 }
0150
0151 template <class T>
0152 template<int N>
0153 inline T constant_half<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0154 {
0155    return static_cast<T>(1) / static_cast<T>(2);
0156 }
0157
0158 template <class T>
0159 template<int M>
0160 inline T constant_euler<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, M>)))
0161 {
0162    BOOST_MATH_STD_USING
0163    //
0164    // This is the method described in:
0165    // "Some New Algorithms for High-Precision Computation of Euler's Constant"
0166    // Richard P Brent and Edwin M McMillan.
0167    // Mathematics of Computation, Volume 34, Number 149, Jan 1980, pages 305-312.
0168    // See equation 17 with p = 2.
0169    //
0170    T n = 3 + (M ? (std::min)(M, tools::digits<T>()) : tools::digits<T>()) / 4;
0171    T lim = M ? ldexp(T(1), 1 - (std::min)(M, tools::digits<T>())) : tools::epsilon<T>();
0172    T lnn = log(n);
0173    T term = 1;
0174    T N = -lnn;
0175    T D = 1;
0176    T Hk = 0;
0177    T one = 1;
0178
0179    for(unsigned k = 1;; ++k)
0180    {
0181       term *= n * n;
0182       term /= k * k;
0183       Hk += one / k;
0184       N += term * (Hk - lnn);
0185       D += term;
0186
0187       if(term < D * lim)
0188          break;
0189    }
0190    return N / D;
0191 }
0192
0193 template <class T>
0194 template<int N>
0195 inline T constant_euler_sqr<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0196 {
0197   BOOST_MATH_STD_USING
0198   return euler<T, policies::policy<policies::digits2<N> > >()
0199      * euler<T, policies::policy<policies::digits2<N> > >();
0200 }
0201
0202 template <class T>
0203 template<int N>
0204 inline T constant_one_div_euler<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0205 {
0206   BOOST_MATH_STD_USING
0207   return static_cast<T>(1)
0208      / euler<T, policies::policy<policies::digits2<N> > >();
0209 }
0210
0211
0212 template <class T>
0213 template<int N>
0214 inline T constant_root_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0215 {
0216    BOOST_MATH_STD_USING
0217    return sqrt(static_cast<T>(2));
0218 }
0219
0220
0221 template <class T>
0222 template<int N>
0223 inline T constant_root_three<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0224 {
0225    BOOST_MATH_STD_USING
0226    return sqrt(static_cast<T>(3));
0227 }
0228
0229 template <class T>
0230 template<int N>
0231 inline T constant_half_root_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0232 {
0233    BOOST_MATH_STD_USING
0234    return sqrt(static_cast<T>(2)) / 2;
0235 }
0236
0237 template <class T>
0238 template<int N>
0239 inline T constant_ln_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0240 {
0241    //
0242    // Although there are good ways to calculate this from scratch, this hooks into
0243    // T's own notion of log(2) which will hopefully be accurate to the full precision
0244    // of T:
0245    //
0246    BOOST_MATH_STD_USING
0247    return log(static_cast<T>(2));
0248 }
0249
0250 template <class T>
0251 template<int N>
0252 inline T constant_ln_ten<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0253 {
0254    BOOST_MATH_STD_USING
0255    return log(static_cast<T>(10));
0256 }
0257
0258 template <class T>
0259 template<int N>
0260 inline T constant_ln_ln_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0261 {
0262    BOOST_MATH_STD_USING
0263    return log(log(static_cast<T>(2)));
0264 }
0265
0266 template <class T>
0267 template<int N>
0268 inline T constant_third<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0269 {
0270    BOOST_MATH_STD_USING
0271    return static_cast<T>(1) / static_cast<T>(3);
0272 }
0273
0274 template <class T>
0275 template<int N>
0276 inline T constant_twothirds<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0277 {
0278    BOOST_MATH_STD_USING
0279    return static_cast<T>(2) / static_cast<T>(3);
0280 }
0281
0282 template <class T>
0283 template<int N>
0284 inline T constant_two_thirds<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0285 {
0286    BOOST_MATH_STD_USING
0287    return static_cast<T>(2) / static_cast<T>(3);
0288 }
0289
0290 template <class T>
0291 template<int N>
0292 inline T constant_three_quarters<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0293 {
0294    BOOST_MATH_STD_USING
0295    return static_cast<T>(3) / static_cast<T>(4);
0296 }
0297
0298 template <class T>
0299 template<int N>
0300 inline T constant_sixth<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0301 {
0302   BOOST_MATH_STD_USING
0303   return static_cast<T>(1) / static_cast<T>(6);
0304 }
0305
0306 // Pi and related constants.
0307 template <class T>
0308 template<int N>
0309 inline T constant_pi_minus_three<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0310 {
0311    return pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(3);
0312 }
0313
0314 template <class T>
0315 template<int N>
0316 inline T constant_four_minus_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0317 {
0318    return static_cast<T>(4) - pi<T, policies::policy<policies::digits2<N> > >();
0319 }
0320
0321 template <class T>
0322 template<int N>
0323 inline T constant_exp_minus_half<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0324 {
0325    BOOST_MATH_STD_USING
0326    return exp(static_cast<T>(-0.5));
0327 }
0328
0329 template <class T>
0330 template<int N>
0331 inline T constant_exp_minus_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0332 {
0333   BOOST_MATH_STD_USING
0334   return exp(static_cast<T>(-1.));
0335 }
0336
0337 template <class T>
0338 template<int N>
0339 inline T constant_one_div_root_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0340 {
0341    return static_cast<T>(1) / root_two<T, policies::policy<policies::digits2<N> > >();
0342 }
0343
0344 template <class T>
0345 template<int N>
0346 inline T constant_one_div_root_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0347 {
0348    return static_cast<T>(1) / root_pi<T, policies::policy<policies::digits2<N> > >();
0349 }
0350
0351 template <class T>
0352 template<int N>
0353 inline T constant_one_div_root_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0354 {
0355    return static_cast<T>(1) / root_two_pi<T, policies::policy<policies::digits2<N> > >();
0356 }
0357
0358 template <class T>
0359 template<int N>
0360 inline T constant_root_one_div_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0361 {
0362    BOOST_MATH_STD_USING
0363    return sqrt(static_cast<T>(1) / pi<T, policies::policy<policies::digits2<N> > >());
0364 }
0365
0366 template <class T>
0367 template<int N>
0368 inline T constant_four_thirds_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0369 {
0370    BOOST_MATH_STD_USING
0371    return pi<T, policies::policy<policies::digits2<N> > >() * static_cast<T>(4) / static_cast<T>(3);
0372 }
0373
0374 template <class T>
0375 template<int N>
0376 inline T constant_half_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0377 {
0378    BOOST_MATH_STD_USING
0379    return pi<T, policies::policy<policies::digits2<N> > >()  / static_cast<T>(2);
0380 }
0381
0382 template <class T>
0383 template<int N>
0384 inline T constant_third_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0385 {
0386    BOOST_MATH_STD_USING
0387    return pi<T, policies::policy<policies::digits2<N> > >()  / static_cast<T>(3);
0388 }
0389
0390 template <class T>
0391 template<int N>
0392 inline T constant_sixth_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0393 {
0394    BOOST_MATH_STD_USING
0395    return pi<T, policies::policy<policies::digits2<N> > >()  / static_cast<T>(6);
0396 }
0397
0398 template <class T>
0399 template<int N>
0400 inline T constant_two_thirds_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0401 {
0402    BOOST_MATH_STD_USING
0403    return pi<T, policies::policy<policies::digits2<N> > >() * static_cast<T>(2) / static_cast<T>(3);
0404 }
0405
0406 template <class T>
0407 template<int N>
0408 inline T constant_three_quarters_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0409 {
0410    BOOST_MATH_STD_USING
0411    return pi<T, policies::policy<policies::digits2<N> > >() * static_cast<T>(3) / static_cast<T>(4);
0412 }
0413
0414 template <class T>
0415 template<int N>
0416 inline T constant_pi_pow_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0417 {
0418    BOOST_MATH_STD_USING
0419    return pow(pi<T, policies::policy<policies::digits2<N> > >(), e<T, policies::policy<policies::digits2<N> > >()); //
0420 }
0421
0422 template <class T>
0423 template<int N>
0424 inline T constant_pi_sqr<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0425 {
0426    BOOST_MATH_STD_USING
0427    return pi<T, policies::policy<policies::digits2<N> > >()
0428    *      pi<T, policies::policy<policies::digits2<N> > >() ; //
0429 }
0430
0431 template <class T>
0432 template<int N>
0433 inline T constant_pi_sqr_div_six<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0434 {
0435    BOOST_MATH_STD_USING
0436    return pi<T, policies::policy<policies::digits2<N> > >()
0437    *      pi<T, policies::policy<policies::digits2<N> > >()
0438    / static_cast<T>(6); //
0439 }
0440
0441 template <class T>
0442 template<int N>
0443 inline T constant_pi_cubed<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0444 {
0445    BOOST_MATH_STD_USING
0446    return pi<T, policies::policy<policies::digits2<N> > >()
0447    *      pi<T, policies::policy<policies::digits2<N> > >()
0448    *      pi<T, policies::policy<policies::digits2<N> > >()
0449    ; //
0450 }
0451
0452 template <class T>
0453 template<int N>
0454 inline T constant_cbrt_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0455 {
0456    BOOST_MATH_STD_USING
0457    return pow(pi<T, policies::policy<policies::digits2<N> > >(), static_cast<T>(1)/ static_cast<T>(3));
0458 }
0459
0460 template <class T>
0461 template<int N>
0462 inline T constant_one_div_cbrt_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0463 {
0464    BOOST_MATH_STD_USING
0465    return static_cast<T>(1)
0466    / pow(pi<T, policies::policy<policies::digits2<N> > >(), static_cast<T>(1)/ static_cast<T>(3));
0467 }
0468
0469 // Euler's e
0470
0471 template <class T>
0472 template<int N>
0473 inline T constant_e_pow_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0474 {
0475    BOOST_MATH_STD_USING
0476    return pow(e<T, policies::policy<policies::digits2<N> > >(), pi<T, policies::policy<policies::digits2<N> > >()); //
0477 }
0478
0479 template <class T>
0480 template<int N>
0481 inline T constant_root_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0482 {
0483    BOOST_MATH_STD_USING
0484    return sqrt(e<T, policies::policy<policies::digits2<N> > >());
0485 }
0486
0487 template <class T>
0488 template<int N>
0489 inline T constant_log10_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0490 {
0491    BOOST_MATH_STD_USING
0492    return log10(e<T, policies::policy<policies::digits2<N> > >());
0493 }
0494
0495 template <class T>
0496 template<int N>
0497 inline T constant_one_div_log10_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0498 {
0499    BOOST_MATH_STD_USING
0500    return  static_cast<T>(1) /
0501      log10(e<T, policies::policy<policies::digits2<N> > >());
0502 }
0503
0504 // Trigonometric
0505
0506 template <class T>
0507 template<int N>
0508 inline T constant_degree<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0509 {
0510    BOOST_MATH_STD_USING
0511    return pi<T, policies::policy<policies::digits2<N> > >()
0512    / static_cast<T>(180)
0513    ; //
0514 }
0515
0516 template <class T>
0517 template<int N>
0518 inline T constant_radian<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0519 {
0520    BOOST_MATH_STD_USING
0521    return static_cast<T>(180)
0522    / pi<T, policies::policy<policies::digits2<N> > >()
0523    ; //
0524 }
0525
0526 template <class T>
0527 template<int N>
0528 inline T constant_sin_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0529 {
0530    BOOST_MATH_STD_USING
0531    return sin(static_cast<T>(1)) ; //
0532 }
0533
0534 template <class T>
0535 template<int N>
0536 inline T constant_cos_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0537 {
0538    BOOST_MATH_STD_USING
0539    return cos(static_cast<T>(1)) ; //
0540 }
0541
0542 template <class T>
0543 template<int N>
0544 inline T constant_sinh_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0545 {
0546    BOOST_MATH_STD_USING
0547    return sinh(static_cast<T>(1)) ; //
0548 }
0549
0550 template <class T>
0551 template<int N>
0552 inline T constant_cosh_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0553 {
0554    BOOST_MATH_STD_USING
0555    return cosh(static_cast<T>(1)) ; //
0556 }
0557
0558 template <class T>
0559 template<int N>
0560 inline T constant_phi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0561 {
0562    BOOST_MATH_STD_USING
0563    return (static_cast<T>(1) + sqrt(static_cast<T>(5)) )/static_cast<T>(2) ; //
0564 }
0565
0566 template <class T>
0567 template<int N>
0568 inline T constant_ln_phi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0569 {
0570    BOOST_MATH_STD_USING
0571    return log((static_cast<T>(1) + sqrt(static_cast<T>(5)) )/static_cast<T>(2) );
0572 }
0573
0574 template <class T>
0575 template<int N>
0576 inline T constant_one_div_ln_phi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0577 {
0578    BOOST_MATH_STD_USING
0579    return static_cast<T>(1) /
0580      log((static_cast<T>(1) + sqrt(static_cast<T>(5)) )/static_cast<T>(2) );
0581 }
0582
0583 // Zeta
0584
0585 template <class T>
0586 template<int N>
0587 inline T constant_zeta_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0588 {
0589    BOOST_MATH_STD_USING
0590
0591      return pi<T, policies::policy<policies::digits2<N> > >()
0592      *  pi<T, policies::policy<policies::digits2<N> > >()
0593      /static_cast<T>(6);
0594 }
0595
0596 template <class T>
0597 template<int N>
0598 inline T constant_zeta_three<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0599 {
0600    // http://mathworld.wolfram.com/AperysConstant.html
0601    // http://en.wikipedia.org/wiki/Mathematical_constant
0602
0603    // http://oeis.org/A002117/constant
0604    //T zeta3("1.20205690315959428539973816151144999076"
0605    // "4986292340498881792271555341838205786313"
0606    // "09018645587360933525814619915");
0607
0608   //"1.202056903159594285399738161511449990, 76498629234049888179227155534183820578631309018645587360933525814619915"  A002117
0609   // 1.202056903159594285399738161511449990, 76498629234049888179227155534183820578631309018645587360933525814619915780, +00);
0610   //"1.2020569031595942 double
0611   // http://www.spaennare.se/SSPROG/ssnum.pdf  // section 11, Algorithm for Apery's constant zeta(3).
0612   // Programs to Calculate some Mathematical Constants to Large Precision, Document Version 1.50
0613
0614   // by Stefan Spannare  September 19, 2007
0615   // zeta(3) = 1/64 * sum
0616    BOOST_MATH_STD_USING
0617    T n_fact=static_cast<T>(1); // build n! for n = 0.
0618    T sum = static_cast<double>(77); // Start with n = 0 case.
0619    // for n = 0, (77/1) /64 = 1.203125
0620    //double lim = std::numeric_limits<double>::epsilon();
0621    T lim = N ? ldexp(T(1), 1 - (std::min)(N, tools::digits<T>())) : tools::epsilon<T>();
0622    for(unsigned int n = 1; n < 40; ++n)
0623    { // three to five decimal digits per term, so 40 should be plenty for 100 decimal digits.
0624       //cout << "n = " << n << endl;
0625       n_fact *= n; // n!
0626       T n_fact_p10 = n_fact * n_fact * n_fact * n_fact * n_fact * n_fact * n_fact * n_fact * n_fact * n_fact; // (n!)^10
0627       T num = ((205 * n * n) + (250 * n) + 77) * n_fact_p10; // 205n^2 + 250n + 77
0628       // int nn = (2 * n + 1);
0629       // T d = factorial(nn); // inline factorial.
0630       T d = 1;
0631       for(unsigned int i = 1; i <= (n+n + 1); ++i) // (2n + 1)
0632       {
0633         d *= i;
0634       }
0635       T den = d * d * d * d * d; // [(2n+1)!]^5
0636       //cout << "den = " << den << endl;
0637       T term = num/den;
0638       if (n % 2 != 0)
0639       { //term *= -1;
0640         sum -= term;
0641       }
0642       else
0643       {
0644         sum += term;
0645       }
0646       //cout << "term = " << term << endl;
0647       //cout << "sum/64  = " << sum/64 << endl;
0648       if(abs(term) < lim)
0649       {
0650          break;
0651       }
0652    }
0653    return sum / 64;
0654 }
0655
0656 template <class T>
0657 template<int N>
0658 inline T constant_catalan<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0659 { // http://oeis.org/A006752/constant
0660   //T c("0.915965594177219015054603514932384110774"
0661  //"149374281672134266498119621763019776254769479356512926115106248574");
0662
0663   // 9.159655941772190150546035149323841107, 74149374281672134266498119621763019776254769479356512926115106248574422619, -01);
0664
0665    // This is equation (entry) 31 from
0666    // http://www-2.cs.cmu.edu/~adamchik/articles/catalan/catalan.htm
0667    // See also http://www.mpfr.org/algorithms.pdf
0668    BOOST_MATH_STD_USING
0669    T k_fact = 1;
0670    T tk_fact = 1;
0671    T sum = 1;
0672    T term;
0673    T lim = N ? ldexp(T(1), 1 - (std::min)(N, tools::digits<T>())) : tools::epsilon<T>();
0674
0675    for(unsigned k = 1;; ++k)
0676    {
0677       k_fact *= k;
0678       tk_fact *= (2 * k) * (2 * k - 1);
0679       term = k_fact * k_fact / (tk_fact * (2 * k + 1) * (2 * k + 1));
0680       sum += term;
0681       if(term < lim)
0682       {
0683          break;
0684       }
0685    }
0686    return boost::math::constants::pi<T, boost::math::policies::policy<> >()
0687       * log(2 + boost::math::constants::root_three<T, boost::math::policies::policy<> >())
0688        / 8
0689       + 3 * sum / 8;
0690 }
0691
0692 namespace khinchin_detail{
0693
0694 template <class T>
0695 T zeta_polynomial_series(T s, T sc, int digits)
0696 {
0697    BOOST_MATH_STD_USING
0698    //
0699    // This is algorithm 3 from:
0700    //
0701    // "An Efficient Algorithm for the Riemann Zeta Function", P. Borwein,
0702    // Canadian Mathematical Society, Conference Proceedings, 2000.
0703    // See: http://www.cecm.sfu.ca/personal/pborwein/PAPERS/P155.pdf
0704    //
0705    BOOST_MATH_STD_USING
0706    int n = (digits * 19) / 53;
0707    T sum = 0;
0708    T two_n = ldexp(T(1), n);
0709    int ej_sign = 1;
0710    for(int j = 0; j < n; ++j)
0711    {
0712       sum += ej_sign * -two_n / pow(T(j + 1), s);
0713       ej_sign = -ej_sign;
0714    }
0715    T ej_sum = 1;
0716    T ej_term = 1;
0717    for(int j = n; j <= 2 * n - 1; ++j)
0718    {
0719       sum += ej_sign * (ej_sum - two_n) / pow(T(j + 1), s);
0720       ej_sign = -ej_sign;
0721       ej_term *= 2 * n - j;
0722       ej_term /= j - n + 1;
0723       ej_sum += ej_term;
0724    }
0725    return -sum / (two_n * (1 - pow(T(2), sc)));
0726 }
0727
0728 template <class T>
0729 T khinchin(int digits)
0730 {
0731    BOOST_MATH_STD_USING
0732    T sum = 0;
0733    T term;
0734    T lim = ldexp(T(1), 1-digits);
0735    T factor = 0;
0736    unsigned last_k = 1;
0737    T num = 1;
0738    for(unsigned n = 1;; ++n)
0739    {
0740       for(unsigned k = last_k; k <= 2 * n - 1; ++k)
0741       {
0742          factor += num / k;
0743          num = -num;
0744       }
0745       last_k = 2 * n;
0746       term = (zeta_polynomial_series(T(2 * n), T(1 - T(2 * n)), digits) - 1) * factor / n;
0747       sum += term;
0748       if(term < lim)
0749          break;
0750    }
0751    return exp(sum / boost::math::constants::ln_two<T, boost::math::policies::policy<> >());
0752 }
0753
0754 }
0755
0756 template <class T>
0757 template<int N>
0758 inline T constant_khinchin<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0759 {
0760    int n = N ? (std::min)(N, tools::digits<T>()) : tools::digits<T>();
0761    return khinchin_detail::khinchin<T>(n);
0762 }
0763
0764 template <class T>
0765 template<int N>
0766 inline T constant_extreme_value_skewness<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0767 {  // N[12 Sqrt[6]  Zeta[3]/Pi^3, 1101]
0768    BOOST_MATH_STD_USING
0769    T ev(12 * sqrt(static_cast<T>(6)) * zeta_three<T, policies::policy<policies::digits2<N> > >()
0770     / pi_cubed<T, policies::policy<policies::digits2<N> > >() );
0771
0772 //T ev(
0773 //"1.1395470994046486574927930193898461120875997958365518247216557100852480077060706857071875468869385150"
0774 //"1894272048688553376986765366075828644841024041679714157616857834895702411080704529137366329462558680"
0775 //"2015498788776135705587959418756809080074611906006528647805347822929577145038743873949415294942796280"
0776 //"0895597703063466053535550338267721294164578901640163603544404938283861127819804918174973533694090594"
0777 //"3094963822672055237678432023017824416203652657301470473548274848068762500300316769691474974950757965"
0778 //"8640779777748741897542093874605477776538884083378029488863880220988107155275203245233994097178778984"
0779 //"3488995668362387892097897322246698071290011857605809901090220903955815127463328974447572119951192970"
0780 //"3684453635456559086126406960279692862247058250100678008419431185138019869693206366891639436908462809"
0781 //"9756051372711251054914491837034685476095423926553367264355374652153595857163724698198860485357368964"
0782 //"3807049634423621246870868566707915720704996296083373077647528285782964567312903914752617978405994377"
0783 //"9064157147206717895272199736902453130842229559980076472936976287378945035706933650987259357729800315");
0784
0785   return ev;
0786 }
0787
0788 namespace detail{
0789 //
0790 // Calculation of the Glaisher constant depends upon calculating the
0791 // derivative of the zeta function at 2, we can then use the relation:
0792 // zeta'(2) = 1/6 pi^2 [euler + ln(2pi)-12ln(A)]
0793 // To get the constant A.
0794 // See equation 45 at http://mathworld.wolfram.com/RiemannZetaFunction.html.
0795 //
0796 // The derivative of the zeta function is computed by direct differentiation
0797 // of the relation:
0798 // (1-2^(1-s))zeta(s) = SUM(n=0, INF){ (-n)^n / (n+1)^s  }
0799 // Which gives us 2 slowly converging but alternating sums to compute,
0800 // for this we use Algorithm 1 from "Convergent Acceleration of Alternating Series",
0801 // Henri Cohen, Fernando Rodriguez Villegas and Don Zagier, Experimental Mathematics 9:1 (1999).
0802 // See http://www.math.utexas.edu/users/villegas/publications/conv-accel.pdf
0803 //
0804 template <class T>
0805 T zeta_series_derivative_2(unsigned digits)
0806 {
0807    // Derivative of the series part, evaluated at 2:
0808    BOOST_MATH_STD_USING
0809    int n = digits * 301 * 13 / 10000;
0810    T d = pow(3 + sqrt(T(8)), n);
0811    d = (d + 1 / d) / 2;
0812    T b = -1;
0813    T c = -d;
0814    T s = 0;
0815    for(int k = 0; k < n; ++k)
0816    {
0817       T a = -log(T(k+1)) / ((k+1) * (k+1));
0818       c = b - c;
0819       s = s + c * a;
0820       b = (k + n) * (k - n) * b / ((k + T(0.5f)) * (k + 1));
0821    }
0822    return s / d;
0823 }
0824
0825 template <class T>
0826 T zeta_series_2(unsigned digits)
0827 {
0828    // Series part of zeta at 2:
0829    BOOST_MATH_STD_USING
0830    int n = digits * 301 * 13 / 10000;
0831    T d = pow(3 + sqrt(T(8)), n);
0832    d = (d + 1 / d) / 2;
0833    T b = -1;
0834    T c = -d;
0835    T s = 0;
0836    for(int k = 0; k < n; ++k)
0837    {
0838       T a = T(1) / ((k + 1) * (k + 1));
0839       c = b - c;
0840       s = s + c * a;
0841       b = (k + n) * (k - n) * b / ((k + T(0.5f)) * (k + 1));
0842    }
0843    return s / d;
0844 }
0845
0846 template <class T>
0847 inline T zeta_series_lead_2()
0848 {
0849    // lead part at 2:
0850    return 2;
0851 }
0852
0853 template <class T>
0854 inline T zeta_series_derivative_lead_2()
0855 {
0856    // derivative of lead part at 2:
0857    return -2 * boost::math::constants::ln_two<T>();
0858 }
0859
0860 template <class T>
0861 inline T zeta_derivative_2(unsigned n)
0862 {
0863    // zeta derivative at 2:
0864    return zeta_series_derivative_2<T>(n) * zeta_series_lead_2<T>()
0865       + zeta_series_derivative_lead_2<T>() * zeta_series_2<T>(n);
0866 }
0867
0868 }  // namespace detail
0869
0870 template <class T>
0871 template<int N>
0872 inline T constant_glaisher<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0873 {
0874
0875    BOOST_MATH_STD_USING
0876    typedef policies::policy<policies::digits2<N> > forwarding_policy;
0877    int n = N ? (std::min)(N, tools::digits<T>()) : tools::digits<T>();
0878    T v = detail::zeta_derivative_2<T>(n);
0879    v *= 6;
0880    v /= boost::math::constants::pi<T, forwarding_policy>() * boost::math::constants::pi<T, forwarding_policy>();
0881    v -= boost::math::constants::euler<T, forwarding_policy>();
0882    v -= log(2 * boost::math::constants::pi<T, forwarding_policy>());
0883    v /= -12;
0884    return exp(v);
0885
0886    /*
0887    // from http://mpmath.googlecode.com/svn/data/glaisher.txt
0888      // 20,000 digits of the Glaisher-Kinkelin constant A = exp(1/2 - zeta'(-1))
0889      // Computed using A = exp((6 (-zeta'(2))/pi^2 + log 2 pi + gamma)/12)
0890    // with Euler-Maclaurin summation for zeta'(2).
0891    T g(
0892    "1.282427129100622636875342568869791727767688927325001192063740021740406308858826"
0893    "46112973649195820237439420646120399000748933157791362775280404159072573861727522"
0894    "14334327143439787335067915257366856907876561146686449997784962754518174312394652"
0895    "76128213808180219264516851546143919901083573730703504903888123418813674978133050"
0896    "93770833682222494115874837348064399978830070125567001286994157705432053927585405"
0897    "81731588155481762970384743250467775147374600031616023046613296342991558095879293"
0898    "36343887288701988953460725233184702489001091776941712153569193674967261270398013"
0899    "52652668868978218897401729375840750167472114895288815996668743164513890306962645"
0900    "59870469543740253099606800842447417554061490189444139386196089129682173528798629"
0901    "88434220366989900606980888785849587494085307347117090132667567503310523405221054"
0902    "14176776156308191919997185237047761312315374135304725819814797451761027540834943"
0903    "14384965234139453373065832325673954957601692256427736926358821692159870775858274"
0904    "69575162841550648585890834128227556209547002918593263079373376942077522290940187");
0905
0906    return g;
0907    */
0908 }
0909
0910 template <class T>
0911 template<int N>
0912 inline T constant_rayleigh_skewness<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0913 {  // 1100 digits of the Rayleigh distribution skewness
0914    // N[2 Sqrt[Pi] (Pi - 3)/((4 - Pi)^(3/2)), 1100]
0915
0916    BOOST_MATH_STD_USING
0917    T rs(2 * root_pi<T, policies::policy<policies::digits2<N> > >()
0918       * pi_minus_three<T, policies::policy<policies::digits2<N> > >()
0919       / pow(four_minus_pi<T, policies::policy<policies::digits2<N> > >(), static_cast<T>(3./2))
0920       );
0921    //   6.31110657818937138191899351544227779844042203134719497658094585692926819617473725459905027032537306794400047264,
0922
0923    //"0.6311106578189371381918993515442277798440422031347194976580945856929268196174737254599050270325373067"
0924    //"9440004726436754739597525250317640394102954301685809920213808351450851396781817932734836994829371322"
0925    //"5797376021347531983451654130317032832308462278373358624120822253764532674177325950686466133508511968"
0926    //"2389168716630349407238090652663422922072397393006683401992961569208109477307776249225072042971818671"
0927    //"4058887072693437217879039875871765635655476241624825389439481561152126886932506682176611183750503553"
0928    //"1218982627032068396407180216351425758181396562859085306247387212297187006230007438534686340210168288"
0929    //"8956816965453815849613622117088096547521391672977226658826566757207615552041767516828171274858145957"
0930    //"6137539156656005855905288420585194082284972984285863898582313048515484073396332610565441264220790791"
0931    //"0194897267890422924599776483890102027823328602965235306539844007677157873140562950510028206251529523"
0932    //"7428049693650605954398446899724157486062545281504433364675815915402937209673727753199567661561209251"
0933    //"4695589950526053470201635372590001578503476490223746511106018091907936826431407434894024396366284848");  ;
0934    return rs;
0935 }
0936
0937 template <class T>
0938 template<int N>
0939 inline T constant_rayleigh_kurtosis_excess<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0940 { // - (6 Pi^2 - 24 Pi + 16)/((Pi - 4)^2)
0941     // Might provide and calculate this using pi_minus_four.
0942    BOOST_MATH_STD_USING
0943    return - (((static_cast<T>(6) * pi<T, policies::policy<policies::digits2<N> > >()
0944         * pi<T, policies::policy<policies::digits2<N> > >())
0945    - (static_cast<T>(24) * pi<T, policies::policy<policies::digits2<N> > >()) + static_cast<T>(16) )
0946    /
0947    ((pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(4))
0948    * (pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(4)))
0949    );
0950 }
0951
0952 template <class T>
0953 template<int N>
0954 inline T constant_rayleigh_kurtosis<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0955 { // 3 - (6 Pi^2 - 24 Pi + 16)/((Pi - 4)^2)
0956   // Might provide and calculate this using pi_minus_four.
0957    BOOST_MATH_STD_USING
0958    return static_cast<T>(3) - (((static_cast<T>(6) * pi<T, policies::policy<policies::digits2<N> > >()
0959         * pi<T, policies::policy<policies::digits2<N> > >())
0960    - (static_cast<T>(24) * pi<T, policies::policy<policies::digits2<N> > >()) + static_cast<T>(16) )
0961    /
0962    ((pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(4))
0963    * (pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(4)))
0964    );
0965 }
0966
0967 template <class T>
0968 template<int N>
0969 inline T constant_log2_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0970 {
0971    return 1 / boost::math::constants::ln_two<T>();
0972 }
0973
0974 template <class T>
0975 template<int N>
0976 inline T constant_quarter_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0977 {
0978    return boost::math::constants::pi<T>() / 4;
0979 }
0980
0981 template <class T>
0982 template<int N>
0983 inline T constant_one_div_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0984 {
0985    return 1 / boost::math::constants::pi<T>();
0986 }
0987
0988 template <class T>
0989 template<int N>
0990 inline T constant_two_div_root_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0991 {
0992    return 2 * boost::math::constants::one_div_root_pi<T>();
0993 }
0994
0995 #if __cplusplus >= 201103L || (defined(_MSC_VER) && _MSC_VER >= 1900)
0996 template <class T>
0997 template<int N>
0998 inline T constant_first_feigenbaum<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
0999 {
1000    // We know the constant to 1018 decimal digits.
1001    // See:  http://www.plouffe.fr/simon/constants/feigenbaum.txt
1002    // Also: https://oeis.org/A006890
1003    // N is in binary digits; so we multiply by log_2(10)
1004
1005    static_assert(N < 3.321*1018, "\nThe first Feigenbaum constant cannot be computed at runtime; it is too expensive. It is known to 1018 decimal digits; you must request less than that.");
1006    T alpha{"4.6692016091029906718532038204662016172581855774757686327456513430041343302113147371386897440239480138171659848551898151344086271420279325223124429888908908599449354632367134115324817142199474556443658237932020095610583305754586176522220703854106467494942849814533917262005687556659523398756038256372256480040951071283890611844702775854285419801113440175002428585382498335715522052236087250291678860362674527213399057131606875345083433934446103706309452019115876972432273589838903794946257251289097948986768334611626889116563123474460575179539122045562472807095202198199094558581946136877445617396074115614074243754435499204869180982648652368438702799649017397793425134723808737136211601860128186102056381818354097598477964173900328936171432159878240789776614391395764037760537119096932066998361984288981837003229412030210655743295550388845849737034727532121925706958414074661841981961006129640161487712944415901405467941800198133253378592493365883070459999938375411726563553016862529032210862320550634510679399023341675"};
1007    return alpha;
1008 }
1009
1010 template <class T>
1011 template<int N>
1012 inline T constant_plastic<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
1013 {
1014    using std::sqrt;
1015    return (cbrt(9-sqrt(T(69))) + cbrt(9+sqrt(T(69))))/cbrt(T(18));
1016 }
1017
1018
1019 template <class T>
1020 template<int N>
1021 inline T constant_gauss<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
1022 {
1023    using std::sqrt;
1024    T a = sqrt(T(2));
1025    T g = 1;
1026    const T scale = sqrt(std::numeric_limits<T>::epsilon())/512;
1027    while (a-g > scale*g)
1028    {
1029       T anp1 = (a + g)/2;
1030       g = sqrt(a*g);
1031       a = anp1;
1032    }
1033
1034    return 2/(a + g);
1035 }
1036
1037 template <class T>
1038 template<int N>
1039 inline T constant_dottie<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
1040 {
1041   // Error analysis: cos(x(1+d)) - x(1+d) = -(sin(x)+1)xd; plug in x = 0.739 gives -1.236d; take d as half an ulp gives the termination criteria we want.
1042   using std::cos;
1043   using std::abs;
1044   using std::sin;
1045   T x{".739085133215160641655312087673873404013411758900757464965680635773284654883547594599376106931766531849801246"};
1046   T residual = cos(x) - x;
1047   do {
1048     x += residual/(sin(x)+1);
1049     residual = cos(x) - x;
1050   } while(abs(residual) > std::numeric_limits<T>::epsilon());
1051   return x;
1052 }
1053
1054
1055 template <class T>
1056 template<int N>
1057 inline T constant_reciprocal_fibonacci<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
1058 {
1059   // Wikipedia says Gosper has deviced a faster algorithm for this, but I read the linked paper and couldn't see it!
1060   // In any case, k bits per iteration is fine, though it would be better to sum from smallest to largest.
1061   // That said, the condition number is unity, so it should be fine.
1062   T x0 = 1;
1063   T x1 = 1;
1064   T sum = 2;
1065   T diff = 1;
1066   while (diff > std::numeric_limits<T>::epsilon()) {
1067     T tmp = x1 + x0;
1068     diff = 1/tmp;
1069     sum += diff;
1070     x0 = x1;
1071     x1 = tmp;
1072   }
1073   return sum;
1074 }
1075
1076 template <class T>
1077 template<int N>
1078 inline T constant_laplace_limit<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
1079 {
1080   // If x is the exact root, then the approximate root is given by x(1+delta).
1081   // Plugging this into the equation for the Laplace limit gives the residual of approximately
1082   // 2.6389delta. Take delta as half an epsilon and give some leeway so we don't get caught in an infinite loop,
1083   // gives a termination condition as 2eps.
1084   using std::abs;
1085   using std::exp;
1086   using std::sqrt;
1087   T x{"0.66274341934918158097474209710925290705623354911502241752039253499097185308651127724965480259895818168"};
1088   T tmp = sqrt(1+x*x);
1089   T etmp = exp(tmp);
1090   T residual = x*exp(tmp) - 1 - tmp;
1091   T df = etmp -x/tmp + etmp*x*x/tmp;
1092   do {
1093     x -= residual/df;
1094     tmp = sqrt(1+x*x);
1095     etmp = exp(tmp);
1096     residual = x*exp(tmp) - 1 - tmp;
1097     df = etmp -x/tmp + etmp*x*x/tmp;
1098   } while(abs(residual) > 2*std::numeric_limits<T>::epsilon());
1099   return x;
1100 }
1101
1102 #endif
1103
1104 }
1105 }
1106 }
1107 } // namespaces
1108
1109 #endif // BOOST_MATH_CALCULATE_CONSTANTS_CONSTANTS_INCLUDED