12

我编写了一个小程序来计算 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它被调用了两次，因此它不是内联的。虽然我不明白整个事情，也无法说出问题所在。

4

2 回答 2

11

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

于 2014-09-17T13:37:13.763 回答
6

浮点数的精度因存储位置而异。如果编译器将一个变量保存在寄存器中，则它作为存储在内存中的变量具有更高的精度。您可以尝试强制将变量存储在内存中，例如：

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。

于 2014-09-17T13:22:15.803 回答