开启优化后浮点数结果不同 - 编译器 bug？

Question

以下代码在Visual Studio 2008上可以正常运行，无论是否进行优化。但是只有在没有优化（O0）的情况下才能在g++上运行。

#include <cstdlib>
#include <iostream>
#include <cmath>

double round(double v, double digit)
{
    double pow = std::pow(10.0, digit);
    double t = v * pow;
    //std::cout << "t:" << t << std::endl;
    double r = std::floor(t + 0.5);
    //std::cout << "r:" << r << std::endl;
    return r / pow;
}

int main(int argc, char *argv[])
{
    std::cout << round(4.45, 1) << std::endl;
    std::cout << round(4.55, 1) << std::endl;
}

输出应该是：

4.5
4.6

但是使用优化的g++（O1 - O3）将输出：

4.5
4.5

如果我在t之前加上volatile关键字，它就能正常工作，那么可能存在某种优化错误吗？

在g++ 4.1.2和4.4.4上进行测试。

这是在ideone上的结果： http://ideone.com/Rz937 我在g++上测试的选项很简单：

g++ -O2 round.cpp

更有趣的是，即使我在Visual Studio 2008中打开了/fp:fast选项，结果仍然正确。

进一步的问题：

我想知道，我是否应该总是打开-ffloat-store选项？

因为我测试的g++版本是随CentOS/Red Hat Linux 5和CentOS/Redhat 6一起提供的。

我在这些平台下编译了许多程序，担心会在我的程序中引起意外的错误。看起来要调查所有我的C++代码和使用的库是否存在这样的问题有点困难。有什么建议吗？

是否有人对为什么即使打开了/fp:fast，Visual Studio 2008仍然可以工作感兴趣？看起来Visual Studio 2008在这个问题上比g++更可靠？

- Bear

---

给所有新的SO用户：这就是你如何提问。+1 - tenfour

ideone使用的是4.3.4版本。http://ideone.com/b8VXg - Daniel A. White

请记住，你的例程很可能无法可靠地处理所有类型的输出。与将双精度浮点数四舍五入为整数不同，这种方法容易受到无法表示所有实数的事实的影响，因此你应该预计会出现更多类似于这样的错误。 - Jakub Wieczorek

对于那些无法重现错误的人：请不要取消注释已注释的调试语句，它们会影响结果。 - n. m.

有趣的是，我在2006年遇到了类似的问题。当时，我被引导到以下错误报告：http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323。 - Throwback1986

Answer 1

Intel x86处理器在内部使用80位扩展精度，而double通常是64位宽。不同的优化级别会影响从CPU获取浮点值保存到内存的频率，从而将其从80位精度舍入为64位精度。

使用-ffloat-store gcc选项可以在不同的优化级别下获得相同的浮点结果。

或者，使用long double类型，在gcc上通常是80位宽，避免从80位精度舍入到64位精度。 man gcc说得很清楚：

   -ffloat-store
       Do not store floating point variables in registers, and inhibit
       other options that might change whether a floating point value is
       taken from a register or memory.

       This option prevents undesirable excess precision on machines such
       as the 68000 where the floating registers (of the 68881) keep more
       precision than a "double" is supposed to have.  Similarly for the
       x86 architecture.  For most programs, the excess precision does
       only good, but a few programs rely on the precise definition of
       IEEE floating point.  Use -ffloat-store for such programs, after
       modifying them to store all pertinent intermediate computations
       into variables.

---

在x86_64编译中，默认情况下编译器会使用SSE寄存器来处理 float 和 double ，因此不会使用扩展精度，也不会出现此问题。 gcc编译器选项-mfpmath 可控制此行为。

- Maxim Egorushkin

---

我认为这就是答案。常数4.55被转换成64位中最接近的二进制表示形式4.54999999999999；乘以10并再次舍入到64位，得到45.5。如果你在一个80位寄存器中保留它而跳过舍入步骤，你得到的是45.4999999999999。 - Mark Ransom

@Bear，调试语句可能会导致某个值从寄存器刷新到内存中。 - Mark Ransom

@Bear，通常情况下，你的应用程序应该从扩展精度中受益，除非它处理极小或极大的值时，64位浮点数可能会发生下溢或上溢并产生“inf”。没有一个好的经验法则，单元测试可以给你一个明确的答案。 - Maxim Egorushkin

一般而言，如果需要完全预测的结果和/或完全按照人类手算得出的结果，则应避免使用浮点数。-ffloat-store可以消除一些不确定性，但这并不是一个万能方案。 - plugwash

@plugwash 如果你必须使用浮点数，那么这个规则就没有帮助。 - Maxim Egorushkin

Answer 2

输出应该是：4.5 4.6 如果您拥有无限精度或使用基于十进制而不是基于二进制的浮点表示的设备，则输出将是这样。但是，大多数计算机使用二进制IEEE浮点标准。正如Maxim Yegorushkin在他的答案中已经指出的那样，问题的一部分是您的计算机在内部使用80位浮点表示。然而，这只是问题的一部分。问题的根源是任何形式为n.nn5的数字都没有一个精确的二进制浮点表示。这些边角情况总是不精确的数字。如果您真的希望舍入能够可靠地舍入这些边角情况，您需要使用一个舍入算法来解决n.n5、n.nn5或n.nnn5等（但不是n.5）始终不精确的事实。找到决定某个输入值向上舍入还是向下舍入的边角情况，并根据与此边角情况的比较返回向上舍入或向下舍入的值。您确实需要注意优化编译器不会将发现的边角情况放在扩展精度寄存器中。请参见How does Excel successfully Rounds Floating numbers even though they are imprecise?以获取此类算法。或者，您可以接受边角情况有时会错误地舍入的事实。

- David Hammen

Answer 3

不同的编译器有不同的优化设置。其中一些更快的优化设置不遵守IEEE 754的严格浮点规则。Visual Studio有一个特定的设置，/fp:strict，/fp:precise，/fp:fast，其中/fp:fast违反了标准上可执行的操作。你可能会发现这个标志是控制此类设置中的优化的。您还可以在GCC中找到类似的设置，可以改变行为。

如果是这种情况，那么编译器之间唯一不同的就是，在更高级别的优化中，GCC默认寻找最快的浮点行为，而Visual Studio不会随着更高的优化级别改变浮点行为。因此，这可能并不一定是实际的错误，而是您没有意识到开启的选项的预期行为。

- Puppy

---

GCC有一个-ffast-math开关，它不会被任何-O优化级别打开，因为引用的话：“对于那些依赖于IEEE或ISO数学规则/规范的程序，它可能导致输出不正确。” - Mat

@Mat：我已经尝试在我的g++ 4.4.3上使用-ffast-math和其他一些方法，但仍无法重现这个问题。 - NPE

很好：使用-ffast-math，在优化级别大于0的情况下，两种情况都可以得到4.5。 - Kerrek SB

（更正：在GCC 4.4.3中，我使用“-O1”和“-O2”得到了“4.5”，但在“-O0”和“-O3”下没有得到。但是，在GCC 4.6.1中，“-O1,2,3”都可以得到“4.5”。） - Kerrek SB

Answer 4

对于那些无法重现错误的人：不要取消注释已注释掉的调试语句，它们会影响结果。这意味着问题与调试语句有关。看起来在输出语句期间将值加载到寄存器中导致了舍入误差，这就是为什么其他人发现你可以通过使用-ffloat-store来解决这个问题。

进一步的问题： 我想知道是否应该始终打开-ffloat-store选项？

说句俏皮话，某些程序员不打开-ffloat-store肯定有原因，否则该选项就不存在（同样，某些程序员打开-ffloat-store也有原因）。我不建议总是打开或关闭它。打开它会防止一些优化，但关闭它允许您获得您正在获取的行为。

但是，一般来说，二进制浮点数（计算机使用的）和十进制浮点数（人们熟悉的）之间存在某些不匹配，这种不匹配可能会导致类似于您获得的行为（要明确的是，您获得的行为不是由这种不匹配引起的，但是类似的行为可以）。问题是，由于在处理浮点数时已经存在一些模糊性，我不能说-ffloat-store会使情况变得更好或更糟。

相反，您可能需要寻找其他解决方案来解决您正在尝试解决的问题（不幸的是，Koenig没有指向实际论文的具体位置，而且我也找不到一个明显的“规范”位置，所以我必须将您发送到Google）。

---

如果您不是为了输出而四舍五入，我可能会考虑使用cmath中的std::modf()和limits中的std::numeric_limits<double>::epsilon()。思考原始的round()函数，我认为用这个函数替换对std::floor(d + .5)的调用会更清晰：

// this still has the same problems as the original rounding function
int round_up(double d)
{
    // return value will be coerced to int, and truncated as expected
    // you can then assign the int to a double, if desired
    return d + 0.5;
}

我认为这表明了以下改进：

// this won't work for negative d ...
// this may still round some numbers up when they should be rounded down
int round_up(double d)
{
    double floor;
    d = std::modf(d, &floor);
    return floor + (d + .5 + std::numeric_limits<double>::epsilon());
}

一个简单的提示：`std::numeric_limits::epsilon()`被定义为“添加到1的最小数，使得得到的数字不等于1。”通常需要使用相对epsilon（即，以某种方式缩放epsilon以考虑您正在处理的不是“1”的数字）。`d`、`.5`和`std::numeric_limits::epsilon()`的总和应该接近1，因此将它们分组意味着`std::numeric_limits::epsilon()`的大小大约适合我们所做的事情。如果有什么不同，`std::numeric_limits::epsilon()`会太大（当所有三个数字的总和小于1时），并可能导致我们将某些数字四舍五入。

---

现在，你应该考虑std::nearbyint()。

- Max Lybbert

---

“相对 epsilon” 称为 1 ulp（最后一位的 1 单位）。x - nextafter(x, INFINITY)与 x 的 1 ulp 有关（但不要使用它；我确定存在边界情况，我只是举了个例子）。epsilon() 的 cppreference 示例演示了将其缩放以获取基于 ULP 的相对误差。 - Peter Cordes

顺便提一下，2016年对于-ffloat-store的答案是：首先不要使用x87。使用SSE2数学运算（64位二进制文件或-mfpmath=sse -msse2用于制作老旧的32位二进制文件），因为SSE/SSE2具有没有额外精度的临时变量。 XMM寄存器中的double和float变量实际上采用IEEE 64位或32位格式。（与x87不同，其中寄存器始终为80位，并且存储到内存会四舍五入为32位或64位。） - Peter Cordes

Answer 5

如果您编译的目标不包括SSE2，那么被接受的答案是正确的。所有现代x86处理器都支持SSE2，因此如果您能利用它，就应该这样做：

-mfpmath=sse -msse2 -ffp-contract=off

让我们来分解一下。

-mfpmath=sse -msse2。使用SSE2寄存器进行四舍五入，比将每个中间结果存储到内存中要快得多。请注意，在x86-64上，这已经是GCC的默认选项。来自GCC wiki：

在支持SSE2的更现代的x86处理器上，指定编译器选项-mfpmath=sse -msse2可确保所有float和double操作都在SSE寄存器中执行并正确舍入。这些选项不影响ABI，因此应尽可能用于预测性数字结果。

-ffp-contract=off。然而，仅控制舍入并不能实现完全匹配。FMA（融合乘加）指令与其非融合对应物相比可以改变舍入行为，因此我们需要禁用它。这是Clang的默认设置，而不是GCC。正如这个答案所解释的那样：

一个FMA只有一个舍入（它实际上保持内部临时乘积结果的无限精度），而ADD + MUL有两个。

通过禁用FMA，我们可以在调试和发布时获得完全匹配的结果，但代价是一些性能（和精度）。我们仍然可以利用SSE和AVX的其他性能优势。

- tmandry

---

在提到FMA指令可能存在问题时，要给予+1的评价。此外，在某些平台上，SSE指令也可能存在缺陷，因此尽可能使用-msse2而不是-msse。 - ManuelAtWork

Answer 6

我进一步研究了这个问题，可以提供更多的细节。首先，在x84_64上，根据gcc的精确表示，4.45和4.55的确切表示如下（使用libquadmath打印最后一位精度）：

float 32:   4.44999980926513671875
double 64:  4.45000000000000017763568394002504646778106689453125
doublex 80: 4.449999999999999999826527652402319290558807551860809326171875
quad 128:   4.45000000000000000000000000000000015407439555097886824447823540679418548304813185723105561919510364532470703125

float 32:   4.55000019073486328125
double 64:  4.54999999999999982236431605997495353221893310546875
doublex 80: 4.550000000000000000173472347597680709441192448139190673828125
quad 128:   4.54999999999999999999999999999999984592560444902113175552176459320581451695186814276894438080489635467529296875

正如Maxim所说，问题是由于FPU寄存器的80位大小引起的。

但为什么在Windows上从未出现过这个问题？在IA-32上，x87 FPU被配置为使用53位（等效于总大小为64位：double）的尾数内部精度。对于Linux和Mac OS，默认精度为64位（等效于总大小为80位：long double）。因此，通过更改FPU的控制字（假设指令序列会触发错误），应该可以在这些不同的平台上可能出现或不出现该问题。该问题已报告给gcc作为bug 323（至少阅读第92条评论！）。

要显示Windows上的尾数精度，可以使用VC ++以32位编译此内容：

#include "stdafx.h"
#include <stdio.h>  
#include <float.h>  

int main(void)
{
    char t[] = { 64, 53, 24, -1 };
    unsigned int cw = _control87(0, 0);
    printf("mantissa is %d bits\n", t[(cw >> 16) & 3]);
}

在 Linux/Cygwin 上：

#include <stdio.h>

int main(int argc, char **argv)
{
    char t[] = { 24, -1, 53, 64 };
    unsigned int cw = 0;
    __asm__ __volatile__ ("fnstcw %0" : "=m" (*&cw));
    printf("mantissa is %d bits\n", t[(cw >> 8) & 3]);
}

请注意，使用gcc可以使用-mpc32/64/80设置FPU精度，但在Cygwin中会被忽略。但请记住，它将修改尾数的大小，而不是指数，这可能导致其他不同行为的问题。

在x86_64架构上，如tmandry所说，使用SSE，因此除非您强制使用旧的x87 FPU进行FP计算（使用-mfpmath=387），或者编译为32位模式（您需要multilib包），否则不会出现问题。我可以在Linux上使用不同的标志和gcc版本组合来复现该问题。

g++-5 -m32 floating.cpp -O1
g++-8 -mfpmath=387 floating.cpp -O1

我在Windows或Cygwin上尝试了几种VC++/gcc/tcc的组合，但从未出现过这个错误。我想生成的指令序列不同。

最后，请注意，防止4.45或4.55出现此问题的一种奇特方法是使用_Decimal32/64/128，但支持非常稀少...我花了很多时间才能使用libdfp进行printf！

- calandoa

Answer 7

就个人而言，我曾经遇到过相反的问题——从GCC到VS。在大多数情况下，我认为最好避免优化。唯一值得进行优化的时候是处理涉及大量浮点数据数组的数值方法。即使是反汇编代码之后，我经常对编译器的选择感到失望。很多时候，使用编译器内置函数或者直接编写汇编语言会更容易。

- cdcdcd