Irmen's blog

Compensated sum (Kahan summation): a+b is niet c

Een rij floating point getallen optellen op de normale manier levert een steeds groter wordende afrondingsfout op. Er bestaat een algorithme dat deze fout behoorlijk verkleint: Kahan Summation.

Ik heb het even geprobeerd in Java en ten opzichte van de absolute precisie die je krijgt door BigDecimal te gebruiken voor de optelling, leverde dat het volgende resultaat:

• 100 random doubles: normal sum error: 7.105427e-15 ; compensated sum error: 0.000000e+00
• 10000 random doubles: normal sum error: -9.094947e-13 ; compensated sum error: 0.000000e+00
• 1 miljoen random doubles: normal sum error: -1.338776e-09 ; compensated sum error: 0.000000e+00

Pas zo'n beetje vanaf 1 miljoen waardes wordt de normale error groter dan 1e-10, dus tenzij je echt absolute precisie wilt zou je gewoon normale optelling kunnen gebruiken... Of je kiest ervoor om sowieso BigDecimal (arbitrary precision decimals) te gebruiken voor dit soort 'kritieke' optellingen. Het is veel langzamer dan doubles optellen, maar als het niet om heel veel waardes gaat is dat wellicht verwaarloosbaar.

De Kahan summation was bijna even snel als de gewone optelling, en leverde in al mijn tests een fout op van nul, dus in speciale gevallen is dit een interessante optie. Let wel op de problemen die je kunt krijgen met optimizing compilers: als die slim genoeg zijn, dan optimaliseren ze het hele algoritme eruit zodat er een normale optelling overblijft. Zie het Wikipedia artikel voor meer info.

Dit is de code voor de compensated sum:

double compensatedSum(double[] numbers)
{
    double s = 0.0;
    double c = 0.0;
    double y;
    double t;
    for (double number : numbers) {
         y = number-c;
         t = s+y;
         c = (t-s) - y;
         s = t;
    }
    return s;
}