|
Boost-Commit : |
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);
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk