Inexact Sum of Floating Point Numbers

sum.h

// Returns s and e, such that s = T(a + b), and e = (a + b) - s, i.e. s + e = a
// + b exactly.
template <class T>
std::pair<T, T> AddPair(const T& a, const T& b) {
  auto s = a + b;
  auto b_v = s - a;
  auto e = (a - (s - b_v)) + (b - b_v);
  // auto a_v = s - b_v;
  // auto b_r = b - b_v;
  // auto a_r = a - a_v;
  // auto e = a_r + b_r;
  return std::make_pair(s, e);
}

template <class T>
constexpr int SeriesLength() {
  return (std::numeric_limits<T>::max_exponent -
          std::numeric_limits<T>::min_exponent) /
             std::numeric_limits<T>::digits +
         2;
}

template <class T>
int SeriesIndex(T v) {
  return (ilogb(v) + std::numeric_limits<T>::digits - 1 -
          std::numeric_limits<T>::min_exponent) /
         std::numeric_limits<T>::digits;
}

template <class T>
void SeriesAddAt(T* vs, T v, int i) {
  assert(v);
  assert(i >= 0);
  assert(i < SeriesLength<T>());
  if (!vs[i]) {
    vs[i] = v;
    return;
  }
  T s, e;
  std::tie(s, e) = AddPair(vs[i], v);
  if (!e) {
    vs[i] = s;
    return;
  }
  auto si = SeriesIndex(s);
  auto ei = SeriesIndex(e);
  assert(si > ei);
  if (ei == i) {
    vs[i] = e;
    SeriesAddAt(vs, s, si);
  } else {
    vs[i] = s;
    SeriesAddAt(vs, e, ei);
  }
}

template <class T>
void SeriesAdd(T* vs, T v) {
  if (!v) {
    return;
  }
  int i = SeriesIndex(v);
  SeriesAddAt(vs, v, i);
}

template <class T>
size_t SeriesNormalize(T* buf) {
  int len = SeriesLength<T>();
  while (len > 1) {
    bool changed = false;
    size_t n = 0;
    T v = buf[0];
    for (size_t i = 1; i < len; ++i) {
      if (buf[i]) {
        auto p = AddPair(buf[i], v);
        changed |= p.first != buf[i];
        if (p.second) {
          buf[n++] = p.second;
        }
        v = p.first;
      }
    }
    if (v) {
      assert(n < len);
      buf[n++] = v;
    }
    len = n;
    if (!changed) {
      break;
    }
  }
  return len;
}

template <class T, class Vs>
T SeriesSum(const Vs& vs) {
  T accs[SeriesLength<T>()] = {0};
  for (const auto& v : vs) {
    SeriesAdd(accs, v);
  }
  size_t n = SeriesNormalize(accs);
  return n > 0 ? accs[n - 1] : 0;
}

sum_test.cc

TEST(Sum) {
  size_t N = 100000;
  std::vector<double> vs(N);
  for (double& x : vs) {
    x = random() / static_cast<double>(RAND_MAX);
  }
  long double p = 0;
  for (double x : vs) {
    p += x;
  }
  for (size_t i = 0; i < 10000; ++i) {
    double x = random() *
               std::exp((random() / static_cast<double>(RAND_MAX) - 0.5) * 200);
    double y = random() *
               std::exp((random() / static_cast<double>(RAND_MAX) - 0.5) * 200);
    double s = x + y;
    double xv = s - x;
    double yv = s - xv;
    vs.push_back(-s);
    vs.push_back(xv);
    vs.push_back(yv);
  }
  std::random_shuffle(vs.begin(), vs.end());
  double s = 0;
  long double l = 0;
  for (double x : vs) {
    s += x;
    l += x;
  }
  CLOG(INFO, "pure long sum=" << p);
  CLOG(INFO, "noise double sum=" << s, "diff=" << s - (double)p);
  CLOG(INFO, "noise long sum=" << l, "diff=" << l - (double)p);
  double e = SeriesSum<double>(vs);
  CLOG(INFO, "noise series sum=" << e, "diff=" << e - (double)p);
}

test_result.log

pure long sum=49986.8133624702753615
noise double sum=-4.9485087854700077e+38 diff=-4.9485087854700077e+38
noise long sum=-4.29587139833043406612e+35 diff=-4.29587139833043406612e+35
noise series sum=49986.813362470275 diff=0