Logo ROOT  
Reference Guide
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 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) {
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 /// Auto-convert to type T
250 operator T() const {
251 return Sum();
252 }
253
254 /// \return The sum used for compensation.
255 T Carry() const {
256 return std::accumulate(std::begin(fCarry), std::end(fCarry), 0.);
257 }
258
259 /// Add `arg` into accumulator. Does not vectorise.
261 Add(arg);
262 return *this;
263 }
264
265 /// Add `arg` into accumulator. Does not vectorise.
266 template<typename U>
268 Add(arg.Sum());
269 fCarry[0] += arg.Carry();
270 return *this;
271 }
272
273 /// Subtract other KahanSum. Does not vectorise.
274 ///
275 /// This is only meaningful when both the sum and carry of each operand are of similar order of magnitude.
277 fSum[0] -= other.Sum();
278 fCarry[0] -= other.Carry();
279 // add zero such that if the summed carry is large enough to be partly represented in the sum,
280 // it will happen
281 Add(T{});
282 return *this;
283 }
284
285 private:
288 };
289
290 } // end namespace Math
291
292} // end namespace ROOT
293
294
295#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
KahanSum & operator+=(const KahanSum< U > &arg)
Add arg into accumulator. Does not vectorise.
Definition: Util.h:267
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< T, N > & operator-=(KahanSum< T, N > const &other)
Subtract other KahanSum.
Definition: Util.h:276
KahanSum(KahanSum< T, M > const &other)
Constructor to create a KahanSum from another KahanSum with a different number of accumulators.
Definition: Util.h:157
T Carry() const
Definition: Util.h:255
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:260
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(Iterator sumBegin, Iterator sumEnd, Iterator carryBegin, Iterator carryEnd)
Initialise the sum with a pre-existing state.
Definition: Util.h:148
RVec< PromoteType< T > > log(const RVec< T > &v)
Definition: RVec.hxx:1765
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.
double T(double x)
Definition: ChebyshevPol.h:34
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
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
fill
Definition: fit1_py.py:6
double epsilon
Definition: triangle.c:618