-- --- -- title: Kahan summation -- --- -- -- Summation of floating-point numbers is vulnerable to roundoff -- error. The [Kahan summation -- algorithm](https://en.wikipedia.org/wiki/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) let 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](https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements) -- 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. let 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] let normal_sum = f64.sum -- > normal_sum [1e100, 1.0, -1e100, 1.0] -- 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 it 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.