c_k 的准确度 = a + ( N + k ) * b答案

【问题标题】：Accuracy of c_k = a + ( N + k ) * bc_k 的准确度 = a + ( N + k ) * b
【发布时间】：2013-05-16 12:57:56
【问题描述】：

a, b 是 32 位浮点值，N 是 32 位整数，k 可以取值 0, 1, 2, ... M。需要计算 c_k = a + ( N + k ) * b ;这些操作需要是 32 位操作（不是双精度）。关注的是准确性——以下哪个更准确？：

I) c_k = a + ( N + k ) * b

II) 首先计算：c_0 = a + N * b
然后通过加法迭代计算c_1、c_2等：
c_1 = c_0 + b;
c_2 = c_1 + b;

【问题讨论】：

我的直觉是选项 II 更好（性能更高），但我不太确定它是否更准确。我想这取决于数据。
链加法是你能做的最糟糕的操作之一，因为最后一个结果中的舍入误差将是链中每个加法的舍入误差之和。使用第一种方式或使用c_i = c_0 + b*i 会更准确。
@PatriciaShanahan 您应该在评论中提交答案。这是这个问题的关键信息。

标签： floating-point floating-accuracy multiplication addition single-precision

【解决方案1】：

链式加法是你能做的最糟糕的操作之一，因为最后一个结果中的舍入误差将是链中每个加法上单个操作舍入误差的净和。使用第一种方式或使用c_i = c_0 + b*i 会更准确。

【讨论】：

【解决方案2】：

由于您似乎并不关心操作的数量，假设 IEEE 754 模型您可以使用 32 位操作完全执行它。
请参阅 Shewchuck 自适应精度浮点算术和快速稳健的几何谓词 - http://www.cs.berkeley.edu/~jrs/papers/robustr.pdf 或 http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps

你定义了两个精确的操作（见论文）

(product,residue) = twoproduct(a,b)
(sum,residue) = twosum(a,b)

然后你必须将 N+k 分解为两个 24 位有效数字，例如

NkH = (N+k) / 256;
NkL = (N+K) % 256;

那么你有两个可能不精确的乘法

( HH , HL ) = twoproduct( NkH , b)
( LH , LL ) = twoproduct( NkL , b)

然后你可以将这些 ( HH , HL ) + ( LH , LL ) + a 相加

这可以通过快速扩展和完全执行（再次参见论文）

(c1,c2,c3,c4,c5) = sort_increasing_magnitude(HH,HL,LH,LL,a)
(s2,s1) = twosum( c2,c1 )
(s3,s2) = twosum( c3,s2 )
(s4,s3) = twosum( c4,s3 )
(s5,s4) = twosum( c5,s4 )

然后您会在 s5 中得到精确舍入的结果，就好像这些操作是使用无限精度算术执行的一样。

【讨论】：