Общий способ обработки сглаженных значений с плавающей запятой

Question

Вчера я отслеживал ошибку в своем проекте, которая - через несколько часов - я сузился до фрагмента кода, который более или менее делал что-то вроде это:

#include <iostream> 
#include <cmath> 
#include <cassert> 

volatile float r = -0.979541123; 
volatile float alpha = 0.375402451; 

int main() 
{ 
    float sx = r * cosf(alpha); // -0.911326 
    float sy = r * sinf(alpha); // -0.359146 
    float ex = r * cosf(alpha); // -0.911326 
    float ey = r * sinf(alpha); // -0.359146 
    float mx = ex - sx;  // should be 0 
    float my = ey - sy;  // should be 0 
    float distance = sqrtf(mx * mx + my * my) * 57.2958f; // should be 0, gives 1.34925e-06 

// std::cout << "sv: {" << sx << ", " << sy << "}" << std::endl; 
// std::cout << "ev: {" << ex << ", " << ey << "}" << std::endl; 
// std::cout << "mv: {" << mx << ", " << my << "}" << std::endl; 
    std::cout << "distance: " << distance << std::endl; 

    assert(distance == 0.f); 
// assert(sx == ex && sy == ey); 
// assert(mx == 0.f && my == 0.f); 
} 

После компиляции и выполнения:

$ g++ -Wall -Wextra -Wshadow -march=native -O2 vfma.cpp && ./a.out 
distance: 1.34925e-06 
a.out: vfma.cpp:23: int main(): Assertion `distance == 0.f' failed. 
Aborted (core dumped) 

с моей точки зрения что-то не так, как я просил 2 вычитаний двух пар поразрядными идентичных (я ожидал получить два нули), затем их возведение в квадрат (два нуля снова) и их объединение (ноль).

Оказалось, что основной причиной проблемы является использование операции с объединенным многократным добавлением, которая где-то вдоль линии делает результат неточным (с моей точки зрения). Вообще-то я ничего не имею против этой оптимизации, так как обещает дать результаты, которые больше точнее, но в этом случае 1.34925e-06 действительно далека от ожидаемого мною 0.

Тестовый чехол очень «хрупкий» - если вы включаете больше отпечатков или больше утверждений, он перестает утверждать, потому что компилятор больше не использует fused-multiply-add. Например, если я раскомментировать все строки:

$ g++ -Wall -Wextra -Wshadow -march=native -O2 vfma.cpp && ./a.out 
sv: {-0.911326, -0.359146} 
ev: {-0.911326, -0.359146} 
mv: {0, 0} 
distance: 0 

Как я считал это ошибка в компиляторе, я сообщил, что, но он был закрыт с объяснением, что это правильное поведение.

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=79436

Так что я задаюсь вопросом - каким образом следует один код таких вычислений, чтобы избежать этой проблемы? Я думал о том, общее решение, но что-то лучше, чем:

mx = ex != sx ? ex - sx : 0.f; 

Я хотел бы, чтобы исправить или улучшить свой код - если есть что-нибудь, чтобы исправить/улучшить - вместо установки -ffp-contract=off для всего моего проекта, fused- multiply-add используется внутри библиотек компилятора в любом случае (я вижу много всего в sinf() и cosf()), поэтому это будет «частичное обходное», а не решение ... Мне также хотелось бы чтобы избежать решений, таких как «не использовать с плавающей точкой» (;

Answer 1

в общем нет: это именно цена, которую вы платите за пользование -ffp-contract=fast (совпадению, именно этот пример, William Kahan notes in the problems with automatic contraction)

Теоретически, если вы использовали C (не C++), а ваш компилятор поддерживал прагмы C-1999 (т. не ССАГПЗ), вы можете использовать

#pragma STDC FP_CONTRACT OFF 
// non-contracted code 
#pragma STDC FP_CONTRACT ON 

Answer 2

Интересно, что благодаря FMA, поплавки тх и мой дает ошибку округления, которая была сделана при умножении г и соз.

fma(r,cos, -r*cos) = theoretical(r*cos) - float(r*cos) 

Так что результат вы получите как-то указывает на то, насколько далеко был вычисленный (Sx, Sy) от теоретического (Sx, Sy), из-за умножения поплавками (но без учета из-за ошибок округления при расчетах сов и син).

Итак, вопрос в том, как ваша программа может полагаться на разницу (ex-sx, ey-sy), которая находится в интервале неопределенности, связанной с округлением с плавающей запятой?