473,218 Members | 1,470 Online
Bytes | Software Development & Data Engineering Community
Post Job

Home Posts Topics Members FAQ

Join Bytes to post your question to a community of 473,218 software developers and data experts.

Implementation of Kahan sum algorithm

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: cl**@plethora.net -- 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.
Oct 3 '06 #1
9 7705
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
Oct 4 '06 #2
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: cl**@plethora.net -- 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: cl**@plethora.net -- 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.
Oct 5 '06 #3
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
es*****@acm-dot-org.invalid
--
comp.lang.c.moderated - moderation address: cl**@plethora.net -- 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.
Oct 5 '06 #4
On 03 Oct 2006 21:44:35 GMT, "asaguiar" <as******@gmail.comwrote:
>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: cl**@plethora.net -- 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.
Oct 5 '06 #5
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: cl**@plethora.net -- 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.
Oct 5 '06 #6
In article <cl****************@plethora.net"asaguiar" <as******@gmail.comwrites:
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: cl**@plethora.net -- 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.
Oct 5 '06 #7
In comp.lang.c.moderated asaguiar <as******@gmail.comwrote:
>
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: cl**@plethora.net -- 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.
Oct 5 '06 #8
dc*****@connx.com 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: cl**@plethora.net -- 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
es*****@acm-dot-org.invalid
--
comp.lang.c.moderated - moderation address: cl**@plethora.net -- 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.
Oct 7 '06 #9
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: cl**@plethora.net -- 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.
Oct 14 '06 #10

This thread has been closed and replies have been disabled. Please start a new discussion.

Similar topics

31
by: Scott Robert Ladd | last post by:
I've posted my revised C++ implementation of the Mersenne Twister at: http://www.coyotegulch.com/libcoyote/TwistedRoad/TwistedRoad.html This is "free-as-in-liberty" and "free-as-in-beer" code....
3
by: Lyn Powell | last post by:
I've translated :-) the algorithm from The Art of Computer Programming volume 3 for AVL balanced tree insertion into C. While I understand the concept behind the algorithm, the actual...
13
by: buda | last post by:
I had some spare time and decided to try to implement the standard library qsort. This is a pretty simpleminded attempt, but it seems to be working. I have to point out that efficiency is not at...
5
by: A_StClaire_ | last post by:
thoughts or criticism anyone? using System; namespace Dijkstra { public class Algorithm { private const int nodes = 10;
8
by: Jef Driesen | last post by:
I'm working on an image segmentation algorithm. An essential part of the algorithm is a graph to keep track of the connectivity between regions. At the moment I have a working implementation, but...
17
by: Soumyadip Rakshit | last post by:
Hi, Could anyone direct me to a very good 1D DCT Implementation in C/C++ My data is in sets of 8 numbers. Thanks a ton, Soumyadip.
35
by: fermineutron | last post by:
For a while now i have been "playing" with a little C program to compute the factorial of large numbers. Currently it takes aboy 1 second per 1000 multiplications, that is 25000P1000 will take...
4
by: thomasau | last post by:
Hi everyone ! Can someone point me to a good and simple implementation of the Dijkstras algorithm using C . Thanks! -Thomas
3
by: yinglcs | last post by:
Hi, Can you please recommend any good books for c++ algorithm/data structure implementation? I am looking for book which has code/ explanation for common algorithm in c++ , e.g. search, sort,...
1
isladogs
by: isladogs | last post by:
The next online meeting of the Access Europe User Group will be on Wednesday 6 Dec 2023 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). In this month's session, Mike...
0
by: veera ravala | last post by:
ServiceNow is a powerful cloud-based platform that offers a wide range of services to help organizations manage their workflows, operations, and IT services more efficiently. At its core, ServiceNow...
0
by: mar23 | last post by:
Here's the situation. I have a form called frmDiceInventory with subform called subfrmDice. The subform's control source is linked to a query called qryDiceInventory. I've been trying to pick up the...
0
by: abbasky | last post by:
### Vandf component communication method one: data sharing ​ Vandf components can achieve data exchange through data sharing, state sharing, events, and other methods. Vandf's data exchange method...
2
by: jimatqsi | last post by:
The boss wants the word "CONFIDENTIAL" overlaying certain reports. He wants it large, slanted across the page, on every page, very light gray, outlined letters, not block letters. I thought Word Art...
2
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 7 Feb 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:30 (7.30PM). In this month's session, the creator of the excellent VBE...
0
by: fareedcanada | last post by:
Hello I am trying to split number on their count. suppose i have 121314151617 (12cnt) then number should be split like 12,13,14,15,16,17 and if 11314151617 (11cnt) then should be split like...
0
by: stefan129 | last post by:
Hey forum members, I'm exploring options for SSL certificates for multiple domains. Has anyone had experience with multi-domain SSL certificates? Any recommendations on reliable providers or specific...
0
by: MeoLessi9 | last post by:
I have VirtualBox installed on Windows 11 and now I would like to install Kali on a virtual machine. However, on the official website, I see two options: "Installer images" and "Virtual machines"....

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.