Implementation of Kahan sum algorithm

Discussion in 'C Programming' started by asaguiar, Oct 3, 2006.

  1. asaguiar

    asaguiar Guest

    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: -- 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.
     
    asaguiar, Oct 3, 2006
    #1
    1. Advertising

  2. asaguiar

    Thad Smith Guest

    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
     
    Thad Smith, Oct 4, 2006
    #2
    1. Advertising

  3. asaguiar

    Guest

    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: -- 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: -- 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, 2006
    #3
  4. asaguiar

    Eric Sosman Guest

    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
    lid
    --
    comp.lang.c.moderated - moderation address: -- 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.
     
    Eric Sosman, Oct 5, 2006
    #4
  5. On 03 Oct 2006 21:44:35 GMT, "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.


    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: -- 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.
     
    Barry Schwarz, Oct 5, 2006
    #5
  6. asaguiar

    Tim Prince Guest

    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: -- 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.
     
    Tim Prince, Oct 5, 2006
    #6
  7. In article <> "asaguiar" <> 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: -- 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.
     
    Dik T. Winter, Oct 5, 2006
    #7
  8. asaguiar

    Guest

    In comp.lang.c.moderated 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.
    > 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: -- 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, 2006
    #8
  9. asaguiar

    Eric Sosman Guest

    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: -- 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
    lid
    --
    comp.lang.c.moderated - moderation address: -- 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.
     
    Eric Sosman, Oct 8, 2006
    #9
  10. asaguiar

    Elan Riesman Guest

    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: -- 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.
     
    Elan Riesman, Oct 15, 2006
    #10
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Rag
    Replies:
    0
    Views:
    672
  2. Ahmed Moustafa
    Replies:
    0
    Views:
    814
    Ahmed Moustafa
    Nov 15, 2003
  3. Bapaiah Katepalli
    Replies:
    1
    Views:
    1,531
    Mike Treseler
    Jun 23, 2006
  4. Luc The Perverse

    [Algorithm] Sum of Primes < 1000000

    Luc The Perverse, Jan 5, 2007, in forum: Java
    Replies:
    50
    Views:
    1,989
  5. Szabolcs Horvát
    Replies:
    25
    Views:
    744
    Rhamphoryncus
    May 9, 2008
Loading...

Share This Page