[Boost-commit] svn:boost r76333 - in sandbox/math_constants: boost/math/constants libs/math/test

Subject: [Boost-commit] svn:boost r76333 - in sandbox/math_constants: boost/math/constants libs/math/test
From: john_at_[hidden]
Date: 2012-01-06 08:17:49

Author: johnmaddock
Date: 2012-01-06 08:17:48 EST (Fri, 06 Jan 2012)
New Revision: 76333
URL: http://svn.boost.org/trac/boost/changeset/76333

Log:
Add calculation for Glaisher's constant.  Plus fix a buglet in constants.hpp.
Text files modified: 
   sandbox/math_constants/boost/math/constants/calculate_constants.hpp |    96 ++++++++++++++++++++++++++++++++++++++++
   sandbox/math_constants/boost/math/constants/constants.hpp           |     4                                         
   sandbox/math_constants/libs/math/test/test_constant_generate.cpp    |     4 -                                       
   3 files changed, 99 insertions(+), 5 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	2012-01-06 08:17:48 EST (Fri, 06 Jan 2012)
@@ -8,6 +8,8 @@
 #ifndef BOOST_MATH_CALCULATE_CONSTANTS_CONSTANTS_INCLUDED
 #define BOOST_MATH_CALCULATE_CONSTANTS_CONSTANTS_INCLUDED
 
+#include <boost/math/special_functions/trunc.hpp>
+
 namespace boost{ namespace math{ namespace constants{ namespace detail{
 
 template <class T>
@@ -800,6 +802,88 @@
   return ev;
 }
 
+namespace detail{
+//
+// Calculation of the Glaisher constant depends upon calculating the
+// derivative of the zeta function at 2, we can then use the relation:
+// zeta'(2) = 1/6 pi^2 [euler + ln(2pi)-12ln(A)]
+// To get the constant A.
+// See equation 45 at http://mathworld.wolfram.com/RiemannZetaFunction.html.
+//
+// The derivative of the zeta function is computed by direct differentiation
+// of the relation:
+// (1-2^(1-s))zeta(s) = SUM(n=0, INF){ (-n)^n / (n+1)^s  } 
+// Which gives us 2 slowly converging but alternating sums to compute,
+// for this we use Algorithm 1 from "Convergent Acceleration of Alternating Series",
+// Henri Cohen, Fernando Rodriguez Villegas and Don Zagier, Experimental Mathematics 9:1.
+// See http://www.math.utexas.edu/users/villegas/publications/conv-accel.pdf
+//
+template <class T>
+T zeta_series_derivative_2()
+{
+   // Derivative of the series part, evaluated at 2:
+   BOOST_MATH_STD_USING
+   int n = boost::math::itrunc((std::numeric_limits<T>::digits10 + 1) * 1.3);
+   T d = pow(3 + sqrt(T(8)), n);
+   d = (d + 1 / d) / 2;
+   T b = -1;
+   T c = -d;
+   T s = 0;
+   for(int k = 0; k < n; ++k)
+   {
+      T a = -log(T(k+1)) / ((k+1) * (k+1));
+      c = b - c;
+      s = s + c * a;
+      b = (k + n) * (k - n) * b / ((k + T(0.5f)) * (k + 1));
+   }
+   return s / d;
+}
+
+template <class T>
+T zeta_series_2()
+{
+   // Series part of zeta at 2:
+   BOOST_MATH_STD_USING
+   int n = boost::math::itrunc((std::numeric_limits<T>::digits10 + 1) * 1.3);
+   T d = pow(3 + sqrt(T(8)), n);
+   d = (d + 1 / d) / 2;
+   T b = -1;
+   T c = -d;
+   T s = 0;
+   for(int k = 0; k < n; ++k)
+   {
+      T a = T(1) / ((k + 1) * (k + 1));
+      c = b - c;
+      s = s + c * a;
+      b = (k + n) * (k - n) * b / ((k + T(0.5f)) * (k + 1));
+   }
+   return s / d;
+}
+
+template <class T>
+inline T zeta_series_lead_2()
+{
+   // lead part at 2:
+   return 2;
+}
+
+template <class T>
+inline T zeta_series_derivative_lead_2()
+{
+   // derivative of lead part at 2:
+   return -2 * boost::math::constants::ln_two<T>();
+}
+
+template <class T>
+inline T zeta_derivative_2()
+{
+   // zeta derivative at 2:
+   return zeta_series_derivative_2<T>() * zeta_series_lead_2<T>()
+      + zeta_series_derivative_lead_2<T>() * zeta_series_2<T>();
+}
+
+}  // detail
+
 template <class T>
 template<int N>
 inline T constant_glaisher<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>))
@@ -810,6 +894,17 @@
   // No good references/algorithms for this found yet.
   // TODO calculate this?
 
+   BOOST_MATH_STD_USING
+   typedef policies::policy<policies::digits2<N> > forwarding_policy;
+   T v = detail::zeta_derivative_2<T>();
+   v *= 6;
+   v /= boost::math::constants::pi<T, forwarding_policy>() * boost::math::constants::pi<T, forwarding_policy>();
+   v -= boost::math::constants::euler<T, forwarding_policy>();
+   v -= log(2 * boost::math::constants::pi<T, forwarding_policy>());
+   v /= -12;
+   return exp(v);
+
+   /*
   T g(
 "1.282427129100622636875342568869791727767688927325001192063740021740406308858826"
 "46112973649195820237439420646120399000748933157791362775280404159072573861727522"
@@ -826,6 +921,7 @@
 "69575162841550648585890834128227556209547002918593263079373376942077522290940187");
 
   return g;
+  */
 }
 
 template <class T>

Modified: sandbox/math_constants/boost/math/constants/constants.hpp
==============================================================================
--- sandbox/math_constants/boost/math/constants/constants.hpp	(original)
+++ sandbox/math_constants/boost/math/constants/constants.hpp	2012-01-06 08:17:48 EST (Fri, 06 Jan 2012)
@@ -165,7 +165,7 @@
    /* The default implementations come next: */ \
    static inline T get_from_string()\
    {\
-      static const T result = detail::convert_from_string<T>(y, boost::is_convertible<const char*, T>());\
+      static const T result = convert_from_string<T>(y, boost::is_convertible<const char*, T>());\
       return result;\
    }\
    /* This one is for very high precision that is none the less known at compile time: */ \
@@ -208,7 +208,7 @@
    namespace long_double_constants{ static const long double name = BOOST_JOIN(x, L); }\
    namespace constants{
 
- BOOST_DEFINE_MATH_CONSTANT(half, 5.000000000000000000000000000000000000e-01, "5.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-01");
+  BOOST_DEFINE_MATH_CONSTANT(half, 5.000000000000000000000000000000000000e-01, "5.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-01");
   BOOST_DEFINE_MATH_CONSTANT(third, 3.333333333333333333333333333333333333e-01, "3.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333e-01");
   BOOST_DEFINE_MATH_CONSTANT(twothirds, 6.666666666666666666666666666666666666e-01, "6.66666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667e-01");
   BOOST_DEFINE_MATH_CONSTANT(two_thirds, 6.666666666666666666666666666666666666e-01, "6.66666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667e-01");

Modified: sandbox/math_constants/libs/math/test/test_constant_generate.cpp
==============================================================================
--- sandbox/math_constants/libs/math/test/test_constant_generate.cpp	(original)
+++ sandbox/math_constants/libs/math/test/test_constant_generate.cpp	2012-01-06 08:17:48 EST (Fri, 06 Jan 2012)
@@ -13,8 +13,6 @@
 #  pragma warning (disable : 4996) //To disable this warning, use -D_SCL_SECURE_NO_WARNINGS.
 #endif
 
-#define USE_CPP_FLOAT
-
 // To add new constants, add a function that calculates the value of the constant to
 // boost/math/constants/calculate_constants.hpp
 
@@ -103,8 +101,8 @@
    BOOST_CONSTANTS_GENERATE(zeta_two);
    BOOST_CONSTANTS_GENERATE(zeta_three);
    BOOST_CONSTANTS_GENERATE(catalan);
-   BOOST_CONSTANTS_GENERATE(khinchin);
    BOOST_CONSTANTS_GENERATE(glaisher);
+   BOOST_CONSTANTS_GENERATE(khinchin);
    BOOST_CONSTANTS_GENERATE(extreme_value_skewness);
    BOOST_CONSTANTS_GENERATE(rayleigh_skewness);
    BOOST_CONSTANTS_GENERATE(rayleigh_kurtosis);