From a334562b2692be16593e9d17db582a4d1662643b Mon Sep 17 00:00:00 2001 From: Xavier Valls Date: Mon, 27 Feb 2017 14:17:26 +0100 Subject: [PATCH 1/6] Add a class to perform the compensate summation Kahan summation algorithm: https://en.wikipedia.org/wiki/Kahan_summation_algorithm --- math/mathcore/inc/Math/Math.h | 38 ++++++++++++++++++++++++++++++++--- 1 file changed, 35 insertions(+), 3 deletions(-) diff --git a/math/mathcore/inc/Math/Math.h b/math/mathcore/inc/Math/Math.h index a5b0645cd5956..a408ba8eb8834 100644 --- a/math/mathcore/inc/Math/Math.h +++ b/math/mathcore/inc/Math/Math.h @@ -48,7 +48,7 @@ /** \namespace ROOT - Namespace for new ROOT classes and functions + Namespace for new ROOT classes and functions */ namespace ROOT { @@ -56,7 +56,7 @@ namespace ROOT { /** \namespace Math Namespace for new Math classes and functions. - See the \ref Math "Math Libraries" page for a detailed description. + See the \ref Math "Math Libraries" page for a detailed description. */ @@ -139,7 +139,39 @@ inline double expm1( double x) { #endif } - + template + class KahanSum { + public: + void Add(const T &x) + { + auto y = x - correction; + auto t = sum + y; + correction = (t - sum) - y; + sum = t; + } + + void Add(const std::vector &elements) + { + for (auto e : elements) + this->Add(e); + } + + static double Accumulate(const std::vector &elements) + { + KahanSum init; + init.Add(elements); + return init.sum; + } + + double Result() + { + return sum; + } + + private: + double sum = 0.; + double correction = 0.; + }; } // end namespace Math } // end namespace ROOT From 4bbf52a39c63884f197adea3cb61b251393f3e06 Mon Sep 17 00:00:00 2001 From: Xavier Valls Date: Mon, 27 Feb 2017 15:52:23 +0100 Subject: [PATCH 2/6] Move KahanSum to Math/Util.h --- math/mathcore/inc/Math/Math.h | 33 --------------------------------- math/mathcore/inc/Math/Util.h | 34 ++++++++++++++++++++++++++++++++++ 2 files changed, 34 insertions(+), 33 deletions(-) diff --git a/math/mathcore/inc/Math/Math.h b/math/mathcore/inc/Math/Math.h index a408ba8eb8834..c31179c704332 100644 --- a/math/mathcore/inc/Math/Math.h +++ b/math/mathcore/inc/Math/Math.h @@ -139,39 +139,6 @@ inline double expm1( double x) { #endif } - template - class KahanSum { - public: - void Add(const T &x) - { - auto y = x - correction; - auto t = sum + y; - correction = (t - sum) - y; - sum = t; - } - - void Add(const std::vector &elements) - { - for (auto e : elements) - this->Add(e); - } - - static double Accumulate(const std::vector &elements) - { - KahanSum init; - init.Add(elements); - return init.sum; - } - - double Result() - { - return sum; - } - - private: - double sum = 0.; - double correction = 0.; - }; } // end namespace Math } // end namespace ROOT diff --git a/math/mathcore/inc/Math/Util.h b/math/mathcore/inc/Math/Util.h index ad803d2af30e5..e1bda26f31c1b 100644 --- a/math/mathcore/inc/Math/Util.h +++ b/math/mathcore/inc/Math/Util.h @@ -18,6 +18,7 @@ #include #include +#include // for defining unused variables in the interfaces @@ -67,6 +68,39 @@ inline double EvalLog(double x) { } // end namespace Util + template + class KahanSum { + public: + void Add(const T &x) + { + auto y = x - correction; + auto t = sum + y; + correction = (t - sum) - y; + sum = t; + } + + void Add(const std::vector &elements) + { + for (auto e : elements) + this->Add(e); + } + + static double Accumulate(const std::vector &elements) + { + KahanSum init; + init.Add(elements); + return init.sum; + } + + double Result() + { + return sum; + } + + private: + double sum = 0.; + double correction = 0.; + }; } // end namespace Math From 5f819428ca47241bbe75afab7759c01091e15770 Mon Sep 17 00:00:00 2001 From: Xavier Valls Date: Mon, 27 Feb 2017 15:53:46 +0100 Subject: [PATCH 3/6] Fix indentation :art: --- math/mathcore/inc/Math/Util.h | 74 +++++++++++++++++------------------ 1 file changed, 37 insertions(+), 37 deletions(-) diff --git a/math/mathcore/inc/Math/Util.h b/math/mathcore/inc/Math/Util.h index e1bda26f31c1b..b5ac3db3952d9 100644 --- a/math/mathcore/inc/Math/Util.h +++ b/math/mathcore/inc/Math/Util.h @@ -31,42 +31,43 @@ namespace ROOT { namespace Math { -/** - namespace defining Utility functions needed by mathcore -*/ -namespace Util { - -/** - Utility function for conversion to strings -*/ -template -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) -inline double EvalLog(double x) { - // evaluate the log -#ifdef __CINT__ - static const double epsilon = 2.*2.2250738585072014e-308; -#else - static const double epsilon = 2.*std::numeric_limits::min(); -#endif - if(x<= epsilon) - return x/epsilon + std::log(epsilon) - 1; - else - return std::log(x); -} - -} // end namespace Util + /** + namespace defining Utility functions needed by mathcore + */ + namespace Util { + + /** + Utility function for conversion to strings + */ + template + 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) + inline double EvalLog(double x) + { + // evaluate the log + #ifdef __CINT__ + static const double epsilon = 2.*2.2250738585072014e-308; + #else + static const double epsilon = 2.*std::numeric_limits::min(); + #endif + if (x <= epsilon) + return x / epsilon + std::log(epsilon) - 1; + else + return std::log(x); + } + + } // end namespace Util template class KahanSum { @@ -107,5 +108,4 @@ inline double EvalLog(double x) { } // end namespace ROOT - #endif /* ROOT_Math_Util */ From ce1c9adfe25da7c33015bc00e40d8179804cf877 Mon Sep 17 00:00:00 2001 From: Xavier Valls Date: Tue, 28 Feb 2017 17:22:31 +0100 Subject: [PATCH 4/6] Add test for the compensate summation --- math/mathcore/test/CMakeLists.txt | 9 +++++---- math/mathcore/test/testKahan.cxx | 17 +++++++++++++++++ 2 files changed, 22 insertions(+), 4 deletions(-) create mode 100644 math/mathcore/test/testKahan.cxx diff --git a/math/mathcore/test/CMakeLists.txt b/math/mathcore/test/CMakeLists.txt index e1d396bbdada0..4cdacd23120ad 100644 --- a/math/mathcore/test/CMakeLists.txt +++ b/math/mathcore/test/CMakeLists.txt @@ -16,19 +16,20 @@ set(TestSource testIntegration.cxx testRootFinder.cxx testSampleQuantiles.cxx - kDTreeTest.cxx + kDTreeTest.cxx testkdTreeBinning.cxx - newKDTreeTest.cxx + newKDTreeTest.cxx binarySearchTime.cxx stdsort.cxx testSpecFuncErf.cxx testSpecFuncGamma.cxx testSpecFuncBeta.cxx - testSpecFuncBetaI.cxx - testSpecFuncSiCi.cxx + testSpecFuncBetaI.cxx + testSpecFuncSiCi.cxx testIntegrationMultiDim.cxx testAnalyticalIntegrals.cxx testTStatistic.cxx + testKahan.cxx fit/testFit.cxx fit/testGraphFit.cxx fit/SparseDataComparer.cxx diff --git a/math/mathcore/test/testKahan.cxx b/math/mathcore/test/testKahan.cxx new file mode 100644 index 0000000000000..493119bc44523 --- /dev/null +++ b/math/mathcore/test/testKahan.cxx @@ -0,0 +1,17 @@ +#include"Math/Util.h" +#include + +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); + if(k.Result()!=result) + return 1; + return 0; +} + +int main(){ + return KahanTest(); +} \ No newline at end of file From 1f0f1e663fa726bd3f6259e16382bb1d8b95d21b Mon Sep 17 00:00:00 2001 From: Xavier Valls Date: Wed, 1 Mar 2017 11:33:30 +0100 Subject: [PATCH 5/6] Fixed the summation for non-real types --- math/mathcore/inc/Math/Util.h | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/math/mathcore/inc/Math/Util.h b/math/mathcore/inc/Math/Util.h index b5ac3db3952d9..17b0fbd865af3 100644 --- a/math/mathcore/inc/Math/Util.h +++ b/math/mathcore/inc/Math/Util.h @@ -86,21 +86,22 @@ namespace ROOT { this->Add(e); } - static double Accumulate(const std::vector &elements) + static T Accumulate(const std::vector &elements) { + KahanSum init; init.Add(elements); return init.sum; } - double Result() + T Result() { return sum; } private: - double sum = 0.; - double correction = 0.; + T sum{}; + T correction{}; }; } // end namespace Math From c3fee5b780c5354f2335754e20756ea7123f6b8b Mon Sep 17 00:00:00 2001 From: Xavier Valls Date: Thu, 6 Apr 2017 09:37:51 +0200 Subject: [PATCH 6/6] Apply clang-format :art: --- math/mathcore/inc/Math/Math.h | 73 ++++++++-------- math/mathcore/inc/Math/Util.h | 141 +++++++++++++++---------------- math/mathcore/test/testKahan.cxx | 10 +-- 3 files changed, 108 insertions(+), 116 deletions(-) diff --git a/math/mathcore/inc/Math/Math.h b/math/mathcore/inc/Math/Math.h index c31179c704332..9c354f235248c 100644 --- a/math/mathcore/inc/Math/Math.h +++ b/math/mathcore/inc/Math/Math.h @@ -53,44 +53,43 @@ namespace ROOT { - /** - \namespace Math - Namespace for new Math classes and functions. - See the \ref Math "Math Libraries" page for a detailed description. - */ - - - namespace Math { - // Enable Vc/VecCore template instantiations to replace std math functions. - // - // Vc declares `std::sqrt(Vc-type)`. To use this for Vc-`SCALAR`s, the call - // to `sqrt()` must only be resolved at the template instantiation time, when - // the Vc headers are guaranteed to be included, and thus its `sqrt()` - // overloads have been declared. - // The trick is to keep sqrt() dependent (on its argument type) by making it - // an unqualified name. The `std::` of `std::sqrt()` makes it a qualified - // name, so the code here has to use `sqrt()`, not `std::sqrt()`. To still - // find `std::sqrt()` we pull `std::sqrt()` into the surrounding namespace. - // - // We don't want to use 'using namespace std' because it would polute the including headers. - using std::atan2; - using std::cos; - using std::cosh; - using std::exp; - using std::floor; - using std::log; - using std::pow; - using std::sin; - using std::sinh; - using std::sqrt; - using std::tan; +/** +\namespace Math +Namespace for new Math classes and functions. +See the \ref Math "Math Libraries" page for a detailed description. +*/ + +namespace Math { +// Enable Vc/VecCore template instantiations to replace std math functions. +// +// Vc declares `std::sqrt(Vc-type)`. To use this for Vc-`SCALAR`s, the call +// to `sqrt()` must only be resolved at the template instantiation time, when +// the Vc headers are guaranteed to be included, and thus its `sqrt()` +// overloads have been declared. +// The trick is to keep sqrt() dependent (on its argument type) by making it +// an unqualified name. The `std::` of `std::sqrt()` makes it a qualified +// name, so the code here has to use `sqrt()`, not `std::sqrt()`. To still +// find `std::sqrt()` we pull `std::sqrt()` into the surrounding namespace. +// +// We don't want to use 'using namespace std' because it would polute the including headers. +using std::atan2; +using std::cos; +using std::cosh; +using std::exp; +using std::floor; +using std::log; +using std::pow; +using std::sin; +using std::sinh; +using std::sqrt; +using std::tan; - /** - Mathematical constants - */ - inline double Pi() - { - return M_PI; +/** + Mathematical constants +*/ +inline double Pi() +{ + return M_PI; } /** diff --git a/math/mathcore/inc/Math/Util.h b/math/mathcore/inc/Math/Util.h index 17b0fbd865af3..48d230880ca34 100644 --- a/math/mathcore/inc/Math/Util.h +++ b/math/mathcore/inc/Math/Util.h @@ -20,7 +20,6 @@ #include #include - // for defining unused variables in the interfaces // and have still them in the documentation #define MATH_UNUSED(var) (void)var @@ -30,79 +29,73 @@ namespace ROOT { namespace Math { - - /** - namespace defining Utility functions needed by mathcore - */ - namespace Util { - - /** - Utility function for conversion to strings - */ - template - 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) - inline double EvalLog(double x) - { - // evaluate the log - #ifdef __CINT__ - static const double epsilon = 2.*2.2250738585072014e-308; - #else - static const double epsilon = 2.*std::numeric_limits::min(); - #endif - if (x <= epsilon) - return x / epsilon + std::log(epsilon) - 1; - else - return std::log(x); - } - - } // end namespace Util - - template - class KahanSum { - public: - void Add(const T &x) - { - auto y = x - correction; - auto t = sum + y; - correction = (t - sum) - y; - sum = t; - } - - void Add(const std::vector &elements) - { - for (auto e : elements) - this->Add(e); - } - - static T Accumulate(const std::vector &elements) - { - - KahanSum init; - init.Add(elements); - return init.sum; - } - - T Result() - { - return sum; - } - - private: - T sum{}; - T correction{}; - }; + /** + namespace defining Utility functions needed by mathcore + */ + namespace Util { + + /** + Utility function for conversion to strings + */ + template + 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) + inline double EvalLog(double x) + { + // evaluate the log +#ifdef __CINT__ + static const double epsilon = 2. * 2.2250738585072014e-308; +#else + static const double epsilon = 2. * std::numeric_limits::min(); +#endif + if (x <= epsilon) + return x / epsilon + std::log(epsilon) - 1; + else + return std::log(x); + } + + } // end namespace Util + + template + class KahanSum { + public: + void Add(const T &x) + { + auto y = x - correction; + auto t = sum + y; + correction = (t - sum) - y; + sum = t; + } + + void Add(const std::vector &elements) + { + for (auto e : elements) this->Add(e); + } + + static T Accumulate(const std::vector &elements) + { + + KahanSum init; + init.Add(elements); + return init.sum; + } + + T Result() { return sum; } + + private: + T sum{}; + T correction{}; + }; } // end namespace Math diff --git a/math/mathcore/test/testKahan.cxx b/math/mathcore/test/testKahan.cxx index 493119bc44523..a2f51ffdb8e69 100644 --- a/math/mathcore/test/testKahan.cxx +++ b/math/mathcore/test/testKahan.cxx @@ -1,4 +1,4 @@ -#include"Math/Util.h" +#include "Math/Util.h" #include int KahanTest() @@ -7,11 +7,11 @@ int KahanTest() ROOT::Math::KahanSum k; k.Add(numbers); auto result = ROOT::Math::KahanSum::Accumulate(numbers); - if(k.Result()!=result) - return 1; + if (k.Result() != result) return 1; return 0; } -int main(){ - return KahanTest(); +int main() +{ + return KahanTest(); } \ No newline at end of file