c - 再次浮点精度

昨天我问了一个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时，我们从较低的精度值中得到了正确的答案。当然，如果较低的精度值可以给出正确的答案，那么较高的精度值也应该给出正确的答案？这与“双舍入”有关吗？

在此先感谢您，这对我来说是一个全新的主题，所以我有点挣扎。

最佳答案

关键是，由于指数范围较大，因此两个数字在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位用于有效位数，而将64位用于x87扩展精度类型，则相加的结果将更接近IEEE754结果0.d9e66553db971 × 2^(-1022)，因此会四舍五入。

实际上，由于x87表示具有较大的指数范围，因此即使在有效位数中有限制的情况下，对于IEEE754次正规数的有效位数，您也拥有比IEEE754表示更多的位数。因此，此处的计算结果在x87中比在IEEE754中具有更高的有效位。

关于c - 再次浮点精度，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/15450918/