Алгоритм Кэхэна: как получить точную разность произведений

Недавно я вернулся к анализу погрешностей чисел с плавающей запятой, чтобы усовершенствовать некоторые детали в следующей редакции книги Physically Based Rendering. Числа с плавающей запятой — интересная область вычислений, полная сюрпризов (хороших и плохих), а также хитрых трюков, позволяющих избавиться от неприятных неожиданностей.

В процессе работы я наткнулся на этот пост на StackOverflow, из которого узнал об изящном алгоритме точного вычисления $$.

Но прежде чем приступать к алгоритму, нужно понять, что же такого хитрого в выражении $$? Возьмём $$, $$, $$ и $$. (Это реальные значения, которые получились у меня во время запуска pbrt.) При 32-битных значениях float получаем: $$ и $$. Выполняем вычитание, и получаем $$. Но если выполнить вычисления с двойной точностью, а в конце преобразовать их во float, то получится $$. Что произошло?

Проблема в том, что значение каждого произведения может сильно выйти за нижнюю границу $$, где расстояние между представимыми значениями с плавающей запятой очень велико — 64. То есть при округлении $$ и $$ по отдельности до ближайшего представимого float, они превращаются в числа, кратные 64. В свою очередь, их разность будет кратной 64, и не останется никакой надежды, что она станет к $$ ближе, чем $$. В нашем случае результат оказался ещё дальше из-за того, как два произведения были округлены в $$. Мы напрямую столкнёмся со старым добрым катастрофическим сокращением¹.

Вот решение получше²:

inline float DifferenceOfProducts(float a, float b, float c, float d) {
    float cd = c * d;
    float err = std::fma(-c, d, cd);
    float dop = std::fma(a, b, -cd);
    return dop + err;
}

DifferenceOfProducts() вычисляет $$ таким образом, что избегает катастрофического сокращения. Впервые эта методика была описана легендарным Уильямом Кэхэном в статье On the Cost of Floating-PointComputation Without Extra-Precise Arithmetic. Стоит заметить, что работы Кэхэна интересно читать в целом, в них есть множество комментариев о текущем состоянии мира плавающих запятых, а также математических и технических рассуждений. Вот один из его выводов:

Те из нас, кто сражался с превратностями арифметики с плавающей запятой и плохо продуманными «оптимизациями» компиляторов, могут справедливо испытывать гордость за победу в этой битве. Но если мы передадим продолжение этой битвы последующим поколениям, это будет противоречить всей сути цивилизации. Наш опыт говорит, что языки программирования и системы разработки являются источниками слишком большого количества хаоса, с которым нам приходится бороться. Слишком от многих ошибок вполне можно обойтись, как и от некоторых привлекательных «оптимизаций», безопасных для целых чисел, но иногда оказывающихся фатальными для чисел с плавающей запятой.

Отдав должное его остроумию, вернёмся к DifferenceOfProducts(): в основе искусности этой функции лежит использование инструкций совмещённого умножения-сложения (fused multiply-add, FMA)³. С математической точки зрения FMA(a,b,c) является $$, поэтому поначалу кажется, что эта операция полезна только в качестве микрооптимизаци: одна инструкция вместо двух. Однако FMA
обладает особым свойством — она округляет только один раз.

В обычном $$ сначала вычисляется $$, а потом это значение, которое в общем случае нельзя представить в формате с плавающей запятой, округляется до ближайшего float. Затем к этому округлённому значению прибавляется $$, и этот результат снова округляется до ближайшего float. FMA реализована таким образом, что округление выполняется только в конце — промежуточное значение $$ сохраняет достаточную точность, поэтому после прибавления к нему $$ готовый результат будет являться ближайшей к истинному значению $$ величиной float.

Разобравшись с FMA, вернёмся к DifferenceOfProducts(). Снова покажу первые две её строки:

    float cd = c * d;
    float err = std::fma(-c, d, cd);

Первая вычисляет округлённое значение $$, а вторая… вычитает $$ из их произведения? Если не знать, как работают FMA, то можно подумать, что err всегда будет равна нулю. Но при работе с FMA вторая строка на самом деле извлекает величину погрешности округления в вычисленном значении $$ и сохраняет её в err. После этого на выходе всё получается очень просто:

    float dop = std::fma(a, b, -cd);
    return dop + err;

Вторая FMA вычисляет разность произведений при помощи FMA, выполняя округление только в самом конце. Следовательно, она устойчива к катастрофическому сокращению, но ей нужно работать с округлённым значением $$. Оператор return «патчит» эту проблему, прибавляя погрешность, выделенную во второй строке. В статье Jeannenrod et al. показано, что результат верен до 1,5 ulps (мер единичной точности), что великолепно: FMA и простые операции с плавающей запятой точны до 0,5 ulps, поэтому алгоритм практически идеален.

Используем новый молоток

Когда начинаешь искать способы применения DifferenceOfProducts(), она оказывается на удивление полезной. Вычисляете дискриминанта квадратного уравнения? Вызывайте DifferenceOfProducts(b, b, 4 * a, c)⁴. Вычисляете определитель матрицы 2x2? Алгоритм решит эту задачу. В реализации следующей версии pbrt я нашёл ему около 80 применений. Из всех них самой любимой является функция векторного произведения. Она всегда была источником проблем, из-за которого приходилось поднимать руки вверх и использовать в реализации числа double, чтобы избежать катастрофического сокращения:

inline Vector3f Cross(const Vector3f &v1, const Vector3f &v2) {
    double v1x = v1.x, v1y = v1.y, v1z = v1.z;
    double v2x = v2.x, v2y = v2.y, v2z = v2.z;
    return Vector3f(v1y * v2z - v1z * v2y,
                    v1z * v2x - v1x * v2z,
                    v1x * v2y - v1y * v2x);
}

А теперь мы можем продолжать работать с float и использовать DifferenceOfProducts().

inline Vector3f Cross(const Vector3f &v1, const Vector3f &v2) {
    return Vector3f(DifferenceOfProducts(v1.y, v2.z, v1.z, v2.y),
                    DifferenceOfProducts(v1.z, v2.x, v1.x, v2.z),
                    DifferenceOfProducts(v1.x, v2.y, v1.y, v2.x));
}

Тот хитрый пример из начала поста на самом деле является частью векторного произведения. На определённом этапе коду pbrt требуется вычислить векторное произведение векторов $$ и $$. При вычислении с помощью float получался бы вектор $$. (Общее правило: если при вычислениях с плавающей запятой, где задействованы большие числа, вы получаете не такие большие целочисленные значения, то это почти верный признак того, что произошло катастрофическое сокращение.)

При двойной точности векторное произведение равно $$. Мы видим, что при float $$ выглядит нормально, $$ уже довольно плох, а $$ становится той катастрофой, которая стала мотивацией к поиску решения. А какие результаты мы получим с DifferenceOfProducts() и значениями float? $$. Значения $$ и $$ соответствуют двойной точности, а $$ слегка смещён — отсюда и взялась лишняя ulp.

А что со скоростью? DifferenceOfProducts() выполняет две FMA, а также умножение и сложение. Наивный алгоритм можно реализовать с одной FMA и одним умножением, что, казалось бы, должно занять в два раза меньше времени. Но на практике оказывается, что после получения значений из регистров это неважно: в синтетическом бенчмарке, проведённом на моём ноутбуке, DifferenceOfProducts() всего в 1,09 раза затратнее наивного алгоритма. Работа с двойной точностью была медленнее в 2,98 раза.

Как только вы узнаёте о возможности катастрофического сокращения, всевозможные невинно выглядящие выражения в коде начинаю казаться подозрительными. Похоже, что DifferenceOfProducts() является хорошим лекарством против большей их части. Она проста в использовании, и у нас нет особых причин не применять её.

Примечания

1. Катастрофическое сокращение не представляет проблемы при вычитании величин с разными знаками или при сложении величин с одинаковым знаком. Однако оно может стать проблемой при сложении значений с разными знаками. То есть суммы следует рассматривать с тем же подозрением, что и разности.
2. В качестве упражнения для читателя я советую написать функцию SumOfProducts(), обеспечивающую защиту от катастрофического сокращения. Если вы хотите усложнить задачу, то объясните, почему в DifferenceOfProducts(), dop + err == dop, ведь знаки a*b и c*d различаются.
3. Инструкция FMA доступна в GPU уже больше десятка лет, а в большинстве ЦП — не менее пяти лет. На ЦП может потребоваться добавить флаги компилятора для их непосредственной генерации при использовании std::fma(); в gcc и clang работает -march=native.
4. В формате чисел с плавающей запятой IEEE умножение на степени двойки выполняется точно, поэтому 4 * a не вызывает никакой ошибки округления, если только не происходит переполнения.