1 files changed, 80 insertions, 0 deletions
diff --git a/lib/deprecated/pure/sums.nim b/lib/deprecated/pure/sums.nim
new file mode 100644
index 000000000..a6ce1b85d
--- /dev/null
+++ b/lib/deprecated/pure/sums.nim
@@ -0,0 +1,80 @@
+#
+#
+#            Nim's Runtime Library
+#        (c) Copyright 2019 b3liever
+#
+#    See the file "copying.txt", included in this
+#    distribution, for details about the copyright.
+
+## Accurate summation functions.
+
+{.deprecated: "use the nimble package `sums` instead.".}
+
+runnableExamples:
+  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 <math.html>`_ for a standard `sum proc <math.html#sum,openArray[T]>`_
+
+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:
+  ## * https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements
+  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.:
+  ## * 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).
+  let n = len(x)
+  if n == 0: T(0) else: sumPairwise(x, 0, n)