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.

Я знаю, что с плавающей точкой следует обращаться с осторожностью. Однако я думал, что операции с плавающей запятой были как бы детерминированными. В этой статье о детерминизме с плавающей запятой говорится, что:

Некоторые из вещей, которые гарантированы, являются результатом сложения, вычитания, умножения, деления и квадратного корня. Результаты этих операций гарантированно будут корректно округлены (более подробно об этом позже), поэтому, если вы укажете одинаковые входные значения, с теми же глобальными настройками и с той же конечной точностью, вам гарантирован тот же результат.

Как вы можете видеть, единственными действиями, которые эта программа выполняет для значений с плавающей запятой, являются сложение, вычитание, умножение и квадратный корень. Если я доверяю цитированной выше статье, считая, что она работает в одном потоке и что я не изменяю режимы округления или что-то еще связанное с плавающей запятой, я думал, что norm(arr) - norm(arr) будет 0, так как я точно те же операции с теми же значениями дважды.

Являются ли мои допущения неправильными, или это случай компилятора не строго соответствует математике с плавающей запятой IEEE? В настоящее время я использую MinGW-W64 GCC 4.9.1 32 бита (я пробовал каждый уровень оптимизации от -O0 до -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 вызывается дважды, поэтому он не встроен. Я все это не понимаю, и не могу понять, в чем проблема.

Ответ 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, но, похоже, это сложная тема...

Ответ 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.