介绍
该答案对以下错误进行了初步检查:
float result = valA * 0.587f + valB * 0.331f;
在这个答案中,浮点格式的值和使用浮点格式计算的表达式将用代码样式表示,如z
or x * y
。数学变量将使用斜体并且不会采用代码样式,如z或x • y。
我假设所有算术都是使用 IEEE-754 基本 32 位二进制浮点数完成的。这种格式通常用于float
类型,尽管一些编程语言实现混合精度,可能double
在评估float
类型表达式时使用或其他精度。我还假设所有算术都是使用四舍五入模式完成的,并与偶数低位的数字相关联。
这种格式为 24 位的有效位,因此最低精度单位 (ULP) 通常是最高有效位值的2 -23倍。这是可表示值之间的步长。例如,对于 [1, 2) 中的值,ULP 为 2 -23。对于 [128, 256) 中的值,ULP 为 2 7 •2 -23 = 2 -16。(对于次正规值,有效数字的位数较少。ULP 的最低值可以是 2 -149。超出最大的有限可表示值,到下一个可表示值的步长是无限的。但是,在这个问题中,只有适度的值涉及值,因此我们可以忽略无穷大。)
使用正确舍入计算任何运算的结果最多与正确答案相差 ½ ULP。也就是说,z = x + y
例如,如果我们计算 ,计算结果z
与精确的数学结果z = x
+最多相差zy
的 ½ ULP 。(虽然z是一个具有无限精度的精确数学结果,但我们使用它的大小来确定它在浮点格式中落入哪个范围,这决定了我们所说的z的 ULP 的含义。)误差最多为 ½ 的原因ULP 是,如果最接近z的两个可表示值是和,我们必须有≤ z ≤ ,并且如果 ½ ULP <z0
z1
z0
z1
z1
− z,然后z − z0
< ½ ULP(因为z1
− z0
= 1 ULP,根据 ULP 的定义。)因此,在选择最接近的可表示值时,我们会选择 和 的更接近者z0
,z1
因此误差永远不会超过 ½ ULP。
如评论中所述valA
,valB
、 和result
位于 [0, 256) 中。
符号分析
当我们开始计算valA * 0.587f + valB * 0.331f
时,之前valA
的valB
操作会出现一些错误。也就是说,理想情况下,使用精确的数学,我们会计算出一些数字A和B,但计算机计算的是valA
和valB
,差异是eA = valA
- A和eB = valB
- B。
理想情况下,我们希望使用精确的数学计算R使得R为A • .587 + B • .331。当我们使用计算机算术时:
0.587f
将从 .587 转换为浮点格式,结果会有一些舍入误差e0,所以结果是0.587f
= 0.587 + e0。
0.331f
将从 .331 转换为浮点格式,结果会有一些舍入误差e1,所以结果是0.331f
= .331 + e1。
valA * 0.587f
将计算出一些错误e2,因此结果将是valA * 0.587f
= valA
• 0.587f
+ e2。
valB * 0.331f
将计算出一些错误e3,因此结果将是valB * 0.331f
= valB
• 0.331f
+ e3。
- 这两个产品将被添加,有一些错误e4,所以结果将是
valA * 0.587f + valB * 0.331f
= valA * 0.587f
++ valB * 0.331f
e4 。
现在我们可以替换表达式,所以:
valA * 0.587f + valB * 0.331f
= ( valA
• 0.587f
+ e2 ) + ( valB
• 0.331f
+ e3 + e4。
valA * 0.587f + valB * 0.331f
= ( valA
• (0.587 + e0 ) + e2 ) + ( valB
• (.331 + e1 ) + e3 ) + e4。
valA * 0.587f + valB * 0.331f
= (( A + eA ) • (0.587 + e0 ) + e2 ) + (( B + eB ) • (.331 + e1 ) + e3 ) + e4。
这样,我们将计算结果 表示valA * 0.587f + valB * 0.331f
为精确的数学表达式(具有不完全已知值的变量), (( A + eA ) • (0.587 + e0 ) + e2 ) + (( B + eB ) • (. 331 + e1 ) + e3 ) + e4。
数值分析
接下来,我们可以对错误设置一些界限。e0和e1很容易,它们的幅度分别为 0.587 和 0.331 的 ½ ULP。.587 在 [½, 1) 中,所以它的 ULP 是 2 -24,而 .331 在 [¼, ½) 中,所以它的 ULP 是 2 -25。所以| e0 | ≤= 2 -25和 | e1 | ≤= 2 -26。
e2和e3的边界取决于 和 的valA * 0.587f
大小valB * 0.331f
。因为val
< 256, valA * 0.587f
< 256,所以它的 ULP 最多为 2 -16,并且 | e2 | ≤ 2 -17。有了valB
,我们可以看到valB * 0.331f
< 128,因此 ULPvalB * 0.331f
最多为 2 -17,并且 | e3 | ≤ 2 -18。
最后,我们在最后添加valA * 0.587f + valB * 0.331f
. 我们假设它小于 256,所以它的 ULP 最多为 2 -16,并且 | e4 | ≤ 2 -17。
查看计算结果的数学表达式 (( A + eA ) • (0.587 + e0 ) + e2 ) + (( B + eB ) • (.331 + e1 ) + e3 ) + e4,我们可以看到当e0、e1、e2、e3和e4具有最大值时,会发生最大可能的错误(除非eA或eB是巨大的负数,我们假设这不是真的)。因此,我们可以替换为这些错误准备的上限:
(( A + eA ) • (0.587 + 2 -25 ) + 2 -17 ) + (( B + eB ) • (.331 + 2 -26 ) + 2 -18 ) + 2 -17。
为了时间的利益,我用 Maple 对此进行了评估。(手动扩展表达式并保留一些因子而不是将系数合并为单个数字可能更有启发性,但我将其留给读者。)结果是:
2462056573/4194304000 • A + 2462056573/4194304000 • eA + 5/262144 + 2776629373/8388608000 • B + 2776629373/8388608000 • eB。
理想的结果是A • .587 + B • .331。当我们从上面减去它时,结果是计算误差的界限:
1/33554432 • A + 2462056573/4194304000 • eA + 5/262144 + 1/67108864 * B + 2776629373/8388608000 • eB。
由于A < 256 和B < 256,我们可以用 256 代替A和B,得到:
1/32768 + 2462056573/4194304000 • eA + 2776629373/8388608000 • eB。
反转一点 Maple 的算术,即:
2 -15 + (.587 + 2 -25 ) • eA + (.331 + 2 -26 ) • eB。
因此,这是 中错误的上限valA * 0.587f + valB * 0.331f
。valA
通过有关和之间关系的附加信息,它可能会进一步减少valB
。此外,将 .587 和 .331 转换为的错误float
是完全已知的,因此应该使用这些错误,而不是我在此答案中用作说明的边界。
还需要建立误差的下限。舍入误差可能是负数,我们必须问 (( A + eA ) • (0.587 + e0 ) + e2 ) + (( B + eB ) • (.331 + e1 ) + e3 )的最低可能值是多少+ e4是。由于我现在没有时间,这留给读者。
附录
e0是 13/1048576000。e1是 1/4194304000。那么误差的上限可以减少到 731/32768000 + 4924113/8388608 • eA + 11106517/33554432 • eB,即:
.731•2 -15 + (.587 + .013•2 -20 ) • eA + (.331 + .001•2 -22 ) • eB。