1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
|
// @(#)root/mathcore:$Id$
// Author: L. Moneta Tue Nov 14 15:44:38 2006
/**********************************************************************
* *
* Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
* *
* *
**********************************************************************/
// Utility functions for all ROOT Math classes
#ifndef ROOT_Math_Util
#define ROOT_Math_Util
#include <string>
#include <sstream>
#include <cmath>
#include <limits>
#include "Types.h"
// for defining unused variables in the interfaces
// and have still them in the documentation
#define MATH_UNUSED(var) (void)var
namespace ROOT {
namespace Math {
/**
namespace defining Utility functions needed by mathcore
*/
namespace Util {
/**
Utility function for conversion to strings
*/
template <class T>
std::string ToString(const T &val)
{
std::ostringstream buf;
buf << val;
std::string ret = buf.str();
return ret;
}
/// safe evaluation of log(x) with a protections against negative or zero argument to the log
/// smooth linear extrapolation below function values smaller than epsilon
/// (better than a simple cut-off)
template<class T>
inline T EvalLog(T x) {
static const T epsilon = T(2.0 * std::numeric_limits<double>::min());
#if !defined(R__HAS_VECCORE)
T logval = x <= epsilon ? x / epsilon + std::log(epsilon) - T(1.0) : std::log(x);
#else
T logval = vecCore::Blend<T>(x <= epsilon, x / epsilon + std::log(epsilon) - T(1.0), std::log(x));
#endif
return logval;
}
} // end namespace Util
///\class KahanSum
/// The Kahan compensate summation algorithm significantly reduces the numerical error in the total obtained
/// by adding a sequence of finite precision floating point numbers.
/// This is done by keeping a separate running compensation (a variable to accumulate small errors).\n
///
/// The intial values of the result and the correction are set to the default value of the type it hass been
/// instantiated with.\n
/// ####Examples:
/// ~~~{.cpp}
/// std::vector<double> numbers = {0.01, 0.001, 0.0001, 0.000001, 0.00000000001};
/// ROOT::Math::KahanSum<double> k;
/// k.Add(numbers);
/// ~~~
/// ~~~{.cpp}
/// auto result = ROOT::Math::KahanSum<double>::Accumulate(numbers);
/// ~~~
template <class T>
class KahanSum {
public:
/// Constructor accepting a initial value for the summation as parameter
KahanSum(const T &initialValue = T{}) : fSum(initialValue) {}
/// Single element accumulated addition.
void Add(const T &x)
{
auto y = x - fCorrection;
auto t = fSum + y;
fCorrection = (t - fSum) - y;
fSum = t;
}
/// Iterate over a datastructure referenced by a pointer and accumulate on the exising result
template <class Iterator>
void Add(const Iterator begin, const Iterator end)
{
static_assert(!std::is_same<decltype(*begin), T>::value,
"Iterator points to an element of the different type than the KahanSum class");
for (auto it = begin; it != end; it++) this->Add(*it);
}
/// Iterate over a datastructure referenced by a pointer and return the result of its accumulation.
/// Can take an initial value as third parameter.
template <class Iterator>
static T Accumulate(const Iterator begin, const Iterator end, const T &initialValue = T{})
{
static_assert(!std::is_same<decltype(*begin), T>::value,
"Iterator points to an element of the different type than the KahanSum class");
KahanSum init(initialValue);
init.Add(begin, end);
return init.fSum;
}
/// Return the result
T Result() { return fSum; }
private:
T fSum{};
T fCorrection{};
};
} // end namespace Math
} // end namespace ROOT
#endif /* ROOT_Math_Util */
|
