Kahan summation algorithm

In numerical analysis, theKahan summation algorithm (also known ascompensated summation) significantly reduces thenumerical error in the total obtained by adding asequence of finiteprecisionfloating point numbers, compared to the obvious approach.

该算法是数值分析中的算法，用来减少浮点数累加和的数值误差，可以在高性能计算对浮点精度要求高的过程中使用。

原理：This is done by keeping a separate running compensation (a variable to accumulate small errors). 例子中所使用的c

例子如下： 

假定我们使用六位有效数字的十进制浮点算法来累加, 第一个值：10000.0，第二个值为： 3.14159，第三个值为：2.71828。那么它们 精确和应该是 10005.85987, 四舍五入为：10005.9. 然而普通的算法，前两个数相加约分后为：10003.1. 然后再和第三个数相加约分后为 10005.8，结果错误。

例子中的计算过程如下（c 初始为0）：

  y = 3.14159 - 0                   y = input[i] - c
  t = 10000.0 + 3.14159
    = 10003.1                       Many digits have been lost!
  c = (10003.1 - 10000.0) - 3.14159 This must be evaluated as written! 
    = 3.10000 - 3.14159             The assimilated part of y recovered, vs. the original full y.
    = -.0415900                     Trailing zeros shown because this is six-digit arithmetic.
sum = 10003.1                       Thus, few digits from input(i) met those of sum.

 

  y = 2.71828 - -.0415900           The shortfall from the previous stage gets included.
    = 2.75987                       It is of a size similar to y: most digits meet.
  t = 10003.1 + 2.75987             But few meet the digits of sum.
    = 10005.85987, rounds to 10005.9
  c = (10005.9 - 10003.1) - 2.75987 This extracts whatever went in.
    = 2.80000 - 2.75987             In this case, too much.
    = .040130                       But no matter, the excess would be subtracted off next time.
sum = 10005.9                       Exact result is 10005.85987, this is correctly rounded to 6 digits.

 

算法伪代码如下：

function KahanSum(input)
    var sum = 0.0
    var c = 0.0          //A running compensation for lost low-order bits.
    for i = 1 to input.length do
        y = input[i] - c    //So far, so good: c is zero.
        t = sum + y         //Alas, sum is big, y small, so low-order digits of y are lost.
        c = (t - sum) - y   //(t - sum) recovers the high-order part of y; subtracting y recovers -(low part of y)
        sum = t             //Algebraically, c should always be zero. Beware eagerly optimising compilers!
        //Next time around, the lost low part will be added to y in a fresh attempt.
    return sum