ceph/src/boost/boost/math/constants/calculate_constants.hpp

   1 //  Copyright John Maddock 2010, 2012.
   2 //  Copyright Paul A. Bristow 2011, 2012.
   3
   4 //  Use, modification and distribution are subject to the
   5 //  Boost Software License, Version 1.0. (See accompanying file
   6 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   7
   8 #ifndef BOOST_MATH_CALCULATE_CONSTANTS_CONSTANTS_INCLUDED
   9 #define BOOST_MATH_CALCULATE_CONSTANTS_CONSTANTS_INCLUDED
  10 #include <type_traits>
  11
  12 namespace boost{ namespace math{ namespace constants{ namespace detail{
  13
  14 template <class T>
  15 template<int N>
  16 inline T constant_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
  17 {
  18    BOOST_MATH_STD_USING
  19
  20    return ldexp(acos(T(0)), 1);
  21
  22    /*
  23    // Although this code works well, it's usually more accurate to just call acos
  24    // and access the number types own representation of PI which is usually calculated
  25    // at slightly higher precision...
  26
  27    T result;
  28    T a = 1;
  29    T b;
  30    T A(a);
  31    T B = 0.5f;
  32    T D = 0.25f;
  33
  34    T lim;
  35    lim = boost::math::tools::epsilon<T>();
  36
  37    unsigned k = 1;
  38
  39    do
  40    {
  41       result = A + B;
  42       result = ldexp(result, -2);
  43       b = sqrt(B);
  44       a += b;
  45       a = ldexp(a, -1);
  46       A = a * a;
  47       B = A - result;
  48       B = ldexp(B, 1);
  49       result = A - B;
  50       bool neg = boost::math::sign(result) < 0;
  51       if(neg)
  52          result = -result;
  53       if(result <= lim)
  54          break;
  55       if(neg)
  56          result = -result;
  57       result = ldexp(result, k - 1);
  58       D -= result;
  59       ++k;
  60       lim = ldexp(lim, 1);
  61    }
  62    while(true);
  63
  64    result = B / D;
  65    return result;
  66    */
  67 }
  68
  69 template <class T>
  70 template<int N>
  71 inline T constant_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
  72 {
  73    return 2 * pi<T, policies::policy<policies::digits2<N> > >();
  74 }
  75
  76 template <class T> // 2 / pi
  77 template<int N>
  78 inline T constant_two_div_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
  79 {
  80    return 2 / pi<T, policies::policy<policies::digits2<N> > >();
  81 }
  82
  83 template <class T>  // sqrt(2/pi)
  84 template <int N>
  85 inline T constant_root_two_div_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
  86 {
  87    BOOST_MATH_STD_USING
  88    return sqrt((2 / pi<T, policies::policy<policies::digits2<N> > >()));
  89 }
  90
  91 template <class T>
  92 template<int N>
  93 inline T constant_one_div_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
  94 {
  95    return 1 / two_pi<T, policies::policy<policies::digits2<N> > >();
  96 }
  97
  98 template <class T>
  99 template<int N>
 100 inline T constant_root_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 101 {
 102    BOOST_MATH_STD_USING
 103    return sqrt(pi<T, policies::policy<policies::digits2<N> > >());
 104 }
 105
 106 template <class T>
 107 template<int N>
 108 inline T constant_root_half_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 109 {
 110    BOOST_MATH_STD_USING
 111    return sqrt(pi<T, policies::policy<policies::digits2<N> > >() / 2);
 112 }
 113
 114 template <class T>
 115 template<int N>
 116 inline T constant_root_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 117 {
 118    BOOST_MATH_STD_USING
 119    return sqrt(two_pi<T, policies::policy<policies::digits2<N> > >());
 120 }
 121
 122 template <class T>
 123 template<int N>
 124 inline T constant_log_root_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 125 {
 126    BOOST_MATH_STD_USING
 127    return log(root_two_pi<T, policies::policy<policies::digits2<N> > >());
 128 }
 129
 130 template <class T>
 131 template<int N>
 132 inline T constant_root_ln_four<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 133 {
 134    BOOST_MATH_STD_USING
 135    return sqrt(log(static_cast<T>(4)));
 136 }
 137
 138 template <class T>
 139 template<int N>
 140 inline T constant_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 141 {
 142    //
 143    // Although we can clearly calculate this from first principles, this hooks into
 144    // T's own notion of e, which hopefully will more accurate than one calculated to
 145    // a few epsilon:
 146    //
 147    BOOST_MATH_STD_USING
 148    return exp(static_cast<T>(1));
 149 }
 150
 151 template <class T>
 152 template<int N>
 153 inline T constant_half<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 154 {
 155    return static_cast<T>(1) / static_cast<T>(2);
 156 }
 157
 158 template <class T>
 159 template<int M>
 160 inline T constant_euler<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, M>)))
 161 {
 162    BOOST_MATH_STD_USING
 163    //
 164    // This is the method described in:
 165    // "Some New Algorithms for High-Precision Computation of Euler's Constant"
 166    // Richard P Brent and Edwin M McMillan.
 167    // Mathematics of Computation, Volume 34, Number 149, Jan 1980, pages 305-312.
 168    // See equation 17 with p = 2.
 169    //
 170    T n = 3 + (M ? (std::min)(M, tools::digits<T>()) : tools::digits<T>()) / 4;
 171    T lim = M ? ldexp(T(1), 1 - (std::min)(M, tools::digits<T>())) : tools::epsilon<T>();
 172    T lnn = log(n);
 173    T term = 1;
 174    T N = -lnn;
 175    T D = 1;
 176    T Hk = 0;
 177    T one = 1;
 178
 179    for(unsigned k = 1;; ++k)
 180    {
 181       term *= n * n;
 182       term /= k * k;
 183       Hk += one / k;
 184       N += term * (Hk - lnn);
 185       D += term;
 186
 187       if(term < D * lim)
 188          break;
 189    }
 190    return N / D;
 191 }
 192
 193 template <class T>
 194 template<int N>
 195 inline T constant_euler_sqr<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 196 {
 197   BOOST_MATH_STD_USING
 198   return euler<T, policies::policy<policies::digits2<N> > >()
 199      * euler<T, policies::policy<policies::digits2<N> > >();
 200 }
 201
 202 template <class T>
 203 template<int N>
 204 inline T constant_one_div_euler<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 205 {
 206   BOOST_MATH_STD_USING
 207   return static_cast<T>(1)
 208      / euler<T, policies::policy<policies::digits2<N> > >();
 209 }
 210
 211
 212 template <class T>
 213 template<int N>
 214 inline T constant_root_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 215 {
 216    BOOST_MATH_STD_USING
 217    return sqrt(static_cast<T>(2));
 218 }
 219
 220
 221 template <class T>
 222 template<int N>
 223 inline T constant_root_three<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 224 {
 225    BOOST_MATH_STD_USING
 226    return sqrt(static_cast<T>(3));
 227 }
 228
 229 template <class T>
 230 template<int N>
 231 inline T constant_half_root_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 232 {
 233    BOOST_MATH_STD_USING
 234    return sqrt(static_cast<T>(2)) / 2;
 235 }
 236
 237 template <class T>
 238 template<int N>
 239 inline T constant_ln_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 240 {
 241    //
 242    // Although there are good ways to calculate this from scratch, this hooks into
 243    // T's own notion of log(2) which will hopefully be accurate to the full precision
 244    // of T:
 245    //
 246    BOOST_MATH_STD_USING
 247    return log(static_cast<T>(2));
 248 }
 249
 250 template <class T>
 251 template<int N>
 252 inline T constant_ln_ten<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 253 {
 254    BOOST_MATH_STD_USING
 255    return log(static_cast<T>(10));
 256 }
 257
 258 template <class T>
 259 template<int N>
 260 inline T constant_ln_ln_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 261 {
 262    BOOST_MATH_STD_USING
 263    return log(log(static_cast<T>(2)));
 264 }
 265
 266 template <class T>
 267 template<int N>
 268 inline T constant_third<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 269 {
 270    BOOST_MATH_STD_USING
 271    return static_cast<T>(1) / static_cast<T>(3);
 272 }
 273
 274 template <class T>
 275 template<int N>
 276 inline T constant_twothirds<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 277 {
 278    BOOST_MATH_STD_USING
 279    return static_cast<T>(2) / static_cast<T>(3);
 280 }
 281
 282 template <class T>
 283 template<int N>
 284 inline T constant_two_thirds<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 285 {
 286    BOOST_MATH_STD_USING
 287    return static_cast<T>(2) / static_cast<T>(3);
 288 }
 289
 290 template <class T>
 291 template<int N>
 292 inline T constant_three_quarters<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 293 {
 294    BOOST_MATH_STD_USING
 295    return static_cast<T>(3) / static_cast<T>(4);
 296 }
 297
 298 template <class T>
 299 template<int N>
 300 inline T constant_sixth<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 301 {
 302   BOOST_MATH_STD_USING
 303   return static_cast<T>(1) / static_cast<T>(6);
 304 }
 305
 306 // Pi and related constants.
 307 template <class T>
 308 template<int N>
 309 inline T constant_pi_minus_three<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 310 {
 311    return pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(3);
 312 }
 313
 314 template <class T>
 315 template<int N>
 316 inline T constant_four_minus_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 317 {
 318    return static_cast<T>(4) - pi<T, policies::policy<policies::digits2<N> > >();
 319 }
 320
 321 template <class T>
 322 template<int N>
 323 inline T constant_exp_minus_half<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 324 {
 325    BOOST_MATH_STD_USING
 326    return exp(static_cast<T>(-0.5));
 327 }
 328
 329 template <class T>
 330 template<int N>
 331 inline T constant_exp_minus_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 332 {
 333   BOOST_MATH_STD_USING
 334   return exp(static_cast<T>(-1.));
 335 }
 336
 337 template <class T>
 338 template<int N>
 339 inline T constant_one_div_root_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 340 {
 341    return static_cast<T>(1) / root_two<T, policies::policy<policies::digits2<N> > >();
 342 }
 343
 344 template <class T>
 345 template<int N>
 346 inline T constant_one_div_root_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 347 {
 348    return static_cast<T>(1) / root_pi<T, policies::policy<policies::digits2<N> > >();
 349 }
 350
 351 template <class T>
 352 template<int N>
 353 inline T constant_one_div_root_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 354 {
 355    return static_cast<T>(1) / root_two_pi<T, policies::policy<policies::digits2<N> > >();
 356 }
 357
 358 template <class T>
 359 template<int N>
 360 inline T constant_root_one_div_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 361 {
 362    BOOST_MATH_STD_USING
 363    return sqrt(static_cast<T>(1) / pi<T, policies::policy<policies::digits2<N> > >());
 364 }
 365
 366 template <class T>
 367 template<int N>
 368 inline T constant_four_thirds_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 369 {
 370    BOOST_MATH_STD_USING
 371    return pi<T, policies::policy<policies::digits2<N> > >() * static_cast<T>(4) / static_cast<T>(3);
 372 }
 373
 374 template <class T>
 375 template<int N>
 376 inline T constant_half_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 377 {
 378    BOOST_MATH_STD_USING
 379    return pi<T, policies::policy<policies::digits2<N> > >()  / static_cast<T>(2);
 380 }
 381
 382 template <class T>
 383 template<int N>
 384 inline T constant_third_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 385 {
 386    BOOST_MATH_STD_USING
 387    return pi<T, policies::policy<policies::digits2<N> > >()  / static_cast<T>(3);
 388 }
 389
 390 template <class T>
 391 template<int N>
 392 inline T constant_sixth_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 393 {
 394    BOOST_MATH_STD_USING
 395    return pi<T, policies::policy<policies::digits2<N> > >()  / static_cast<T>(6);
 396 }
 397
 398 template <class T>
 399 template<int N>
 400 inline T constant_two_thirds_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 401 {
 402    BOOST_MATH_STD_USING
 403    return pi<T, policies::policy<policies::digits2<N> > >() * static_cast<T>(2) / static_cast<T>(3);
 404 }
 405
 406 template <class T>
 407 template<int N>
 408 inline T constant_three_quarters_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 409 {
 410    BOOST_MATH_STD_USING
 411    return pi<T, policies::policy<policies::digits2<N> > >() * static_cast<T>(3) / static_cast<T>(4);
 412 }
 413
 414 template <class T>
 415 template<int N>
 416 inline T constant_pi_pow_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 417 {
 418    BOOST_MATH_STD_USING
 419    return pow(pi<T, policies::policy<policies::digits2<N> > >(), e<T, policies::policy<policies::digits2<N> > >()); //
 420 }
 421
 422 template <class T>
 423 template<int N>
 424 inline T constant_pi_sqr<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 425 {
 426    BOOST_MATH_STD_USING
 427    return pi<T, policies::policy<policies::digits2<N> > >()
 428    *      pi<T, policies::policy<policies::digits2<N> > >() ; //
 429 }
 430
 431 template <class T>
 432 template<int N>
 433 inline T constant_pi_sqr_div_six<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 434 {
 435    BOOST_MATH_STD_USING
 436    return pi<T, policies::policy<policies::digits2<N> > >()
 437    *      pi<T, policies::policy<policies::digits2<N> > >()
 438    / static_cast<T>(6); //
 439 }
 440
 441 template <class T>
 442 template<int N>
 443 inline T constant_pi_cubed<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 444 {
 445    BOOST_MATH_STD_USING
 446    return pi<T, policies::policy<policies::digits2<N> > >()
 447    *      pi<T, policies::policy<policies::digits2<N> > >()
 448    *      pi<T, policies::policy<policies::digits2<N> > >()
 449    ; //
 450 }
 451
 452 template <class T>
 453 template<int N>
 454 inline T constant_cbrt_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 455 {
 456    BOOST_MATH_STD_USING
 457    return pow(pi<T, policies::policy<policies::digits2<N> > >(), static_cast<T>(1)/ static_cast<T>(3));
 458 }
 459
 460 template <class T>
 461 template<int N>
 462 inline T constant_one_div_cbrt_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 463 {
 464    BOOST_MATH_STD_USING
 465    return static_cast<T>(1)
 466    / pow(pi<T, policies::policy<policies::digits2<N> > >(), static_cast<T>(1)/ static_cast<T>(3));
 467 }
 468
 469 // Euler's e
 470
 471 template <class T>
 472 template<int N>
 473 inline T constant_e_pow_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 474 {
 475    BOOST_MATH_STD_USING
 476    return pow(e<T, policies::policy<policies::digits2<N> > >(), pi<T, policies::policy<policies::digits2<N> > >()); //
 477 }
 478
 479 template <class T>
 480 template<int N>
 481 inline T constant_root_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 482 {
 483    BOOST_MATH_STD_USING
 484    return sqrt(e<T, policies::policy<policies::digits2<N> > >());
 485 }
 486
 487 template <class T>
 488 template<int N>
 489 inline T constant_log10_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 490 {
 491    BOOST_MATH_STD_USING
 492    return log10(e<T, policies::policy<policies::digits2<N> > >());
 493 }
 494
 495 template <class T>
 496 template<int N>
 497 inline T constant_one_div_log10_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 498 {
 499    BOOST_MATH_STD_USING
 500    return  static_cast<T>(1) /
 501      log10(e<T, policies::policy<policies::digits2<N> > >());
 502 }
 503
 504 // Trigonometric
 505
 506 template <class T>
 507 template<int N>
 508 inline T constant_degree<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 509 {
 510    BOOST_MATH_STD_USING
 511    return pi<T, policies::policy<policies::digits2<N> > >()
 512    / static_cast<T>(180)
 513    ; //
 514 }
 515
 516 template <class T>
 517 template<int N>
 518 inline T constant_radian<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 519 {
 520    BOOST_MATH_STD_USING
 521    return static_cast<T>(180)
 522    / pi<T, policies::policy<policies::digits2<N> > >()
 523    ; //
 524 }
 525
 526 template <class T>
 527 template<int N>
 528 inline T constant_sin_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 529 {
 530    BOOST_MATH_STD_USING
 531    return sin(static_cast<T>(1)) ; //
 532 }
 533
 534 template <class T>
 535 template<int N>
 536 inline T constant_cos_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 537 {
 538    BOOST_MATH_STD_USING
 539    return cos(static_cast<T>(1)) ; //
 540 }
 541
 542 template <class T>
 543 template<int N>
 544 inline T constant_sinh_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 545 {
 546    BOOST_MATH_STD_USING
 547    return sinh(static_cast<T>(1)) ; //
 548 }
 549
 550 template <class T>
 551 template<int N>
 552 inline T constant_cosh_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 553 {
 554    BOOST_MATH_STD_USING
 555    return cosh(static_cast<T>(1)) ; //
 556 }
 557
 558 template <class T>
 559 template<int N>
 560 inline T constant_phi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 561 {
 562    BOOST_MATH_STD_USING
 563    return (static_cast<T>(1) + sqrt(static_cast<T>(5)) )/static_cast<T>(2) ; //
 564 }
 565
 566 template <class T>
 567 template<int N>
 568 inline T constant_ln_phi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 569 {
 570    BOOST_MATH_STD_USING
 571    return log((static_cast<T>(1) + sqrt(static_cast<T>(5)) )/static_cast<T>(2) );
 572 }
 573
 574 template <class T>
 575 template<int N>
 576 inline T constant_one_div_ln_phi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 577 {
 578    BOOST_MATH_STD_USING
 579    return static_cast<T>(1) /
 580      log((static_cast<T>(1) + sqrt(static_cast<T>(5)) )/static_cast<T>(2) );
 581 }
 582
 583 // Zeta
 584
 585 template <class T>
 586 template<int N>
 587 inline T constant_zeta_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 588 {
 589    BOOST_MATH_STD_USING
 590
 591      return pi<T, policies::policy<policies::digits2<N> > >()
 592      *  pi<T, policies::policy<policies::digits2<N> > >()
 593      /static_cast<T>(6);
 594 }
 595
 596 template <class T>
 597 template<int N>
 598 inline T constant_zeta_three<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 599 {
 600    // http://mathworld.wolfram.com/AperysConstant.html
 601    // http://en.wikipedia.org/wiki/Mathematical_constant
 602
 603    // http://oeis.org/A002117/constant
 604    //T zeta3("1.20205690315959428539973816151144999076"
 605    // "4986292340498881792271555341838205786313"
 606    // "09018645587360933525814619915");
 607
 608   //"1.202056903159594285399738161511449990, 76498629234049888179227155534183820578631309018645587360933525814619915"  A002117
 609   // 1.202056903159594285399738161511449990, 76498629234049888179227155534183820578631309018645587360933525814619915780, +00);
 610   //"1.2020569031595942 double
 611   // http://www.spaennare.se/SSPROG/ssnum.pdf  // section 11, Algorithm for Apery's constant zeta(3).
 612   // Programs to Calculate some Mathematical Constants to Large Precision, Document Version 1.50
 613
 614   // by Stefan Spannare  September 19, 2007
 615   // zeta(3) = 1/64 * sum
 616    BOOST_MATH_STD_USING
 617    T n_fact=static_cast<T>(1); // build n! for n = 0.
 618    T sum = static_cast<double>(77); // Start with n = 0 case.
 619    // for n = 0, (77/1) /64 = 1.203125
 620    //double lim = std::numeric_limits<double>::epsilon();
 621    T lim = N ? ldexp(T(1), 1 - (std::min)(N, tools::digits<T>())) : tools::epsilon<T>();
 622    for(unsigned int n = 1; n < 40; ++n)
 623    { // three to five decimal digits per term, so 40 should be plenty for 100 decimal digits.
 624       //cout << "n = " << n << endl;
 625       n_fact *= n; // n!
 626       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
 627       T num = ((205 * n * n) + (250 * n) + 77) * n_fact_p10; // 205n^2 + 250n + 77
 628       // int nn = (2 * n + 1);
 629       // T d = factorial(nn); // inline factorial.
 630       T d = 1;
 631       for(unsigned int i = 1; i <= (n+n + 1); ++i) // (2n + 1)
 632       {
 633         d *= i;
 634       }
 635       T den = d * d * d * d * d; // [(2n+1)!]^5
 636       //cout << "den = " << den << endl;
 637       T term = num/den;
 638       if (n % 2 != 0)
 639       { //term *= -1;
 640         sum -= term;
 641       }
 642       else
 643       {
 644         sum += term;
 645       }
 646       //cout << "term = " << term << endl;
 647       //cout << "sum/64  = " << sum/64 << endl;
 648       if(abs(term) < lim)
 649       {
 650          break;
 651       }
 652    }
 653    return sum / 64;
 654 }
 655
 656 template <class T>
 657 template<int N>
 658 inline T constant_catalan<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 659 { // http://oeis.org/A006752/constant
 660   //T c("0.915965594177219015054603514932384110774"
 661  //"149374281672134266498119621763019776254769479356512926115106248574");
 662
 663   // 9.159655941772190150546035149323841107, 74149374281672134266498119621763019776254769479356512926115106248574422619, -01);
 664
 665    // This is equation (entry) 31 from
 666    // http://www-2.cs.cmu.edu/~adamchik/articles/catalan/catalan.htm
 667    // See also http://www.mpfr.org/algorithms.pdf
 668    BOOST_MATH_STD_USING
 669    T k_fact = 1;
 670    T tk_fact = 1;
 671    T sum = 1;
 672    T term;
 673    T lim = N ? ldexp(T(1), 1 - (std::min)(N, tools::digits<T>())) : tools::epsilon<T>();
 674
 675    for(unsigned k = 1;; ++k)
 676    {
 677       k_fact *= k;
 678       tk_fact *= (2 * k) * (2 * k - 1);
 679       term = k_fact * k_fact / (tk_fact * (2 * k + 1) * (2 * k + 1));
 680       sum += term;
 681       if(term < lim)
 682       {
 683          break;
 684       }
 685    }
 686    return boost::math::constants::pi<T, boost::math::policies::policy<> >()
 687       * log(2 + boost::math::constants::root_three<T, boost::math::policies::policy<> >())
 688        / 8
 689       + 3 * sum / 8;
 690 }
 691
 692 namespace khinchin_detail{
 693
 694 template <class T>
 695 T zeta_polynomial_series(T s, T sc, int digits)
 696 {
 697    BOOST_MATH_STD_USING
 698    //
 699    // This is algorithm 3 from:
 700    //
 701    // "An Efficient Algorithm for the Riemann Zeta Function", P. Borwein,
 702    // Canadian Mathematical Society, Conference Proceedings, 2000.
 703    // See: http://www.cecm.sfu.ca/personal/pborwein/PAPERS/P155.pdf
 704    //
 705    BOOST_MATH_STD_USING
 706    int n = (digits * 19) / 53;
 707    T sum = 0;
 708    T two_n = ldexp(T(1), n);
 709    int ej_sign = 1;
 710    for(int j = 0; j < n; ++j)
 711    {
 712       sum += ej_sign * -two_n / pow(T(j + 1), s);
 713       ej_sign = -ej_sign;
 714    }
 715    T ej_sum = 1;
 716    T ej_term = 1;
 717    for(int j = n; j <= 2 * n - 1; ++j)
 718    {
 719       sum += ej_sign * (ej_sum - two_n) / pow(T(j + 1), s);
 720       ej_sign = -ej_sign;
 721       ej_term *= 2 * n - j;
 722       ej_term /= j - n + 1;
 723       ej_sum += ej_term;
 724    }
 725    return -sum / (two_n * (1 - pow(T(2), sc)));
 726 }
 727
 728 template <class T>
 729 T khinchin(int digits)
 730 {
 731    BOOST_MATH_STD_USING
 732    T sum = 0;
 733    T term;
 734    T lim = ldexp(T(1), 1-digits);
 735    T factor = 0;
 736    unsigned last_k = 1;
 737    T num = 1;
 738    for(unsigned n = 1;; ++n)
 739    {
 740       for(unsigned k = last_k; k <= 2 * n - 1; ++k)
 741       {
 742          factor += num / k;
 743          num = -num;
 744       }
 745       last_k = 2 * n;
 746       term = (zeta_polynomial_series(T(2 * n), T(1 - T(2 * n)), digits) - 1) * factor / n;
 747       sum += term;
 748       if(term < lim)
 749          break;
 750    }
 751    return exp(sum / boost::math::constants::ln_two<T, boost::math::policies::policy<> >());
 752 }
 753
 754 }
 755
 756 template <class T>
 757 template<int N>
 758 inline T constant_khinchin<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 759 {
 760    int n = N ? (std::min)(N, tools::digits<T>()) : tools::digits<T>();
 761    return khinchin_detail::khinchin<T>(n);
 762 }
 763
 764 template <class T>
 765 template<int N>
 766 inline T constant_extreme_value_skewness<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 767 {  // N[12 Sqrt[6]  Zeta[3]/Pi^3, 1101]
 768    BOOST_MATH_STD_USING
 769    T ev(12 * sqrt(static_cast<T>(6)) * zeta_three<T, policies::policy<policies::digits2<N> > >()
 770     / pi_cubed<T, policies::policy<policies::digits2<N> > >() );
 771
 772 //T ev(
 773 //"1.1395470994046486574927930193898461120875997958365518247216557100852480077060706857071875468869385150"
 774 //"1894272048688553376986765366075828644841024041679714157616857834895702411080704529137366329462558680"
 775 //"2015498788776135705587959418756809080074611906006528647805347822929577145038743873949415294942796280"
 776 //"0895597703063466053535550338267721294164578901640163603544404938283861127819804918174973533694090594"
 777 //"3094963822672055237678432023017824416203652657301470473548274848068762500300316769691474974950757965"
 778 //"8640779777748741897542093874605477776538884083378029488863880220988107155275203245233994097178778984"
 779 //"3488995668362387892097897322246698071290011857605809901090220903955815127463328974447572119951192970"
 780 //"3684453635456559086126406960279692862247058250100678008419431185138019869693206366891639436908462809"
 781 //"9756051372711251054914491837034685476095423926553367264355374652153595857163724698198860485357368964"
 782 //"3807049634423621246870868566707915720704996296083373077647528285782964567312903914752617978405994377"
 783 //"9064157147206717895272199736902453130842229559980076472936976287378945035706933650987259357729800315");
 784
 785   return ev;
 786 }
 787
 788 namespace detail{
 789 //
 790 // Calculation of the Glaisher constant depends upon calculating the
 791 // derivative of the zeta function at 2, we can then use the relation:
 792 // zeta'(2) = 1/6 pi^2 [euler + ln(2pi)-12ln(A)]
 793 // To get the constant A.
 794 // See equation 45 at http://mathworld.wolfram.com/RiemannZetaFunction.html.
 795 //
 796 // The derivative of the zeta function is computed by direct differentiation
 797 // of the relation:
 798 // (1-2^(1-s))zeta(s) = SUM(n=0, INF){ (-n)^n / (n+1)^s  }
 799 // Which gives us 2 slowly converging but alternating sums to compute,
 800 // for this we use Algorithm 1 from "Convergent Acceleration of Alternating Series",
 801 // Henri Cohen, Fernando Rodriguez Villegas and Don Zagier, Experimental Mathematics 9:1 (1999).
 802 // See http://www.math.utexas.edu/users/villegas/publications/conv-accel.pdf
 803 //
 804 template <class T>
 805 T zeta_series_derivative_2(unsigned digits)
 806 {
 807    // Derivative of the series part, evaluated at 2:
 808    BOOST_MATH_STD_USING
 809    int n = digits * 301 * 13 / 10000;
 810    T d = pow(3 + sqrt(T(8)), n);
 811    d = (d + 1 / d) / 2;
 812    T b = -1;
 813    T c = -d;
 814    T s = 0;
 815    for(int k = 0; k < n; ++k)
 816    {
 817       T a = -log(T(k+1)) / ((k+1) * (k+1));
 818       c = b - c;
 819       s = s + c * a;
 820       b = (k + n) * (k - n) * b / ((k + T(0.5f)) * (k + 1));
 821    }
 822    return s / d;
 823 }
 824
 825 template <class T>
 826 T zeta_series_2(unsigned digits)
 827 {
 828    // Series part of zeta at 2:
 829    BOOST_MATH_STD_USING
 830    int n = digits * 301 * 13 / 10000;
 831    T d = pow(3 + sqrt(T(8)), n);
 832    d = (d + 1 / d) / 2;
 833    T b = -1;
 834    T c = -d;
 835    T s = 0;
 836    for(int k = 0; k < n; ++k)
 837    {
 838       T a = T(1) / ((k + 1) * (k + 1));
 839       c = b - c;
 840       s = s + c * a;
 841       b = (k + n) * (k - n) * b / ((k + T(0.5f)) * (k + 1));
 842    }
 843    return s / d;
 844 }
 845
 846 template <class T>
 847 inline T zeta_series_lead_2()
 848 {
 849    // lead part at 2:
 850    return 2;
 851 }
 852
 853 template <class T>
 854 inline T zeta_series_derivative_lead_2()
 855 {
 856    // derivative of lead part at 2:
 857    return -2 * boost::math::constants::ln_two<T>();
 858 }
 859
 860 template <class T>
 861 inline T zeta_derivative_2(unsigned n)
 862 {
 863    // zeta derivative at 2:
 864    return zeta_series_derivative_2<T>(n) * zeta_series_lead_2<T>()
 865       + zeta_series_derivative_lead_2<T>() * zeta_series_2<T>(n);
 866 }
 867
 868 }  // namespace detail
 869
 870 template <class T>
 871 template<int N>
 872 inline T constant_glaisher<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 873 {
 874
 875    BOOST_MATH_STD_USING
 876    typedef policies::policy<policies::digits2<N> > forwarding_policy;
 877    int n = N ? (std::min)(N, tools::digits<T>()) : tools::digits<T>();
 878    T v = detail::zeta_derivative_2<T>(n);
 879    v *= 6;
 880    v /= boost::math::constants::pi<T, forwarding_policy>() * boost::math::constants::pi<T, forwarding_policy>();
 881    v -= boost::math::constants::euler<T, forwarding_policy>();
 882    v -= log(2 * boost::math::constants::pi<T, forwarding_policy>());
 883    v /= -12;
 884    return exp(v);
 885
 886    /*
 887    // from http://mpmath.googlecode.com/svn/data/glaisher.txt
 888      // 20,000 digits of the Glaisher-Kinkelin constant A = exp(1/2 - zeta'(-1))
 889      // Computed using A = exp((6 (-zeta'(2))/pi^2 + log 2 pi + gamma)/12)
 890    // with Euler-Maclaurin summation for zeta'(2).
 891    T g(
 892    "1.282427129100622636875342568869791727767688927325001192063740021740406308858826"
 893    "46112973649195820237439420646120399000748933157791362775280404159072573861727522"
 894    "14334327143439787335067915257366856907876561146686449997784962754518174312394652"
 895    "76128213808180219264516851546143919901083573730703504903888123418813674978133050"
 896    "93770833682222494115874837348064399978830070125567001286994157705432053927585405"
 897    "81731588155481762970384743250467775147374600031616023046613296342991558095879293"
 898    "36343887288701988953460725233184702489001091776941712153569193674967261270398013"
 899    "52652668868978218897401729375840750167472114895288815996668743164513890306962645"
 900    "59870469543740253099606800842447417554061490189444139386196089129682173528798629"
 901    "88434220366989900606980888785849587494085307347117090132667567503310523405221054"
 902    "14176776156308191919997185237047761312315374135304725819814797451761027540834943"
 903    "14384965234139453373065832325673954957601692256427736926358821692159870775858274"
 904    "69575162841550648585890834128227556209547002918593263079373376942077522290940187");
 905
 906    return g;
 907    */
 908 }
 909
 910 template <class T>
 911 template<int N>
 912 inline T constant_rayleigh_skewness<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 913 {  // 1100 digits of the Rayleigh distribution skewness
 914    // N[2 Sqrt[Pi] (Pi - 3)/((4 - Pi)^(3/2)), 1100]
 915
 916    BOOST_MATH_STD_USING
 917    T rs(2 * root_pi<T, policies::policy<policies::digits2<N> > >()
 918       * pi_minus_three<T, policies::policy<policies::digits2<N> > >()
 919       / pow(four_minus_pi<T, policies::policy<policies::digits2<N> > >(), static_cast<T>(3./2))
 920       );
 921    //   6.31110657818937138191899351544227779844042203134719497658094585692926819617473725459905027032537306794400047264,
 922
 923    //"0.6311106578189371381918993515442277798440422031347194976580945856929268196174737254599050270325373067"
 924    //"9440004726436754739597525250317640394102954301685809920213808351450851396781817932734836994829371322"
 925    //"5797376021347531983451654130317032832308462278373358624120822253764532674177325950686466133508511968"
 926    //"2389168716630349407238090652663422922072397393006683401992961569208109477307776249225072042971818671"
 927    //"4058887072693437217879039875871765635655476241624825389439481561152126886932506682176611183750503553"
 928    //"1218982627032068396407180216351425758181396562859085306247387212297187006230007438534686340210168288"
 929    //"8956816965453815849613622117088096547521391672977226658826566757207615552041767516828171274858145957"
 930    //"6137539156656005855905288420585194082284972984285863898582313048515484073396332610565441264220790791"
 931    //"0194897267890422924599776483890102027823328602965235306539844007677157873140562950510028206251529523"
 932    //"7428049693650605954398446899724157486062545281504433364675815915402937209673727753199567661561209251"
 933    //"4695589950526053470201635372590001578503476490223746511106018091907936826431407434894024396366284848");  ;
 934    return rs;
 935 }
 936
 937 template <class T>
 938 template<int N>
 939 inline T constant_rayleigh_kurtosis_excess<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 940 { // - (6 Pi^2 - 24 Pi + 16)/((Pi - 4)^2)
 941     // Might provide and calculate this using pi_minus_four.
 942    BOOST_MATH_STD_USING
 943    return - (((static_cast<T>(6) * pi<T, policies::policy<policies::digits2<N> > >()
 944         * pi<T, policies::policy<policies::digits2<N> > >())
 945    - (static_cast<T>(24) * pi<T, policies::policy<policies::digits2<N> > >()) + static_cast<T>(16) )
 946    /
 947    ((pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(4))
 948    * (pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(4)))
 949    );
 950 }
 951
 952 template <class T>
 953 template<int N>
 954 inline T constant_rayleigh_kurtosis<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 955 { // 3 - (6 Pi^2 - 24 Pi + 16)/((Pi - 4)^2)
 956   // Might provide and calculate this using pi_minus_four.
 957    BOOST_MATH_STD_USING
 958    return static_cast<T>(3) - (((static_cast<T>(6) * pi<T, policies::policy<policies::digits2<N> > >()
 959         * pi<T, policies::policy<policies::digits2<N> > >())
 960    - (static_cast<T>(24) * pi<T, policies::policy<policies::digits2<N> > >()) + static_cast<T>(16) )
 961    /
 962    ((pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(4))
 963    * (pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(4)))
 964    );
 965 }
 966
 967 template <class T>
 968 template<int N>
 969 inline T constant_log2_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 970 {
 971    return 1 / boost::math::constants::ln_two<T>();
 972 }
 973
 974 template <class T>
 975 template<int N>
 976 inline T constant_quarter_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 977 {
 978    return boost::math::constants::pi<T>() / 4;
 979 }
 980
 981 template <class T>
 982 template<int N>
 983 inline T constant_one_div_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 984 {
 985    return 1 / boost::math::constants::pi<T>();
 986 }
 987
 988 template <class T>
 989 template<int N>
 990 inline T constant_two_div_root_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 991 {
 992    return 2 * boost::math::constants::one_div_root_pi<T>();
 993 }
 994
 995 #if __cplusplus >= 201103L || (defined(_MSC_VER) && _MSC_VER >= 1900)
 996 template <class T>
 997 template<int N>
 998 inline T constant_first_feigenbaum<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
 999 {
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::cbrt;
1015    using std::sqrt;
1016    return (cbrt(9-sqrt(T(69))) + cbrt(9+sqrt(T(69))))/cbrt(T(18));
1017 }
1018
1019
1020 template <class T>
1021 template<int N>
1022 inline T constant_gauss<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
1023 {
1024    using std::sqrt;
1025    T a = sqrt(T(2));
1026    T g = 1;
1027    const T scale = sqrt(std::numeric_limits<T>::epsilon())/512;
1028    while (a-g > scale*g)
1029    {
1030       T anp1 = (a + g)/2;
1031       g = sqrt(a*g);
1032       a = anp1;
1033    }
1034
1035    return 2/(a + g);
1036 }
1037
1038 template <class T>
1039 template<int N>
1040 inline T constant_dottie<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
1041 {
1042   // 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.
1043   using std::cos;
1044   using std::abs;
1045   using std::sin;
1046   T x{".739085133215160641655312087673873404013411758900757464965680635773284654883547594599376106931766531849801246"};
1047   T residual = cos(x) - x;
1048   do {
1049     x += residual/(sin(x)+1);
1050     residual = cos(x) - x;
1051   } while(abs(residual) > std::numeric_limits<T>::epsilon());
1052   return x;
1053 }
1054
1055
1056 template <class T>
1057 template<int N>
1058 inline T constant_reciprocal_fibonacci<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
1059 {
1060   // Wikipedia says Gosper has deviced a faster algorithm for this, but I read the linked paper and couldn't see it!
1061   // In any case, k bits per iteration is fine, though it would be better to sum from smallest to largest.
1062   // That said, the condition number is unity, so it should be fine.
1063   T x0 = 1;
1064   T x1 = 1;
1065   T sum = 2;
1066   T diff = 1;
1067   while (diff > std::numeric_limits<T>::epsilon()) {
1068     T tmp = x1 + x0;
1069     diff = 1/tmp;
1070     sum += diff;
1071     x0 = x1;
1072     x1 = tmp;
1073   }
1074   return sum;
1075 }
1076
1077 template <class T>
1078 template<int N>
1079 inline T constant_laplace_limit<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
1080 {
1081   // If x is the exact root, then the approximate root is given by x(1+delta).
1082   // Plugging this into the equation for the Laplace limit gives the residual of approximately
1083   // 2.6389delta. Take delta as half an epsilon and give some leeway so we don't get caught in an infinite loop,
1084   // gives a termination condition as 2eps.
1085   using std::abs;
1086   using std::exp;
1087   using std::sqrt;
1088   T x{"0.66274341934918158097474209710925290705623354911502241752039253499097185308651127724965480259895818168"};
1089   T tmp = sqrt(1+x*x);
1090   T etmp = exp(tmp);
1091   T residual = x*exp(tmp) - 1 - tmp;
1092   T df = etmp -x/tmp + etmp*x*x/tmp;
1093   do {
1094     x -= residual/df;
1095     tmp = sqrt(1+x*x);
1096     etmp = exp(tmp);
1097     residual = x*exp(tmp) - 1 - tmp;
1098     df = etmp -x/tmp + etmp*x*x/tmp;
1099   } while(abs(residual) > 2*std::numeric_limits<T>::epsilon());
1100   return x;
1101 }
1102
1103 #endif
1104
1105 }
1106 }
1107 }
1108 } // namespaces
1109
1110 #endif // BOOST_MATH_CALCULATE_CONSTANTS_CONSTANTS_INCLUDED