Source file: kahan.fut

Kahan summation

Summation of floating-point numbers is vulnerable to roundoff error. The Kahan summation algorithm improves the accuracy over naive summation by tracking an additional error term, which collects small values that would otherwise be lost.

type kahan = (f64, f64)

def kahan_add ((s1, c1) : kahan) ((s2, c2) : kahan) : kahan =
  let s = s1 + s2
  let d  = if f64.abs s1 >= f64.abs s2 then (s1 - s) + s2
           else (s2 - s) + s1
  let c = (c1 + c2) + d
  in (s, c)

The above actually incorporates a tweak by Neumaier that better handles the case where the term to be added is larger than the running sum. We can pack it up as a map-reduce composition for summing an entire array.

def kahan_sum (xs: []f64) : f64 =
  let (s,c) = reduce kahan_add (0,0) (map (\x -> (x,0)) xs)
  in s + c

And it really is more accurate than normal summation, as easily shown with some pathological examples:

> kahan_sum [1e100, 1.0, -1e100, 1.0]

2.0f64

def normal_sum = f64.sum

> normal_sum [1e100, 1.0, -1e100, 1.0]

1.0f64

Is kahan_add actually associative, as required by reduce? No, it’s not, but neither is floating-point addition in the first place. Passing a non-associative operator to reduce doesn’t mean the result can be arbitrary, it just means the result is nondeterministic, based on how the compiler or runtime decides to “parenthesise” the application of the operator. We’re already implicitly accepting this when we parallelise a floating-point summation, so we should be fine with it for Kahan summation as well.