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;
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;
链式加法是您可以做的最糟糕的操作之一,因为最后一个结果中的舍入误差将是链中每个加法上单个操作舍入误差的净和。使用第一种方式或使用c_i = c_0 + b*i
.
由于您似乎并不关心操作的数量,假设 IEEE 754 模型您可以使用 32 位操作完全执行它。
请参阅 Shewchuck 自适应精度浮点算术和快速稳健的几何谓词 - http://www.cs.berkeley.edu/~jrs/papers/robusr.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 中获得精确舍入的结果,就好像这些操作是使用无限精度算术执行的一样。