Accurate summation functions.
Example:
import std/sums import std/math template `~=`(x, y: float): bool = abs(x - y) < 1e-4 let n = 1_000_000 first = 1e10 small = 0.1 var data = @[first] for _ in 1 .. n: data.add(small) let result = first + small * n.float doAssert abs(sum(data) - result) > 0.3 doAssert sumKbn(data) ~= result doAssert sumPairs(data) ~= result
See also
- math module for a standard sum proc
Procs
func sumKbn[T](x: openArray[T]): T
-
Kahan-Babuška-Neumaier summation: O(1) error growth, at the expense of a considerable increase in computational cost.
See:
Source Edit func sumPairs[T](x: openArray[T]): T
-
Pairwise (cascade) summation of x[i0:i0+n-1], with O(log n) error growth (vs O(n) for a simple loop) with negligible performance cost if the base case is large enough.
See, e.g.:
- https://en.wikipedia.org/wiki/Pairwise_summation
- Higham, Nicholas J. (1993), "The accuracy of floating point summation", SIAM Journal on Scientific Computing 14 (4): 783–799.
In fact, the root-mean-square error growth, assuming random roundoff errors, is only O(sqrt(log n)), which is nearly indistinguishable from O(1) in practice. See:
- Manfred Tasche and Hansmartin Zeuner, Handbook of Analytic-Computational Methods in Applied Mathematics (2000).