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 
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  /// 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:
263  T fSum[N];
265  };
266 
267  } // end namespace Math
268 
269 } // end namespace ROOT
270 
271 
272 #endif /* ROOT_Math_Util */
n
const Int_t n
Definition: legend1.C:16
fit1_py.fill
fill
Definition: fit1_py.py:6
ROOT::Math::KahanSum::Result
T Result() const
Definition: Util.h:222
ROOT::Math::KahanSum::Add
void Add(const Container_t &inputs)
Fill from a container that supports index access.
Definition: Util.h:172
log
double log(double)
N
#define N
x
Double_t x[n]
Definition: legend1.C:17
ROOT::Math::KahanSum::AddIndexed
void AddIndexed(T input, std::size_t index)
Add input to the sum.
Definition: Util.h:208
ROOT::Math::KahanSum::operator-=
KahanSum< T, N > & operator-=(KahanSum< T, N > const &other)
Subtract other KahanSum.
Definition: Util.h:253
ROOT::Math::KahanSum::Accumulate
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
ROOT::Math::KahanSum::Add
void Add(Iterator begin, Iterator end)
Accumulate from a range denoted by iterators.
Definition: Util.h:157
ROOT::MacOSX::Util
Definition: CocoaUtils.h:22
ROOT::Math::KahanSum::operator+=
KahanSum< T, N > & operator+=(T arg)
Add arg into accumulator. Does not vectorise.
Definition: Util.h:237
ROOT::Math::KahanSum::fCarry
T fCarry[N]
Definition: Util.h:264
epsilon
REAL epsilon
Definition: triangle.c:617
ROOT::Math::KahanSum::operator+=
KahanSum & operator+=(const KahanSum< U > &arg)
Add arg into accumulator. Does not vectorise.
Definition: Util.h:244
ROOT::Math::Util::ToString
std::string ToString(const T &val)
Utility function for conversion to strings.
Definition: Util.h:50
y
Double_t y[n]
Definition: legend1.C:17
Types.h
ROOT::Math::KahanSum::Add
void Add(T x)
Single-element accumulation. Will not vectorise.
Definition: Util.h:142
ROOT::Math::Util::EvalLog
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
ROOT::Math::KahanSum::Sum
T Sum() const
Definition: Util.h:217
ROOT::Math::KahanSum::fSum
T fSum[N]
Definition: Util.h:263
ROOT::Math::Chebyshev::T
double T(double x)
Definition: ChebyshevPol.h:34
type
int type
Definition: TGX11.cxx:121
ROOT::Math::KahanSum
The Kahan summation is a compensated summation algorithm, which significantly reduces numerical error...
Definition: Util.h:122
ROOT::Math::KahanSum::KahanSum
KahanSum(T initialValue=T{})
Initialise the sum.
Definition: Util.h:126
ROOT::Math::KahanSum::Carry
T Carry() const
Definition: Util.h:232
ROOT
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: EExecutionPolicy.hxx:4
Math
Namespace for new Math classes and functions.
ROOT::Math::KahanSum::KahanSum
KahanSum(KahanSum< T, M > const &other)
Constructor to create a KahanSum from another KahanSum with a different number of accumulators.
Definition: Util.h:134