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

Question

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

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

sqrt(1977061516825203605555216616167125005658976571589721139027150498657494589171970335387417823661417383745964289845929120708819092392090053015474001800648403714048.0) = 44464159913633855548904943164666890000299422761159637702558734139742800916250624.0

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

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

sqrt(1977061516825203605555216616167125005658976571589721139027150498657494589171970335387417823661417383745964289845929120708819092392090053015474001800648403714048.0) =44464159913633852501611468455197640079591886932526256694498106717014555047373210.truncated

Видно, что 17-е число слева отличается, что означает ошибку типа:

3047293474709469249920707535828633381008060627422728245868877413.0

Вопрос 1. Допускается ли это огромное количество ошибок?

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

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

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

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

Answer 1

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

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

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

Вопрос 1. Допускается ли это огромное количество ошибок?

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

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

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

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

Answer 2

В двоичном формате первые 58 бит ответа с произвольной точностью: 1011111111111111111111110101010111111111111111111011010001 ...

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

10111111111111111111111101010101011110111010101111101111111 * 100означает, что двойное значение правильно округлено до 53 значащих бит и находится в пределах 1/2 ULP.(То, что ошибка «большая», только потому, что само число большое).