]>
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 | } | |
92f5a8d4 | 38 | while((abs(factor * result) < abs(next_term)) && --counter); |
7c673cae FG |
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 | } | |
92f5a8d4 TL |
88 | // |
89 | // Checked summation: | |
90 | // | |
91 | template <class Functor, class U, class V> | |
92 | inline typename Functor::result_type checked_sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms, const V& init_value, V& norm) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()())) | |
93 | { | |
94 | BOOST_MATH_STD_USING | |
95 | ||
96 | typedef typename Functor::result_type result_type; | |
97 | ||
98 | boost::uintmax_t counter = max_terms; | |
99 | ||
100 | result_type result = init_value; | |
101 | result_type next_term; | |
102 | do { | |
103 | next_term = func(); | |
104 | result += next_term; | |
105 | norm += fabs(next_term); | |
106 | } while ((abs(factor * result) < abs(next_term)) && --counter); | |
107 | ||
108 | // set max_terms to the actual number of terms of the series evaluated: | |
109 | max_terms = max_terms - counter; | |
110 | ||
111 | return result; | |
112 | } | |
113 | ||
7c673cae FG |
114 | |
115 | // | |
116 | // Algorithm kahan_sum_series invokes Functor func until the N'th | |
117 | // term is too small to have any effect on the total, the terms | |
118 | // are added using the Kahan summation method. | |
119 | // | |
120 | // CAUTION: Optimizing compilers combined with extended-precision | |
121 | // machine registers conspire to render this algorithm partly broken: | |
122 | // double rounding of intermediate terms (first to a long double machine | |
123 | // register, and then to a double result) cause the rounding error computed | |
124 | // by the algorithm to be off by up to 1ulp. However this occurs rarely, and | |
125 | // in any case the result is still much better than a naive summation. | |
126 | // | |
127 | template <class Functor> | |
128 | 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>()())) | |
129 | { | |
130 | BOOST_MATH_STD_USING | |
131 | ||
132 | typedef typename Functor::result_type result_type; | |
133 | ||
134 | result_type factor = pow(result_type(2), bits); | |
135 | result_type result = func(); | |
136 | result_type next_term, y, t; | |
137 | result_type carry = 0; | |
138 | do{ | |
139 | next_term = func(); | |
140 | y = next_term - carry; | |
141 | t = result + y; | |
142 | carry = t - result; | |
143 | carry -= y; | |
144 | result = t; | |
145 | } | |
146 | while(fabs(result) < fabs(factor * next_term)); | |
147 | return result; | |
148 | } | |
149 | ||
150 | template <class Functor> | |
151 | 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>()())) | |
152 | { | |
153 | BOOST_MATH_STD_USING | |
154 | ||
155 | typedef typename Functor::result_type result_type; | |
156 | ||
157 | boost::uintmax_t counter = max_terms; | |
158 | ||
159 | result_type factor = ldexp(result_type(1), bits); | |
160 | result_type result = func(); | |
161 | result_type next_term, y, t; | |
162 | result_type carry = 0; | |
163 | do{ | |
164 | next_term = func(); | |
165 | y = next_term - carry; | |
166 | t = result + y; | |
167 | carry = t - result; | |
168 | carry -= y; | |
169 | result = t; | |
170 | } | |
171 | while((fabs(result) < fabs(factor * next_term)) && --counter); | |
172 | ||
173 | // set max_terms to the actual number of terms of the series evaluated: | |
174 | max_terms = max_terms - counter; | |
175 | ||
176 | return result; | |
177 | } | |
178 | ||
179 | } // namespace tools | |
180 | } // namespace math | |
181 | } // namespace boost | |
182 | ||
183 | #endif // BOOST_MATH_TOOLS_SERIES_INCLUDED | |
184 |