精确计算倍数之和

引言

许多算法都需要计算倍数之和。它非常重要，以至于已被添加到IEEE 浮点算术标准（IEEE Std 754-2008，第 9.4 节）中的推荐函数列表中。需要进行此类计算的操作示例包括

• 复数乘法和除法
• 向量长度
• 向量点积

不幸的是，执行这些操作的最常见方法可能会由于中间结果的舍入而导致不可接受的错误。

无误差转换

精确乘法

在fma()操作可用的情况下（如 IEEE Std 754-2008 所要求，并且存在于许多现代 CPU 中），并且假设没有溢出/下溢，则可以使用两个浮点运算来计算a*b的精确值

2MultFMA(a, b):
  s = a * b
  if s is non-finite
    return (s, 0)
  t = fma(a, b, -s)
  return (s,t)

精确结果由s+t给出。

• 如果没有检查，如果a*b溢出，则s将等于 +infinity，而t将等于 -infinity，并且总和将为 NaN
• 如果a*b下溢，则结果将不精确，但将给出最佳近似值

精确加法

2Sum(a, b):
  s = RN(a + b)
  if s is non-finite
    return (s, 0)
  b’ = RN(s − a)
  a’ = RN(s – b’)
  δb = RN(b – b’)
  δa = RN(a – a’)
  t = RN(δa + δb)
  return (s, t)

其中RN()表示“四舍五入到最近的偶数”。

精确结果由s+t给出。

• 如果没有检查，如果a+b溢出，则结果将为 NaN
• 如果a+b下溢，结果仍然精确

多个分量的求和

对两个以上分量进行无误差求和非常困难。 此处介绍的算法是补偿加法的一种变体，其中所有运算都是精确的。

给定一个包含 N 个浮点数的向量，我们可以按如下方式求和

for i = 1 to N-1
    (a[i], a[i-1]) = 2Sum(a[i], a[i-1])

这将向量a无误差地转换为一种形式，其中a[N-1]包含总和的近似值，而a[0] .. a[N-2]包含求和运算的残差。

重复此操作会“扫描”残差朝向总和，从而减小残差的大小并提高结果的准确性。 实际上，重复该算法一次(K = 2)通常就足够了。

Rump、Ogita 和 Oishi 分析了该算法，结果表明，在最坏的情况下，误差（所有残差的总和）收敛到略大于总和的 0.5 ulp 的值。 因此，该算法不能保证正确的舍入，但在大多数情况下，结果都是正确舍入的。

准确性

Muller 等人对各种求和算法进行了分析。 除了正式的分析之外，还给出了两个示例加法的结果

对于 N=5000 和 binary32 算术，Xi = RN(cos(i)) 的各种方法的误差

方法	ulps 中的误差
按升序添加值（需要排序！）	18.90625
按降序添加值（需要排序！）	16.90625
补偿加法 (Kahan)	6.90625
双重补偿 (Priest)（需要排序！）	0.09375
Rump、Ogita 和 Oishi	0.09375

对于 N = 10 和 binary32 算术，Xi = RN(1/i) 的各种方法的误差

方法	ulps 中的误差
按升序添加值（需要排序！）	6.860
按降序添加值（需要排序！）	738.900
补偿加法 (Kahan)	0.137
双重补偿 (Priest)（需要排序！）	0.137
Rump、Ogita 和 Oishi	0.137

正如我们所看到的，唯一始终给出准确结果（误差 < 0.5 ulp）的算法是 Priest 的和 Rump-Ogita-Oishi 的。 但是，Priest 的方法要求对操作数进行预排序。

结论

一般操作方法是

• 将每个乘法转换为一对浮点数，它们的总和是结果的无误差表示。
• 使用 Rump-Ogita-Oishi K 折加法算法对这 2*N 个数的向量求和，得到的精度略高于正确舍入结果的 0.5 ulp。

简单方法需要 N 次乘法和 N-1 次加法才能提供准确度可疑的结果。 此处描述的方法需要 N 次乘法、N 次fma()运算和6*K*(2*N-1)次加法/减法，但可以提供高度准确的结果。

参考文献

1. IEEE 浮点算术标准（Std 754-2008）
2. O. Møller. Quasi double-precision in floating-point addition. BIT, 5:37–50, 1965.
3. T. J. Dekker. A floating-point technique for extending the available precision. Numerische Mathematik, 18(3):224–242, 1971.
4. D. Knuth. The Art of Computer Programming, volume 2. Addison-Wesley, Reading, MA, 3rd edition, 1998.
5. S. M. Rump, T. Ogita, and S. Oishi. Accurate floating-point summation part I: Faithful rounding. SIAM Journal on Scientific Computing, 31(1):189–224, 2008.
6. Muller et al., Handbook of Floating-Point Arithmetic, Birkh¨auser, 2010