Pasted

#include <stdio.h>
/* using printf */
 
 
#define KAHAN_ACCUMULATE(type, sum, compensation, added) \
  do {                           \
  type _y = (added) - compensation;          \
  type _t = sum + _y;                    \
  compensation = (_t - sum) - _y;            \
  sum = _t;                      \
  } while(0);
 
struct Accumulator
{
  double sum, compensation;
};
 
void accumulator_init(struct Accumulator* self)
{
  self->sum = self->compensation = 0.;
}
 
double accumulator_add(struct Accumulator* self, double added)
{
  KAHAN_ACCUMULATE(double, self->sum, self->compensation, added);
  return self->sum;
}
 
double accumulator_sum(const struct Accumulator* self)
{
  return self->sum;
}
 
struct AccumulatorF
{
  float sum, compensation;
};
 
void accumulatorf_init(struct AccumulatorF* self)
{
  self->sum = self->compensation = 0.f;
}
 
float accumulatorf_add(struct AccumulatorF* self, float added)
{
  KAHAN_ACCUMULATE(float, self->sum, self->compensation, added);
  return self->sum;
}
 
float accumulatorf_sum(const struct AccumulatorF* self)
{
  return self->sum;
}
 
double test_accumulation(long start, long end, long step)
{
  long i;
  struct Accumulator adder;
  accumulator_init(&adder);
  for(i = start; i != end; i += step)
    accumulator_add(&adder, 1. / (((double) i) * i));
  return accumulator_sum(&adder);
}
 
float test_accumulationf(long start, long end, long step)
{
  long i;
  struct AccumulatorF adder;
  accumulatorf_init(&adder);
  for(i = start; i != end; i += step)
    accumulatorf_add(&adder, 1.f / (((float) i) * i));
  return accumulatorf_sum(&adder);
}
 
/*
 * Demonstration of the Kahan summation algorithm
 */
int main()
{
  printf("%.8f\n", (double) test_accumulationf(1, 100000001, 1));
  printf("%.8f\n", (double) test_accumulationf(100000000, 0, -1));
  printf("%.8f\n", test_accumulation(1, 100000001, 1));
  printf("%.8f\n", test_accumulation(100000000, 0, -1));
  return 0;
}