]> git.proxmox.com Git - ceph.git/blame - ceph/src/boost/libs/math/include/boost/math/tools/series.hpp
bump version to 12.2.2-pve1
[ceph.git] / ceph / src / boost / libs / math / include / boost / math / tools / series.hpp
CommitLineData
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
18namespace boost{ namespace math{ namespace tools{
19
20//
21// Simple series summation come first:
22//
23template <class Functor, class U, class V>
24inline 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
46template <class Functor, class U>
47inline 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
53template <class Functor, class U>
54inline 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
62template <class Functor>
63inline 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
72template <class Functor>
73inline 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
81template <class Functor, class U>
82inline 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//
101template <class Functor>
102inline 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
124template <class Functor>
125inline 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