Недавно я вернулся к анализу погрешностей чисел с плавающей запятой, чтобы усовершенствовать некоторые детали в следующей редакции книги Physically Based Rendering. Числа с плавающей запятой — интересная область вычислений, полная сюрпризов (хороших и плохих), а также хитрых трюков, позволяющих избавиться от неприятных неожиданностей.
В процессе работы я наткнулся на этот пост на StackOverflow, из которого узнал об изящном алгоритме точного вычисления
.
Но прежде чем приступать к алгоритму, нужно понять, что же такого хитрого в выражении
? Возьмём
,
,
и
. (Это реальные значения, которые получились у меня во время запуска pbrt.) При 32-битных значениях float получаем:
и
. Выполняем вычитание, и получаем
. Но если выполнить вычисления с двойной точностью, а в конце преобразовать их во float, то получится
. Что произошло?
Проблема в том, что значение каждого произведения может сильно выйти за нижнюю границу
, где расстояние между представимыми значениями с плавающей запятой очень велико — 64. То есть при округлении
и
по отдельности до ближайшего представимого float, они превращаются в числа, кратные 64. В свою очередь, их разность будет кратной 64, и не останется никакой надежды, что она станет к
ближе, чем
. В нашем случае результат оказался ещё дальше из-за того, как два произведения были округлены в
. Мы напрямую столкнёмся со старым добрым катастрофическим сокращением1.
Читать полностью »