#
#
# Nim's Runtime Library
# (c) Copyright 2019 b3liever
#
# See the file "copying.txt", included in this
# distribution, for details about the copyright.
## Fast sumation functions.
func sumKbn*[T](x: openArray[T]): T =
## Kahan (compensated) summation: O(1) error growth, at the expense
## of a considerable increase in computational expense.
if len(x) == 0: return
var sum = x[0]
var c = T(0)
for i in 1 ..< len(x):
let xi = x[i]
let t = sum + xi
if abs(sum) >= abs(xi):
c += (sum - t) + xi
else:
c += (xi - t) + sum
sum = t
result = sum + c
func sumPairwise[T](x: openArray[T], i0, n: int): T =
if n < 128:
result = x[i0]
for i in i0 + 1 ..< i0 + n:
result += x[i]
else:
let n2 = n div 2
result = sumPairwise(x, i0, n2) + sumPairwise(x, i0 + n2, n - n2)
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.:
## * http://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).
##
let n = len(x)
if n == 0: T(0) else: sumPairwise(x, 0, n)
when isMainModule:
from math import pow
var epsilon = 1.0
while 1.0 + epsilon != 1.0:
epsilon /= 2.0
let data = @[1.0, epsilon, -epsilon]
assert sumKbn(data) == 1.0
assert sumPairs(data) != 1.0 # known to fail
assert (1.0 + epsilon) - epsilon != 1.0
var tc1: seq[float]
for n in 1 .. 1000:
tc1.add 1.0 / n.float
assert sumKbn(tc1) == 7.485470860550345
assert sumPairs(tc1) == 7.485470860550345
var tc2: seq[float]
for n in 1 .. 1000:
tc2.add pow(-1.0, n.float) / n.float
assert sumKbn(tc2) == -0.6926474305598203
assert sumPairs(tc2) == -0.6926474305598204