c++ - MinGW GCC 4.9.1 和浮点确定性

Question

我编写了一个小程序来计算 3 坐标 vector 的欧几里得范数。这里是:

#include <array>
#include <cmath>
#include <iostream>

template<typename T, std::size_t N>
auto norm(const std::array<T, N>& arr)
    -> T
{
    T res{};
    for (auto value: arr)
    {
        res += value * value;
    }
    return std::sqrt(res);
}

int main()
{
    std::array<double, 3u> arr = { 4.0, -2.0, 6.0 };
    std::cout << norm(arr) - norm(arr) << '\n';
}

在我的电脑上，它打印出 -1.12323e-016。

---

我知道浮点类型应该小心处理。但是，我认为浮点运算至少在某种程度上是确定性的。 This article关于浮点确定性指出:

Some of the things that are guaranteed are the results of addition, subtraction, multiplication, division, and square root. The results of these operations are guaranteed to be the exact result correctly rounded (more on that later) so if you supply the same input value(s), with the same global settings, and the same destination precision you are guaranteed the same result.

如您所见，该程序对浮点值执行的唯一操作是加法、减法、乘法和平方根。如果我相信我上面引用的文章，考虑到它在单个线程中运行并且我没有更改舍入模式或任何其他与浮点相关的内容，我认为 norm(arr) - norm(arr) 将是 0 因为我对相同的值执行完全相同的操作两次。

我的假设是错误的，还是编译器在 IEEE float 学方面不严格符合？我目前正在使用 MinGW-W64 GCC 4.9.1 32 位(我尝试了从 -O0 到 -O3 的每个优化级别)。使用 MinGW-W64 GCC 4.8.x 的同一程序显示 0，这是我所期望的。

---

编辑:我反汇编了代码。我不会发布整个生成的程序集，因为它太大了。但是，我相信相关部分在这里:

call    ___main
fldl    LC0
fstpl   -32(%ebp)
fldl    LC1
fstpl   -24(%ebp)
fldl    LC2
fstpl   -16(%ebp)
leal    -32(%ebp), %eax
movl    %eax, (%esp)
call    __Z4normIdLj3EET_RKSt5arrayIS0_XT0_EE
fstpl   -48(%ebp)
leal    -32(%ebp), %eax
movl    %eax, (%esp)
call    __Z4normIdLj3EET_RKSt5arrayIS0_XT0_EE
fsubrl  -48(%ebp)
fstpl   (%esp)
movl    $__ZSt4cout, %ecx
call    __ZNSolsEd
subl    $8, %esp
movl    $10, 4(%esp)
movl    %eax, (%esp)
call    __ZStlsISt11char_traitsIcEERSt13basic_ostreamIcT_ES5_c
movl    $0, %eax
movl    -4(%ebp), %ecx
.cfi_def_cfa 1, 0
leave

如您所见，__Z4normIdLj3EET_RKSt5arrayIS0_XT0_EE 被调用了两次，因此它不是内联的。不过我不明白整件事，也说不出是什么问题。

Answer 1

正如@MatthiasB 所指出的，这似乎是 gcc 将 80 位浮点值临时存储到 64 位寄存器/内存位置的问题。考虑以下仍然重现问题的简化程序:

#include <cmath>
#include <iostream>

double norm() {
    double res = 4.0 * 4.0 + (-2.0 * -2.0) + (6.0 * 6.0);
    return std::sqrt(res);
}

int main() {
    std::cout << norm() - norm() << '\n';
    return 0;
}

主要部分norm() - norm()的汇编代码如下(使用32位mingw 4.8.0编译器)

...
call    __Z4normv     ; call norm()
fstpl   -16(%ebp)     ; store result (80 bit) in temporary (64 bit!)
call    __Z4normv     ; call norm() again
fsubrl  -16(%ebp)     ; subtract result (80 bit) from temporary (64 bit!)
...

基本上，我认为这是一个 gcc 错误，但它似乎是 complicated topic ...