williamr@2: //  (C) Copyright John Maddock 2005.
williamr@2: //  Use, modification and distribution are subject to the
williamr@2: //  Boost Software License, Version 1.0. (See accompanying file
williamr@2: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
williamr@2: 
williamr@2: #ifndef BOOST_MATH_SERIES_INCLUDED
williamr@2: #define BOOST_MATH_SERIES_INCLUDED
williamr@2: 
williamr@2: #include <cmath>
williamr@2: 
williamr@2: #ifdef BOOST_NO_STDC_NAMESPACE
williamr@2: namespace std{ using ::pow; using ::fabs; }
williamr@2: #endif
williamr@2: 
williamr@2: 
williamr@2: namespace boost{ namespace math{ namespace detail{
williamr@2: 
williamr@2: //
williamr@2: // Algorithm kahan_sum_series invokes Functor func until the N'th
williamr@2: // term is too small to have any effect on the total, the terms 
williamr@2: // are added using the Kahan summation method.
williamr@2: //
williamr@2: // CAUTION: Optimizing compilers combined with extended-precision
williamr@2: // machine registers conspire to render this algorithm partly broken:
williamr@2: // double rounding of intermediate terms (first to a long double machine
williamr@2: // register, and then to a double result) cause the rounding error computed
williamr@2: // by the algorithm to be off by up to 1ulp.  However this occurs rarely, and 
williamr@2: // in any case the result is still much better than a naive summation.
williamr@2: //
williamr@2: template <class Functor>
williamr@2: typename Functor::result_type kahan_sum_series(Functor& func, int bits)
williamr@2: {
williamr@2:    typedef typename Functor::result_type result_type;
williamr@2:    result_type factor = std::pow(result_type(2), bits);
williamr@2:    result_type result = func();
williamr@2:    result_type next_term, y, t;
williamr@2:    result_type carry = 0;
williamr@2:    do{
williamr@2:       next_term = func();
williamr@2:       y = next_term - carry;
williamr@2:       t = result + y;
williamr@2:       carry = t - result;
williamr@2:       carry -= y;
williamr@2:       result = t;
williamr@2:    }
williamr@2:    while(std::fabs(result) < std::fabs(factor * next_term));
williamr@2:    return result;
williamr@2: }
williamr@2: 
williamr@2: } } } // namespaces
williamr@2: 
williamr@2: #endif // BOOST_MATH_SERIES_INCLUDED