]> git.proxmox.com Git - ceph.git/blob - ceph/src/boost/boost/math/tools/series.hpp
import new upstream nautilus stable release 14.2.8
[ceph.git] / ceph / src / boost / boost / math / tools / series.hpp
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((abs(factor * result) < abs(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 // 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
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