Add more summation functions

lib/std/sums.nim

@@ -53,13 +53,111 @@ func sumPairs*[T](x: openArray[T]): T =
   let n = len(x)
   if n == 0: T(0) else: sumPairwise(x, 0, n)
 
+func fastTwoSum*[T: SomeFloat](a, b: T): (T, T) =
+  ## Deker's algorithm
+  ## pre-condition: |a| >= |b|
+  ## you must swap a and b if pre-condition is not satisfied
+  ## s + r = a + b exactly (s is result[0] and r is result[1])
+  ## s is the nearest FP number of a + B
+  assert(abs(a) >= abs(b))
+  result[0] = a + b
+  result[1] = b - (result[0] - a)
+
+func twoSum*[T](a, b: T): (T, T) =
+  ## Møller-Knuth's algorithm.
+  ## Improve Deker's algorithm, no branch needed
+  ## More operations, but still cheap.
+  result[0] = a + b
+  let z = result[0] - a
+  result[1] = (a - (result[0] - z)) + (b - z)
+
+func sum2*[T: SomeFloat](v: openArray[T]): T =
+  ## sum an array v using twoSum function
+  if len(v) == 0:
+    return 0.0
+
+  var s = v[0]
+  var e: float
+  for i in 1..<v.len-1:
+    let sum = twoSum(s, v[i])
+    s = sum[0]
+    e += sum[1]
+  return s + e
+
+func sumShewchuck_add[T: SomeFloat](v: openArray[T]): seq[T] =
+  for x in v:
+    var x = x
+    var i = 0
+    for y in result:
+      let sum = twoSum(x, y)
+      let hi = sum[0]
+      let lo = sum[1]
+      if lo != 0.0:
+        result[i] = lo
+        i.inc
+      x = hi
+    setLen(result, i + 1)
+    result[i] = x
+
+func sumShewchuck_total[T: SomeFloat](partials: openArray[T]): T =
+  var hi = 0.0
+  if len(partials) > 0:
+    var n = len(partials)
+    dec(n)
+    hi = partials[n]
+    var lo = 0.0
+    while n > 0:
+      var x = hi
+      dec(n)
+      var y = partials[n]
+      let sum = twoSum(x, y)
+      hi = sum[0]
+      lo = sum[1]
+      if lo != 0.0:
+        break
+      if (n > 0 and
+          (
+            (lo < 0.0 and partials[n - 1] < 0.0) or
+            (lo > 0.0 and partials[n - 1] > 0.0)
+          )
+        ):
+        y = lo * 2.0
+        x = hi + y
+        var yr = x - hi
+        if y == yr:
+          hi = x
+  result = hi
+
+func sumShewchuck*[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
+  ##
+  ## See also:
+  ## - https://docs.python.org/3/library/math.html#math.fsum
+  ## - https://code.activestate.com/recipes/393090/
+  ##
+  ## Reference:
+  ## Shewchuk, JR. (1996) Adaptive Precision Floating-Point Arithmetic and \
+  ## Fast Robust GeometricPredicates.
+  ## http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps
+  sumShewchuck_total(sumShewchuck_add(x))
+
+func fsum*[T: SomeFloat](x: openArray[T]): T =
+  ## Alias of shewchuckSum function as python fsum
+  sumShewchuck(x)
 
 runnableExamples:
   static:
     block:
       const data = [1, 2, 3, 4, 5, 6, 7, 8, 9]
       doAssert sumKbn(data) == 45
       doAssert sumPairs(data) == 45
+      const dataFloat = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
+      doAssert sumPairs(dataFloat) != 1.0 # 0.9999...
+      doAssert sumKbn(dataFloat) == 1.0
+      doAssert sumShewchuck(dataFloat) == 1.0
 
 
 when isMainModule:
@@ -70,17 +168,20 @@ when isMainModule:
     epsilon /= 2.0
   let data = @[1.0, epsilon, -epsilon]
   assert sumKbn(data) == 1.0
+  assert sumShewchuck(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 sumShewchuck(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 sumShewchuck(tc2) == -0.6926474305598203
   assert sumPairs(tc2) == -0.6926474305598204