[Boost-commit] svn:boost r56503 - in trunk: boost/math/policies boost/math/special_functions boost/math/tools libs/math/test

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);