Numerically stable summation

It is known that adding up a large number of small floating-point values can lead to significant rounding errors unless a smart summation algorithm is applied. In fact, such a problem can easily be triggered in DAPHNE at the moment.

The following script generates a matrix of 100 M single-precision floating-point values uniformly distributed in [0, 1]. The sum should be roughly 0.5 * 100 M = 50 M = 5.0e+07. However, the output (on my system) is 1.67772e+07.

X = rand(10^4, 10^4, as.f32(0.0), as.f32(1.0), 1, -1);
print(sum(X));

This task is to implement a numerically stable summation in the aggAll-kernel (src/runtime/local/kernels/AggAll.h). The Kahan summation algorithm could be a candidate.

One could extend the task to check in how far other aggregates (e.g., column-wise (AggCol.h), row-wise (AggRow.h), and cumulative (AggCum.h)) or other aggregation functions (e.g., mean, var, stddev) or other values types (f64) are affected by similar problems.