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