Implementation of Kahan sum algorithm

[asaguiar]

Hi,

Given the pseudocode explanation of the Kahan algorithm at
http://en.wikipedia.org/wiki/Kahan_summation_algorithm, I tried to
implement it in C. It is supposed to minimize the effect of base
conversion errors after repeated sum operations.
However, the effect is null. The rounding error persists without
change.
My implementation is:

double kahanSum(double input, double tosum, double times)
{
double c=0.0, sum=input,y,t;
int count;
for(count=0; count<times; count++)
{
y=tosum-c;
t=sum+y;
c=(t-sum)-y;
sum=t;
}
return(sum);
}

Could you, please, point where I went wrong?

Thanks.

Alexandre Aguiar, MD

 

[Thad Smith]

Given the pseudocode explanation of the Kahan algorithm at
http://en.wikipedia.org/wiki/Kahan_summation_algorithm, I tried to
implement it in C. It is supposed to minimize the effect of base
conversion errors after repeated sum operations.

No, it minimizes summation roundoff or truncation errors due to unequal
scale factors, not conversion errors. There is no base conversion being
done in the summation loop.

However, the effect is null. The rounding error persists without
change.
My implementation is:

double kahanSum(double input, double tosum, double times)
{
double c=0.0, sum=input,y,t;
int count;
for(count=0; count<times; count++)
{
y=tosum-c;
t=sum+y;
c=(t-sum)-y;
sum=t;
}
return(sum);
}

Could you, please, point where I went wrong?

That seems correct to me, although it is a little strange to use a
floating point value for the counter. Here is a full implementation of
a test program. This demonstrates, on my system, a definite retention
of precision in the summation.

/* Test of kahanSum */
#include <stdio.h>

double kahanSum(double input, double tosum, double times)
{
double c=0.0, sum=input,y,t;
int count;
for(count=0; count<times; count++)
{
y=tosum-c;
t=sum+y;
c=(t-sum)-y;
sum=t;
}
return(sum);
}

#define PI 3.1415926535897932384626

int main(void) {
double s;
long int i;

s = 0.0;
for (i=0; i < 10000; i++) {
s += PI;
}

printf ("(double)pi = %.18g\n", PI);
printf ("10000 additions of pi = %.18g, error = %.6g\n",
s, 10000.*PI-s);
printf ("10000 Kahan additions of pi = %.18g, error = %.6g\n",
kahanSum(0.0,PI, 10000), 10000.*PI-kahanSum(0.0,PI, 10000));
return 0;
}

 

[dcorbit]

Hi,

Given the pseudocode explanation of the Kahan algorithm at
http://en.wikipedia.org/wiki/Kahan_summation_algorithm, I tried to
implement it in C. It is supposed to minimize the effect of base
conversion errors after repeated sum operations.
However, the effect is null. The rounding error persists without
change.
My implementation is:

double kahanSum(double input, double tosum, double times)
{
double c=0.0, sum=input,y,t;
int count;
for(count=0; count<times; count++)
{
y=tosum-c;
t=sum+y;
c=(t-sum)-y;
sum=t;
}
return(sum);
}

Could you, please, point where I went wrong?

Thanks.

Alexandre Aguiar, MD

typedef struct KahanAdder_t {
double sum_; /* The current working sum. */
double carry_; /* The carry from the previous
operation */
double temp_; /* A temporary used to capture residual
bits of precision */
double y_; /* A temporary used to capture residual
bits of precision */
} KahanAdder_t;

/* After each add operation, the carry will contain the additional */
/* bits that could be left out from the previous operation. */
void add(const double datum, KahanAdder_t * kp)
{
kp->y_ = datum - kp->carry_;
kp->temp_ = kp->sum_ + kp->y_;
kp->carry_ = (kp->temp_ - kp->sum_) - kp->y_;
kp->sum_ = kp->temp_;
}

#ifdef UNIT_TEST
#include <stdio.h>
#include <stdlib.h>
int main(void)
{
KahanAdder_t k = {0};
double d;
double standard_sum = 0;
size_t s;

for (s = 0; s < 10000; s++) {
d = rand() / (rand() * 1.0 + 1.0);
add(d, &k);
standard_sum += d;
}
printf("Standard sum = %20.15f, Kahan sum = %20.15f\n",
standard_sum, k.sum_);
return 0;
}
#endif


If you can use C++, here is a Kahan Adder as a template:
http://cap.connx.com/tournament_software/Kahan.Hpp

 

[Eric Sosman]

Hi,

Given the pseudocode explanation of the Kahan algorithm at
http://en.wikipedia.org/wiki/Kahan_summation_algorithm, I tried to
implement it in C. It is supposed to minimize the effect of base
conversion errors after repeated sum operations.
However, the effect is null. The rounding error persists without
change.
My implementation is:

double kahanSum(double input, double tosum, double times)
{
double c=0.0, sum=input,y,t;
int count;
for(count=0; count<times; count++)
{
y=tosum-c;
t=sum+y;
c=(t-sum)-y;
sum=t;
}
return(sum);
}

Could you, please, point where I went wrong?

It would help if you described what this function was
intended to do, and why you think the result is "wrong."
For what it's worth, I think you may have confused `sum'
and `tosum' in one or perhaps two places -- but I can't be
sure, as I don't know what the function is supposed to do.

 

[Barry Schwarz]

Hi,

Given the pseudocode explanation of the Kahan algorithm at
http://en.wikipedia.org/wiki/Kahan_summation_algorithm, I tried to
implement it in C. It is supposed to minimize the effect of base
conversion errors after repeated sum operations.
However, the effect is null. The rounding error persists without
change.

What rounding error? What was your input to the function? What
output were you expecting? What output did you get?

My implementation is:

double kahanSum(double input, double tosum, double times)

The algorithm in wikipedia is for summing several different values.
You are summing the value of tosum repeatedly plus the value of input
once. That's OK if it's deliberate.

{
double c=0.0, sum=input,y,t;
int count;
for(count=0; count<times; count++)

This loop iterates one time too many compared to the algorithm.

{
y=tosum-c;
t=sum+y;

The algorithm, or at least the sample, assumes that this computation
is rounded.

c=(t-sum)-y;
sum=t;
}
return(sum);
}

Did you turn off optimization for your compiler based on the warnings
in the article?


Remove del for email

 

[Tim Prince]

Hi,

Given the pseudocode explanation of the Kahan algorithm at
http://en.wikipedia.org/wiki/Kahan_summation_algorithm, I tried to

double kahanSum(double input, double tosum, double times)
{
double c=0.0, sum=input,y,t;
int count;
for(count=0; count<times; count++)
{
y=tosum-c;
t=sum+y;
c=(t-sum)-y;
sum=t;
}
return(sum);
}

Could you, please, point where I went wrong?

Did you miss the purpose of the algorithm, to sum the elements of an
array, while minimizing the numerical roundoff effect of the order in
which they are taken? It's a bit strange that you choose double as the
type of the number of elements in tosum[], but totally disabling that
you use just one element, rather than tosum[count], as shown in your
reference. Beyond that, you don't reveal whether you took into account
the caution in that reference about using a strict standard mode to obey
the parentheses, or the necessity mentioned elsewhere of ensuring that t
is rounded to the target precision.
Depending on your choice of OS and compiler, either requirement may be
satisfied by default, or may require an appropriate option. Needless to
say, it also implies disabling any parallel batching optimization.

 

[Dik T. Winter]

> Given the pseudocode explanation of the Kahan algorithm at
> http://en.wikipedia.org/wiki/Kahan_summation_algorithm, I tried to
> implement it in C. It is supposed to minimize the effect of base
> conversion errors after repeated sum operations.
> However, the effect is null. The rounding error persists without
> change.

See the paragraph on that page starting with:
"The algorithm's execution can also be affected by..."

 

[lawrence.jones]

Given the pseudocode explanation of the Kahan algorithm at
http://en.wikipedia.org/wiki/Kahan_summation_algorithm, I tried to
implement it in C. It is supposed to minimize the effect of base
conversion errors after repeated sum operations.
However, the effect is null. The rounding error persists without
change. [...]
Could you, please, point where I went wrong?

Most likely, you didn't read the stuff at the bottom of that page under
"Caution! Beware Optimising Compilers!".

-Larry Jones

My dreams are getting way too literal. -- Calvin

 

[Eric Sosman]

typedef struct KahanAdder_t {
double sum_; /* The current working sum. */
double carry_; /* The carry from the previous
operation */
double temp_; /* A temporary used to capture residual
bits of precision */
double y_; /* A temporary used to capture residual
bits of precision */

Why store the two temporaries in the struct instead of just
using a pair of `auto' variables in add()?

} KahanAdder_t;

/* After each add operation, the carry will contain the additional */
/* bits that could be left out from the previous operation. */
void add(const double datum, KahanAdder_t * kp)
{
kp->y_ = datum - kp->carry_;
kp->temp_ = kp->sum_ + kp->y_;
kp->carry_ = (kp->temp_ - kp->sum_) - kp->y_;
kp->sum_ = kp->temp_;
}

#ifdef UNIT_TEST
#include <stdio.h>
#include <stdlib.h>
int main(void)
{
KahanAdder_t k = {0};
double d;
double standard_sum = 0;
size_t s;

for (s = 0; s < 10000; s++) {
d = rand() / (rand() * 1.0 + 1.0);
add(d, &k);
standard_sum += d;
}
printf("Standard sum = %20.15f, Kahan sum = %20.15f\n",
standard_sum, k.sum_);
return 0;
}
#endif

It might be better to compute a sum whose value is more
predictable, and subject to independent verification. With
the test above you can observe that the naive and Kahan sums
are different and by how much, but you can't tell whether the
Kahan sum is "significantly better." (What *should* the sum
have turned out to be?)

To test my version (Java, so I won't post it here), I used
the harmonic series 1/1 + 1/2 + 1/3 + ... But now that I
think of it, it might have been better to compute a geometric
series 1 + r + r^2 + ..., for two reasons: First, it forces the
algorithm to add small numbers to large numbers, where the naive
approach is most liable to lose precision and thus where Kahan's
method shines. Second, the sum of the first n terms is easy to
compute (if pow() is trustworthy) as (r^(n+1) - 1) / (r - 1), so
there's a way to assess the accuracy of the answer.

 

[Elan Riesman]

Dear Dr. Aguiar,

Your implementation is correct. The problem could be any or all of
three things: the inputs, the way the result is printed, or optimizing
away the roundoff errors in what I call simpleSum below. The following
implementation of simpleSum enforces rounding after each step. Try
printf ("%.18g %.18g", simpleSum (1e15, .1, 1000), kahanSum (1e15,
..1, 1000)). The %.18g displays the full precision of a double. With
this example, I found simpleSum erroneously added 125 where kahanSum
added 100. Hope this helps.

Sincerely,
Elan Riesman
Acton, Massachusetts



void addDouble (double *c, double *a, double *b) {*c = *a + *b;}

double simpleSum(double input, double tosum, double times)
{
while (times-- > 0) addDouble (&input, &input, &tosum);
return input;
}