Kahan summation

Kahan summation is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

The Kahan summation algorithm is a method of summing a series of numbers represented in a limited precision in a way that minimises the loss of precision in the result.

The task is to follow the previously linked Wikipedia articles algorithm and its worked example by:

1. Do all arithmetic in decimal to a precision of six digits.
2. Write a function/method/procedure/subroutine/... to perform Kahan summation on an ordered collection of numbers, (such as a list of numbers).
3. Create the three numbers a, b, c equal to 10000.0, 3.14159, 2.71828 respectively.
4. Show that the simple left-to-right summation, equivalent to (a + b) + c gives an answer of 10005.8
5. Show that the Kahan function applied to the sequence of values a, b, c results in the more precise answer of 10005.9

Show all output on this page.

Go

<lang go>package main

import "fmt"

func kahan(s ...float64) float64 {

   var tot, c float64
   for _, x := range s {
       y := x - c
       t := tot + y
       c = (t - tot) - y
       tot = t
   }
   return tot

}

func main() {

   a := 1e15
   b := 3.14159
   fmt.Println("Left associative:", a+b+b+b+b+b+b+b+b+b+b)
   fmt.Println("Kahan summation: ", kahan(a, b, b, b, b, b, b, b, b, b, b))

}</lang>

Output:

Left associative: 1.0000000000000312e+15
Kahan summation:  1.0000000000000314e+15

Python

<lang python>>>> from decimal import * >>> >>> getcontext().prec = 6 >>> >>> def kahansum(input):

   summ = c = 0
   for num in input:
       y = num - c
       t = summ + y
       c = (t - summ) - y
       summ = t
   return summ

>>> a, b, c = [Decimal(n) for n in '10000.0 3.14159 2.71828'.split()] >>> a, b, c (Decimal('10000.0'), Decimal('3.14159'), Decimal('2.71828')) >>> >>> (a + b) + c Decimal('10005.8') >>> kahansum([a, b, c]) Decimal('10005.9') >>> >>> >>> sum([a, b, c]) Decimal('10005.8') >>> # it seems Python's sum() doesn't make use of this technique. >>> >>> # More info on the current Decimal context: >>> getcontext() Context(prec=6, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[Inexact, Rounded], traps=[InvalidOperation, DivisionByZero, Overflow]) >>> >>> >>> ## Lets try the simple summation with more precision for comparison >>> getcontext().prec = 20 >>> (a + b) + c Decimal('10005.85987') >>> </lang>

PHP

Interactive mode

<lang php>$ php -a Interactive mode enabled

php > function kahansum($input) { php { $sum = $c = 0; php { foreach($input as $item) { php { $y = $item + $c; php { $t = $sum + $y; php { $c = ($t - $sum) - $y; php { $sum = $t; php { } php { return $sum; php { } php > $input = array(10000.0, 3.14159, 2.71828); php > list($a, $b, $c) = $input; php > echo $sumabc."\n"; 10005.85987 php > echo kahansum($input)."\n"; 10005.85987 php > echo array_sum($input)."\n"; 10005.85987 php > // does not seem to be any difference in PHP 5.5 </lang>

Script

<lang php><?php function kahansum($input) {

   $sum = $c = 0;                                                              
   foreach($input as $item) {                                                  
       $y = $item + $c;                                                        
       $t = $sum + $y;                                                         
       $c = ($t - $sum) - $y;                                                  
       $sum = $t;                                                              
   }                                                                           
   return $sum;                                                                

}

$input = array(10000.0, 3.14159, 2.71828); list($a, $b, $c) = $input; $sumabc = $a + $b + $c; echo $sumabc."\n"; echo kahansum($input)."\n"; echo array_sum($input)."\n"; </lang>

Script output: <lang bash> $ php kahansum.php 10005.85987 10005.85987 10005.85987 </lang>

R

<lang r> > # A straightforward implementation of the algorithm > kahansum <- function(x) { + ks <- 0 + c <- 0 + for(i in 1:length(x)) { + y <- x[i] - c + kt <- ks + y + c = (kt - ks) - y + ks = kt + } + ks + } > # make R display more digits in the result > options(digits=10) > a <- 10000.0 > b <- 3.14159 > c <- 2.71828 > a + b + c [1] 10005.85987 > # no problem there > input <- c(10000.0, 3.14159, 2.71828) > kahansum(input) [1] 10005.85987 > # for completeness, let's compare with the standard sum() function > sum(input) [1] 10005.85987 > # seems like R does the "right thing" > # bonus: > # let's use the Dr. Dobbs example (http://www.drdobbs.com/floating-point-summation/184403224) > input <- c(1.0, 1.0e10, -1.0e10) > kahansum(input) [1] 1 > sum(input) [1] 1 > # again, R does it right </lang>