浮点运算的准确性

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编号.前五张照片显示的数字是我预期的两倍.令我困惑的是,接下来的五个版画显示不同的数字.如果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双精度计算.

不幸的是,"尽可能接近"并不总是足够接近.为了模拟双精度,历史80位浮点寄存器的有效位置可以限制其宽度,但它们始终保留指数的整个范围.这种差异尤其表现在操纵非正规(如你的y)时.

怎么了

我的解释是,in printf("%23.16e\n", 1.6*y);,1.6*y被计算为80位递减有效数和全指数(因此是正常数),然后转换为IEEE 754双精度(导致非正规),然后打印.

另一方面,in printf("%23.16e\n", x + 1.6*y);,x + 1.6*y用所有80位有效数字和全指数数字计算(同样所有中间结果都是正常数字),然后转换为IEEE 754双精度,然后打印.

这可以解释为什么1.6*y打印相同2.0*y但添加时效果不同x.打印的数字是双精度非正规.添加的数字是x80位减少有效数和全指数正常数(不同).

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

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

David Monniaux撰写的这篇文章包含了您可能需要的所有细节以及更多细节.

删除不需要的行为

要摆脱这个问题(如果你认为它是一个),找到告诉你的编译器为浮点生成SSE2指令的标志.这些实现了单精度和双精度的IEEE 754语义.