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

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

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

Алгоритм стал широко известен благодаря реализации в компьютерной игре Quake III Arena, в своё время налаживавшей красивое освещение даже полупрограммно, без использования Transform & Lighting.

Мотивация

Поле нормалей: а) для призмы (угловатый объект); б) для низкополигонального цилиндра (сглаженный объект)^[4]

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

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

Алгоритм

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

1. Трактуя 32-битное дробное число как целое, провести операцию y₀ = 5F3759DF₁₆ − (x >> 1), где >> — битовый сдвиг вправо. Результат снова трактуется как 32-битное дробное число.
2. Для уточнения можно провести одну итерацию метода Ньютона: y₁ = y₀(1,5 − 0,5xy₀²).

Реализация алгоритма в том виде, в котором он был впервые опубликован в исходном коде Quake III^[5]:

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;
}

Эта реализация ориентирована на 32-битную архитектуру процессоров x86, в которой размеры float и long одинаковы и равны 4 байтам. Корректная по меркам современного Си реализация, с учётом возможных оптимизаций и кроссплатформенности:

#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;
}

В C++ для устранения неопределённого поведения при преобразовании числа с плавающей запятой в целочисленное алгоритм можно реализовать следующим образом (используются возможности стандартов C++20 и C++23):

import std;

constexpr std::float32_t Q_rsqrt(std::float32_t number) noexcept
{
  const auto y = std::bit_cast<std::float32_t>(
    0x5f3759df - (std::bit_cast<std::uint32_t>(number) >> 1));
  return y * (1.5f32 - (number * 0.5f32 * y * y));
}

История

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

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

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

Приближённый обратный корень не раз реализовывали аппаратно: pfrsqrt (AMD 3DNow!, 1998), rsqrtss (Intel SSE, 1999), последний значительно точнее, 0,01 %^[12]. А функция Transform & Lighting (на ПК — nVidia GeForce 256, 1999) перенесла расчёт освещения на видеоадаптер. Так что программистам на ПК уже с 2000-х годов алгоритм не нужен. Тем не менее, он остаётся полезным в проектах для встраиваемых систем или систем с ограниченными ресурсами^[13].

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

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

Битовое представление 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	${\displaystyle =(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₂^[a].

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

В примере выше целочисленное представление равняется 0x3E20.0000, и оно раскладывается так: знаковое поле 0, смещённый порядок^[a] 011.1110.0₂=124, несмещённый 124−127=−3, мантисса (вместе с неявной единицей) 1,01₂=1,25, и дробное значение 1,25·2⁻³=0,15625.

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

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

${\displaystyle \log _{2}(1+m_{x})\approx m_{x}+\sigma }$. Приведены крайние случаи — σ = 0 и 0,086.

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

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

${\displaystyle I_{x}\equiv L(e_{x}+B+m_{x})\approx L\log _{2}x+L(B-\sigma )}$

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

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

${\displaystyle -{\tfrac {1}{2}}L\log _{2}x\approx I_{y}-L(B-\sigma )}$ 2️⃣

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

${\displaystyle -{\tfrac {1}{2}}[I_{x}-L(B-\sigma )]\approx I_{y}-L(B-\sigma )}$

${\displaystyle y\approx I^{-1}\left[{\tfrac {3}{2}}L(B-\sigma )-{\tfrac {1}{2}}I_{x}\right]\equiv I^{-1}\left(Q-{\tfrac {1}{2}}I_{x}\right)}$

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

Магическая константа Q

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

${\displaystyle e_{Q}=\left\lfloor {\tfrac {Q}{L}}-B\right\rfloor =\left\lfloor {\tfrac {3L}{2L}}(B-\sigma )-B\right\rfloor =\left\lfloor {\tfrac {B-3\sigma }{2}}\right\rfloor =63}$.

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

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

а в двоичной записи^[a] — 0|101.1111.0|01₁… (1 — неявная единица; 0,5 пришли из ${\displaystyle {\tfrac {3B}{2}}}$; маленькая единица соответствует диапазону [1,375; 1,5) и потому крайне вероятна, но не гарантирована нашими прикидочными расчётами.)

Первое (кусочно-линейное) приближение быстрого обратного квадратного корня (c = 1,43)

Через константу c можно вычислить, чему равняется первое кусочно-линейное приближение^[15] (в источнике используется не сама мантисса, а её явная часть ${\displaystyle m_{Q}\equiv t=c-1}$):

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

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

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

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

После одного шага метода Ньютона результат получается довольно точный (+0 % −0,18 %)^[1][2], что для целей компьютерной графики более чем подходит (1⁄256 ≈ 0,39 %). Такая погрешность сохраняется на всём диапазоне нормированных дробных чисел. Два шага дают точность в 5 цифр^[1], после четырёх достигается погрешность double.

Метод Ньютона может испортить монотонность. Однако как компьютерный перебор, так и аналитические выкладки говорят, что монотонность остаётся.

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

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

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

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

Чех Ян Ка́длец двоичным поиском, а затем перебором в окрестности найденного улучшил алгоритм^[16]. Его метод даёт в 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);
}

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

${\displaystyle y_{n+1}=y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}=y_{n}\left({\frac {3+xy_{n}^{2}}{1+3xy_{n}^{2}}}\right),}$

где ${\displaystyle xy_{n}^{2}}$ нужно рассчитать всего один раз и сохранить во временной переменой.

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

Комментарии

1. Здесь и далее точки — границы полубайтов, вертикальные линии — границы полей компьютерного дробного.
2. Здесь и далее ≡ означает «равно по определению».
3. Здесь стрелка ↔ означает объяснённую выше биекцию двоичного представления целого числа и двоичного представления числа с плавающей запятой в формате IEEE 754.
4. Неинтуитивное округление теоретического c=1,432450084… до 1,4324500 на самом деле двойное: сначала к ближайшему двоичному (0x3FB75A86 ≈ 1,432450056), а потом к самому круглому десятичному — единица младшего разряда равняется 1,19·10⁻⁷, и 1,19⁄2 > 0,56, так что 1,43245 — самое круглое, преобразующееся в 0x3FB75A86.