HDK
UT_Accumulator.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
7  * Accumulators for scalar and vector types. These are used to collect the
8  * sum of incrementally accumulated quantities. Using the Neumaier's
9  * extension to Kahan's summation method, these should vastly limit the
10  * truncation error of floating-points.
11  */
12
13 #ifndef __UT_Accumulator__
14 #define __UT_Accumulator__
15
16 #include "UT_FixedVector.h"
17 #include <SYS/SYS_Types.h>
18
19
20 namespace UT
21 {
22
23 template <typename T>
25 {
26 public:
27  explicit Accumulator(const T &init_val = T(0.0)) :
28  mySum(init_val), myCompensation(T(0.0)) { }
29
31  T &operator+=(const T& rhs)
32  {
33
34  T t = mySum + rhs;
35  if (SYSabs(mySum) >= SYSabs(rhs))
36  myCompensation += (mySum - t) + rhs;
37  else
38  myCompensation += (rhs - t) + mySum;
39  mySum = t;
40  return mySum;
41  }
42
43  T sum() const { return mySum + myCompensation; }
44
45 private:
46  T mySum;
47  T myCompensation;
48 };
49
50 template <typename T, int S>
52 {
53 public:
55
56  explicit Accumulator(const V &init_val = V(T(0.0))) :
57  mySum(init_val), myCompensation(T(0.0)) { }
58
60  V &operator+=(const V& rhs)
61  {
62  for (int i = 0, ie = S; i < ie; i++)
63  {
64  T t = mySum(i) + rhs(i);
65  if (SYSabs(mySum(i)) >= SYSabs(rhs(i)))
66  myCompensation(i) += (mySum(i) - t) + rhs(i);
67  else
68  myCompensation(i) += (rhs(i) - t) + mySum(i);
69  mySum(i) = t;
70  }
71  return mySum;
72  }
73
74  V sum() const { return mySum + myCompensation; }
75
76 private:
77  V mySum;
78  V myCompensation;
79 };
80
81
82 template <typename T>
83 class Accumulator<std::complex<T>>
84 {
85 public:
86  using C = std::complex<T>;
87
88  explicit Accumulator(const C &init_val = C(0.0)) :
89  mySum(init_val), myCompensation(C(0.0)) { }
90
92  C &operator+=(const C& rhs)
93  {
94  T t = mySum.real() + rhs.real();
95  if (SYSabs(mySum.real()) >= SYSabs(rhs.real()))
96  myCompensation.real(myCompensation.real()
97  + (mySum.real() - t) + rhs.real());
98  else
99  myCompensation.real(myCompensation.real()
100  + (rhs.real() - t) + mySum.real());
101  mySum.real(t);
102
103  t = mySum.imag() + rhs.imag();
104  if (SYSabs(mySum.imag()) >= SYSabs(rhs.imag()))
105  myCompensation.imag(myCompensation.imag()
106  + (mySum.imag() - t) + rhs.imag());
107  else
108  myCompensation.imag(myCompensation.imag()
109  + (rhs.imag() - t) + mySum.imag());
110  mySum.imag(t);
111
112  return mySum;
113  }
114
115  C sum() const { return mySum + myCompensation; }
116
117 private:
118  C mySum;
119  C myCompensation;
120 };
121
122
123 } // namespace UT
124
125 template <typename T>
127
128
129 #endif
#define SYSabs(a)
Definition: SYS_Math.h:1540
Accumulator(const T &init_val=T(0.0))
SYS_FORCE_INLINE C & operator+=(const C &rhs)
SYS_FORCE_INLINE V & operator+=(const V &rhs)
#define SYS_FORCE_INLINE
Definition: SYS_Inline.h:45
Accumulator(const C &init_val=C(0.0))
GLdouble t