Точность с плавающей точкой IEEE-754: насколько допустима ошибка?

Я работаю над переносом функции sqrt (для 64-разрядных удвоений) из fdlibm в инструмент проверки модели Я использую в данный момент (cbmc).
В рамках моих работ я много читал о стандарте ieee-754, но, думаю, я не понял гарантии точности для основных операций (включая sqrt).

Тестирование моего порта fdlibm sqrt, я получил следующий расчет с sqrt в 64-битном двойном:

sqrt(1977061516825203605555216616167125005658976571589721139027150498657494589171970335387417823661417383745964289845929120708819092392090053015474001800648403714048.0) = 44464159913633855548904943164666890000299422761159637702558734139742800916250624.0

(этот случай нарушил простое пост-условие в моем тесте относительно точности, я не уверен, если это пост-условие возможно с IEEE-754)

Для сравнения несколько инструментов с несколькими точками вычисляли что-то вроде:

sqrt(1977061516825203605555216616167125005658976571589721139027150498657494589171970335387417823661417383745964289845929120708819092392090053015474001800648403714048.0) =44464159913633852501611468455197640079591886932526256694498106717014555047373210.truncated

Можно видеть, что 17-е число слева отличается, что означает ошибку, например:

3047293474709469249920707535828633381008060627422728245868877413.0

Вопрос 1: разрешено ли это огромное количество ошибок?

В стандарте говорится, что каждая базовая операция (+, -, *,/, sqrt) должна быть в пределах 0,5 ulps, что означает, что она должна быть равна математически точному результату, округленному до ближайшего fp-представления (wiki говорит что в некоторых библиотеках только 1 ulp, но это не так важно на данный момент).

Вопрос 2: Означает ли это, что каждая базовая операция должна иметь ошибку < 2.220446e-16 с 64-битным удвоением (машинный-эпсилон)?

Я подсчитал то же самое с Linux-системой x86-32 (glibc/eglibc) и получил тот же результат, что и с fdlibm, что позволяет мне думать, что:

a: Я сделал что-то неправильно (но как: printf был кандидатом, но я не знаю, может ли это быть причиной).
b: ошибка/точность распространена в этих библиотеках.

Ответ 1

Стандарт IEEE-754 требует, чтобы так называемые "основные операции" (включая добавление, умножение, деление и квадратный корень) были правильно округлены. Это означает, что существует однозначный разрешенный ответ, и он является самым близким представляемым числом с плавающей запятой к так называемому "бесконечно точному" результату операции.

В двойной точности числа имеют 53 двоичных цифры точности, поэтому правильный ответ - это точный ответ, округленный до 53 значащих цифр. Поскольку Рик Риган показал в своем ответе, это именно то, что вы получили.

Ответы на ваши вопросы:

Вопрос 1: разрешено ли это огромное количество ошибок?

Да, но довольно ошибочно назвать эту ошибку "огромной". Дело в том, что нет значения двойной точности, которое может быть возвращено, что будет иметь меньшую ошибку.

Вопрос 2: Означает ли это, что каждая базовая операция должна иметь ошибку < 2.220446e-16 с 64-битным удвоением (машинный-эпсилон)?

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

Вопрос 3: Какой результат вы получите с помощью своего оборудования x86 и gcc + libc?

Тот же ответ, который вы сделали, потому что sqrt правильно округлен на любой разумной платформе.

Ответ 2

В двоичном выражении первые 58 бит произвольного точного ответа составляют 1011111111111111111111110101010101111111110111111011010001...

53-битное значение двойного значения

10111111111111111111111101010101011111111111111110111

Это означает, что двойное значение правильно округлено до 53 значащих бит и находится в пределах 1/2 ULP. (Что ошибка "большая" - это только потому, что само число является большим).