1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
|
#
#
# 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-Babuška-Neumaier 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)
runnableExamples:
static:
block:
const data = [1, 2, 3, 4, 5, 6, 7, 8, 9]
doAssert sumKbn(data) == 45
doAssert sumPairs(data) == 45
|