[math] make KahanSum::operator+= and -= consistent

math/mathcore/inc/Math/Util.h

@@ -262,31 +262,62 @@ namespace ROOT {
          return *this;
        }
 
-       /// Add `arg` into accumulator. Does not vectorise.
-       template<typename U>
-       KahanSum& operator+=(const KahanSum<U>& arg) {
-         Add(arg.Sum());
-         fCarry[0] += arg.Carry();
+       /// Add other KahanSum into accumulator. Does not vectorise.
+       ///
+       /// Based on KahanIncrement from:
+       /// Y. Tian, S. Tatikonda and B. Reinwald, "Scalable and Numerically Stable Descriptive Statistics in SystemML," 2012 IEEE 28th International Conference on Data Engineering, 2012, pp. 1351-1359, doi: 10.1109/ICDE.2012.12.
+       /// Note that while Tian et al. add the carry in the first step, we subtract
+       /// the carry, in accordance with the Add(Indexed) implementation(s) above.
+       /// This is purely an implementation choice that has no impact on performance.
+       template<typename U, unsigned int M>
+       KahanSum<T, N>& operator+=(const KahanSum<U, M>& other) {
+         U corrected_arg_sum = other.Sum() - (fCarry[0] + other.Carry());
+         U sum = fSum[0] + corrected_arg_sum;
+         U correction = (sum - fSum[0]) - corrected_arg_sum;
+         fSum[0] = sum;
+         fCarry[0] = correction;
          return *this;
        }
 
        /// Subtract other KahanSum. Does not vectorise.
        ///
-       /// This is only meaningful when both the sum and carry of each operand are of similar order of magnitude.
-       KahanSum<T, N>& operator-=(KahanSum<T, N> const& other) {
-         fSum[0]   -= other.Sum();
-         fCarry[0] -= other.Carry();
-         // add zero such that if the summed carry is large enough to be partly represented in the sum,
-         // it will happen
-         Add(T{});
+       /// Based on KahanIncrement from: Tian et al., 2012 (see operator+= documentation).
+       template<typename U, unsigned int M>
+       KahanSum<T, N>& operator-=(KahanSum<U, M> const& other) {
+         U corrected_arg_sum = -other.Sum() - (fCarry[0] - other.Carry());
+         U sum = fSum[0] + corrected_arg_sum;
+         U correction = (sum - fSum[0]) - corrected_arg_sum;
+         fSum[0] = sum;
+         fCarry[0] = correction;
          return *this;
        }
 
+       KahanSum<T, N> operator-()
+       {
+          return {-this->fSum[0], -this->fCarry[0]};
+       }
+
      private:
        T fSum[N];
        T fCarry[N];
    };
 
+   /// Add two non-vectorized KahanSums.
+   template<typename T, unsigned int N, typename U, unsigned int M>
+   KahanSum<T, N> operator+(const KahanSum<T, N>& left, const KahanSum<U, M>& right) {
+      KahanSum<T, N> sum(left);
+      sum += right;
+      return sum;
+   }
+
+   /// Subtract two non-vectorized KahanSums.
+   template<typename T, unsigned int N, typename U, unsigned int M>
+   KahanSum<T, N> operator-(const KahanSum<T, N>& left, const KahanSum<U, M>& right) {
+      KahanSum<T, N> sum(left);
+      sum -= right;
+      return sum;
+   }
+
    } // end namespace Math
 
 } // end namespace ROOT