JuliaLang/julia — pull request #4039 by stevengj (2013-08-13 03:28:53+00:00)
Title: RFC: use pairwise summation for sum, cumsum, and cumprod

This patch modifies the `sum` and `cumsum` functions (and `cumprod`) to use [pairwise summation](http://en.wikipedia.org/wiki/Pairwise_summation) for summing `AbstractArray`s, in order to achieve much better accuracy at a negligible computational cost.

Pairwise summation recursively divides the array in half, sums the halves recursively, and then adds the two sums. As long as the base case is large enough (here, n=128 seemed to suffice), the overhead of the recursion is negligible compared to naive summation (a simple loop). The advantage of this algorithm is that achieves O(sqrt(log n)) mean error growth, versus O(sqrt(n)) for naive summation, which is almost indistinguishable from the O(1) error growth of Kahan compensated summation.

For example:
```
A = rand(10^8)
es = sum_kbn(A)
(oldsum(A) - es)/es, (newsum(A) - es)/es
```
where `oldsum` and `newsum` are the old and new implementations, respectively, gives `(-1.2233732622777777e-13,0.0)` on my machine in a typical run: the old `sum` loses three significant digits, and the new sum loses none. On the other hand, their execution time is nearly identical:
```
@time oldsum(A);
@time newsum(A);
```
gives
```
elapsed time: 0.116988841 seconds (105816 bytes allocated)
elapsed time: 0.110502884 seconds (64 bytes allocated)
```