Неправильный результат для atan2 с glibc / libm и float32

Question

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

Чтобы смоделировать поведение и доказать правильность моих формул и кода, я перенес соответствующую математическую частьиз встроенного ПО для инструментария GCC в Linux (gcc 6.3.0, libc6 2.24), дважды проверяя, что float32 используется повсеместно и что не используется переключатель компилятора, который может снизить точность или совместимость стандартов математических операций;Примечательно, что нет ни одного -ffast-math или его друзей.

Теперь оказалось, что я получаю неожиданные результаты для небольшого набора входных параметров.Я разыскал проблему и пришел к выводу, что libm вычисляет неверный результат для arctan (если быть точным: atan2) для очень небольшого набора входных параметров.

Например,, если у меня есть

#include <math.h>

#define C_RAD2DEG (57.29577951308f)

int main(void)
{
  float f_Temp = C_RAD2DEG * atan2f(0.713114202f, 0.665558934f);
}

f_Temp вычисляется как 46.9755516f, где правильный результат будет 46.975548972f.

Обратите внимание, что я в целом осведомлен о проблемахс различными типами данных с плавающей запятой, ошибками округления и т. п.

Однако мне кажется, что показанная выше ошибка слишком велика на порядок даже при низкой точности float32, и, к сожалению,для последующих расчетов эта ошибка слишком велика.

Кроме того, проблема связана только с очень небольшим подмножеством возможных входных параметров для функции atan2.

Может кто-нибудь пожалуйстакратко объясните, является ли это ошибкой в ​​libm или это просто из-за неточности float32 и большого количества последовательных операций, необходимых для вычисления atan2?

Answer 1

Число, которое вы сообщаете как наблюдаемый результат, 46.9755516f, соответствует значению float 46,975551605224609375.

Число, которое вы сообщаете как ожидаемый результат, 46.975548972f, соответствует floatзначение 46,97554779052734375.

Это смежные float значения, означающие, что они отличаются на 1 единицу наименьшей точности (ULP).(Их разность составляет 3.814697265625e-06, что является значением младшего значащего бита в значении float, и когда старший значащий бит имеет значение 32, как это имеет место для чисел около 47.) Это наименьшая возможная величина, на которуюfloat может измениться в этом масштабе.

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

Даже если используемый вами код libc дал правильно округленный результат, преобразование его из радиан в градусы приводит к еще двум ошибкам округления (преобразование 180 / π в aпредставительное значение и умножение на него).Неразумно ожидать, что конечным результатом будет float, ближайший к идеальному математическому результату;Вы должны ожидать несколько ошибок ULP.