Accumulator.hpp
Go to the documentation of this file.
1 
10 #if !defined(GEOGRAPHICLIB_ACCUMULATOR_HPP)
11 #define GEOGRAPHICLIB_ACCUMULATOR_HPP 1
12 
14 
15 namespace GeographicLib {
16 
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 T fastsum(T u, T v, T& t) {
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  // Accumulate starting at least significant end
56  y = Math::sum(y, _t, u);
57  _s = Math::sum(y, _s, _t);
58  // Start is _s, _t decreasing and non-adjacent. Sum is now (s + t + u)
59  // exactly with s, t, u non-adjacent and in decreasing order (except for
60  // possible zeros). The following code tries to normalize the result.
61  // Ideally, we want _s = round(s+t+u) and _u = round(s+t+u - _s). The
62  // following does an approximate job (and maintains the decreasing
63  // non-adjacent property). Here are two "failures" using 3-bit floats:
64  //
65  // Case 1: _s is not equal to round(s+t+u) -- off by 1 ulp
66  // [12, -1] - 8 -> [4, 0, -1] -> [4, -1] = 3 should be [3, 0] = 3
67  //
68  // Case 2: _s+_t is not as close to s+t+u as it shold be
69  // [64, 5] + 4 -> [64, 8, 1] -> [64, 8] = 72 (off by 1)
70  // should be [80, -7] = 73 (exact)
71  //
72  // "Fixing" these problems is probably not worth the expense. The
73  // representation inevitably leads to small errors in the accumulated
74  // values. The additional errors illustrated here amount to 1 ulp of the
75  // less significant word during each addition to the Accumulator and an
76  // additional possible error of 1 ulp in the reported sum.
77  //
78  // Incidentally, the "ideal" representation described above is not
79  // canonical, because _s = round(_s + _t) may not be true. For example,
80  // with 3-bit floats:
81  //
82  // [128, 16] + 1 -> [160, -16] -- 160 = round(145).
83  // But [160, 0] - 16 -> [128, 16] -- 128 = round(144).
84  //
85  if (_s == 0) // This implies t == 0,
86  _s = u; // so result is u
87  else
88  _t += u; // otherwise just accumulate u to t.
89  }
90  T Sum(T y) const {
91  Accumulator a(*this);
92  a.Add(y);
93  return a._s;
94  }
95  public:
102  Accumulator(T y = T(0)) : _s(y), _t(0) {
103  GEOGRAPHICLIB_STATIC_ASSERT(!std::numeric_limits<T>::is_integer,
104  "Accumulator type is not floating point");
105  }
111  Accumulator& operator=(T y) { _s = y; _t = 0; return *this; }
117  T operator()() const { return _s; }
125  T operator()(T y) const { return Sum(y); }
131  Accumulator& operator+=(T y) { Add(y); return *this; }
137  Accumulator& operator-=(T y) { Add(-y); return *this; }
145  Accumulator& operator*=(int n) { _s *= n; _t *= n; return *this; }
153  T d = _s; _s *= y;
154  d = Math::fma(y, d, -_s); // the error in the first multiplication
155  _t = Math::fma(y, _t, d); // add error to the second term
156  return *this;
157  }
161  bool operator==(T y) const { return _s == y; }
165  bool operator!=(T y) const { return _s != y; }
169  bool operator<(T y) const { return _s < y; }
173  bool operator<=(T y) const { return _s <= y; }
177  bool operator>(T y) const { return _s > y; }
181  bool operator>=(T y) const { return _s >= y; }
182  };
183 
184 } // namespace GeographicLib
185 
186 #endif // GEOGRAPHICLIB_ACCUMULATOR_HPP
static T sum(T u, T v, T &t)
Definition: Math.hpp:399
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:91
Scalar * y
Accumulator & operator=(T y)
int n
#define GEOGRAPHICLIB_VOLATILE
Definition: Math.hpp:84
bool operator!=(T y) const
An accumulator for sums.
Definition: Accumulator.hpp:40
Accumulator & operator*=(T y)
bool operator<=(T y) const
Array< int, Dynamic, 1 > v
Eigen::Triplet< double > T
Namespace for GeographicLib.
RealScalar s
Accumulator & operator+=(T y)
Accumulator & operator-=(T y)
bool operator==(T y) const
static T fastsum(T u, T v, T &t)
Definition: Accumulator.hpp:46
Header for GeographicLib::Constants class.
bool operator>=(T y) const
static T fma(T x, T y, T z)
Definition: Math.hpp:369
Point2 t(10, 10)
Accumulator & operator*=(int n)


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:33:53