跳转至

Kahan 求和

引入

Kahan 求和 算法，又名補償求和或進位求和算法，是一個用來 降低有限精度浮點數序列累加值誤差 的算法。它主要通過保持一個單獨變量用來累積誤差（常用變量名為 \(c\)）來完成的。

該算法主要由 William Kahan 於 1960s 發現。因為 Ivo Babuška 也曾獨立提出了一個類似的算法，Kahan 求和算法又名為 Kahan–Babuška 求和算法。

舍入誤差

在計算機程序中，我們需要用有限位數對實數做近似表示，如今的大多數計算機都使用 IEEE-754 規定的浮點數來作為這個近似表示。對於 \(\frac{1}{3}\)，由於我們不能在有限位數內對它進行精準表示，因此在使用 IEEE-754 表示法時，必須四捨五入一部分數值（truncate）。這種 舍入誤差（Rounding off error）是浮點計算的一個特徵。

在浮點加法計算中，交換律（commutativity）成立，但結合律（associativity）不成立。也就是説，\(a+b = b+a\) 但 \((a+b)+c \neq a+(b+c)\)。因此在浮點序列加法計算中，我們可以從左到右一個個累加，也可以在原有順序上，將他們兩兩分成一對。第二種算法會相對較慢並需要更多內存，也常被一些語言的特定求和函數使用，但相對結果更準確。

為了得到更準確的浮點累加結果，我們需要使用 Kahan 求和算法。

在計算 \(S_{new}=S_{old}+a\)（\(a\) 為浮點序列的一個數值）時，定義實際計算加入 \(S\) 的值為 \(a_{eff}=S_{new}-S_{old}\), 如果 \(a_{eff}\) 比 \(a\) 大，則證明有向上舍入誤差；如果 \(a_{eff}\) 比 \(a\) 小，則證明有向下舍入誤差。則舍入誤差定義為 \(E_{roundoff} = a_{eff} - a\)。那麼用來糾正這部分舍入誤差的值就為 \(a-a_{eff}\), 即 \(E_{roundoff}\) 的負值。定義 \(c\) 是對丟失的低位進行運算補償的變量，就可以得到 \(c_{new} = c_{old} + (a - a_{eff})\)。

過程

Kahan 求和算法主要通過一個單獨變量用來累積誤差。如下方參考代碼所示，\(sum\) 為最終返回的累加結果。\(c\) 是對丟失的低位進行運算補償的變量（其被捨去的部分），也是 Kahan 求和算法中的必要變量。

因為 \(sum\) 大，\(y\) 小，所以 \(y\) 的低位數丟失。\((t - sum)\) 抵消了 \(y\) 的高階部分，減去 \(y\) 則會恢復負值（\(y\) 的低價部分）。因此代數值中 \(c\) 始終為零。在下一輪迭代中，丟失的低位部分會被更新添加到 \(y\)。

實現

參考代碼 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
float kahanSum(vector<float> nums) {
  float sum = 0.0f;
  float c = 0.0f;
  for (auto num : nums) {
    float y = num - c;
    float t = sum + y;
    c = (t - sum) - y;
    sum = t;
  }
  return sum;
}

習題

在 OI 中，Kahan 求和主要作為輔助工具存在，為計算結果提供誤差更小的值。

例題 CodeForces Contest 800 Problem A. Voltage Keepsake

有 \(n\) 個同時使用的設備。第 \(i\) 個設備每秒使用 \(a_{i}\) 單位的功率。這種用法是連續的。也就是説，在 \(\lambda\) 秒內，設備將使用 \(\lambda \times a_{i}\) 單位的功率。第 \(i\) 個設備當前存儲了 \(b_{i}\) 單位的電力。所有設備都可以存儲任意數量的電量。有一個可以插入任何單個設備的充電器。充電器每秒會為設備增加 \(p\) 個單位的電量。這種充電是連續的。也就是説，如果將設備插入 \(\lambda\) 秒，它將獲得 \(\lambda \times p\) 單位的功率。我們可以在任意時間單位內（包括實數）切換哪個設備正在充電（切換所需時間忽略不計）。求其中一個設備達到 \(0\) 單位功率前，可以使用這些設備的最長時間。

例題 CodeForces Contest 504 Problem B. Misha and Permutations Summation

定義數字 \(0, 1, \cdots, (n - 1)\) 的兩個排列 \(p\) 和 \(q\) 的和為 \(Perm((Ord(p)+Ord(q))\bmod n!)\)，其中 \(Perm(x)\) 是數字 \(0, 1, \cdots, (n-1)\) 的第 \(x\) 個字典排列（從零開始計數），\(Ord(p)\) 是字典序排列 \(p\) 的個數。例如，\(Perm(0) = (0, 1, \cdots , n - 2, n - 1)\)，\(Perm(n! - 1) = (n - 1, n-2,\cdots, 1,0))\)。Misha 有兩個排列 \(p\) 和 \(q\)，找到它們的總和。

編程語言的求和

Python 的標準庫指定了精確舍入求和的 fsum 函數可用於返回可迭代對象中值的準確浮點總和，它通過使用 Shewchuk 算法跟蹤多箇中間部分和來避免精度損失。

Julia 語言中，sum 函數的默認實現是成對求和，以獲得高精度和良好的性能。同時外部庫函數 sum_kbn 為需要更高精度的情況提供了 Neumaier 變體的實現，具體可見 KahanSummation.jl。

參考資料與註釋

  1. Kahan_summation_algorithm - Wikipedia
  2. Kahan summation - Rosetta Code
  3. VK Cup Round 2 + Codeforces Round 409 Announcement
  4. Rounding off errors in Java - GeeksforGeeks
本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用