c - 再次浮点精度

Question

昨天我问了一个question关于为什么我在浮点运算中失去准确性。我收到了关于中间结果如何保存在 x87 寄存器中的答案。这很有帮助，但一些细节仍然在逃避我。这是我在上一个问题中提供的程序的变体，我在 Debug模式下使用 VC++ 2010 Express。

int main()
{
    double x = 1.8939201459282359e-308; /* subnormal number */
    double tiny = 4.9406564584124654e-324; /* smallest IEEE double */
    double scale = 1.6;
    double temp = scale*tiny;
    printf("%23.16e\n", x + temp);
    printf("%23.16e\n", x + scale*tiny);
}

这输出

1.8939201459282369e-308
1.8939201459282364e-308

根据 IEEE 标准，第一个值是正确的。为 scale 变量赋予 2.0 的值可为两种计算提供正确的值。我知道第一次计算中的 temp 是一个低于正常值的值，因此会失去精度。我还了解到 scale*tiny 的值保存在 x87 寄存器中，该寄存器具有更大的指数范围，因此该值比 temp 更精确。我不明白的是，当将值添加到 x 时，我们从较低精度的值中得到了正确的答案。当然，如果精度较低的值可以给出正确答案，那么精度较高的值也应该给出正确答案吗？这与“双舍入”有关吗？

提前致谢，这对我来说是一个全新的主题，所以我有点挣扎。

Answer 1

关键是由于指数范围较大，这两个数字在 x87 表示中不是次正规的。

在 IEEE754 表示中，

x    = 0.d9e66553db96f × 2^(-1022)
tiny = 0.0000000000001 × 2^(-1022)

但在 x87 表示中，

x    = 1.b3cccaa7b72de × 2^(-1023)
tiny = 1.0000000000000 × 2^(-1074)

现在，当在 IEEE754 表示中计算 1.6*tiny 时，它会四舍五入为 0.0000000000002 × 2^(-1022)，因为这是最接近的可表示数字数学结果。将其添加到 x 结果是

  0.d9e66553db96f × 2^(-1022)
+ 0.0000000000002 × 2^(-1022)
-----------------------------
  0.d9e66553db971 × 2^(-1022)

但是在 x87 表示中，1.6*tiny 变成了

1.999999999999a × 2^(-1074)

什么时候添加

  1.b3cccaa7b72de × 2^(-1023)
+ 0.0000000000003333333333334 × 2^(-1023)
-----------------------------------------
  1.b3cccaa7b72e1333333333334 × 2^(-1023)

四舍五入到 53 位有效位的结果是

  1.b3cccaa7b72e1 × 2^(-1023)

尾数为 1 的最后一位。如果将其转换为 IEEE754 表示形式(尾数最多可以有 52 位，因为它是次正规数)，因为它恰好位于两个相邻可表示数之间的一半0.d9e66553db970 × 2^(-1022) 和 0.d9e66553db971 × 2^(-1022) 默认情况下四舍五入到尾数为零的最后一位.

请注意，如果 FPU 未配置为仅使用 53 位作为尾数，而是 x87 扩展精度类型的全 64 位，则加法的结果将更接近 IEEE754 结果 0.d9e66553db971 × 2^(-1022) 并因此四舍五入。

实际上，由于 x87 表示具有更大的指数范围，因此与 IEEE754 表示相比，IEEE754 次正规数的尾数有更多的位数，即使尾数中的位数有限。因此，这里的计算结果在 x87 中比在 IEEE754 中多了一位有效位。