GNU C++中程序使用浮点数的奇怪行为

看看这个程序：

#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）。

and 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 编译可以解决问题。

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