Velocity Reviews - Computer Hardware Reviews

Velocity Reviews > Newsgroups > Programming > C Programming > Implementation of Kahan sum algorithm

Reply
Thread Tools

Implementation of Kahan sum algorithm

 
 
asaguiar
Guest
Posts: n/a
 
      10-03-2006
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
--
comp.lang.c.moderated - moderation address: http://www.velocityreviews.com/forums/(E-Mail Removed) -- you must
have an appropriate newsgroups line in your header for your mail to be seen,
or the newsgroup name in square brackets in the subject line. Sorry.
 
Reply With Quote
 
 
 
 
Thad Smith
Guest
Posts: n/a
 
      10-04-2006
asaguiar wrote:

> 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;
}

--
Thad
 
Reply With Quote
 
 
 
 
dcorbit@connx.com
Guest
Posts: n/a
 
      10-05-2006
asaguiar wrote:
> 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
> --
> comp.lang.c.moderated - moderation address: (E-Mail Removed) -- you must
> have an appropriate newsgroups line in your header for your mail to be seen,
> or the newsgroup name in square brackets in the subject line. Sorry.


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
--
comp.lang.c.moderated - moderation address: (E-Mail Removed) -- you must
have an appropriate newsgroups line in your header for your mail to be seen,
or the newsgroup name in square brackets in the subject line. Sorry.
 
Reply With Quote
 
Eric Sosman
Guest
Posts: n/a
 
      10-05-2006
asaguiar wrote:
> 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.

--
Eric Sosman
(E-Mail Removed)lid
--
comp.lang.c.moderated - moderation address: (E-Mail Removed) -- you must
have an appropriate newsgroups line in your header for your mail to be seen,
or the newsgroup name in square brackets in the subject line. Sorry.
 
Reply With Quote
 
Barry Schwarz
Guest
Posts: n/a
 
      10-05-2006
On 03 Oct 2006 21:44:35 GMT, "asaguiar" <(E-Mail Removed)> wrote:

>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
--
comp.lang.c.moderated - moderation address: (E-Mail Removed) -- you must
have an appropriate newsgroups line in your header for your mail to be seen,
or the newsgroup name in square brackets in the subject line. Sorry.
 
Reply With Quote
 
Tim Prince
Guest
Posts: n/a
 
      10-05-2006
asaguiar wrote:
> 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.
--
comp.lang.c.moderated - moderation address: (E-Mail Removed) -- you must
have an appropriate newsgroups line in your header for your mail to be seen,
or the newsgroup name in square brackets in the subject line. Sorry.
 
Reply With Quote
 
Dik T. Winter
Guest
Posts: n/a
 
      10-05-2006
In article <(E-Mail Removed)> "asaguiar" <(E-Mail Removed)> writes:
> 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..."
--
dik t. winter, cwi, kruislaan 413, 1098 sj amsterdam, nederland, +31205924131
home: bovenover 215, 1025 jn amsterdam, nederland; http://www.cwi.nl/~dik/
--
comp.lang.c.moderated - moderation address: (E-Mail Removed) -- you must
have an appropriate newsgroups line in your header for your mail to be seen,
or the newsgroup name in square brackets in the subject line. Sorry.
 
Reply With Quote
 
lawrence.jones@ugs.com
Guest
Posts: n/a
 
      10-05-2006
In comp.lang.c.moderated asaguiar <(E-Mail Removed)> wrote:
>
> 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
--
comp.lang.c.moderated - moderation address: (E-Mail Removed) -- you must
have an appropriate newsgroups line in your header for your mail to be seen,
or the newsgroup name in square brackets in the subject line. Sorry.
 
Reply With Quote
 
Eric Sosman
Guest
Posts: n/a
 
      10-07-2006
(E-Mail Removed) wrote:
> asaguiar wrote:
>
>>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
>>--
>>comp.lang.c.moderated - moderation address: (E-Mail Removed) -- you must
>>have an appropriate newsgroups line in your header for your mail to be seen,
>>or the newsgroup name in square brackets in the subject line. Sorry.

>
>
> 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.

--
Eric Sosman
(E-Mail Removed)lid
--
comp.lang.c.moderated - moderation address: (E-Mail Removed) -- you must
have an appropriate newsgroups line in your header for your mail to be seen,
or the newsgroup name in square brackets in the subject line. Sorry.
 
Reply With Quote
 
Elan Riesman
Guest
Posts: n/a
 
      10-14-2006
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;
}

> 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?

--
comp.lang.c.moderated - moderation address: (E-Mail Removed) -- you must
have an appropriate newsgroups line in your header for your mail to be seen,
or the newsgroup name in square brackets in the subject line. Sorry.
 
Reply With Quote
 
 
 
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are Off


Similar Threads
Thread Thread Starter Forum Replies Last Post
Feature suggestion: sum() ought to use a compensated summation algorithm Szabolcs Horvát Python 25 05-09-2008 04:12 PM
[Algorithm] Sum of Primes < 1000000 Luc The Perverse Java 50 01-18-2007 11:18 PM
Filtered Back Projection Algorithm (FBP Algorithm) Bapaiah Katepalli VHDL 1 06-23-2006 04:50 PM
Hardware implementation of Safer+ algorithm blocks 'e', 'l' Rag VHDL 0 01-30-2006 08:36 PM
Key generation algorithm and Cipher algorithm Ahmed Moustafa Java 0 11-15-2003 06:35 AM



Advertisments