Accurate summation functions
Example:
import 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 assert abs(sum(data) - result) > 0.3 assert sumKbn(data) =~ result assert sumKbk(data) =~ result assert sumPairs(data) =~ result assert sumShewchuk(data) =~ result
See also
- math module for a standard sum proc
Procs
func sumKbn[T: SomeFloat](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 sumKbk[T: SomeFloat](x: openArray[T]): T
-
Kahan-Babuška-Klein variant, a second-order "iterative Kahan–Babuška algorithm".
See:
Source Edit func sumPairs[T: SomeFloat](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).
func sumShewchuk[T: SomeFloat](x: openArray[T]): T
-
Shewchuk's summation Full precision sum of values in iterable. Returns the value of the sum, rounded to the nearest representable floating-point number using the round-half-to-even rule
Reference:
- Shewchuk, JR. (1996) Adaptive Precision Floating-Point Arithmetic and Fast Robust GeometricPredicates.