Nim/lib/std/sums.nim

#
#
#            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.

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)