diff --git a/math/mathcore/inc/Math/Util.h b/math/mathcore/inc/Math/Util.h index 48d230880ca34..77dc0f6ab8560 100644 --- a/math/mathcore/inc/Math/Util.h +++ b/math/mathcore/inc/Math/Util.h @@ -18,7 +18,6 @@ #include #include -#include // for defining unused variables in the interfaces // and have still them in the documentation @@ -66,35 +65,64 @@ namespace ROOT { } // 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 numbers = {0.01, 0.001, 0.0001, 0.000001, 0.00000000001}; + /// ROOT::Math::KahanSum k; + /// k.Add(numbers); + /// ~~~ + /// ~~~{.cpp} + /// auto result = ROOT::Math::KahanSum::Accumulate(numbers); + /// ~~~ template 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 - correction; - auto t = sum + y; - correction = (t - sum) - y; - sum = t; + auto y = x - fCorrection; + auto t = fSum + y; + fCorrection = (t - fSum) - y; + fSum = t; } - void Add(const std::vector &elements) + /// Iterate over a datastructure referenced by a pointer and accumulate on the exising result + template + void Add(const Iterator begin, const Iterator end) { - for (auto e : elements) this->Add(e); + static_assert(!std::is_same::value, + "Iterator points to an element of the different type than the KahanSum class"); + for (auto it = begin; it != end; it++) this->Add(*it); } - static T Accumulate(const std::vector &elements) + /// Iterate over a datastructure referenced by a pointer and return the result of its accumulation. + /// Can take an initial value as third parameter. + template + static T Accumulate(const Iterator begin, const Iterator end, const T &initialValue = T{}) { - - KahanSum init; - init.Add(elements); - return init.sum; + static_assert(!std::is_same::value, + "Iterator points to an element of the different type than the KahanSum class"); + KahanSum init(initialValue); + init.Add(begin, end); + return init.fSum; } - T Result() { return sum; } + /// Return the result + T Result() { return fSum; } private: - T sum{}; - T correction{}; + T fSum{}; + T fCorrection{}; }; } // end namespace Math diff --git a/math/mathcore/test/testKahan.cxx b/math/mathcore/test/testKahan.cxx index a2f51ffdb8e69..3521954b24da3 100644 --- a/math/mathcore/test/testKahan.cxx +++ b/math/mathcore/test/testKahan.cxx @@ -5,9 +5,17 @@ int KahanTest() { std::vector numbers = {0.01, 0.001, 0.0001, 0.000001, 0.00000000001}; ROOT::Math::KahanSum k; - k.Add(numbers); - auto result = ROOT::Math::KahanSum::Accumulate(numbers); + k.Add(numbers.begin(), numbers.end()); + auto result = ROOT::Math::KahanSum::Accumulate(numbers.begin(), numbers.end()); if (k.Result() != result) return 1; + + ROOT::Math::KahanSum k2; + ROOT::Math::KahanSum k3(1); + k2.Add(1); + k2.Add(numbers.begin(), numbers.end()); + k3.Add(numbers.begin(), numbers.end()); + if (k2.Result() != k3.Result()) return 2; + return 0; }