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