Pasted

#include <iostream>
 
/*
 * Demonstration of the Kahan summation algorithm
 */
 
int main()
{
  float s = 0.;
  float c = 0.;
  for(long i = 1; i <= 100000000; ++i) {
    float cur = 1.f / i / i;
    float y = cur - c;
    float t = s + y;
    c = (t - s) - y;
    s = t;
  }
  std::cout << s << std::endl;
  s = c = 0.;
  for(long i = 100000000; i >= 1; --i) {
    float cur = 1.f / i / i;
    float y = cur - c;
    float t = s + y;
    c = (t - s) - y;
    s = t;
  }
  std::cout << s << std::endl;
}
 
/*
 * You can use this template to use it easier
 */
 
template<class Real>
class Accumulator
{
  Real _sum, _compensation;
public:
  Accumulator()
    : _sum(0.), _compensation(0.)
  {}
  Real sum() const
  { return _sum; }
  Real operator()(Real added)
  {
    Real y = added - _compensation;
    Real t = _sum + y;
    _compensation = (t - _sum) - y;
    _sum = t;
    return _sum;
  }
  void reset(Real sum = 0.)
  {
    _sum = sum;
    _compensation = 0.;
  }
};
 
void demonstration()
{
  Accumulator<float> add;
  for(long i = 1; i <= 100000000; ++i)
    add(1.f / i / i);
  float sum = add.sum();
}