KahanSum.h

Go to the documentation of this file.

#ifndef AMPGEN_KAHANSUM_H
#define AMPGEN_KAHANSUM_H
 
#ifdef __FAST_MATH__
#error Compensated summation is unsafe with -ffast-math (/fp:fast)
#endif
#include "AmpGen/simd/utils.h"
 
namespace AmpGen {
  template <typename T> struct KahanSum {
    T sum = 0.;
    T cor = 0.;
    KahanSum &operator+=(const T &var) {
      T t = sum + var;
      if constexpr(utils::is_vector_type<T>::value) {
        cor += select(abs(sum) >= abs(var), (sum - t) + var, (var - t) + sum);
      } else {
        cor += std::abs(sum) >= std::abs(var) ? (sum - t) + var : (var - t) + sum;
      }
      sum = t;
      return *this;
    }
    T operator()() const { return sum + cor; }
  };
  template <typename T> KahanSum<T> operator+(const KahanSum<T> &l, const T &var) {
    KahanSum<T> rt;
    T y = var - l.cor;
    T t = l.sum + y;
    rt.cor = (t - l.sum) - y;
    rt.sum = t;
    return rt;
  }
 
  template <typename T> T KahanBabushkaKleinSum(const std::vector<T> &container) {
    T sum = 0.0;
    T cs = 0.0;
    T ccs = 0.0;
    T c = 0.0;
    T cc = 0.0;
 
    for(auto &input : container) {
      T t = sum + input;
      if(std::fabs(sum) >= std::fabs(input)) {
        c = (sum - t) + input;
      } else
        c = (input - t) + sum;
      sum = t;
      t = cs + c;
      if(std::fabs(cs) >= std::fabs(c))
        cc = (cs - t) + c;
      else
        cc = (c - t) + cs;
      cs = t;
      ccs = ccs + cc;
    }
    return sum + cs + ccs;
  }
}
 
#endif