Chris@16: //  (C) Copyright John Maddock 2005-2006.
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_TOOLS_SERIES_INCLUDED
Chris@16: #define BOOST_MATH_TOOLS_SERIES_INCLUDED
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #pragma once
Chris@16: #endif
Chris@16: 
Chris@16: #include <boost/config/no_tr1/cmath.hpp>
Chris@16: #include <boost/cstdint.hpp>
Chris@16: #include <boost/limits.hpp>
Chris@16: #include <boost/math/tools/config.hpp>
Chris@16: 
Chris@16: namespace boost{ namespace math{ namespace tools{
Chris@16: 
Chris@16: //
Chris@16: // Simple series summation come first:
Chris@16: //
Chris@16: template <class Functor, class U, class V>
Chris@16: inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms, const V& init_value)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    typedef typename Functor::result_type result_type;
Chris@16: 
Chris@16:    boost::uintmax_t counter = max_terms;
Chris@16: 
Chris@16:    result_type result = init_value;
Chris@16:    result_type next_term;
Chris@16:    do{
Chris@16:       next_term = func();
Chris@16:       result += next_term;
Chris@16:    }
Chris@16:    while((fabs(factor * result) < fabs(next_term)) && --counter);
Chris@16: 
Chris@16:    // set max_terms to the actual number of terms of the series evaluated:
Chris@16:    max_terms = max_terms - counter;
Chris@16: 
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class Functor, class U>
Chris@16: inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms)
Chris@16: {
Chris@16:    typename Functor::result_type init_value = 0;
Chris@16:    return sum_series(func, factor, max_terms, init_value);
Chris@16: }
Chris@16: 
Chris@16: template <class Functor, class U>
Chris@16: inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms, const U& init_value)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    typedef typename Functor::result_type result_type;
Chris@16:    result_type factor = ldexp(result_type(1), 1 - bits);
Chris@16:    return sum_series(func, factor, max_terms, init_value);
Chris@16: }
Chris@16: 
Chris@16: template <class Functor>
Chris@16: inline typename Functor::result_type sum_series(Functor& func, int bits)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    typedef typename Functor::result_type result_type;
Chris@16:    boost::uintmax_t iters = (std::numeric_limits<boost::uintmax_t>::max)();
Chris@16:    result_type init_val = 0;
Chris@16:    return sum_series(func, bits, iters, init_val);
Chris@16: }
Chris@16: 
Chris@16: template <class Functor>
Chris@16: inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    typedef typename Functor::result_type result_type;
Chris@16:    result_type init_val = 0;
Chris@16:    return sum_series(func, bits, max_terms, init_val);
Chris@16: }
Chris@16: 
Chris@16: template <class Functor, class U>
Chris@16: inline typename Functor::result_type sum_series(Functor& func, int bits, const U& init_value)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    boost::uintmax_t iters = (std::numeric_limits<boost::uintmax_t>::max)();
Chris@16:    return sum_series(func, bits, iters, init_value);
Chris@16: }
Chris@16: 
Chris@16: //
Chris@16: // Algorithm kahan_sum_series invokes Functor func until the N'th
Chris@16: // term is too small to have any effect on the total, the terms
Chris@16: // are added using the Kahan summation method.
Chris@16: //
Chris@16: // CAUTION: Optimizing compilers combined with extended-precision
Chris@16: // machine registers conspire to render this algorithm partly broken:
Chris@16: // double rounding of intermediate terms (first to a long double machine
Chris@16: // register, and then to a double result) cause the rounding error computed
Chris@16: // by the algorithm to be off by up to 1ulp.  However this occurs rarely, and
Chris@16: // in any case the result is still much better than a naive summation.
Chris@16: //
Chris@16: template <class Functor>
Chris@16: inline typename Functor::result_type kahan_sum_series(Functor& func, int bits)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    typedef typename Functor::result_type result_type;
Chris@16: 
Chris@16:    result_type factor = pow(result_type(2), bits);
Chris@16:    result_type result = func();
Chris@16:    result_type next_term, y, t;
Chris@16:    result_type carry = 0;
Chris@16:    do{
Chris@16:       next_term = func();
Chris@16:       y = next_term - carry;
Chris@16:       t = result + y;
Chris@16:       carry = t - result;
Chris@16:       carry -= y;
Chris@16:       result = t;
Chris@16:    }
Chris@16:    while(fabs(result) < fabs(factor * next_term));
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class Functor>
Chris@16: inline typename Functor::result_type kahan_sum_series(Functor& func, int bits, boost::uintmax_t& max_terms)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    typedef typename Functor::result_type result_type;
Chris@16: 
Chris@16:    boost::uintmax_t counter = max_terms;
Chris@16: 
Chris@16:    result_type factor = ldexp(result_type(1), bits);
Chris@16:    result_type result = func();
Chris@16:    result_type next_term, y, t;
Chris@16:    result_type carry = 0;
Chris@16:    do{
Chris@16:       next_term = func();
Chris@16:       y = next_term - carry;
Chris@16:       t = result + y;
Chris@16:       carry = t - result;
Chris@16:       carry -= y;
Chris@16:       result = t;
Chris@16:    }
Chris@16:    while((fabs(result) < fabs(factor * next_term)) && --counter);
Chris@16: 
Chris@16:    // set max_terms to the actual number of terms of the series evaluated:
Chris@16:    max_terms = max_terms - counter;
Chris@16: 
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: } // namespace tools
Chris@16: } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #endif // BOOST_MATH_TOOLS_SERIES_INCLUDED
Chris@16: