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 /// Constructor to create a KahanSum from another KahanSum with a different number of accumulators
133 template <unsigned int M>
134 KahanSum(KahanSum<T,M> const& other) {
135 fSum[0] = other.Sum();
136 fCarry[0] = other.Carry();
137 std::fill(std::begin(fSum)+1, std::end(fSum), 0.);
138 std::fill(std::begin(fCarry)+1, std::end(fCarry), 0.);
139 }
140
141 /// Single-element accumulation. Will not vectorise.
142 void Add(T x) {
143 auto y = x - fCarry[0];
144 auto t = fSum[0] + y;
145 fCarry[0] = (t - fSum[0]) - y;
146 fSum[0] = t;
147 }
148
149
150 /// Accumulate from a range denoted by iterators.
151 ///
152 /// This function will auto-vectorise with random-access iterators.
153 /// \param[in] begin Beginning of a range. Needs to be a random access iterator for automatic
154 /// vectorisation, because a contiguous block of memory needs to be read.
155 /// \param[in] end End of the range.
156 template <class Iterator>
157 void Add(Iterator begin, Iterator end) {
158 static_assert(std::is_floating_point<
159 typename std::remove_reference<decltype(*begin)>::type>::value,
160 "Iterator needs to point to floating-point values.");
161 const std::size_t n = std::distance(begin, end);
162
163 for (std::size_t i=0; i<n; ++i) {
164 AddIndexed(*(begin++), i);
165 }
166 }
167
168
169 /// Fill from a container that supports index access.
170 /// \param[in] inputs Container with index access such as std::vector or array.
171 template<class Container_t>
172 void Add(const Container_t& inputs) {
173 static_assert(std::is_floating_point<typename Container_t::value_type>::value,
174 "Container does not hold floating-point values.");
175 for (std::size_t i=0; i < inputs.size(); ++i) {
176 AddIndexed(inputs[i], i);
177 }
178 }
179
180
181 /// Iterate over a range and return an instance of a KahanSum.
182 ///
183 /// See Add(Iterator,Iterator) for details.
184 /// \param[in] begin Beginning of a range.
185 /// \param[in] end End of the range.
186 /// \param[in] initialValue Optional initial value.
187 template <class Iterator>
188 static KahanSum<T, N> Accumulate(Iterator begin, Iterator end,
189 T initialValue = T{}) {
190 KahanSum<T, N> theSum(initialValue);
191 theSum.Add(begin, end);
192
193 return theSum;
194 }
195
196
197 /// Add `input` to the sum.
198 ///
199 /// Particularly helpful when filling from a for loop.
200 /// This function can be inlined and auto-vectorised if
201 /// the index parameter is used to enumerate *consecutive* fills.
202 /// Use Add() or Accumulate() when no index is available.
203 /// \param[in] input Value to accumulate.
204 /// \param[in] index Index of the value. Determines internal accumulator that this
205 /// value is added to. Make sure that consecutive fills have consecutive indices
206 /// to make a loop auto-vectorisable. The actual value of the index does not matter,
207 /// as long as it is consecutive.
208 void AddIndexed(T input, std::size_t index) {
209 const unsigned int i = index % N;
210 const T y = input - fCarry[i];
211 const T t = fSum[i] + y;
212 fCarry[i] = (t - fSum[i]) - y;
213 fSum[i] = t;
214 }
215
216 /// \return Compensated sum.
217 T Sum() const {
218 return std::accumulate(std::begin(fSum), std::end(fSum), 0.);
219 }
220
221 /// \return Compensated sum.
222 T Result() const {
223 return Sum();
224 }
225
226 /// Auto-convert to type T
227 operator T() const {
228 return Sum();
229 }
230
231 /// \return The sum used for compensation.
232 T Carry() const {
233 return std::accumulate(std::begin(fCarry), std::end(fCarry), 0.);
234 }
235
236 /// Add `arg` into accumulator. Does not vectorise.
238 Add(arg);
239 return *this;
240 }
241
242 /// Add `arg` into accumulator. Does not vectorise.
243 template<typename U>
245 Add(arg.Sum());
246 fCarry[0] += arg.Carry();
247 return *this;
248 }
249
250 /// Subtract other KahanSum. Does not vectorise.
251 ///
252 /// This is only meaningfull when both the sum and carry of each operand are of similar order of magnitude.
254 fSum[0] -= other.Sum();
255 fCarry[0] -= other.Carry();
256 // add zero such that if the summed carry is large enough to be partly represented in the sum,
257 // it will happen
258 Add(T{});
259 return *this;
260 }
261
262 private:
265 };
266
267 } // end namespace Math
268
269} // end namespace ROOT
270
271
272#endif /* ROOT_Math_Util */
#define N
int type
Definition: TGX11.cxx:121
double log(double)
The Kahan summation is a compensated summation algorithm, which significantly reduces numerical error...
Definition: Util.h:122
T Sum() const
Definition: Util.h:217
KahanSum & operator+=(const KahanSum< U > &arg)
Add arg into accumulator. Does not vectorise.
Definition: Util.h:244
T Result() const
Definition: Util.h:222
void Add(Iterator begin, Iterator end)
Accumulate from a range denoted by iterators.
Definition: Util.h:157
KahanSum< T, N > & operator-=(KahanSum< T, N > const &other)
Subtract other KahanSum.
Definition: Util.h:253
void Add(const Container_t &inputs)
Fill from a container that supports index access.
Definition: Util.h:172
KahanSum< T, N > & operator+=(T arg)
Add arg into accumulator. Does not vectorise.
Definition: Util.h:237
KahanSum(KahanSum< T, M > const &other)
Constructor to create a KahanSum from another KahanSum with a different number of accumulators.
Definition: Util.h:134
T Carry() const
Definition: Util.h:232
void AddIndexed(T input, std::size_t index)
Add input to the sum.
Definition: Util.h:208
KahanSum(T initialValue=T{})
Initialise the sum.
Definition: Util.h:126
void Add(T x)
Single-element accumulation. Will not vectorise.
Definition: Util.h:142
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:188
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
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: RNumpyDS.hxx:30
fill
Definition: fit1_py.py:6
REAL epsilon
Definition: triangle.c:618