MinGW GCC 4.9.1和浮点确定性 [英] MinGW GCC 4.9.1 and floating-point determinism

问题描述

我写了一个小程序来计算3坐标矢量的欧几里得范数.在这里:

I wrote a small program to compute the Euclidean norm of a 3-coordinate vector. Here it is:

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

On my computer, it prints -1.12323e-016.

我知道浮点类型应谨慎处理.但是，我认为浮点运算至少在某种程度上是确定性的. 这篇关于浮点确定性的文章指出:

I know that floating point types should be handled with care. However, I thought that floating-point operations were at least somehow deterministic. This article about floating-point determinism states that:

可以保证的某些事情是加，减，乘，除和平方根的结果.保证这些运算的结果是正确舍入的精确结果(稍后会详细介绍)，因此，如果您提供相同的输入值，相同的全局设置和相同的目标精度，则可以保证得到相同的结果.

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，因为我确实做了对相同的值执行两次相同的操作.

As you can see, the only operations that this program does on floating-point values are addition, subtraction, multiplication and square root. If I trust the article I quoted above, considering that it runs in a single thread and that I do not change the rounding modes or anything else floating-point related, I thought that norm(arr) - norm(arr) would be 0 since I do exactly the same operations on the same values twice.

我的假设是否正确，或者这是编译器在IEEE浮点数学方面不严格符合的情况吗?我当前使用的是MinGW-W64 GCC 4.9.1 32位(我尝试了从-O0到-O3的每个优化级别).与MinGW-W64 GCC 4.8.x相同的程序显示0，这正是我所期望的.

Are my assumptions wrong, or is this a case of compiler not strictly conformant with regard to IEEE floating-point math? I am currently using MinGW-W64 GCC 4.9.1 32 bits (I tried every optimization level from -O0 to -O3). The same program with MinGW-W64 GCC 4.8.x displayed 0, which is what I would have expected.

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

I disassembled the code. I won't post the whole generated assembly because it is too big. However, I believe that the relevant part is here:

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被调用了两次，因此没有内联.虽然我不了解整个事情，也无法说出问题所在.

As you can see, __Z4normIdLj3EET_RKSt5arrayIS0_XT0_EE is called twice, therefore, it is not inlined. I don't understand the whole thing though, and can't tell what is the problem.

推荐答案

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

As @MatthiasB pointed out, this seems to be an issue of gcc temporarily storing an 80 bit floating point value into a 64 bit register/memory location. Consider the following simplified program which still reproduces the issue:

#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编译器)

The assembly code of the essential part norm() - norm() looks like this (using 32 bit mingw 4.8.0 compiler)

...
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错误，但这似乎是一个复杂的主题 ...

Essentially, I would consider this a gcc bug, but it seems to be a complicated topic ...

这篇关于MinGW GCC 4.9.1和浮点确定性的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！