Быстрый обратный квадратный корень

При расчёте освещения OpenArena (свободный порт Quake III: Arena) вычисляет углы падения и отражения через быстрый обратный квадратный корень. Обратите внимание на кожух оружия — при очень низкой детализации (8 четырёхугольников) игра делает вид, что он криволинейный.

Бы́стрый обра́тный квадра́тный ко́рень (также быстрый InvSqrt() или 0x5F3759DF по используемой «магической» константе) — это приближённый алгоритм вычисления обратного квадратного корня $y=\frac{1}{\sqrt{x}}$ для положительных 32-битных чисел с плавающей запятой. Алгоритм использует целочисленные операции «вычесть» и «битовый сдвиг», а также дробные «вычесть» и «умножить» — без медленных операций «разделить» и «квадратный корень». Несмотря на «хакерство» на битовом уровне, приближение монотонно и непрерывно: близкие аргументы дают близкий результат. Точности (менее 0,2 % в меньшую сторону и никогда — в большую) не хватает для настоящих численных расчётов и даже для нормирования матриц поворота в трёхмерной графике, однако вполне достаточно для маловажных эффектов вроде освещения и теней.

Алгоритм

Алгоритм принимает 32-битное число с плавающей запятой (одинарной точности в формате IEEE 754) в качестве исходных данных и производит над ним следующие операции:

1. Трактуя 32-битное дробное число как целое, провести операцию Шаблон:Nobr, где >> — битовый сдвиг вправо. Результат снова трактуется как 32-битное дробное число.
2. Для уточнения можно провести одну итерацию метода Ньютона: Шаблон:Nobr.

Реализация из Quake III:

float Q_rsqrt( float number )
{
	long i;
	float x2, y;
	const float threehalfs = 1.5F;

	x2 = number * 0.5F;
	y  = number;
	i  = * ( long * ) &y;                       // страшное хакерство на битовом уровне
	i  = 0x5f3759df - ( i >> 1 );               // что за чёрт?
	y  = * ( float * ) &i;
	y  = y * ( threehalfs - ( x2 * y * y ) );   // 1-я итерация
//	y  = y * ( threehalfs - ( x2 * y * y ) );   // 2-я итерация, можно убрать

	return y;
}

Эта реализация считает, что Шаблон:Cpp по длине равен Шаблон:Cpp, и использует для преобразования указатели (может ошибочно сработать оптимизация «если изменился Шаблон:Cpp, ни один Шаблон:Cpp не менялся»; на GCC при компиляции в «выпуск» срабатывает предупреждение). По комментариям видно, что Джон Кармак, выкладывая игру в открытый доступ, не понял, что там делается.

Корректная по меркам современного Си реализация, с учётом возможных оптимизаций и кроссплатформенности:

#include <stdint.h>

float Q_rsqrt( float number )
{	
	const float x2 = number * 0.5F;
	const float threehalfs = 1.5F;

	union {
		float f;
		uint32_t i;
	} conv = {number}; // такая инициализация присвоит поле «f»
	conv.i = 0x5f3759df - ( conv.i >> 1 );
	conv.f *= threehalfs - x2 * conv.f * conv.f;
	return conv.f;
}

На Си++20 можно использовать новую функцию Шаблон:Cpp.

#include <bit>
#include <limits>
#include <cstdint>

constexpr float Q_rsqrt(float number) noexcept
{
	static_assert(std::numeric_limits<float>::is_iec559); // Проверка совместимости целевой машины

	float const y = std::bit_cast<float>(
		0x5f3759df - (std::bit_cast<std::uint32_t>(number) >> 1));
	return y * (1.5f - (number * 0.5f * y * y));
}

GCC и Clang (-std=c++20 -mx32 -O3) дают одинаковый машинный код для всех трёх вариантов и близкий — друг относительно друга. У MSVC (/std:c++20 /O2) третья функция незначительно отличается от первых двух.

История

Саму идею приближения дробного числа целым для вычисления $\sqrt{x}$ придумали Уильям Кэхэн и К. Ын в 1986. До этой идеи добрались Грег Уолш, Клив Моулер и Гэри Таролли, работавшие тогда в Ardent Computer. Грегу Уолшу и приписывается знаменитая константа 0x5F3759DF.

Таролли, перешедший в 3dfx, принёс алгоритм туда, где он и применялся для расчёта углов падения и отражения света в трёхмерной графике. Джим Блинн, специалист по 3D-графике, переизобрёл метод в 1997 году с более простой константой 1,5. Более сложный табличный метод, который считает до 4 знаков (0,01 %), найден при дизассемблировании игры Interstate ’76 (1997).

Однако данный метод не появлялся на общедоступных форумах, таких как Usenet, до 2002—2003-х годов. Метод обнаружили в Quake III: Arena, опубликованном в 2005, и приписали авторство Джону Кармаку. Тот предположил, что его в id Software принёс Майкл Абраш, специалист по графике, или Терье Матисен, специалист по ассемблеру; другие ссылались на Брайана Хука, выходца из 3dfx. Изучение вопроса показало, что код имел более глубокие корни как в аппаратной, так и в программной сферах компьютерной графики. Исправления и изменения производились как Silicon Graphics, так и 3dfx Interactive, при этом самая ранняя известная версия написана Гэри Таролли для SGI Indigo.

С выходом в свет в 1998 году набора инструкций 3DNow! в процессорах фирмы AMD появилась ассемблерная инструкция PFRSQRT для быстрого приближенного вычисления обратного квадратного корня. Версия для double бессмысленна — точность вычислений не увеличится — потому её не добавили. В 2000 году в SSE2 добавили функцию RSQRTSS более точную, чем данный алгоритм (0,04 % против 0,2 %). Потому данный метод больше не имеет смысла в современных компьютерах (все x64-процессоры поддерживают SSE2, и все настольные видеоускорители поддерживают Transform & Lighting), зато важен как дань истории и в более ограниченных машинах.

Анализ и погрешность

Преобразование «дробное ↔ целое»

Битовое представление 4-байтового дробного числа в формате IEEE 754 выглядит так:

Знак
	Порядок								Мантисса
0	0	1	1	1	1	1	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	$=(1+2^{-2})\cdot 2^{-3}=0,15625$
31				24				23				16				15				8				7				0

Имеем дело только с положительными числами (знаковый бит равен нулю), не денормализованными, не ∞ и не NaN. Такие числа в стандартном виде записываются как 1,mmmm₂·2^e. Часть 1,mmmm называется мантиссой, e — порядком. Головную единицу не хранят (неявная единица), так что величину 0,mmmm назовём явной частью мантиссы. Кроме того, у машинных дробных чисел смещённый порядок: 2⁰ записывается как 011.1111.1₂Шаблон:Efn.

На положительных числах биекция «дробное ↔ целое» (ниже обозначенная как $I_x$) непрерывна как кусочно-линейная функция и монотонна. Отсюда сразу же можно заявить, что быстрый обратный корень, как комбинация непрерывных функций, непрерывен. А первая его часть — сдвиг-вычитание — к тому же монотонна и кусочно-линейна. Биекция сложна, но почти «бесплатна»: в зависимости от архитектуры процессора и соглашений вызова, нужно или ничего не делать, или переместить число из дробного регистра в целочисленный.

В примере выше целочисленное представление равняется 0x3E20.0000, и оно раскладывается так: знаковое поле 0, смещённый порядокШаблон:Efn 011.1110.0₂=124, несмещённый 124−127=−3, мантисса (вместе с неявной единицей) 1,01₂=1,25, и дробное значение Шаблон:Nobr.

Обозначим $m_x\in [0,1)$ явную часть мантиссы числа $x$, $e_x\in \mathbb{Z}$ — несмещённый порядок, $L=2^{23}$ — разрядность мантиссы, $B=127$ — смещение порядка. ЧислоШаблон:Efn $x\equiv 2^{e_x}(1+m_x)$ будет иметь целочисленное представление $I_x\equiv L(e_x+B+m_x)$. Можно выписать и обратное преобразование: $e_x=\lfloor \frac{I_x}{L}-B\rfloor$, $m_x=\left\{\frac{I_x}{L}\right\}$, ведь первое целое, а второе от 0 до 1.

Первое приближение

${\log }_2(1+m_x)\approx m_x+\sigma$. Приведены крайние случаи — Шаблон:Nobr и 0,086.

Поскольку ${\log }_{2}1=0$ и ${\log }_{2}2=1$, нелинейную функцию «логарифм» можно приблизить линейной ${\log }_2(1+m_x)\approx m_x+\sigma$, где $\sigma$ — параметр, используемый для настройки точности приближения. Этот параметр варьируется от 0 (формула точна при $m_x=0$ и $1$) до 0,086 (точна в одной точке, $m_x=0,443$).

Аргумент $x$, записанный в линейно-логарифмической разрядной сетке компьютерных дробных, можно приблизить логарифмической сеткой как ${\log }_{2}x\equiv e_x+{\log }_2(1+m_x)\approx e_x+m_x+\sigma$. Перегруппируем $e_x+m_x\approx {\log }_{2}x-\sigma$, тогда целочисленное представление числа $x$ можно приблизить как

$I_x\equiv L(e_x+B+m_x)\approx L{\log }_{2}x+L(B-\sigma )$

Соответственно, $L{\log }_{2}x\approx I_x-L(B-\sigma )$. 1️⃣

Проделаем это же для $y=\frac{1}{\sqrt{x}}$ (соответственно ${\log }_{2}y=-\frac{1}{2}{\log }_{2}x$)

$-\frac{1}{2}L{\log }_{2}x\approx I_y-L(B-\sigma )$ 2️⃣

Соединив 1️⃣ и 2️⃣, получаем

$-\frac{1}{2}[I_x-L(B-\sigma )]\approx I_y-L(B-\sigma )$

$y\approx I^{-1}[\frac{3}{2}L(B-\sigma )-\frac{1}{2}{I}_x]\equiv I^{-1}(Q-\frac{1}{2}{I}_x)$

Это и есть формула первого приближения.

Магическая константа Шаблон:Math

Магическая константа $Q\equiv \frac{3}{2}L(B-\sigma )$ находится в пространстве компьютерных целых, но её дробное представление $I^{-1}(Q)$ также важно для исследователей. Его несмещённый порядок

$e_Q=\lfloor \frac{Q}{L}-B\rfloor =\lfloor \frac{3L}{2L}(B-\sigma )-B\rfloor =\lfloor \frac{B-3\sigma }{2}\rfloor$.

Поскольку $\sigma <\frac{1}{6}$, то независимо от чётности числа Шаблон:Math

$e_Q=\lfloor \frac{B-1}{2}\rfloor =63$

Смещение порядка Шаблон:Math нечётное, и полная мантисса числа $I^{-1}(Q)$ равняется

$c\equiv 1+m_Q=1+\left\{\frac{Q}{L}\right\}=1+0,5-\frac{3}{2}\sigma \in (1,37;1,5$),

а в двоичной записиШаблон:Efn — 0|101.1111.0|01₁… (1 — неявная единица; 0,5 пришли из порядка; маленькая единица соответствует диапазону [1,375; 1,5) и потому крайне вероятна, но не гарантирована нашими прикидочными расчётами.)

Первое (кусочно-линейное) приближение быстрого обратного квадратного корня Шаблон:Nobr

Через константу Шаблон:Math можно вычислить, чему равняется первое кусочно-линейное приближение (в источнике используется не сама мантисса, а её явная часть $m_Q\equiv t=c-1$):

• Для $x\in [0,5;c-0,5)$: $y_{01}=-x+t+\frac{3}{2}=-x+c+\frac{1}{2}$;
• Для $x\in [c-0,5;1)$: $y_{02}=-\frac{1}{2}x+\frac{1}{2}t+\frac{5}{4}=-\frac{1}{2}x+\frac{1}{2}c+\frac{3}{4}$;
• Для $x\in [1;2)$: $y_{03}=-\frac{1}{4}x+\frac{1}{2}t+1=-\frac{1}{4}x+\frac{1}{2}c+\frac{1}{2}$.

На бо́льших или меньших $x$ результат пропорционально меняется: при учетверении $x$ результат уменьшается ровно вдвое.

Метод Ньютона

Метод Ньютона даёт $f(y)=\frac{1}{y^2}-x$, $f^{\prime }(y)=-\frac{2}{y^3}$, и $y_{n+1}=y_n-\frac{f(y_n)}{f^{\prime }(y_n)}=\frac{y_n(3-x{y}_n^2)}{2}=y_n(1,5-0,5x{y}_n^2)$. Функция $f(y)$ убывает и выпукла вниз, на таких функциях метод Ньютона подбирается к истинному значению слева — потому алгоритм всегда занижает ответ.

После одного шага метода Ньютона результат получается довольно точный (Шаблон:Nobr), что для целей компьютерной графики более чем подходит (Шаблон:Nobr). Такая погрешность сохраняется на всём диапазоне нормированных дробных чисел. Два шага дают точность в 5 цифр, после четырёх достигается погрешность double.

Метод Ньютона не гарантирует монотонности, но компьютерный перебор показывает, что монотонность всё-таки есть. Шаблон:Начало скрытого блока

#include <iostream>

union FloatInt {
    float asFloat;
    int32_t asInt;
};

int floatToInt(float x)
{
    FloatInt r;
    r.asFloat = x;
    return r.asInt;
}

float intToFloat(int x)
{
    FloatInt r;
    r.asInt = x;
    return r.asFloat;
}

float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;

    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;                       // evil floating point bit level hacking
    i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
    y  = * ( float * ) &i;                      // i don't know, what the fuck!
    y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration

    return y;
}

int main()
{
    int iStart = floatToInt(1.0);
    int iEnd = floatToInt(4.0);
    std::cout << "Numbers to go: " << iEnd - iStart << std::endl;
    int nProblems = 0;
    float oldResult = std::numeric_limits<float>::infinity();

    for (int i = iStart; i <= iEnd; ++i) {
        float x = intToFloat(i);
        float result = Q_rsqrt(x);
        if (result > oldResult) {
            std::cout << "Found a problem on " << x << std::endl;
            ++nProblems;
        }
    }
    std::cout << "Total problems: " << nProblems << std::endl;

    return 0;
}

Существуют аналогичные алгоритмы для других степеней, например, квадратного или кубического корня.

Мотивация

Поле нормалей: Шаблон:Nobr призмы (угловатый объект); Шаблон:Nobr низкополигонального цилиндра (шесть рёбер сглажены)

«Прямое» наложение освещения на трёхмерную модель, даже высокополигональную, даже с учётом закона Ламберта и других сложных формул отражения и рассеивания, сразу же выдаст полигональный вид — зритель увидит разницу в освещении по рёбрам многогранника. Иногда так и нужно — если предмет действительно угловатый. А для криволинейных предметов поступают так: в трёхмерной программе указывают, острое ребро или сглаженное. В зависимости от этого ещё при экспорте модели в игру по углам треугольников вычисляют нормаль единичной длины к криволинейной поверхности. При анимации и поворотах игра преобразует эти нормали вместе с остальными трёхмерными данными; при наложении освещения — интерполирует по всему треугольнику и нормализует (доводит до единичной длины).

Чтобы нормализовать вектор, надо разделить все три его компонента на длину. Или, что лучше, умножить их на величину, обратную длине: $(x^{\prime },y^{\prime },z^{\prime })=(x,y,z)\frac{1}{\sqrt{x^2+y^2+z^2}}$. За секунду должны вычисляться миллионы этих корней. До того как было создано специальное аппаратное обеспечение для обработки трансформаций и освещения, программное обеспечение вычислений могло быть медленным. В частности, в начале 1990-х, когда код был разработан, большинство вычислений с плавающей запятой отставало по производительности от операций с целыми числами.

Quake III Arena использует алгоритм быстрого обратного квадратного корня для ускорения обработки графики центральным процессором, но с тех пор алгоритм уже был реализован в некоторых специализированных аппаратных вершинных шейдерах, используя специальные программируемые матрицы (FPGA).

Даже на компьютерах 2010-х годов, в зависимости от загрузки дробного сопроцессора, скорость может быть втрое-вчетверо выше, чем с использованием стандартных функций.

Дальнейшие улучшения

При желании можно перебалансировать погрешность, умножив коэффициенты 1,5 и 0,5 на 1,0009, чтобы метод давал симметрично ±0,09 % — так поступили в игре Interstate ’76, которая также делает итерацию метода Ньютона.

Константа Уолша 0x5F3759DF ↔Шаблон:Efn 1,4324301·2⁶³ оказалась очень хорошей. Крис Ломонт и Мэттью Робертсон незначительно уменьшили предельную относительную погрешность, отыскав перебором константуШаблон:Efn 0x5F375A86 ↔ 1,4324500·2⁶³, а для double — 0x5FE6EB50C7B537A9. Правда, для double алгоритм бессмысленный (не даёт выигрыша в точности по сравнению с float). Константу Ломонта удалось получить и аналитически (Шаблон:Nobr)Шаблон:Efn, но расчёты довольно сложны. Крайний случай — константа Блинна 1,5 — даёт без перебалансировок и улучшений около −0,6 %.

Чех Ян Ка́длец двоичным поиском, а затем перебором в окрестности найденного улучшил алгоритм. Его метод даёт в 1,3 раза меньшую симметричную погрешность — не ±0,09 %, а ±0,065.

float inv_sqrt(float x)
{ union { float f; uint32 u; } y = {x};
  y.u = 0x5F1FFFF9ul - (y.u >> 1);
  return 0.703952253f * y.f * (2.38924456f - x * y.f * y.f);
}

Вместо метода Ньютона можно использовать метод Галлея, в данной задаче эквивалентный методу Ньютона для уравнения $f(y)=\frac{1}{y^{1/2}}-x{y}^{3/2}=0$. Он точнее одного шага метода Ньютона, но не дотягивает до двух и требует деление:

$y_{n+1}=y_n-\frac{f(y_n)}{f^{\prime }(y_n)}=y_n\left(\frac{3+x{y}_n^2}{1+3x{y}_n^2}\right),$

где $x{y}_n^2$ нужно рассчитать всего один раз и сохранить во временной переменой.

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

Ссылки

• C. Lomont, Fast inverse square root, Technical Report, 2003.
• A Brief History of InvSqrt by Matthew Robertson
• 0x5f3759df, further investigations into accuracy and generalizability of the algorithm by Christian Plesner Hansen