GeographicLib
1.38
Main Page
Related Pages
Namespaces
Classes
Files
File List
File Members
All
Classes
Namespaces
Files
Functions
Variables
Typedefs
Enumerations
Enumerator
Friends
Macros
Pages
include
GeographicLib
Accumulator.hpp
Go to the documentation of this file.
1
/**
2
* \file Accumulator.hpp
3
* \brief Header for GeographicLib::Accumulator class
4
*
5
* Copyright (c) Charles Karney (2010-2011) <charles@karney.com> and licensed
6
* under the MIT/X11 License. For more information, see
7
* http://geographiclib.sourceforge.net/
8
**********************************************************************/
9
10
#if !defined(GEOGRAPHICLIB_ACCUMULATOR_HPP)
11
#define GEOGRAPHICLIB_ACCUMULATOR_HPP 1
12
13
#include <
GeographicLib/Constants.hpp
>
14
15
namespace
GeographicLib {
16
17
/**
18
* \brief An accumulator for sums
19
*
20
* This allow many numbers of floating point type \e T to be added together
21
* with twice the normal precision. Thus if \e T is double, the effective
22
* precision of the sum is 106 bits or about 32 decimal places.
23
*
24
* The implementation follows J. R. Shewchuk,
25
* <a href="http://dx.doi.org/10.1007/PL00009321"> Adaptive Precision
26
* Floating-Point Arithmetic and Fast Robust Geometric Predicates</a>,
27
* Discrete & Computational Geometry 18(3) 305--363 (1997).
28
*
29
* Approximate timings (summing a vector<double>)
30
* - double: 2ns
31
* - Accumulator<double>: 23ns
32
*
33
* In the documentation of the member functions, \e sum stands for the value
34
* currently held in the accumulator.
35
*
36
* Example of use:
37
* \include example-Accumulator.cpp
38
**********************************************************************/
39
template
<
typename
T = Math::real>
40
class
GEOGRAPHICLIB_EXPORT
Accumulator
{
41
private
:
42
// _s + _t accumulators for the sum.
43
T _s, _t;
44
// Same as Math::sum, but requires abs(u) >= abs(v). This isn't currently
45
// used.
46
static
inline
T fastsum(T u, T v, T& t) {
47
GEOGRAPHICLIB_VOLATILE
T s = u + v;
48
GEOGRAPHICLIB_VOLATILE
T vp = s - u;
49
t = v - vp;
50
return
s;
51
}
52
void
Add(T y) {
53
// Here's Shewchuk's solution...
54
T u;
// hold exact sum as [s, t, u]
55
y =
Math::sum
(y, _t, u);
// Accumulate starting at least significant end
56
_s =
Math::sum
(y, _s, _t);
57
// Start is _s, _t decreasing and non-adjacent. Sum is now (s + t + u)
58
// exactly with s, t, u non-adjacent and in decreasing order (except for
59
// possible zeros). The following code tries to normalize the result.
60
// Ideally, we want _s = round(s+t+u) and _u = round(s+t+u - _s). The
61
// following does an approximate job (and maintains the decreasing
62
// non-adjacent property). Here are two "failures" using 3-bit floats:
63
//
64
// Case 1: _s is not equal to round(s+t+u) -- off by 1 ulp
65
// [12, -1] - 8 -> [4, 0, -1] -> [4, -1] = 3 should be [3, 0] = 3
66
//
67
// Case 2: _s+_t is not as close to s+t+u as it shold be
68
// [64, 5] + 4 -> [64, 8, 1] -> [64, 8] = 72 (off by 1)
69
// should be [80, -7] = 73 (exact)
70
//
71
// "Fixing" these problems is probably not worth the expense. The
72
// representation inevitably leads to small errors in the accumulated
73
// values. The additional errors illustrated here amount to 1 ulp of the
74
// less significant word during each addition to the Accumulator and an
75
// additional possible error of 1 ulp in the reported sum.
76
//
77
// Incidentally, the "ideal" representation described above is not
78
// canonical, because _s = round(_s + _t) may not be true. For example,
79
// with 3-bit floats:
80
//
81
// [128, 16] + 1 -> [160, -16] -- 160 = round(145).
82
// But [160, 0] - 16 -> [128, 16] -- 128 = round(144).
83
//
84
if
(_s == 0)
// This implies t == 0,
85
_s = u;
// so result is u
86
else
87
_t += u;
// otherwise just accumulate u to t.
88
}
89
T Sum(T y)
const
{
90
Accumulator
a(*
this
);
91
a.Add(y);
92
return
a._s;
93
}
94
public
:
95
/**
96
* Construct from a \e T. This is not declared explicit, so that you can
97
* write <code>Accumulator<double> a = 5;</code>.
98
*
99
* @param[in] y set \e sum = \e y.
100
**********************************************************************/
101
Accumulator
(T y = T(0)) : _s(y), _t(0) {
102
GEOGRAPHICLIB_STATIC_ASSERT(!std::numeric_limits<T>::is_integer,
103
"Accumulator type is not floating point"
);
104
}
105
/**
106
* Set the accumulator to a number.
107
*
108
* @param[in] y set \e sum = \e y.
109
**********************************************************************/
110
Accumulator
&
operator=
(T y) { _s = y; _t = 0;
return
*
this
; }
111
/**
112
* Return the value held in the accumulator.
113
*
114
* @return \e sum.
115
**********************************************************************/
116
T
operator()
()
const
{
return
_s; }
117
/**
118
* Return the result of adding a number to \e sum (but don't change \e sum).
119
*
120
* @param[in] y the number to be added to the sum.
121
* @return \e sum + \e y.
122
**********************************************************************/
123
T
operator()
(T y)
const
{
return
Sum(y); }
124
/**
125
* Add a number to the accumulator.
126
*
127
* @param[in] y set \e sum += \e y.
128
**********************************************************************/
129
Accumulator
&
operator+=
(T y) { Add(y);
return
*
this
; }
130
/**
131
* Subtract a number from the accumulator.
132
*
133
* @param[in] y set \e sum -= \e y.
134
**********************************************************************/
135
Accumulator
&
operator-=
(T y) { Add(-y);
return
*
this
; }
136
/**
137
* Multiply accumulator by an integer. To avoid loss of accuracy, use only
138
* integers such that \e n × \e T is exactly representable as a \e T
139
* (i.e., ± powers of two). Use \e n = −1 to negate \e sum.
140
*
141
* @param[in] n set \e sum *= \e n.
142
**********************************************************************/
143
Accumulator
&
operator*=
(
int
n) { _s *= n; _t *= n;
return
*
this
; }
144
/**
145
* Test equality of an Accumulator with a number.
146
**********************************************************************/
147
bool
operator==
(T y)
const
{
return
_s == y; }
148
/**
149
* Test inequality of an Accumulator with a number.
150
**********************************************************************/
151
bool
operator!=
(T y)
const
{
return
_s != y; }
152
/**
153
* Less operator on an Accumulator and a number.
154
**********************************************************************/
155
bool
operator<
(T y)
const
{
return
_s < y; }
156
/**
157
* Less or equal operator on an Accumulator and a number.
158
**********************************************************************/
159
bool
operator<=
(T y)
const
{
return
_s <= y; }
160
/**
161
* Greater operator on an Accumulator and a number.
162
**********************************************************************/
163
bool
operator>
(T y)
const
{
return
_s > y; }
164
/**
165
* Greater or equal operator on an Accumulator and a number.
166
**********************************************************************/
167
bool
operator>=
(T y)
const
{
return
_s >= y; }
168
};
169
170
}
// namespace GeographicLib
171
172
#endif // GEOGRAPHICLIB_ACCUMULATOR_HPP
Generated by
1.8.3.1