GeographicLib  1.37
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 
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>
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 &times; \e T is exactly representable as a \e T
139  * (i.e., &plusmn; powers of two). Use \e n = &minus;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
static T sum(T u, T v, T &t)
Definition: Math.hpp:378
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:70
bool operator==(T y) const
Accumulator & operator=(T y)
bool operator!=(T y) const
#define GEOGRAPHICLIB_VOLATILE
Definition: Math.hpp:83
An accumulator for sums.
Definition: Accumulator.hpp:40
bool operator>=(T y) const
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
Accumulator & operator+=(T y)
Accumulator & operator-=(T y)
bool operator>(T y) const
Header for GeographicLib::Constants class.
bool operator<=(T y) const
bool operator<(T y) const
Accumulator & operator*=(int n)