1 // (C) Copyright John Maddock 2005-2006.
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
6 #ifndef BOOST_MATH_TOOLS_SERIES_INCLUDED
7 #define BOOST_MATH_TOOLS_SERIES_INCLUDED
13 #include <boost/config/no_tr1/cmath.hpp>
14 #include <boost/cstdint.hpp>
15 #include <boost/limits.hpp>
16 #include <boost/math/tools/config.hpp>
18 namespace boost{ namespace math{ namespace tools{
21 // Simple series summation come first:
23 template <class Functor, class U, class V>
24 inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms, const V& init_value)
28 typedef typename Functor::result_type result_type;
30 boost::uintmax_t counter = max_terms;
32 result_type result = init_value;
33 result_type next_term;
38 while((fabs(factor * result) < fabs(next_term)) && --counter);
40 // set max_terms to the actual number of terms of the series evaluated:
41 max_terms = max_terms - counter;
46 template <class Functor, class U>
47 inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms)
49 typename Functor::result_type init_value = 0;
50 return sum_series(func, factor, max_terms, init_value);
53 template <class Functor, class U>
54 inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms, const U& init_value)
57 typedef typename Functor::result_type result_type;
58 result_type factor = ldexp(result_type(1), 1 - bits);
59 return sum_series(func, factor, max_terms, init_value);
62 template <class Functor>
63 inline typename Functor::result_type sum_series(Functor& func, int bits)
66 typedef typename Functor::result_type result_type;
67 boost::uintmax_t iters = (std::numeric_limits<boost::uintmax_t>::max)();
68 result_type init_val = 0;
69 return sum_series(func, bits, iters, init_val);
72 template <class Functor>
73 inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms)
76 typedef typename Functor::result_type result_type;
77 result_type init_val = 0;
78 return sum_series(func, bits, max_terms, init_val);
81 template <class Functor, class U>
82 inline typename Functor::result_type sum_series(Functor& func, int bits, const U& init_value)
85 boost::uintmax_t iters = (std::numeric_limits<boost::uintmax_t>::max)();
86 return sum_series(func, bits, iters, init_value);
90 // Algorithm kahan_sum_series invokes Functor func until the N'th
91 // term is too small to have any effect on the total, the terms
92 // are added using the Kahan summation method.
94 // CAUTION: Optimizing compilers combined with extended-precision
95 // machine registers conspire to render this algorithm partly broken:
96 // double rounding of intermediate terms (first to a long double machine
97 // register, and then to a double result) cause the rounding error computed
98 // by the algorithm to be off by up to 1ulp. However this occurs rarely, and
99 // in any case the result is still much better than a naive summation.
101 template <class Functor>
102 inline typename Functor::result_type kahan_sum_series(Functor& func, int bits)
106 typedef typename Functor::result_type result_type;
108 result_type factor = pow(result_type(2), bits);
109 result_type result = func();
110 result_type next_term, y, t;
111 result_type carry = 0;
114 y = next_term - carry;
120 while(fabs(result) < fabs(factor * next_term));
124 template <class Functor>
125 inline typename Functor::result_type kahan_sum_series(Functor& func, int bits, boost::uintmax_t& max_terms)
129 typedef typename Functor::result_type result_type;
131 boost::uintmax_t counter = max_terms;
133 result_type factor = ldexp(result_type(1), bits);
134 result_type result = func();
135 result_type next_term, y, t;
136 result_type carry = 0;
139 y = next_term - carry;
145 while((fabs(result) < fabs(factor * next_term)) && --counter);
147 // set max_terms to the actual number of terms of the series evaluated:
148 max_terms = max_terms - counter;
157 #endif // BOOST_MATH_TOOLS_SERIES_INCLUDED