c - 浮点运算的精度

Question

我无法理解这个程序的输出

int main()
{
    double x = 1.8939201459282359e-308;
    double y = 4.9406564584124654e-324;
    printf("%23.16e\n", 1.6*y);
    printf("%23.16e\n", 1.7*y);
    printf("%23.16e\n", 1.8*y);
    printf("%23.16e\n", 1.9*y);
    printf("%23.16e\n", 2.0*y);
    printf("%23.16e\n", x + 1.6*y);
    printf("%23.16e\n", x + 1.7*y);
    printf("%23.16e\n", x + 1.8*y);
    printf("%23.16e\n", x + 1.9*y);
    printf("%23.16e\n", x + 2.0*y);
}

输出是

9.8813129168249309e-324
9.8813129168249309e-324
9.8813129168249309e-324
9.8813129168249309e-324
9.8813129168249309e-324
1.8939201459282364e-308
1.8939201459282364e-308
1.8939201459282369e-308
1.8939201459282369e-308
1.8939201459282369e-308

我正在使用 IEEE 算法。变量 y 包含可能的最小 IEEE 编号。前五张打印品显示的数字是我预期的两倍 y 。令我困惑的是接下来的五张打印品显示不同的数字。如果 1.6*y 与 2.0*y 相同，那么 x + 1.6*y 如何不同于 x + 2.0* y?

Answer 1

一言以蔽之

你说你的编译器是Visual C++ 2010 Express。我无权访问此编译器，但我知道它会生成最初将 x87 CPU 配置为使用 53 位精度的程序，以便尽可能接近地模拟 IEEE 754 double 计算。

不幸的是，“尽可能接近”并不总是足够接近。出于模拟 double 的目的，历史上的 80 位浮点寄存器可以限制有效数的宽度，但它们始终保留指数的完整范围。这种差异在处理非正规化(如您的 y)时尤其明显。

发生了什么

我的解释是，在 printf("%23.16e\n", 1.6*y); 中，1.6*y 被计算为 80 位缩减-significand full-exponent number(因此它是一个普通数)，然后转换为 IEEE 754 double (导致非正规数)，然后打印。

另一方面，在 printf("%23.16e\n", x + 1.6*y); 中，x + 1.6*y 是用所有80 位缩减尾数全指数(所有中间结果都是正常数)，然后转换为 IEEE 754 double ，然后打印。

这可以解释为什么 1.6*y 打印与 2.0*y 相同，但添加到 x 时具有不同的效果。打印的数字是 double 非正规数字。添加到 x 上的数字是一个 80 位缩减有效全指数正规数(不是同一个)。

生成 x87 指令时其他编译器会发生什么

其他编译器，如 GCC，不会配置 x87 FPU 来操作 53 位尾数。这可能会产生相同的结果(在这种情况下，x + 1.6*y 将使用所有 80 位完整有效数字完整指数进行计算，然后转换为 double 以打印或存储在内存中) .在这种情况下，问题会更加明显(您不需要涉及非正规数或无限数来注意到差异)。

这article由 David Monniaux 撰写，包含您可能希望了解的所有详细信息以及更多内容。

移除不需要的行为

要解决这个问题(如果您认为它是一个问题)，请找到告诉您的编译器为浮点生成 SSE2 指令的标志。它们完全实现了单精度和 double 的 IEEE 754 语义。