Boost logo

Boost-Commit :

Subject: [Boost-commit] svn:boost r56503 - in trunk: boost/math/policies boost/math/special_functions boost/math/tools libs/math/test
From: john_at_[hidden]
Date: 2009-10-01 12:03:44


Author: johnmaddock
Date: 2009-10-01 12:03:42 EDT (Thu, 01 Oct 2009)
New Revision: 56503
URL: http://svn.boost.org/trac/boost/changeset/56503

Log:
Another round of performance tweaks for issue #3408.
These should make our igamma implementation comparable in performance to the dcdflib FORTRAN routine - at least as far as MSVC is concerned.
Text files modified:
   trunk/boost/math/policies/policy.hpp | 68 +++++++++++++++++++++++
   trunk/boost/math/special_functions/bessel.hpp | 8 +-
   trunk/boost/math/special_functions/beta.hpp | 24 ++++----
   trunk/boost/math/special_functions/erf.hpp | 4
   trunk/boost/math/special_functions/expint.hpp | 6 +-
   trunk/boost/math/special_functions/expm1.hpp | 4
   trunk/boost/math/special_functions/gamma.hpp | 22 +++---
   trunk/boost/math/special_functions/log1p.hpp | 8 +-
   trunk/boost/math/special_functions/zeta.hpp | 2
   trunk/boost/math/tools/fraction.hpp | 116 +++++++++++++++++----------------------
   trunk/boost/math/tools/series.hpp | 88 ++++++++++++------------------
   trunk/libs/math/test/test_zeta.cpp | 4
   12 files changed, 196 insertions(+), 158 deletions(-)

Modified: trunk/boost/math/policies/policy.hpp
==============================================================================
--- trunk/boost/math/policies/policy.hpp (original)
+++ trunk/boost/math/policies/policy.hpp 2009-10-01 12:03:42 EDT (Thu, 01 Oct 2009)
@@ -34,6 +34,8 @@
 
 template <class T>
 int digits(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T));
+template <class T>
+T epsilon(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T));
 
 }
 
@@ -854,6 +856,72 @@
 
 namespace detail{
 
+template <class T, class Digits, class Small, class Default>
+struct series_factor_calc
+{
+ static const T get()
+ {
+ return ldexp(T(1.0), 1 - Digits::value);
+ }
+};
+
+template <class T, class Digits>
+struct series_factor_calc<T, Digits, mpl::true_, mpl::true_>
+{
+ static const T get()
+ {
+ return boost::math::tools::epsilon<T>();
+ }
+};
+template <class T, class Digits>
+struct series_factor_calc<T, Digits, mpl::true_, mpl::false_>
+{
+ static const T get()
+ {
+ static const boost::uintmax_t v = static_cast<boost::uintmax_t>(1u) << (Digits::value - 1);
+ return 1 / static_cast<T>(v);
+ }
+};
+template <class T, class Digits>
+struct series_factor_calc<T, Digits, mpl::false_, mpl::true_>
+{
+ static const T get()
+ {
+ return boost::math::tools::epsilon<T>();
+ }
+};
+
+template <class T, class Policy>
+inline T get_epsilon_imp(mpl::true_ const&)
+{
+#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
+ BOOST_STATIC_ASSERT( ::std::numeric_limits<T>::is_specialized);
+#else
+ BOOST_ASSERT(::std::numeric_limits<T>::is_specialized);
+#endif
+ typedef typename boost::math::policies::precision<T, Policy>::type p_t;
+ typedef mpl::bool_<p_t::value <= std::numeric_limits<boost::uintmax_t>::digits> is_small_int;
+ typedef mpl::bool_<p_t::value >= std::numeric_limits<T>::digits> is_default_value;
+ return series_factor_calc<T, p_t, is_small_int, is_default_value>::get();
+}
+
+template <class T, class Policy>
+inline T get_epsilon_imp(mpl::false_ const&)
+{
+ return tools::epsilon<T>();
+}
+
+} // namespace detail
+
+template <class T, class Policy>
+inline T get_epsilon(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T))
+{
+ typedef mpl::bool_< std::numeric_limits<T>::is_specialized > tag_type;
+ return detail::get_epsilon_imp<T, Policy>(tag_type());
+}
+
+namespace detail{
+
 template <class A1,
           class A2,
           class A3,

Modified: trunk/boost/math/special_functions/bessel.hpp
==============================================================================
--- trunk/boost/math/special_functions/bessel.hpp (original)
+++ trunk/boost/math/special_functions/bessel.hpp 2009-10-01 12:03:42 EDT (Thu, 01 Oct 2009)
@@ -94,9 +94,9 @@
    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
    T zero = 0;
- T result = boost::math::tools::sum_series(s, boost::math::policies::digits<T, Policy>(), max_iter, zero);
+ T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
 #else
- T result = boost::math::tools::sum_series(s, boost::math::policies::digits<T, Policy>(), max_iter);
+ T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
 #endif
    policies::check_series_iterations("boost::math::bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
    return result;
@@ -110,9 +110,9 @@
    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
    T zero = 0;
- T result = boost::math::tools::sum_series(s, boost::math::policies::digits<T, Policy>(), max_iter, zero);
+ T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
 #else
- T result = boost::math::tools::sum_series(s, boost::math::policies::digits<T, Policy>(), max_iter);
+ T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
 #endif
    policies::check_series_iterations("boost::math::sph_bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
    return result * sqrt(constants::pi<T>() / 4);

Modified: trunk/boost/math/special_functions/beta.hpp
==============================================================================
--- trunk/boost/math/special_functions/beta.hpp (original)
+++ trunk/boost/math/special_functions/beta.hpp 2009-10-01 12:03:42 EDT (Thu, 01 Oct 2009)
@@ -161,11 +161,11 @@
 
    // calculate the fraction parts:
    T sa = detail::lower_gamma_series(a, la, pol) / a;
- sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::digits<T, Policy>());
+ sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::get_epsilon<T, Policy>());
    T sb = detail::lower_gamma_series(b, lb, pol) / b;
- sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::digits<T, Policy>());
+ sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::get_epsilon<T, Policy>());
    T sc = detail::lower_gamma_series(c, lc, pol) / c;
- sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::digits<T, Policy>());
+ sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::get_epsilon<T, Policy>());
 
    // and the exponent part:
    result = exp(lc - la - lb) * pow(la/lc, a) * pow(lb/lc, b);
@@ -398,11 +398,11 @@
    T lc = a + b + 5;
    // gamma function partials:
    T sa = detail::lower_gamma_series(a, la, pol) / a;
- sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::digits<T, Policy>());
+ sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::get_epsilon<T, Policy>());
    T sb = detail::lower_gamma_series(b, lb, pol) / b;
- sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::digits<T, Policy>());
+ sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::get_epsilon<T, Policy>());
    T sc = detail::lower_gamma_series(c, lc, pol) / c;
- sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::digits<T, Policy>());
+ sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::get_epsilon<T, Policy>());
    // gamma function powers combined with incomplete beta powers:
 
    T b1 = (x * lc) / la;
@@ -502,7 +502,7 @@
       return s0; // Safeguard: series can't cope with denorms.
    ibeta_series_t<T> s(a, b, x, result);
    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
- result = boost::math::tools::sum_series(s, boost::math::policies::digits<T, Policy>(), max_iter, s0);
+ result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0);
    policies::check_series_iterations("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (with lanczos)", max_iter, pol);
    return result;
 }
@@ -532,11 +532,11 @@
 
       // calculate the gamma parts:
       T sa = detail::lower_gamma_series(a, la, pol) / a;
- sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::digits<T, Policy>());
+ sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::get_epsilon<T, Policy>());
       T sb = detail::lower_gamma_series(b, lb, pol) / b;
- sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::digits<T, Policy>());
+ sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::get_epsilon<T, Policy>());
       T sc = detail::lower_gamma_series(c, lc, pol) / c;
- sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::digits<T, Policy>());
+ sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::get_epsilon<T, Policy>());
 
       // and their combined power-terms:
       T b1 = (x * lc) / la;
@@ -582,7 +582,7 @@
       return s0; // Safeguard: series can't cope with denorms.
    ibeta_series_t<T> s(a, b, x, result);
    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
- result = boost::math::tools::sum_series(s, boost::math::policies::digits<T, Policy>(), max_iter, s0);
+ result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0);
    policies::check_series_iterations("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (without lanczos)", max_iter, pol);
    return result;
 }
@@ -634,7 +634,7 @@
       return result;
 
    ibeta_fraction2_t<T> f(a, b, x);
- T fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::digits<T, Policy>());
+ T fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon<T, Policy>());
    return result / fract;
 }
 //

Modified: trunk/boost/math/special_functions/erf.hpp
==============================================================================
--- trunk/boost/math/special_functions/erf.hpp (original)
+++ trunk/boost/math/special_functions/erf.hpp 2009-10-01 12:03:42 EDT (Thu, 01 Oct 2009)
@@ -133,7 +133,7 @@
    {
       detail::erf_asympt_series_t<T> s(z);
       boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
- result = boost::math::tools::sum_series(s, policies::digits<T, Policy>(), max_iter, 1);
+ result = boost::math::tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, 1);
       policies::check_series_iterations("boost::math::erf<%1%>(%1%, %1%)", max_iter, pol);
    }
    else
@@ -160,7 +160,7 @@
          invert = !invert;
          result = z * exp(-x);
          result /= sqrt(boost::math::constants::pi<T>());
- result *= upper_gamma_fraction(T(0.5f), x, policies::digits<T, Policy>());
+ result *= upper_gamma_fraction(T(0.5f), x, policies::get_epsilon<T, Policy>());
       }
    }
    if(invert)

Modified: trunk/boost/math/special_functions/expint.hpp
==============================================================================
--- trunk/boost/math/special_functions/expint.hpp (original)
+++ trunk/boost/math/special_functions/expint.hpp 2009-10-01 12:03:42 EDT (Thu, 01 Oct 2009)
@@ -365,7 +365,7 @@
    expint_fraction<T> f(n, z);
    T result = tools::continued_fraction_b(
       f,
- tools::digits<T>(),
+ boost::math::policies::get_epsilon<T, Policy>(),
       max_iter);
    policies::check_series_iterations("boost::math::expint_continued_fraction<%1%>(unsigned,%1%)", max_iter, pol);
    BOOST_MATH_INSTRUMENT_VARIABLE(result)
@@ -422,7 +422,7 @@
    BOOST_MATH_INSTRUMENT_VARIABLE(result)
 
    expint_series<T> s(k, z, x_k, denom, fact);
- result = tools::sum_series(s, policies::digits<T, Policy>(), max_iter, result);
+ result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, result);
    policies::check_series_iterations("boost::math::expint_series<%1%>(unsigned,%1%)", max_iter, pol);
    BOOST_MATH_INSTRUMENT_VARIABLE(result)
    BOOST_MATH_INSTRUMENT_VARIABLE(max_iter)
@@ -495,7 +495,7 @@
    result += constants::euler<T>();
    expint_i_series<T> s(z);
    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
- result = tools::sum_series(s, policies::digits<T, Policy>(), max_iter, result);
+ result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, result);
    policies::check_series_iterations("boost::math::expint_i_series<%1%>(%1%)", max_iter, pol);
    return result;
 }

Modified: trunk/boost/math/special_functions/expm1.hpp
==============================================================================
--- trunk/boost/math/special_functions/expm1.hpp (original)
+++ trunk/boost/math/special_functions/expm1.hpp 2009-10-01 12:03:42 EDT (Thu, 01 Oct 2009)
@@ -90,10 +90,10 @@
    detail::expm1_series<T> s(x);
    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) && !BOOST_WORKAROUND(__EDG_VERSION__, <= 245)
- T result = tools::sum_series(s, policies::digits<T, Policy>(), max_iter);
+ T result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter);
 #else
    T zero = 0;
- T result = tools::sum_series(s, policies::digits<T, Policy>(), max_iter, zero);
+ T result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, zero);
 #endif
    policies::check_series_iterations("boost::math::expm1<%1%>(%1%)", max_iter, pol);
    return result;

Modified: trunk/boost/math/special_functions/gamma.hpp
==============================================================================
--- trunk/boost/math/special_functions/gamma.hpp (original)
+++ trunk/boost/math/special_functions/gamma.hpp 2009-10-01 12:03:42 EDT (Thu, 01 Oct 2009)
@@ -296,13 +296,13 @@
 };
 
 template <class T>
-inline T upper_gamma_fraction(T a, T z, int bits)
+inline T upper_gamma_fraction(T a, T z, T eps)
 {
    // Multiply result by z^a * e^-z to get the full
    // upper incomplete integral. Divide by tgamma(z)
    // to normalise.
    upper_incomplete_gamma_fract<T> f(a, z);
- return 1 / (z - a + 1 + boost::math::tools::continued_fraction_a(f, bits));
+ return 1 / (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps));
 }
 
 template <class T>
@@ -331,8 +331,8 @@
    // to get the normalised value.
    lower_incomplete_gamma_series<T> s(a, z);
    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
- int bits = policies::digits<T, Policy>();
- T result = boost::math::tools::sum_series(s, bits, max_iter, init_value);
+ T factor = policies::get_epsilon<T, Policy>();
+ T result = boost::math::tools::sum_series(s, factor, max_iter, init_value);
    policies::check_series_iterations("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
    return result;
 }
@@ -382,7 +382,7 @@
       BOOST_MATH_INSTRUMENT_CODE(prefix);
       T sum = detail::lower_gamma_series(z, z, pol) / z;
       BOOST_MATH_INSTRUMENT_CODE(sum);
- sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::digits<T, Policy>());
+ sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::get_epsilon<T, Policy>());
       BOOST_MATH_INSTRUMENT_CODE(sum);
       if(fabs(tools::max_value<T>() / prefix) < fabs(sum))
          return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
@@ -421,7 +421,7 @@
       T limit = (std::max)(z+1, T(10));
       T prefix = z * log(limit) - limit;
       T sum = detail::lower_gamma_series(z, limit, pol) / z;
- sum += detail::upper_gamma_fraction(z, limit, ::boost::math::policies::digits<T, Policy>());
+ sum += detail::upper_gamma_fraction(z, limit, ::boost::math::policies::get_epsilon<T, Policy>());
       result = log(sum) + prefix;
    }
    if(sign)
@@ -686,7 +686,7 @@
 
    T limit = (std::max)(T(10), a);
    T sum = detail::lower_gamma_series(a, limit, pol) / a;
- sum += detail::upper_gamma_fraction(a, limit, ::boost::math::policies::digits<T, Policy>());
+ sum += detail::upper_gamma_fraction(a, limit, ::boost::math::policies::get_epsilon<T, Policy>());
 
    if(a < 10)
    {
@@ -747,7 +747,7 @@
    if(pderivative)
       *pderivative = p / (*pgam * exp(x));
    T init_value = invert ? *pgam : 0;
- result = -p * tools::sum_series(s, boost::math::policies::digits<T, Policy>(), max_iter, (init_value - result) / p);
+ result = -p * tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, (init_value - result) / p);
    policies::check_series_iterations("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
    if(invert)
       result = -result;
@@ -1011,7 +1011,7 @@
          if(p_derivative)
             *p_derivative = result;
          if(result != 0)
- result *= upper_gamma_fraction(a, x, policies::digits<T, Policy>());
+ result *= upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>());
          break;
       }
    case 5:
@@ -1126,9 +1126,9 @@
    }
    prefix *= pow(constants::e<T>() / zd, delta);
    T sum = detail::lower_gamma_series(z, z, pol) / z;
- sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::digits<T, Policy>());
+ sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::get_epsilon<T, Policy>());
    T sumd = detail::lower_gamma_series(zd, zd, pol) / zd;
- sumd += detail::upper_gamma_fraction(zd, zd, ::boost::math::policies::digits<T, Policy>());
+ sumd += detail::upper_gamma_fraction(zd, zd, ::boost::math::policies::get_epsilon<T, Policy>());
    sum /= sumd;
    if(fabs(tools::max_value<T>() / prefix) < fabs(sum))
       return policies::raise_overflow_error<T>("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Result of tgamma is too large to represent.", pol);

Modified: trunk/boost/math/special_functions/log1p.hpp
==============================================================================
--- trunk/boost/math/special_functions/log1p.hpp (original)
+++ trunk/boost/math/special_functions/log1p.hpp 2009-10-01 12:03:42 EDT (Thu, 01 Oct 2009)
@@ -95,10 +95,10 @@
    detail::log1p_series<result_type> s(x);
    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) && !BOOST_WORKAROUND(__EDG_VERSION__, <= 245)
- result_type result = tools::sum_series(s, policies::digits<result_type, Policy>(), max_iter);
+ result_type result = tools::sum_series(s, policies::get_epsilon<result_type, Policy>(), max_iter);
 #else
    result_type zero = 0;
- result_type result = tools::sum_series(s, policies::digits<result_type, Policy>(), max_iter, zero);
+ result_type result = tools::sum_series(s, policies::get_epsilon<result_type, Policy>(), max_iter, zero);
 #endif
    policies::check_series_iterations(function, max_iter, pol);
    return result;
@@ -443,9 +443,9 @@
    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
    T zero = 0;
- T result = boost::math::tools::sum_series(s, policies::digits<T, Policy>(), max_iter, zero);
+ T result = boost::math::tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, zero);
 #else
- T result = boost::math::tools::sum_series(s, policies::digits<T, Policy>(), max_iter);
+ T result = boost::math::tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter);
 #endif
    policies::check_series_iterations(function, max_iter, pol);
    return result;

Modified: trunk/boost/math/special_functions/zeta.hpp
==============================================================================
--- trunk/boost/math/special_functions/zeta.hpp (original)
+++ trunk/boost/math/special_functions/zeta.hpp 2009-10-01 12:03:42 EDT (Thu, 01 Oct 2009)
@@ -113,7 +113,7 @@
    zeta_series2<T> f(s);
    T result = tools::sum_series(
       f,
- tools::digits<T>(),
+ policies::get_epsilon<T, Policy>(),
       max_iter);
    policies::check_series_iterations("boost::math::zeta_series2<%1%>(%1%)", max_iter, pol);
    return result;

Modified: trunk/boost/math/tools/fraction.hpp
==============================================================================
--- trunk/boost/math/tools/fraction.hpp (original)
+++ trunk/boost/math/tools/fraction.hpp 2009-10-01 12:03:42 EDT (Thu, 01 Oct 2009)
@@ -84,8 +84,8 @@
 //
 // Note that the first a0 returned by generator Gen is disarded.
 //
-template <class Gen>
-typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits)
+template <class Gen, class U>
+inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor, boost::uintmax_t& max_terms)
 {
    BOOST_MATH_STD_USING // ADL of std names
 
@@ -93,7 +93,6 @@
    typedef typename traits::result_type result_type;
    typedef typename traits::value_type value_type;
 
- result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
    result_type tiny = tools::min_value<result_type>();
 
    value_type v = g();
@@ -105,6 +104,8 @@
    C = f;
    D = 0;
 
+ boost::uintmax_t counter(max_terms);
+
    do{
       v = g();
       D = traits::b(v) + traits::a(v) * D;
@@ -116,50 +117,43 @@
       D = 1/D;
       delta = C*D;
       f = f * delta;
- }while(fabs(delta - 1) > factor);
+ }while((fabs(delta - 1) > factor) && --counter);
+
+ max_terms = max_terms - counter;
 
    return f;
 }
 
+template <class Gen, class U>
+inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor)
+{
+ boost::uintmax_t max_terms = (std::numeric_limits<boost::uintmax_t>::max)();
+ return continued_fraction_b(g, factor, max_terms);
+}
+
 template <class Gen>
-typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits, boost::uintmax_t& max_terms)
+inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits)
 {
    BOOST_MATH_STD_USING // ADL of std names
 
    typedef detail::fraction_traits<Gen> traits;
    typedef typename traits::result_type result_type;
- typedef typename traits::value_type value_type;
 
    result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
- result_type tiny = tools::min_value<result_type>();
-
- value_type v = g();
-
- result_type f, C, D, delta;
- f = traits::b(v);
- if(f == 0)
- f = tiny;
- C = f;
- D = 0;
-
- boost::uintmax_t counter(max_terms);
+ boost::uintmax_t max_terms = (std::numeric_limits<boost::uintmax_t>::max)();
+ return continued_fraction_b(g, factor, max_terms);
+}
 
- do{
- v = g();
- D = traits::b(v) + traits::a(v) * D;
- if(D == 0)
- D = tiny;
- C = traits::b(v) + traits::a(v) / C;
- if(C == 0)
- C = tiny;
- D = 1/D;
- delta = C*D;
- f = f * delta;
- }while((fabs(delta - 1) > factor) && --counter);
+template <class Gen>
+inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits, boost::uintmax_t& max_terms)
+{
+ BOOST_MATH_STD_USING // ADL of std names
 
- max_terms = max_terms - counter;
+ typedef detail::fraction_traits<Gen> traits;
+ typedef typename traits::result_type result_type;
 
- return f;
+ result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
+ return continued_fraction_b(g, factor, max_terms);
 }
 
 //
@@ -176,8 +170,8 @@
 //
 // Note that the first a1 and b1 returned by generator Gen are both used.
 //
-template <class Gen>
-typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits)
+template <class Gen, class U>
+inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor, boost::uintmax_t& max_terms)
 {
    BOOST_MATH_STD_USING // ADL of std names
 
@@ -185,7 +179,6 @@
    typedef typename traits::result_type result_type;
    typedef typename traits::value_type value_type;
 
- result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
    result_type tiny = tools::min_value<result_type>();
 
    value_type v = g();
@@ -198,6 +191,8 @@
    C = f;
    D = 0;
 
+ boost::uintmax_t counter(max_terms);
+
    do{
       v = g();
       D = traits::b(v) + traits::a(v) * D;
@@ -209,51 +204,44 @@
       D = 1/D;
       delta = C*D;
       f = f * delta;
- }while(fabs(delta - 1) > factor);
+ }while((fabs(delta - 1) > factor) && --counter);
+
+ max_terms = max_terms - counter;
 
    return a0/f;
 }
 
+template <class Gen, class U>
+inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor)
+{
+ boost::uintmax_t max_iter = (std::numeric_limits<boost::uintmax_t>::max)();
+ return continued_fraction_a(g, factor, max_iter);
+}
+
 template <class Gen>
-typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits, boost::uintmax_t& max_terms)
+inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits)
 {
    BOOST_MATH_STD_USING // ADL of std names
 
    typedef detail::fraction_traits<Gen> traits;
    typedef typename traits::result_type result_type;
- typedef typename traits::value_type value_type;
 
    result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
- result_type tiny = tools::min_value<result_type>();
-
- value_type v = g();
-
- result_type f, C, D, delta, a0;
- f = traits::b(v);
- a0 = traits::a(v);
- if(f == 0)
- f = tiny;
- C = f;
- D = 0;
+ boost::uintmax_t max_iter = (std::numeric_limits<boost::uintmax_t>::max)();
 
- boost::uintmax_t counter(max_terms);
+ return continued_fraction_a(g, factor, max_iter);
+}
 
- do{
- v = g();
- D = traits::b(v) + traits::a(v) * D;
- if(D == 0)
- D = tiny;
- C = traits::b(v) + traits::a(v) / C;
- if(C == 0)
- C = tiny;
- D = 1/D;
- delta = C*D;
- f = f * delta;
- }while((fabs(delta - 1) > factor) && --counter);
+template <class Gen>
+inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits, boost::uintmax_t& max_terms)
+{
+ BOOST_MATH_STD_USING // ADL of std names
 
- max_terms = max_terms - counter;
+ typedef detail::fraction_traits<Gen> traits;
+ typedef typename traits::result_type result_type;
 
- return a0/f;
+ result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
+ return continued_fraction_a(g, factor, max_terms);
 }
 
 } // namespace tools

Modified: trunk/boost/math/tools/series.hpp
==============================================================================
--- trunk/boost/math/tools/series.hpp (original)
+++ trunk/boost/math/tools/series.hpp 2009-10-01 12:03:42 EDT (Thu, 01 Oct 2009)
@@ -12,6 +12,7 @@
 
 #include <boost/config/no_tr1/cmath.hpp>
 #include <boost/cstdint.hpp>
+#include <boost/limits.hpp>
 #include <boost/math/tools/config.hpp>
 
 namespace boost{ namespace math{ namespace tools{
@@ -19,26 +20,8 @@
 //
 // Simple series summation come first:
 //
-template <class Functor>
-inline typename Functor::result_type sum_series(Functor& func, int bits)
-{
- BOOST_MATH_STD_USING
-
- typedef typename Functor::result_type result_type;
-
- result_type factor = pow(result_type(2), bits);
- result_type result = func();
- result_type next_term;
- do{
- next_term = func();
- result += next_term;
- }
- while(fabs(result) < fabs(factor * next_term));
- return result;
-}
-
-template <class Functor>
-inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms)
+template <class Functor, class U, class V>
+inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms, const V& init_value)
 {
    BOOST_MATH_STD_USING
 
@@ -46,14 +29,13 @@
 
    boost::uintmax_t counter = max_terms;
 
- result_type factor = ldexp(result_type(1), bits);
- result_type result = func();
+ result_type result = init_value;
    result_type next_term;
    do{
       next_term = func();
       result += next_term;
    }
- while((fabs(result) < fabs(factor * next_term)) && --counter);
+ while((fabs(factor * result) < fabs(next_term)) && --counter);
 
    // set max_terms to the actual number of terms of the series evaluated:
    max_terms = max_terms - counter;
@@ -62,46 +44,46 @@
 }
 
 template <class Functor, class U>
-inline typename Functor::result_type sum_series(Functor& func, int bits, U init_value)
+inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms)
 {
- BOOST_MATH_STD_USING
-
- typedef typename Functor::result_type result_type;
-
- result_type factor = ldexp(result_type(1), bits);
- result_type result = static_cast<result_type>(init_value);
- result_type next_term;
- do{
- next_term = func();
- result += next_term;
- }
- while(fabs(result) < fabs(factor * next_term));
-
- return result;
+ typename Functor::result_type init_value = 0;
+ return sum_series(func, factor, max_terms, init_value);
 }
 
 template <class Functor, class U>
-inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms, U init_value)
+inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms, const U& init_value)
 {
    BOOST_MATH_STD_USING
-
    typedef typename Functor::result_type result_type;
+ result_type factor = ldexp(result_type(1), 1 - bits);
+ return sum_series(func, factor, max_terms, init_value);
+}
 
- boost::uintmax_t counter = max_terms;
-
- result_type factor = ldexp(result_type(1), bits);
- result_type result = init_value;
- result_type next_term;
- do{
- next_term = func();
- result += next_term;
- }
- while((fabs(result) < fabs(factor * next_term)) && --counter);
+template <class Functor>
+inline typename Functor::result_type sum_series(Functor& func, int bits)
+{
+ BOOST_MATH_STD_USING
+ typedef typename Functor::result_type result_type;
+ boost::uintmax_t iters = (std::numeric_limits<boost::uintmax_t>::max)();
+ result_type init_val = 0;
+ return sum_series(func, bits, iters, init_val);
+}
 
- // set max_terms to the actual number of terms of the series evaluated:
- max_terms = max_terms - counter;
+template <class Functor>
+inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms)
+{
+ BOOST_MATH_STD_USING
+ typedef typename Functor::result_type result_type;
+ result_type init_val = 0;
+ return sum_series(func, bits, max_terms, init_val);
+}
 
- return result;
+template <class Functor, class U>
+inline typename Functor::result_type sum_series(Functor& func, int bits, const U& init_value)
+{
+ BOOST_MATH_STD_USING
+ boost::uintmax_t iters = (std::numeric_limits<boost::uintmax_t>::max)();
+ return sum_series(func, bits, iters, init_value);
 }
 
 //

Modified: trunk/libs/math/test/test_zeta.cpp
==============================================================================
--- trunk/libs/math/test/test_zeta.cpp (original)
+++ trunk/libs/math/test/test_zeta.cpp 2009-10-01 12:03:42 EDT (Thu, 01 Oct 2009)
@@ -88,7 +88,7 @@
       ".*", // platform
       "real_concept", // test type(s)
       ".*", // test data group
- ".*", 10, 5); // test function
+ ".*", 16, 5); // test function
 
    std::cout << "Tests run with " << BOOST_COMPILER << ", "
       << BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl;
@@ -176,7 +176,7 @@
    BOOST_CHECK_CLOSE(::boost::math::zeta(static_cast<T>(4.5)), static_cast<T>(1.05470751076145426402296728896028011727249383295625173068468L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::zeta(static_cast<T>(6.5)), static_cast<T>(1.01200589988852479610078491680478352908773213619144808841031L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::zeta(static_cast<T>(7.5)), static_cast<T>(1.00582672753652280770224164440459408011782510096320822989663L), tolerance);
- BOOST_CHECK_CLOSE(::boost::math::zeta(static_cast<T>(8.125)), static_cast<T>(1.0037305205308161603183307711439385250181080293472L), tolerance);
+ BOOST_CHECK_CLOSE(::boost::math::zeta(static_cast<T>(8.125)), static_cast<T>(1.0037305205308161603183307711439385250181080293472L), 2 * tolerance);
    BOOST_CHECK_CLOSE(::boost::math::zeta(static_cast<T>(16.125)), static_cast<T>(1.0000140128224754088474783648500235958510030511915L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::zeta(static_cast<T>(0)), static_cast<T>(-0.5L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::zeta(static_cast<T>(-0.125)), static_cast<T>(-0.39906966894504503550986928301421235400280637468895L), tolerance);


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