Add a class to perform compensate summation

math/mathcore/inc/Math/Math.h

@@ -48,49 +48,48 @@
 
 /**
    \namespace ROOT
-   Namespace for new ROOT classes and functions 
+   Namespace for new ROOT classes and functions
  */
 
 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;
    }
 
    /**
@@ -139,7 +138,6 @@ inline double expm1( double x) {
 #endif
 }
 
-
    } // end namespace Math
 
 } // end namespace ROOT

math/mathcore/inc/Math/Util.h

@@ -18,7 +18,7 @@
 
 #include <cmath>
 #include <limits>
-
+#include <vector>
 
 // for defining unused variables in the interfaces
 //  and have still them in the documentation
@@ -29,49 +29,77 @@ 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)
-inline double EvalLog(double x) {
+   /**
+      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)
+   inline double EvalLog(double x)
+   {
    // evaluate the log
 #ifdef __CINT__
-   static const double epsilon = 2.*2.2250738585072014e-308;
+      static const double epsilon = 2. * 2.2250738585072014e-308;
 #else
-   static const double epsilon = 2.*std::numeric_limits<double>::min();
+      static const double epsilon = 2. * std::numeric_limits<double>::min();
 #endif
-   if(x<= epsilon)
-      return x/epsilon + std::log(epsilon) - 1;
-   else
-      return std::log(x);
-}
-
-}  // end namespace Util
-
+      if (x <= epsilon)
+         return x / epsilon + std::log(epsilon) - 1;
+      else
+         return std::log(x);
+   }
+
+   } // end namespace Util
+
+   template <class T>
+   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<T> &elements)
+      {
+         for (auto e : elements) this->Add(e);
+      }
+
+      static T Accumulate(const std::vector<T> &elements)
+      {
+
+         KahanSum init;
+         init.Add(elements);
+         return init.sum;
+      }
+
+      T Result() { return sum; }
+
+   private:
+      T sum{};
+      T correction{};
+   };
 
    } // end namespace Math
 
 } // end namespace ROOT
 
 
-
 #endif /* ROOT_Math_Util */

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

math/mathcore/test/testKahan.cxx

@@ -0,0 +1,17 @@
+#include "Math/Util.h"
+#include <vector>
+
+int KahanTest()
+{
+   std::vector<double> numbers = {0.01, 0.001, 0.0001, 0.000001, 0.00000000001};
+   ROOT::Math::KahanSum<double> k;
+   k.Add(numbers);
+   auto result = ROOT::Math::KahanSum<double>::Accumulate(numbers);
+   if (k.Result() != result) return 1;
+   return 0;
+}
+
+int main()
+{
+   return KahanTest();
+}