]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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) | |
5 | ||
6 | #ifndef BOOST_MATH_TOOLS_SERIES_INCLUDED | |
7 | #define BOOST_MATH_TOOLS_SERIES_INCLUDED | |
8 | ||
9 | #ifdef _MSC_VER | |
10 | #pragma once | |
11 | #endif | |
12 | ||
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> | |
17 | ||
18 | namespace boost{ namespace math{ namespace tools{ | |
19 | ||
20 | // | |
21 | // Simple series summation come first: | |
22 | // | |
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) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()())) | |
25 | { | |
26 | BOOST_MATH_STD_USING | |
27 | ||
28 | typedef typename Functor::result_type result_type; | |
29 | ||
30 | boost::uintmax_t counter = max_terms; | |
31 | ||
32 | result_type result = init_value; | |
33 | result_type next_term; | |
34 | do{ | |
35 | next_term = func(); | |
36 | result += next_term; | |
37 | } | |
38 | while((fabs(factor * result) < fabs(next_term)) && --counter); | |
39 | ||
40 | // set max_terms to the actual number of terms of the series evaluated: | |
41 | max_terms = max_terms - counter; | |
42 | ||
43 | return result; | |
44 | } | |
45 | ||
46 | template <class Functor, class U> | |
47 | inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()())) | |
48 | { | |
49 | typename Functor::result_type init_value = 0; | |
50 | return sum_series(func, factor, max_terms, init_value); | |
51 | } | |
52 | ||
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) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()())) | |
55 | { | |
56 | BOOST_MATH_STD_USING | |
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); | |
60 | } | |
61 | ||
62 | template <class Functor> | |
63 | inline typename Functor::result_type sum_series(Functor& func, int bits) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()())) | |
64 | { | |
65 | BOOST_MATH_STD_USING | |
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); | |
70 | } | |
71 | ||
72 | template <class Functor> | |
73 | inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()())) | |
74 | { | |
75 | BOOST_MATH_STD_USING | |
76 | typedef typename Functor::result_type result_type; | |
77 | result_type init_val = 0; | |
78 | return sum_series(func, bits, max_terms, init_val); | |
79 | } | |
80 | ||
81 | template <class Functor, class U> | |
82 | inline typename Functor::result_type sum_series(Functor& func, int bits, const U& init_value) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()())) | |
83 | { | |
84 | BOOST_MATH_STD_USING | |
85 | boost::uintmax_t iters = (std::numeric_limits<boost::uintmax_t>::max)(); | |
86 | return sum_series(func, bits, iters, init_value); | |
87 | } | |
88 | ||
89 | // | |
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. | |
93 | // | |
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. | |
100 | // | |
101 | template <class Functor> | |
102 | inline typename Functor::result_type kahan_sum_series(Functor& func, int bits) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()())) | |
103 | { | |
104 | BOOST_MATH_STD_USING | |
105 | ||
106 | typedef typename Functor::result_type result_type; | |
107 | ||
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; | |
112 | do{ | |
113 | next_term = func(); | |
114 | y = next_term - carry; | |
115 | t = result + y; | |
116 | carry = t - result; | |
117 | carry -= y; | |
118 | result = t; | |
119 | } | |
120 | while(fabs(result) < fabs(factor * next_term)); | |
121 | return result; | |
122 | } | |
123 | ||
124 | template <class Functor> | |
125 | inline typename Functor::result_type kahan_sum_series(Functor& func, int bits, boost::uintmax_t& max_terms) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()())) | |
126 | { | |
127 | BOOST_MATH_STD_USING | |
128 | ||
129 | typedef typename Functor::result_type result_type; | |
130 | ||
131 | boost::uintmax_t counter = max_terms; | |
132 | ||
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; | |
137 | do{ | |
138 | next_term = func(); | |
139 | y = next_term - carry; | |
140 | t = result + y; | |
141 | carry = t - result; | |
142 | carry -= y; | |
143 | result = t; | |
144 | } | |
145 | while((fabs(result) < fabs(factor * next_term)) && --counter); | |
146 | ||
147 | // set max_terms to the actual number of terms of the series evaluated: | |
148 | max_terms = max_terms - counter; | |
149 | ||
150 | return result; | |
151 | } | |
152 | ||
153 | } // namespace tools | |
154 | } // namespace math | |
155 | } // namespace boost | |
156 | ||
157 | #endif // BOOST_MATH_TOOLS_SERIES_INCLUDED | |
158 |