1110 lines · plain
1// Copyright John Maddock 2010, 2012.2// Copyright Paul A. Bristow 2011, 2012.3 4// Use, modification and distribution are subject to the5// Boost Software License, Version 1.0. (See accompanying file6// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)7 8#ifndef BOOST_MATH_CALCULATE_CONSTANTS_CONSTANTS_INCLUDED9#define BOOST_MATH_CALCULATE_CONSTANTS_CONSTANTS_INCLUDED10#include <type_traits>11 12namespace boost{ namespace math{ namespace constants{ namespace detail{13 14template <class T>15template<int N>16inline T constant_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))17{18 BOOST_MATH_STD_USING19 20 return ldexp(acos(T(0)), 1);21 22 /*23 // Although this code works well, it's usually more accurate to just call acos24 // and access the number types own representation of PI which is usually calculated25 // 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 do40 {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 69template <class T>70template<int N>71inline 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 76template <class T> // 2 / pi77template<int N>78inline 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 83template <class T> // sqrt(2/pi)84template <int N>85inline T constant_root_two_div_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))86{87 BOOST_MATH_STD_USING88 return sqrt((2 / pi<T, policies::policy<policies::digits2<N> > >()));89}90 91template <class T>92template<int N>93inline 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 98template <class T>99template<int N>100inline T constant_root_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))101{102 BOOST_MATH_STD_USING103 return sqrt(pi<T, policies::policy<policies::digits2<N> > >());104}105 106template <class T>107template<int N>108inline T constant_root_half_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))109{110 BOOST_MATH_STD_USING111 return sqrt(pi<T, policies::policy<policies::digits2<N> > >() / 2);112}113 114template <class T>115template<int N>116inline T constant_root_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))117{118 BOOST_MATH_STD_USING119 return sqrt(two_pi<T, policies::policy<policies::digits2<N> > >());120}121 122template <class T>123template<int N>124inline T constant_log_root_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))125{126 BOOST_MATH_STD_USING127 return log(root_two_pi<T, policies::policy<policies::digits2<N> > >());128}129 130template <class T>131template<int N>132inline T constant_root_ln_four<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))133{134 BOOST_MATH_STD_USING135 return sqrt(log(static_cast<T>(4)));136}137 138template <class T>139template<int N>140inline 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 into144 // T's own notion of e, which hopefully will more accurate than one calculated to145 // a few epsilon:146 //147 BOOST_MATH_STD_USING148 return exp(static_cast<T>(1));149}150 151template <class T>152template<int N>153inline 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 158template <class T>159template<int M>160inline T constant_euler<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, M>)))161{162 BOOST_MATH_STD_USING163 //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 193template <class T>194template<int N>195inline T constant_euler_sqr<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))196{197 BOOST_MATH_STD_USING198 return euler<T, policies::policy<policies::digits2<N> > >()199 * euler<T, policies::policy<policies::digits2<N> > >();200}201 202template <class T>203template<int N>204inline T constant_one_div_euler<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))205{206 BOOST_MATH_STD_USING207 return static_cast<T>(1)208 / euler<T, policies::policy<policies::digits2<N> > >();209}210 211 212template <class T>213template<int N>214inline T constant_root_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))215{216 BOOST_MATH_STD_USING217 return sqrt(static_cast<T>(2));218}219 220 221template <class T>222template<int N>223inline T constant_root_three<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))224{225 BOOST_MATH_STD_USING226 return sqrt(static_cast<T>(3));227}228 229template <class T>230template<int N>231inline T constant_half_root_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))232{233 BOOST_MATH_STD_USING234 return sqrt(static_cast<T>(2)) / 2;235}236 237template <class T>238template<int N>239inline 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 into243 // T's own notion of log(2) which will hopefully be accurate to the full precision244 // of T:245 //246 BOOST_MATH_STD_USING247 return log(static_cast<T>(2));248}249 250template <class T>251template<int N>252inline T constant_ln_ten<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))253{254 BOOST_MATH_STD_USING255 return log(static_cast<T>(10));256}257 258template <class T>259template<int N>260inline T constant_ln_ln_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))261{262 BOOST_MATH_STD_USING263 return log(log(static_cast<T>(2)));264}265 266template <class T>267template<int N>268inline T constant_third<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))269{270 BOOST_MATH_STD_USING271 return static_cast<T>(1) / static_cast<T>(3);272}273 274template <class T>275template<int N>276inline T constant_twothirds<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))277{278 BOOST_MATH_STD_USING279 return static_cast<T>(2) / static_cast<T>(3);280}281 282template <class T>283template<int N>284inline T constant_two_thirds<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))285{286 BOOST_MATH_STD_USING287 return static_cast<T>(2) / static_cast<T>(3);288}289 290template <class T>291template<int N>292inline T constant_three_quarters<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))293{294 BOOST_MATH_STD_USING295 return static_cast<T>(3) / static_cast<T>(4);296}297 298template <class T>299template<int N>300inline T constant_sixth<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))301{302 BOOST_MATH_STD_USING303 return static_cast<T>(1) / static_cast<T>(6);304}305 306// Pi and related constants.307template <class T>308template<int N>309inline 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 314template <class T>315template<int N>316inline 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 321template <class T>322template<int N>323inline T constant_exp_minus_half<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))324{325 BOOST_MATH_STD_USING326 return exp(static_cast<T>(-0.5));327}328 329template <class T>330template<int N>331inline T constant_exp_minus_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))332{333 BOOST_MATH_STD_USING334 return exp(static_cast<T>(-1.));335}336 337template <class T>338template<int N>339inline 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 344template <class T>345template<int N>346inline 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 351template <class T>352template<int N>353inline 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 358template <class T>359template<int N>360inline T constant_root_one_div_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))361{362 BOOST_MATH_STD_USING363 return sqrt(static_cast<T>(1) / pi<T, policies::policy<policies::digits2<N> > >());364}365 366template <class T>367template<int N>368inline T constant_four_thirds_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))369{370 BOOST_MATH_STD_USING371 return pi<T, policies::policy<policies::digits2<N> > >() * static_cast<T>(4) / static_cast<T>(3);372}373 374template <class T>375template<int N>376inline T constant_half_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))377{378 BOOST_MATH_STD_USING379 return pi<T, policies::policy<policies::digits2<N> > >() / static_cast<T>(2);380}381 382template <class T>383template<int N>384inline T constant_third_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))385{386 BOOST_MATH_STD_USING387 return pi<T, policies::policy<policies::digits2<N> > >() / static_cast<T>(3);388}389 390template <class T>391template<int N>392inline T constant_sixth_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))393{394 BOOST_MATH_STD_USING395 return pi<T, policies::policy<policies::digits2<N> > >() / static_cast<T>(6);396}397 398template <class T>399template<int N>400inline T constant_two_thirds_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))401{402 BOOST_MATH_STD_USING403 return pi<T, policies::policy<policies::digits2<N> > >() * static_cast<T>(2) / static_cast<T>(3);404}405 406template <class T>407template<int N>408inline T constant_three_quarters_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))409{410 BOOST_MATH_STD_USING411 return pi<T, policies::policy<policies::digits2<N> > >() * static_cast<T>(3) / static_cast<T>(4);412}413 414template <class T>415template<int N>416inline T constant_pi_pow_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))417{418 BOOST_MATH_STD_USING419 return pow(pi<T, policies::policy<policies::digits2<N> > >(), e<T, policies::policy<policies::digits2<N> > >()); //420}421 422template <class T>423template<int N>424inline T constant_pi_sqr<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))425{426 BOOST_MATH_STD_USING427 return pi<T, policies::policy<policies::digits2<N> > >()428 * pi<T, policies::policy<policies::digits2<N> > >() ; //429}430 431template <class T>432template<int N>433inline T constant_pi_sqr_div_six<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))434{435 BOOST_MATH_STD_USING436 return pi<T, policies::policy<policies::digits2<N> > >()437 * pi<T, policies::policy<policies::digits2<N> > >()438 / static_cast<T>(6); //439}440 441template <class T>442template<int N>443inline T constant_pi_cubed<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))444{445 BOOST_MATH_STD_USING446 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 452template <class T>453template<int N>454inline T constant_cbrt_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))455{456 BOOST_MATH_STD_USING457 return pow(pi<T, policies::policy<policies::digits2<N> > >(), static_cast<T>(1)/ static_cast<T>(3));458}459 460template <class T>461template<int N>462inline T constant_one_div_cbrt_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))463{464 BOOST_MATH_STD_USING465 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 e470 471template <class T>472template<int N>473inline T constant_e_pow_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))474{475 BOOST_MATH_STD_USING476 return pow(e<T, policies::policy<policies::digits2<N> > >(), pi<T, policies::policy<policies::digits2<N> > >()); //477}478 479template <class T>480template<int N>481inline T constant_root_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))482{483 BOOST_MATH_STD_USING484 return sqrt(e<T, policies::policy<policies::digits2<N> > >());485}486 487template <class T>488template<int N>489inline T constant_log10_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))490{491 BOOST_MATH_STD_USING492 return log10(e<T, policies::policy<policies::digits2<N> > >());493}494 495template <class T>496template<int N>497inline T constant_one_div_log10_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))498{499 BOOST_MATH_STD_USING500 return static_cast<T>(1) /501 log10(e<T, policies::policy<policies::digits2<N> > >());502}503 504// Trigonometric505 506template <class T>507template<int N>508inline T constant_degree<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))509{510 BOOST_MATH_STD_USING511 return pi<T, policies::policy<policies::digits2<N> > >()512 / static_cast<T>(180)513 ; //514}515 516template <class T>517template<int N>518inline T constant_radian<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))519{520 BOOST_MATH_STD_USING521 return static_cast<T>(180)522 / pi<T, policies::policy<policies::digits2<N> > >()523 ; //524}525 526template <class T>527template<int N>528inline T constant_sin_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))529{530 BOOST_MATH_STD_USING531 return sin(static_cast<T>(1)) ; //532}533 534template <class T>535template<int N>536inline T constant_cos_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))537{538 BOOST_MATH_STD_USING539 return cos(static_cast<T>(1)) ; //540}541 542template <class T>543template<int N>544inline T constant_sinh_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))545{546 BOOST_MATH_STD_USING547 return sinh(static_cast<T>(1)) ; //548}549 550template <class T>551template<int N>552inline T constant_cosh_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))553{554 BOOST_MATH_STD_USING555 return cosh(static_cast<T>(1)) ; //556}557 558template <class T>559template<int N>560inline T constant_phi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))561{562 BOOST_MATH_STD_USING563 return (static_cast<T>(1) + sqrt(static_cast<T>(5)) )/static_cast<T>(2) ; //564}565 566template <class T>567template<int N>568inline T constant_ln_phi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))569{570 BOOST_MATH_STD_USING571 return log((static_cast<T>(1) + sqrt(static_cast<T>(5)) )/static_cast<T>(2) );572}573 574template <class T>575template<int N>576inline T constant_one_div_ln_phi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))577{578 BOOST_MATH_STD_USING579 return static_cast<T>(1) /580 log((static_cast<T>(1) + sqrt(static_cast<T>(5)) )/static_cast<T>(2) );581}582 583// Zeta584 585template <class T>586template<int N>587inline T constant_zeta_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))588{589 BOOST_MATH_STD_USING590 591 return pi<T, policies::policy<policies::digits2<N> > >()592 * pi<T, policies::policy<policies::digits2<N> > >()593 /static_cast<T>(6);594}595 596template <class T>597template<int N>598inline T constant_zeta_three<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))599{600 // http://mathworld.wolfram.com/AperysConstant.html601 // http://en.wikipedia.org/wiki/Mathematical_constant602 603 // http://oeis.org/A002117/constant604 //T zeta3("1.20205690315959428539973816151144999076"605 // "4986292340498881792271555341838205786313"606 // "09018645587360933525814619915");607 608 //"1.202056903159594285399738161511449990, 76498629234049888179227155534183820578631309018645587360933525814619915" A002117609 // 1.202056903159594285399738161511449990, 76498629234049888179227155534183820578631309018645587360933525814619915780, +00);610 //"1.2020569031595942 double611 // 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.50613 614 // by Stefan Spannare September 19, 2007615 // zeta(3) = 1/64 * sum616 BOOST_MATH_STD_USING617 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.203125620 //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!)^10627 T num = ((205 * n * n) + (250 * n) + 77) * n_fact_p10; // 205n^2 + 250n + 77628 // 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)!]^5636 //cout << "den = " << den << endl;637 T term = num/den;638 if (n % 2 != 0)639 { //term *= -1;640 sum -= term;641 }642 else643 {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 656template <class T>657template<int N>658inline T constant_catalan<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))659{ // http://oeis.org/A006752/constant660 //T c("0.915965594177219015054603514932384110774"661 //"149374281672134266498119621763019776254769479356512926115106248574");662 663 // 9.159655941772190150546035149323841107, 74149374281672134266498119621763019776254769479356512926115106248574422619, -01);664 665 // This is equation (entry) 31 from666 // http://www-2.cs.cmu.edu/~adamchik/articles/catalan/catalan.htm667 // See also http://www.mpfr.org/algorithms.pdf668 BOOST_MATH_STD_USING669 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 / 8689 + 3 * sum / 8;690}691 692namespace khinchin_detail{693 694template <class T>695T zeta_polynomial_series(T s, T sc, int digits)696{697 BOOST_MATH_STD_USING698 //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.pdf704 //705 BOOST_MATH_STD_USING706 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 728template <class T>729T khinchin(int digits)730{731 BOOST_MATH_STD_USING732 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 756template <class T>757template<int N>758inline 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 764template <class T>765template<int N>766inline 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_USING769 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 788namespace detail{789//790// Calculation of the Glaisher constant depends upon calculating the791// 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 differentiation797// 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.pdf803//804template <class T>805T zeta_series_derivative_2(unsigned digits)806{807 // Derivative of the series part, evaluated at 2:808 BOOST_MATH_STD_USING809 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 825template <class T>826T zeta_series_2(unsigned digits)827{828 // Series part of zeta at 2:829 BOOST_MATH_STD_USING830 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 846template <class T>847inline T zeta_series_lead_2()848{849 // lead part at 2:850 return 2;851}852 853template <class T>854inline 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 860template <class T>861inline 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 detail869 870template <class T>871template<int N>872inline T constant_glaisher<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))873{874 875 BOOST_MATH_STD_USING876 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.txt888 // 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 910template <class T>911template<int N>912inline T constant_rayleigh_skewness<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))913{ // 1100 digits of the Rayleigh distribution skewness914 // N[2 Sqrt[Pi] (Pi - 3)/((4 - Pi)^(3/2)), 1100]915 916 BOOST_MATH_STD_USING917 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 937template <class T>938template<int N>939inline 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_USING943 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 952template <class T>953template<int N>954inline 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_USING958 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 967template <class T>968template<int N>969inline 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 974template <class T>975template<int N>976inline 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 981template <class T>982template<int N>983inline 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 988template <class T>989template<int N>990inline 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)996template <class T>997template<int N>998inline 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.txt1002 // Also: https://oeis.org/A0068901003 // 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 1010template <class T>1011template<int N>1012inline 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 1019template <class T>1020template<int N>1021inline 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 1037template <class T>1038template<int N>1039inline 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 1055template <class T>1056template<int N>1057inline 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 1076template <class T>1077template<int N>1078inline 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 approximately1082 // 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#endif1103 1104}1105}1106}1107} // namespaces1108 1109#endif // BOOST_MATH_CALCULATE_CONSTANTS_CONSTANTS_INCLUDED1110