| // (C) Copyright John Maddock 2005-2006. |
| // 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_TOOLS_SERIES_INCLUDED |
| #define BOOST_MATH_TOOLS_SERIES_INCLUDED |
| |
| #ifdef _MSC_VER |
| #pragma once |
| #endif |
| |
| #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{ |
| |
| // |
| // Simple series summation come first: |
| // |
| 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 |
| |
| typedef typename Functor::result_type result_type; |
| |
| boost::uintmax_t counter = max_terms; |
| |
| result_type result = init_value; |
| result_type next_term; |
| do{ |
| next_term = func(); |
| result += next_term; |
| } |
| 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; |
| |
| return result; |
| } |
| |
| template <class Functor, class U> |
| inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms) |
| { |
| 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, 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); |
| } |
| |
| 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); |
| } |
| |
| 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); |
| } |
| |
| 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); |
| } |
| |
| // |
| // 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> |
| inline typename Functor::result_type kahan_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, 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(fabs(result) < fabs(factor * next_term)); |
| return result; |
| } |
| |
| template <class Functor> |
| inline typename Functor::result_type kahan_sum_series(Functor& func, int bits, boost::uintmax_t& max_terms) |
| { |
| BOOST_MATH_STD_USING |
| |
| typedef typename Functor::result_type result_type; |
| |
| boost::uintmax_t counter = max_terms; |
| |
| result_type factor = ldexp(result_type(1), 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((fabs(result) < fabs(factor * next_term)) && --counter); |
| |
| // set max_terms to the actual number of terms of the series evaluated: |
| max_terms = max_terms - counter; |
| |
| return result; |
| } |
| |
| } // namespace tools |
| } // namespace math |
| } // namespace boost |
| |
| #endif // BOOST_MATH_TOOLS_SERIES_INCLUDED |
| |