Logo ROOT   master
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 
37 namespace 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:
234  T fSum[N];
236  };
237 
238  } // end namespace Math
239 
240 } // end namespace ROOT
241 
242 
243 #endif /* ROOT_Math_Util */
void AddIndexed(T input, std::size_t index)
Add input to the sum.
Definition: Util.h:199
Returns the available number of logical cores.
Definition: StringConv.hxx:21
double T(double x)
Definition: ChebyshevPol.h:34
fill
Definition: fit1_py.py:6
T Carry() const
Definition: Util.h:223
#define N
void Add(T x)
Single-element accumulation. Will not vectorise.
Definition: Util.h:133
void Add(Iterator begin, Iterator end)
Accumulate from a range denoted by iterators.
Definition: Util.h:148
Double_t x[n]
Definition: legend1.C:17
The Kahan summation is a compensated summation algorithm, which significantly reduces numerical error...
Definition: Util.h:122
KahanSum(T initialValue=T{})
Initialise the sum.
Definition: Util.h:126
T Sum() const
Definition: Util.h:208
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
T Result() const
Definition: Util.h:213
REAL epsilon
Definition: triangle.c:617
KahanSum< T, N > & operator+=(T arg)
Add arg into accumulator. Does not vectorise.
Definition: Util.h:228
int type
Definition: TGX11.cxx:120
Double_t y[n]
Definition: legend1.C:17
void Add(const Container_t &inputs)
Fill from a container that supports index access.
Definition: Util.h:163
Namespace for new Math classes and functions.
std::string ToString(const T &val)
Utility function for conversion to strings.
Definition: Util.h:50
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
const Int_t n
Definition: legend1.C:16
double log(double)