MinGW GCC 4.9.1 和浮点确定性

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

#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.

我知道应该小心处理浮点类型。然而，我认为浮点运算至少在某种程度上是确定性的。本文 http://randomascii.wordpress.com/2013/07/16/floating-point-determinism/关于浮点决定论指出：

一些被保证的东西是加法、减法、乘法、除法和平方根的结果。这些操作的结果保证是正确舍入的精确结果（稍后会详细介绍），因此，如果您提供相同的输入值、相同的全局设置和相同的目标精度，则可以保证获得相同的结果。

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

我的假设是错误的，还是编译器不严格遵守 IEEE 浮点数学的情况？我目前正在使用 MinGW-W64 GCC 4.9.1 32 位（我尝试了从-O0 to -O3）。显示与 MinGW-W64 GCC 4.8.x 相同的程序0，这正是我所期望的。

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

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被调用两次，因此，它不是内联的。但我不明白整个事情，也不知道问题出在哪里。

正如@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 bug，但它似乎是一个复杂的话题 https://gcc.gnu.org/bugzilla/show_bug.cgi?id=12331 ...