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 /// Single-element accumulation. Will not vectorise.
133 void Add(T x) {
134 auto y = x - fCarry[0];
135 auto t = fSum[0] + y;
136 fCarry[0] = (t - fSum[0]) - y;
137 fSum[0] = t;
138 }
139
140
141 /// Accumulate from a range denoted by iterators.
142 ///
143 /// This function will auto-vectorise with random-access iterators.
144 /// \param[in] begin Beginning of a range. Needs to be a random access iterator for automatic
145 /// vectorisation, because a contiguous block of memory needs to be read.
146 /// \param[in] end End of the range.
147 template <class Iterator>
148 void Add(Iterator begin, Iterator end) {
149 static_assert(std::is_floating_point<
150 typename std::remove_reference<decltype(*begin)>::type>::value,
151 "Iterator needs to point to floating-point values.");
152 const std::size_t n = std::distance(begin, end);
153
154 for (std::size_t i=0; i<n; ++i) {
155 AddIndexed(*(begin++), i);
156 }
157 }
158
159
160 /// Fill from a container that supports index access.
161 /// \param[in] inputs Container with index access such as std::vector or array.
162 template<class Container_t>
163 void Add(const Container_t& inputs) {
164 static_assert(std::is_floating_point<typename Container_t::value_type>::value,
165 "Container does not hold floating-point values.");
166 for (std::size_t i=0; i < inputs.size(); ++i) {
167 AddIndexed(inputs[i], i);
168 }
169 }
170
171
172 /// Iterate over a range and return an instance of a KahanSum.
173 ///
174 /// See Add(Iterator,Iterator) for details.
175 /// \param[in] begin Beginning of a range.
176 /// \param[in] end End of the range.
177 /// \param[in] initialValue Optional initial value.
178 template <class Iterator>
179 static KahanSum<T, N> Accumulate(Iterator begin, Iterator end,
180 T initialValue = T{}) {
181 KahanSum<T, N> theSum(initialValue);
182 theSum.Add(begin, end);
183
184 return theSum;
185 }
186
187
188 /// Add `input` to the sum.
189 ///
190 /// Particularly helpful when filling from a for loop.
191 /// This function can be inlined and auto-vectorised if
192 /// the index parameter is used to enumerate *consecutive* fills.
193 /// Use Add() or Accumulate() when no index is available.
194 /// \param[in] input Value to accumulate.
195 /// \param[in] index Index of the value. Determines internal accumulator that this
196 /// value is added to. Make sure that consecutive fills have consecutive indices
197 /// to make a loop auto-vectorisable. The actual value of the index does not matter,
198 /// as long as it is consecutive.
199 void AddIndexed(T input, std::size_t index) {
200 const unsigned int i = index % N;
201 const T y = input - fCarry[i];
202 const T t = fSum[i] + y;
203 fCarry[i] = (t - fSum[i]) - y;
204 fSum[i] = t;
205 }
206
207 /// \return Compensated sum.
208 T Sum() const {
209 return std::accumulate(std::begin(fSum), std::end(fSum), 0.);
210 }
211
212 /// \return Compensated sum.
213 T Result() const {
214 return Sum();
215 }
216
217 /// Auto-convert to type T
218 operator T() const {
219 return Sum();
220 }
221
222 /// \return The sum used for compensation.
223 T Carry() const {
224 return std::accumulate(std::begin(fCarry), std::end(fCarry), 0.);
225 }
226
227 /// Add `arg` into accumulator. Does not vectorise.
229 Add(arg);
230 return *this;
231 }
232
233 private:
236 };
237
238 } // end namespace Math
239
240} // end namespace ROOT
241
242
243#endif /* ROOT_Math_Util */
#define N
int type
Definition: TGX11.cxx:120
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:208
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:179
T Result() const
Definition: Util.h:213
void Add(Iterator begin, Iterator end)
Accumulate from a range denoted by iterators.
Definition: Util.h:148
void Add(const Container_t &inputs)
Fill from a container that supports index access.
Definition: Util.h:163
T Carry() const
Definition: Util.h:223
void AddIndexed(T input, std::size_t index)
Add input to the sum.
Definition: Util.h:199
KahanSum< T, N > & operator+=(T arg)
Add arg into accumulator. Does not vectorise.
Definition: Util.h:228
KahanSum(T initialValue=T{})
Initialise the sum.
Definition: Util.h:126
void Add(T x)
Single-element accumulation. Will not vectorise.
Definition: Util.h:133
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
VSD Structures.
Definition: StringConv.hxx:21
fill
Definition: fit1_py.py:6
REAL epsilon
Definition: triangle.c:617