Improvements and new interfaces in KahanSum

math/mathcore/inc/Math/Util.h

@@ -18,7 +18,6 @@
 
 #include <cmath>
 #include <limits>
-#include <vector>
 
 // 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<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 - 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<T> &elements)
+      /// 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)
       {
-         for (auto e : elements) this->Add(e);
+         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);
       }
 
-      static T Accumulate(const std::vector<T> &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 <class Iterator>
+      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<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;
       }
 
-      T Result() { return sum; }
+      /// Return the result
+      T Result() { return fSum; }
 
    private:
-      T sum{};
-      T correction{};
+      T fSum{};
+      T fCorrection{};
    };
 
    } // end namespace Math

math/mathcore/test/testKahan.cxx

@@ -5,9 +5,17 @@ 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);
+   k.Add(numbers.begin(), numbers.end());
+   auto result = ROOT::Math::KahanSum<double>::Accumulate(numbers.begin(), numbers.end());
    if (k.Result() != result) return 1;
+
+   ROOT::Math::KahanSum<double> k2;
+   ROOT::Math::KahanSum<double> 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;
 }