[Boost-commit] svn:boost r75965 - sandbox/math_constants/boost/math/constants

Subject: [Boost-commit] svn:boost r75965 - sandbox/math_constants/boost/math/constants
From: pbristow_at_[hidden]
Date: 2011-12-15 12:08:02

Author: pbristow
Date: 2011-12-15 12:08:01 EST (Thu, 15 Dec 2011)
New Revision: 75965
URL: http://svn.boost.org/trac/boost/changeset/75965

Log:
added calculation for catalan and zeta3(3) Apery.  more TODO
Text files modified: 
   sandbox/math_constants/boost/math/constants/calculate_constants.hpp |   102 ++++++++++++++++++++++++++++++++++++--- 
   1 files changed, 92 insertions(+), 10 deletions(-)

Modified: sandbox/math_constants/boost/math/constants/calculate_constants.hpp
==============================================================================
--- sandbox/math_constants/boost/math/constants/calculate_constants.hpp	(original)
+++ sandbox/math_constants/boost/math/constants/calculate_constants.hpp	2011-12-15 12:08:01 EST (Thu, 15 Dec 2011)
@@ -498,6 +498,8 @@
 );
  //    T gamma("0.57721566490153286060651209008240243104215933593992");
      return gamma;
+       // TODO calculate this - JM has code.
+
 }
 
 template <class T, int N> 
@@ -523,6 +525,7 @@
 template <class T, int N> 
 inline T calculate_cf10(const mpl::int_<N>&BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T))
 {
+    // TODO calculate this.
    BOOST_MATH_STD_USING
      T cf10("1.030640834100712935881776094116936840925");
      return cf10;
@@ -543,26 +546,101 @@
 template <class T, int N> 
 inline T calculate_zeta_three(const mpl::int_<N>&BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T))
 {
+   // http://mathworld.wolfram.com/AperysConstant.html
+   // http://en.wikipedia.org/wiki/Mathematical_constant
+
    // http://oeis.org/A002117/constant
-    T zeta3("1.20205690315959428539973816151144999076"
-     "4986292340498881792271555341838205786313"
-     "09018645587360933525814619915");
-     return zeta3;
+   //T zeta3("1.20205690315959428539973816151144999076"
+   // "4986292340498881792271555341838205786313"
+   // "09018645587360933525814619915");
+
+  //"1.202056903159594285399738161511449990, 76498629234049888179227155534183820578631309018645587360933525814619915"  A002117
+  // 1.202056903159594285399738161511449990, 76498629234049888179227155534183820578631309018645587360933525814619915780, +00);
+  //"1.2020569031595942 double
+  // http://www.spaennare.se/SSPROG/ssnum.pdf  // section 11, Algorithmfor Apery’s constant zeta(3).
+  // Programs to Calculate some Mathematical Constants to Large Precision, Document Version 1.50
+ 
+  // by Stefan Spannare  September 19, 2007
+  // zeta(3) = 1/64 * sum 
+   T n_fact=static_cast<T>(1); // build n! for n = 0.
+   T sum = static_cast<double>(77); // Start with n = 0 case.
+   // for n = 0, (77/1) /64 = 1.203125
+   //double lim = std::numeric_limits<double>::epsilon();
+   T lim = boost::math::tools::epsilon<T>();
+   for(unsigned int n = 1; n < 40; ++n)
+   { // three to five decimal digits per term, so 40 should be plenty for 100 decimal digits.
+      //cout << "n = " << n << endl;
+      n_fact *= n; // n!
+      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
+      T num = ((205 * n * n) + (250 * n) + 77) * n_fact_p10; // 205n^2 + 250n + 77
+      // int nn = (2 * n + 1);
+      // T d = factorial(nn); // inline factorial.
+      T d = 1;
+      for(unsigned int i = 1; i <= (n+n + 1); ++i) // (2n + 1)
+      {
+        d *= i;
+      }
+      T den = d * d * d * d * d; // [(2n+1)!]^5
+      //cout << "den = " << den << endl;
+      T term = num/den;
+      if (n % 2 != 0)
+      { //term *= -1;
+        sum -= term;
+      }
+      else
+      {
+        sum += term;
+      }
+      //cout << "term = " << term << endl;
+      //cout << "sum/64  = " << sum/64 << endl;
+      if(abs(term) < lim)
+      {
+         break;
+      }
+   }
+   return sum / 64;
 }
 
 template <class T, int N> 
 inline T calculate_catalan(const mpl::int_<N>&BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T))
 { // http://oeis.org/A006752/constant
-  T c("0.915965594177219015054603514932384110774"
- "1493742816721342664981196217630197762547"
- "69479356512926115106248574");
- return c;
+  //T c("0.915965594177219015054603514932384110774"
+ //"149374281672134266498119621763019776254769479356512926115106248574");
+
+  // 9.159655941772190150546035149323841107, 74149374281672134266498119621763019776254769479356512926115106248574422619, -01);
+
+   // This is equation (entry) 31 from
+   // http://www-2.cs.cmu.edu/~adamchik/articles/catalan/catalan.htm
+   // See also http://www.mpfr.org/algorithms.pdf
+   T k_fact = 1;
+   T tk_fact = 1;
+   T sum = 1;
+   T term;
+   T lim = boost::math::tools::epsilon<T>();
+
+   for(unsigned k = 1;; ++k)
+   {
+      k_fact *= k;
+      tk_fact *= (2 * k) * (2 * k - 1);
+      term = k_fact * k_fact / (tk_fact * (2 * k + 1) * (2 * k + 1));
+      sum += term;
+      if(term < lim)
+      {
+         break;
+      }
+   }
+   return boost::math::constants::pi<T, boost::math::policies::policy<> >()
+      * log(2 + boost::math::constants::root_three<T, boost::math::policies::policy<> >())
+       / 8
+      + 3 * sum / 8;
 }
 
 
 template <class T, int N> 
 inline T calculate_khinchin(const mpl::int_<N>&BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T))
 { // from e_float constants.cpp
+  // http://mathworld.wolfram.com/KhinchinsConstant.html
+  // TODO calculate this.
   T k( "2.6854520010653064453097148354817956938203822939944629530511523455572188595371520028011411749318476979"
 "9515346590528809008289767771641096305179253348325966838185231542133211949962603932852204481940961806"
 "8664166428930847788062036073705350103367263357728904990427070272345170262523702354581068631850103237"
@@ -581,6 +659,7 @@
 inline T calculate_extreme_value_skewness(const mpl::int_<N>&BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T))
 { // from e_float constants.cpp
   // Mathematica: N[12 Sqrt[6]  Zeta[3]/Pi^3, 1101]
+  // TODO Calculate this.
 
 T ev(
 "1.1395470994046486574927930193898461120875997958365518247216557100852480077060706857071875468869385150"
@@ -602,8 +681,9 @@
 inline T calculate_glaisher(const mpl::int_<N>&BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T))
 { // from http://mpmath.googlecode.com/svn/data/glaisher.txt
   // 20,000 digits of the Glaisher-Kinkelin constant A = exp(1/2 - zeta'(-1))
-  // Computed using A = exp((6 (-zeta'(2))/pi^2 + log 2 pi + gamma)/12)
+  // Computed using A = exp((6 (-zeta(2))/pi^2 + log 2 pi + gamma)/12)
   // with Euler-Maclaurin summation for zeta'(2).
+  // TODO calculate this.
 
   T g(
 "1.282427129100622636875342568869791727767688927325001192063740021740406308858826"
@@ -628,7 +708,7 @@
 {  // From e_float
   // 1100 digits of the Rayleigh distribution skewness
   // Mathematica: N[2 Sqrt[Pi] (Pi - 3)/((4 - Pi)^(3/2)), 1100]
-
+  // TODO Calculate this using pi_minus_three etc.
 
 T rs(  
   "0.6311106578189371381918993515442277798440422031347194976580945856929268196174737254599050270325373067"
@@ -648,6 +728,7 @@
 template <class T, int N> 
 inline T calculate_rayleigh_kurtosis_excess(const mpl::int_<N>&BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T))
 { // - (6 Pi^2 - 24 Pi + 16)/((Pi - 4)^2)
+    // TODO Calculate this using pi_minus_four etc.
    return - (((static_cast<T>(6) * pi<T, policies::policy<policies::digits2<N> > >()
         * pi<T, policies::policy<policies::digits2<N> > >())
    - (static_cast<T>(24) * pi<T, policies::policy<policies::digits2<N> > >()) + static_cast<T>(16) )
@@ -660,6 +741,7 @@
 template <class T, int N> 
 inline T calculate_rayleigh_kurtosis(const mpl::int_<N>&BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T))
 { // 3 - (6 Pi^2 - 24 Pi + 16)/((Pi - 4)^2)
+   // TODO Calculate this using pi_minus_four etc.
    return static_cast<T>(3) - (((static_cast<T>(6) * pi<T, policies::policy<policies::digits2<N> > >()
         * pi<T, policies::policy<policies::digits2<N> > >())
    - (static_cast<T>(24) * pi<T, policies::policy<policies::digits2<N> > >()) + static_cast<T>(16) )