Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #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;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement