//  (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