MinGW GCC 4.9.1和浮点确定性

Question

我写了一个小程序来计算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.

---

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

保证的一些事情是加法,减法,乘法,除法和平方根的结果.这些操作的结果保证是正确舍入的精确结果(稍后会详细说明),因此如果您提供相同的输入值,具有相同的全局设置和相同的目标精度,则可以保证相同的结果.

如您所见,此程序对浮点值执行的唯一操作是加法,减法,乘法和平方根.如果我相信我上面引用的文章,考虑到它在单个线程中运行并且我不改变舍入模式或其他浮点相关的东西,我认为那norm(arr) - norm(arr)将是0因为我对相同的值执行完全相同的操作两次.

我的假设是错误的,或者这是编译器在IEEE浮点数学方面不严格符合的情况？我目前使用的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 bug,但它似乎是一个复杂的话题 ......

Answer 2

浮点数的精度存在差异,具体取决于存储位置.如果编译器将一个变量保存在寄存器中,则它作为存储在内存中的变量具有更高的精度.您可以尝试强制将变量存储在内存中,例如:

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

这样可以得到0的预期结果.