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