1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
|
// (C) Copyright John Maddock 2005-2006.
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_MATH_TOOLS_SERIES_INCLUDED
#define BOOST_MATH_TOOLS_SERIES_INCLUDED
#ifdef _MSC_VER
#pragma once
#endif
#include <boost/config/no_tr1/cmath.hpp>
#include <boost/cstdint.hpp>
#include <boost/limits.hpp>
#include <boost/math/tools/config.hpp>
namespace boost{ namespace math{ namespace tools{
//
// Simple series summation come first:
//
template <class Functor, class U, class V>
inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms, const V& init_value)
{
BOOST_MATH_STD_USING
typedef typename Functor::result_type result_type;
boost::uintmax_t counter = max_terms;
result_type result = init_value;
result_type next_term;
do{
next_term = func();
result += next_term;
}
while((fabs(factor * result) < fabs(next_term)) && --counter);
// set max_terms to the actual number of terms of the series evaluated:
max_terms = max_terms - counter;
return result;
}
template <class Functor, class U>
inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms)
{
typename Functor::result_type init_value = 0;
return sum_series(func, factor, max_terms, init_value);
}
template <class Functor, class U>
inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms, const U& init_value)
{
BOOST_MATH_STD_USING
typedef typename Functor::result_type result_type;
result_type factor = ldexp(result_type(1), 1 - bits);
return sum_series(func, factor, max_terms, init_value);
}
template <class Functor>
inline typename Functor::result_type sum_series(Functor& func, int bits)
{
BOOST_MATH_STD_USING
typedef typename Functor::result_type result_type;
boost::uintmax_t iters = (std::numeric_limits<boost::uintmax_t>::max)();
result_type init_val = 0;
return sum_series(func, bits, iters, init_val);
}
template <class Functor>
inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms)
{
BOOST_MATH_STD_USING
typedef typename Functor::result_type result_type;
result_type init_val = 0;
return sum_series(func, bits, max_terms, init_val);
}
template <class Functor, class U>
inline typename Functor::result_type sum_series(Functor& func, int bits, const U& init_value)
{
BOOST_MATH_STD_USING
boost::uintmax_t iters = (std::numeric_limits<boost::uintmax_t>::max)();
return sum_series(func, bits, iters, init_value);
}
//
// Algorithm kahan_sum_series invokes Functor func until the N'th
// term is too small to have any effect on the total, the terms
// are added using the Kahan summation method.
//
// CAUTION: Optimizing compilers combined with extended-precision
// machine registers conspire to render this algorithm partly broken:
// double rounding of intermediate terms (first to a long double machine
// register, and then to a double result) cause the rounding error computed
// by the algorithm to be off by up to 1ulp. However this occurs rarely, and
// in any case the result is still much better than a naive summation.
//
template <class Functor>
inline typename Functor::result_type kahan_sum_series(Functor& func, int bits)
{
BOOST_MATH_STD_USING
typedef typename Functor::result_type result_type;
result_type factor = pow(result_type(2), bits);
result_type result = func();
result_type next_term, y, t;
result_type carry = 0;
do{
next_term = func();
y = next_term - carry;
t = result + y;
carry = t - result;
carry -= y;
result = t;
}
while(fabs(result) < fabs(factor * next_term));
return result;
}
template <class Functor>
inline typename Functor::result_type kahan_sum_series(Functor& func, int bits, boost::uintmax_t& max_terms)
{
BOOST_MATH_STD_USING
typedef typename Functor::result_type result_type;
boost::uintmax_t counter = max_terms;
result_type factor = ldexp(result_type(1), bits);
result_type result = func();
result_type next_term, y, t;
result_type carry = 0;
do{
next_term = func();
y = next_term - carry;
t = result + y;
carry = t - result;
carry -= y;
result = t;
}
while((fabs(result) < fabs(factor * next_term)) && --counter);
// set max_terms to the actual number of terms of the series evaluated:
max_terms = max_terms - counter;
return result;
}
} // namespace tools
} // namespace math
} // namespace boost
#endif // BOOST_MATH_TOOLS_SERIES_INCLUDED
|
