Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
Util.h
Go to the documentation of this file.
1// @(#)root/mathcore:$Id$
2// Author: L. Moneta Tue Nov 14 15:44:38 2006
3
4/**********************************************************************
5 * *
6 * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
7 * *
8 * *
9 **********************************************************************/
10
11// Utility functions for all ROOT Math classes
12
13#ifndef ROOT_Math_Util
14#define ROOT_Math_Util
15
16#include <string>
17#include <sstream>
18
19#include <cmath>
20#include <limits>
21#include <numeric>
22
23
24// This can be protected against by defining ROOT_Math_VecTypes
25// This is only used for the R__HAS_VECCORE define
26// and a single VecCore function in EvalLog
27#ifndef ROOT_Math_VecTypes
28#include "Types.h"
29#endif
30
31
32// for defining unused variables in the interfaces
33// and have still them in the documentation
34#define MATH_UNUSED(var) (void)var
35
36
37namespace ROOT {
38
39 namespace Math {
40
41 /**
42 namespace defining Utility functions needed by mathcore
43 */
44 namespace Util {
45
46 /**
47 Utility function for conversion to strings
48 */
49 template <class T>
50 std::string ToString(const T &val)
51 {
52 std::ostringstream buf;
53 buf << val;
54
55 std::string ret = buf.str();
56 return ret;
57 }
58
59 /// safe evaluation of log(x) with a protections against negative or zero argument to the log
60 /// smooth linear extrapolation below function values smaller than epsilon
61 /// (better than a simple cut-off)
62
63 template<class T>
64 inline T EvalLog(T x) {
65 static const T epsilon = T(2.0 * std::numeric_limits<double>::min());
66#ifdef R__HAS_VECCORE
67 T logval = vecCore::Blend<T>(x <= epsilon, x / epsilon + std::log(epsilon) - T(1.0), std::log(x));
68#else
69 T logval = x <= epsilon ? x / epsilon + std::log(epsilon) - T(1.0) : std::log(x);
70#endif
71 return logval;
72 }
73
74 } // end namespace Util
75
76 /// \class KahanSum
77 /// The Kahan summation is a compensated summation algorithm, which significantly reduces numerical errors
78 /// when adding a sequence of finite-precision floating point numbers.
79 /// This is done by keeping a separate running compensation (a variable to accumulate small errors).
80 ///
81 /// ### Auto-vectorisable accumulation
82 /// This class can internally use multiple accumulators (template parameter `N`).
83 /// When filled from a collection that supports index access from a *contiguous* block of memory,
84 /// compilers such as gcc, clang and icc can auto-vectorise the accumulation. This happens by cycling
85 /// through the internal accumulators based on the value of "`index % N`", so `N` accumulators can be filled from a block
86 /// of `N` numbers in a single instruction.
87 ///
88 /// The usage of multiple accumulators might slightly increase the precision in comparison to the single-accumulator version
89 /// with `N = 1`.
90 /// This depends on the order and magnitude of the numbers being accumulated. Therefore, in rare cases, the accumulation
91 /// result can change *in dependence of N*, even when the data are identical.
92 /// The magnitude of such differences is well below the precision of the floating point type, and will therefore mostly show
93 /// in the compensation sum(see Carry()). Increasing the number of accumulators therefore only makes sense to
94 /// speed up the accumulation, but not to increase precision.
95 ///
96 /// \param T The type of the values to be accumulated.
97 /// \param N Number of accumulators. Defaults to 1. Ideal values are the widths of a vector register on the relevant architecture.
98 /// Depending on the instruction set, good values are:
99 /// - AVX2-float: 8
100 /// - AVX2-double: 4
101 /// - AVX512-float: 16
102 /// - AVX512-double: 8
103 ///
104 /// ### Examples
105 ///
106 /// ~~~{.cpp}
107 /// std::vector<double> numbers(1000);
108 /// for (std::size_t i=0; i<1000; ++i) {
109 /// numbers[i] = rand();
110 /// }
111 ///
112 /// ROOT::Math::KahanSum<double, 4> k;
113 /// k.Add(numbers.begin(), numbers.end());
114 /// // or
115 /// k.Add(numbers);
116 /// ~~~
117 /// ~~~{.cpp}
118 /// double offset = 10.;
119 /// auto result = ROOT::Math::KahanSum<double, 4>::Accumulate(numbers.begin(), numbers.end(), offset);
120 /// ~~~
121 template<typename T = double, unsigned int N = 1>
122 class KahanSum {
123 public:
124 /// Initialise the sum.
125 /// \param[in] initialValue Initialise with this value. Defaults to 0.
126 explicit KahanSum(T initialValue = T{}) {
127 fSum[0] = initialValue;
128 std::fill(std::begin(fSum)+1, std::end(fSum), 0.);
129 std::fill(std::begin(fCarry), std::end(fCarry), 0.);
130 }
131
132 /// Initialise with a sum value and a carry value.
133 /// \param[in] initialSumValue Initialise the sum with this value.
134 /// \param[in] initialCarryValue Initialise the carry with this value.
135 KahanSum(T initialSumValue, T initialCarryValue) {
136 fSum[0] = initialSumValue;
137 fCarry[0] = initialCarryValue;
138 std::fill(std::begin(fSum)+1, std::end(fSum), 0.);
139 std::fill(std::begin(fCarry)+1, std::end(fCarry), 0.);
140 }
141
142 /// Initialise the sum with a pre-existing state.
143 /// \param[in] sumBegin Begin of iterator range with values to initialise the sum with.
144 /// \param[in] sumEnd End of iterator range with values to initialise the sum with.
145 /// \param[in] carryBegin Begin of iterator range with values to initialise the carry with.
146 /// \param[in] carryEnd End of iterator range with values to initialise the carry with.
147 template<class Iterator>
148 KahanSum(Iterator sumBegin, Iterator sumEnd, Iterator carryBegin, Iterator carryEnd) {
149 assert(std::distance(sumBegin, sumEnd) == N);
150 assert(std::distance(carryBegin, carryEnd) == N);
151 std::copy(sumBegin, sumEnd, std::begin(fSum));
152 std::copy(carryBegin, carryEnd, std::begin(fCarry));
153 }
154
155 /// Constructor to create a KahanSum from another KahanSum with a different number of accumulators
156 template <unsigned int M>
157 KahanSum(KahanSum<T,M> const& other) {
158 fSum[0] = other.Sum();
159 fCarry[0] = other.Carry();
160 std::fill(std::begin(fSum)+1, std::end(fSum), 0.);
161 std::fill(std::begin(fCarry)+1, std::end(fCarry), 0.);
162 }
163
164 /// Single-element accumulation. Will not vectorise.
165 void Add(T x) {
166 auto y = x - fCarry[0];
167 auto t = fSum[0] + y;
168 fCarry[0] = (t - fSum[0]) - y;
169 fSum[0] = t;
170 }
171
172
173 /// Accumulate from a range denoted by iterators.
174 ///
175 /// This function will auto-vectorise with random-access iterators.
176 /// \param[in] begin Beginning of a range. Needs to be a random access iterator for automatic
177 /// vectorisation, because a contiguous block of memory needs to be read.
178 /// \param[in] end End of the range.
179 template <class Iterator>
180 void Add(Iterator begin, Iterator end) {
181 static_assert(std::is_floating_point<
182 typename std::remove_reference<decltype(*begin)>::type>::value,
183 "Iterator needs to point to floating-point values.");
184 const std::size_t n = std::distance(begin, end);
185
186 for (std::size_t i=0; i<n; ++i) {
187 AddIndexed(*(begin++), i);
188 }
189 }
190
191
192 /// Fill from a container that supports index access.
193 /// \param[in] inputs Container with index access such as std::vector or array.
194 template<class Container_t>
195 void Add(const Container_t& inputs) {
196 static_assert(std::is_floating_point<typename Container_t::value_type>::value,
197 "Container does not hold floating-point values.");
198 for (std::size_t i=0; i < inputs.size(); ++i) {
199 AddIndexed(inputs[i], i);
200 }
201 }
202
203
204 /// Iterate over a range and return an instance of a KahanSum.
205 ///
206 /// See Add(Iterator,Iterator) for details.
207 /// \param[in] begin Beginning of a range.
208 /// \param[in] end End of the range.
209 /// \param[in] initialValue Optional initial value.
210 template <class Iterator>
211 static KahanSum<T, N> Accumulate(Iterator begin, Iterator end,
212 T initialValue = T{}) {
213 KahanSum<T, N> theSum(initialValue);
214 theSum.Add(begin, end);
215
216 return theSum;
217 }
218
219
220 /// Add `input` to the sum.
221 ///
222 /// Particularly helpful when filling from a for loop.
223 /// This function can be inlined and auto-vectorised if
224 /// the index parameter is used to enumerate *consecutive* fills.
225 /// Use Add() or Accumulate() when no index is available.
226 /// \param[in] input Value to accumulate.
227 /// \param[in] index Index of the value. Determines internal accumulator that this
228 /// value is added to. Make sure that consecutive fills have consecutive indices
229 /// to make a loop auto-vectorisable. The actual value of the index does not matter,
230 /// as long as it is consecutive.
231 void AddIndexed(T input, std::size_t index) {
232 const unsigned int i = index % N;
233 const T y = input - fCarry[i];
234 const T t = fSum[i] + y;
235 fCarry[i] = (t - fSum[i]) - y;
236 fSum[i] = t;
237 }
238
239 /// \return Compensated sum.
240 T Sum() const {
241 return std::accumulate(std::begin(fSum), std::end(fSum), 0.);
242 }
243
244 /// \return Compensated sum.
245 T Result() const {
246 return Sum();
247 }
248
249 /// \return The sum used for compensation.
250 T Carry() const {
251 return std::accumulate(std::begin(fCarry), std::end(fCarry), 0.);
252 }
253
254 /// Add `arg` into accumulator. Does not vectorise.
256 Add(arg);
257 return *this;
258 }
259
260 /// Add other KahanSum into accumulator. Does not vectorise.
261 ///
262 /// Based on KahanIncrement from:
263 /// Y. Tian, S. Tatikonda and B. Reinwald, "Scalable and Numerically Stable Descriptive Statistics in SystemML," 2012 IEEE 28th International Conference on Data Engineering, 2012, pp. 1351-1359, doi: 10.1109/ICDE.2012.12.
264 /// Note that while Tian et al. add the carry in the first step, we subtract
265 /// the carry, in accordance with the Add(Indexed) implementation(s) above.
266 /// This is purely an implementation choice that has no impact on performance.
267 ///
268 /// \note Take care when using += (and -=) to add other KahanSums into a zero-initialized
269 /// KahanSum. The operator behaves correctly in this case, but the result may be slightly
270 /// off if you expect 0 + x to yield exactly x (where 0 is the zero-initialized KahanSum
271 /// and x another KahanSum). In particular, x's carry term may get lost. This doesn't
272 /// just happen with zero-initialized KahanSums; see the SubtractWithABitTooSmallCarry
273 /// test case in the testKahan unittest for other examples. This behavior is internally
274 /// consistent: the carry also gets lost if you switch the operands and it also happens with
275 /// other KahanSum operators.
276 template<typename U, unsigned int M>
278 U corrected_arg_sum = other.Sum() - (fCarry[0] + other.Carry());
279 U sum = fSum[0] + corrected_arg_sum;
280 U correction = (sum - fSum[0]) - corrected_arg_sum;
281 fSum[0] = sum;
282 fCarry[0] = correction;
283 return *this;
284 }
285
286 /// Subtract other KahanSum. Does not vectorise.
287 ///
288 /// Based on KahanIncrement from: Tian et al., 2012 (see operator+= documentation).
289 template<typename U, unsigned int M>
291 U corrected_arg_sum = -other.Sum() - (fCarry[0] - other.Carry());
292 U sum = fSum[0] + corrected_arg_sum;
293 U correction = (sum - fSum[0]) - corrected_arg_sum;
294 fSum[0] = sum;
295 fCarry[0] = correction;
296 return *this;
297 }
298
300 {
301 return {-this->fSum[0], -this->fCarry[0]};
302 }
303
304 template<typename U, unsigned int M>
305 bool operator ==(KahanSum<U, M> const& other) const {
306 return (this->Sum() == other.Sum()) && (this->Carry() == other.Carry());
307 }
308
309 template<typename U, unsigned int M>
310 bool operator !=(KahanSum<U, M> const& other) const {
311 return !(*this == other);
312 }
313
314 private:
315 T fSum[N];
317 };
318
319 /// Add two non-vectorized KahanSums.
320 template<typename T, unsigned int N, typename U, unsigned int M>
322 KahanSum<T, N> sum(left);
323 sum += right;
324 return sum;
325 }
326
327 /// Subtract two non-vectorized KahanSums.
328 template<typename T, unsigned int N, typename U, unsigned int M>
330 KahanSum<T, N> sum(left);
331 sum -= right;
332 return sum;
333 }
334
335 } // end namespace Math
336
337} // end namespace ROOT
338
339
340#endif /* ROOT_Math_Util */
#define N
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void input
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
The Kahan summation is a compensated summation algorithm, which significantly reduces numerical error...
Definition Util.h:122
T Sum() const
Definition Util.h:240
bool operator==(KahanSum< U, M > const &other) const
Definition Util.h:305
KahanSum< T, N > operator-()
Definition Util.h:299
static KahanSum< T, N > Accumulate(Iterator begin, Iterator end, T initialValue=T{})
Iterate over a range and return an instance of a KahanSum.
Definition Util.h:211
T Result() const
Definition Util.h:245
void Add(Iterator begin, Iterator end)
Accumulate from a range denoted by iterators.
Definition Util.h:180
void Add(const Container_t &inputs)
Fill from a container that supports index access.
Definition Util.h:195
KahanSum(KahanSum< T, M > const &other)
Constructor to create a KahanSum from another KahanSum with a different number of accumulators.
Definition Util.h:157
bool operator!=(KahanSum< U, M > const &other) const
Definition Util.h:310
T Carry() const
Definition Util.h:250
void AddIndexed(T input, std::size_t index)
Add input to the sum.
Definition Util.h:231
KahanSum< T, N > & operator+=(T arg)
Add arg into accumulator. Does not vectorise.
Definition Util.h:255
KahanSum< T, N > & operator+=(const KahanSum< U, M > &other)
Add other KahanSum into accumulator.
Definition Util.h:277
KahanSum(T initialValue=T{})
Initialise the sum.
Definition Util.h:126
KahanSum(T initialSumValue, T initialCarryValue)
Initialise with a sum value and a carry value.
Definition Util.h:135
void Add(T x)
Single-element accumulation. Will not vectorise.
Definition Util.h:165
KahanSum< T, N > & operator-=(KahanSum< U, M > const &other)
Subtract other KahanSum.
Definition Util.h:290
KahanSum(Iterator sumBegin, Iterator sumEnd, Iterator carryBegin, Iterator carryEnd)
Initialise the sum with a pre-existing state.
Definition Util.h:148
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
Namespace for new Math classes and functions.
T EvalLog(T x)
safe evaluation of log(x) with a protections against negative or zero argument to the log smooth line...
Definition Util.h:64
std::string ToString(const T &val)
Utility function for conversion to strings.
Definition Util.h:50
DisplacementVector2D< CoordSystem1, U > operator+(DisplacementVector2D< CoordSystem1, U > v1, const DisplacementVector2D< CoordSystem2, U > &v2)
Addition of DisplacementVector2D vectors.
DisplacementVector2D< CoordSystem1, U > operator-(DisplacementVector2D< CoordSystem1, U > v1, DisplacementVector2D< CoordSystem2, U > const &v2)
Difference between two DisplacementVector2D vectors.
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345
double epsilon
Definition triangle.c:618