Add summation algorithms (#9284)

1 files changed, 78 insertions, 0 deletions

diff --git a/lib/std/sums.nim b/lib/std/sums.nim
new file mode 100644
index 000000000..dc9d8d507
--- /dev/null
+++ b/lib/std/sums.nim
@@ -0,0 +1,78 @@
+#
+#
+#            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