c++ - GNU C++ 中程序的奇怪行为，使用 float

Question

看看这个程序:

#include <iostream>
#include <cmath>

using namespace std;

typedef pair<int, int> coords;

double dist(coords a, coords b)
{
    return sqrt((a.first - b.first) * (a.first - b.first) +
              (a.second - b.second) * (a.second - b.second));
}

int main()
{
    coords A = make_pair(1, 0);
    coords B = make_pair(0, 1);
    coords C = make_pair(-1, 0);
    coords D = make_pair(0, -1);

    cerr.precision(20);
    cerr << dist(A, B) + dist(C, D) << endl;
    cerr << dist(A, D) + dist(B, C) << endl;

    if(dist(A, B) + dist(C, D) > dist(A, D) + dist(B, C))
    {
        cerr << "*" << endl;
    }
    return 0;
}

函数 dist 返回两点之间的距离。A、B、C、D 是正方形的角。

应该是 dist(A, B) == dist(B, C) == dist(C, D) == dist(D, A) == sqrt(2)。

并且 dist(A, B) + dist(C, D) == dist(A, D) + dist(B, C) == 2 * sqrt(2)

我正在使用 GNU/Linux、i586、GCC 4.8.2。

我编译这个程序并运行:

$ g++ test.cpp ; ./a.out 
2.8284271247461902909
2.8284271247461902909
*

我们看到，该程序输出 dist(A, B) + dist(C, D) 和 dist(A, D) + dist(B, C) 的相等值，但条件 dist(A, B) + dist (C, D) > dist(A, D) + dist(B, C) 为真!

当我用 -O2 编译时，它看起来不错:

$ g++ test.cpp -O2 ; ./a.out 
2.8284271247461902909
2.8284271247461902909

我认为，它是一个 gcc-bug，它与浮点运算精度没有直接关系，因为在这种情况下值 dist(A, B) + dist(C, D) 和 dist(A, D ) + dist(B, C) 必须相等(每个都是 sqrt(2) + sqrt(2))。

当我更改函数 dist 时:

double dist(coords a, coords b)
{
    double x = sqrt((a.first - b.first) * (a.first - b.first) + (a.second - b.second) * (a.second - b.second));
    return x;
}

程序运行正确。所以问题不在于数字的浮点表示，而在于 gcc 代码。

编辑:

32 位编译器的简化示例:

#include <iostream>
#include <cmath>
using namespace std;
int main()
{
    if (sqrt(2) != sqrt(2))
    {
        cout << "Unequal" << endl;
    }
    else
    {
        cout << "Equal" << endl;
    }
    return 0;
}

未经优化运行:

$ g++ test.cpp ; ./a.out 
Unequal

使用 -ffloat-store 运行:

$ g++ test.cpp -ffloat-store ; ./a.out 
Equal

解决方案:

可能，它在 GCC#323 中“不是错误”:https://gcc.gnu.org/bugzilla/show_bug.cgi?id=323

用-ffloat-store编译解决问题。

Answer 1

这种看似奇怪的行为是由于旧的 x87 浮点单元的工作方式:它使用 80 位 long double 类型作为其寄存器格式，精度为 64 位，而临时 double 的长度为 64 位，精度为 53 位。发生的事情是编译器将 sqrt(2) 结果中的 1 个溢出到内存中(因为 sqrt 返回一个 double，这四舍五入到 53 - 那种类型的位有效数)，这样 FP 寄存器堆栈就可以为下一次调用 sqrt(2) 清除。然后它将从内存中加载的 53 位精度值与从另一个 sqrt(2) 调用返回的未舍入的 64 位精度值进行比较，然后它们出现不同是因为它们的舍入方式不同，正如您从这个汇编程序输出中看到的那样(我的注释，使用您的第二个代码片段，为了清楚起见，将 2s 更改为 2.0s 和 -Wall -O0 -m32 -mfpmath=387 -march= i586 -fno-builtin 用于 Godbolt 上的编译标志:

main:
    # Prologue
    leal    4(%esp), %ecx
    andl    $-16, %esp
    pushl   -4(%ecx)
    pushl   %ebp
    movl    %esp, %ebp
    pushl   %ecx
    subl    $20, %esp
    # Argument push (2.0)
    subl    $8, %esp
    movl    $0, %eax
    movl    $1073741824, %edx
    pushl   %edx
    pushl   %eax
    # sqrt(2.0)
    call    sqrt
    # Return value spill
    addl    $16, %esp
    fstpl   -16(%ebp)
    # Argument push (2.0)
    subl    $8, %esp
    movl    $0, %eax
    movl    $1073741824, %edx
    pushl   %edx
    pushl   %eax
    # sqrt(2.0)
    call    sqrt
    addl    $16, %esp
    # Comparison -- see how one arg is loaded from a spill slot while the other is
    # coming from the ST(0) i387 register
    fldl    -16(%ebp)
    fucompp
    fnstsw  %ax
    # Status word interpretation
    andb    $69, %ah
    xorb    $64, %ah
    setne   %al
    testb   %al, %al
    # The branch was here, but it and the code below was snipped for brevity's sake

结论:x87 是个怪人。 -mfpmath=sse 是针对此行为的最终修复——它将使 GCC 将 FLT_EVAL_METHOD 定义为 0，因为 SSE(2) 浮点支持是单/双只要。 -ffloat-store 开关也适用于此程序，但不建议将其作为通用解决方法——由于额外的溢出/填充，它会使您的程序变慢，并且不会在所有情况下都有效。当然，使用 64 位 CPU/OS/编译器组合也可以解决这个问题，因为 x86-64 ABI 默认使用 SSE2 进行浮点运算。