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

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

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

Алгоритм

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

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

Реализация из 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;                       // evil floating point bit level hacking
	i  = 0x5f3759df - ( i >> 1 );               // what the fuck? 
	y  = * ( float * ) &i;
	y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//	y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

	return y;
}

Эта реализация считает, что float по длине равен long , и использует для преобразования указатели (может ошибочно сработать оптимизация «если изменился float , ни один long не менялся»; на 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}; // member 'f' set to value of 'number'.
	conv.i = 0x5f3759df - ( conv.i >> 1 );
	conv.f *= threehalfs - x2 * conv.f * conv.f;
	return conv.f;
}

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

#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 . До этой идеи добрались Грег Уолш, Клив Моулер и Гэри Таролли, работавшие тогда в . Грегу Уолшу и приписывается знаменитая константа 0x5F3759DF.

Таролли, перешедший в 3dfx , принёс алгоритм туда, где он и применялся для расчёта углов падения и отражения света в трёхмерной графике. Джим Блинн, специалист по 3D-графике, переизобрёл метод в 1997 году с более простой константой . Более сложный табличный метод, который считает до 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, и все настольные видеоускорители поддерживают ), зато важен как дань истории и в более ограниченных машинах .

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

Битовое представление 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

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

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

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

Например, двоичное представление 16-ричного целого числа 0x5F3759DF есть 0|101.1111.0|011.0111.0101.1001.1101.1111 ₂ (точки — границы полубайтов, вертикальные линии — границы полей компьютерного дробного). Порядок 101 1111 0 ₂ равен 190 ₁₀ , после вычитания смещения 127 ₁₀ получаем показатель степени 63 ₁₀ . Явная часть мантиссы 01 101 110 101 100 111 011 111 ₂ после добавления неявной ведущей единицы превращается в 1,011 011 101 011 001 110 111 11 ₂ = 1,432 430 148… ₁₀ . С учётом реальной точности компьютерных дробных 0x5F3759DF ↔ 1,4324301 ₁₀ ·2 ⁶³ .

Обозначим $m_{x}\in [0,1)$ явную часть мантиссы числа $x$ , $e_{x}\in \mathbb {Z}$ — несмещённый порядок, $L=2^{23}$ — разрядность мантиссы, $B=127$ — смещение порядка. Число $x\equiv 2^{e_{x}}(1+m_{x})$ , записанное в линейно-логарифмической разрядной сетке компьютерных дробных, можно приблизить логарифмической сеткой как $\log _{2}x\equiv e_{x}+\log _{2}(1+m_{x})\approx e_{x}+m_{x}+\sigma$ , где $\sigma$ — параметр, используемый для настройки точности приближения. Этот параметр варьируется от 0 (формула точна при $m_{x}=0$ и $1$ ) до 0,086 (точна в одной точке, $m_{x}=0{,}443$ )

Воспользовавшись этим приближением, целочисленное представление числа $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={\tfrac {1}{\sqrt {x}}}$ (соответственно $\log _{2}y=-{\tfrac {1}{2}}\log _{2}x$ )

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

2️⃣

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

I_{y}\approx {\tfrac {3}{2}}L(B-\sigma )-{\tfrac {1}{2}}I_{x}

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

Магическая константа ${\tfrac {3}{2}}L(B-\sigma )$ , с учётом границ $\sigma$ , в арифметике дробных чисел будет иметь несмещённый порядок $\left\lfloor {\tfrac {3L}{2L}}(B-\sigma )\right\rfloor -B=63$ и мантиссу $c=1+0{,}5-{\tfrac {3}{2}}\sigma \in (1{,}37;1{,}5$ ), а в двоичной записи — 0|101.1111.0|01 ₁ … (1 — неявная единица; 0,5 пришли из порядка; маленькая единица соответствует диапазону [1,375; 1,5) и потому крайне вероятна, но не гарантирована нашими прикидочными расчётами.)

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

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

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

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

Метод Ньютона даёт $f(y)={\frac {1}{y^{2}}}-x$ , $f'(y)=-{\frac {2}{y^{3}}}$ , и $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})$ . Функция $f(y)$ убывает и выпукла вниз, на таких функциях метод Ньютона подбирается к истинному значению слева — потому алгоритм всегда занижает ответ.

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

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

Исходный текст (C++)

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

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

Мотивация

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

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

Чтобы нормализовать вектор, надо разделить все три его компонента на длину. Или, что лучше, умножить их на величину, обратную длине: $(x',y',z')=(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 ↔ 1,4324301·2 ⁶³ . Перебором Крис Ломонт и Мэттью Робертсон выяснили , что наилучшая по предельной относительной погрешности константа для float — 0x5F375A86 ↔ 1,4324500·2 ⁶³ , для double — 0x5FE6EB50C7B537A9. Правда, для double алгоритм бессмысленный (не даёт выигрыша в точности по сравнению с float) . Константу Ломонта удалось получить и аналитически ( c = 1,432450084790142642179 ) , но расчёты довольно сложны .

Чех Ян Ка́длец двоичным поиском , а затем перебором в окрестности найденного улучшил алгоритм . Его метод даёт в 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}}}-xy^{3/2}=0$ . Он точнее одного шага метода Ньютона, но не дотягивает до двух и требует деление :

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),

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

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

Здесь стрелка ↔ означает объяснённую выше биекцию двоичного представления целого числа и двоичного представления числа с плавающей запятой в формате IEEE 754 .
↑ Если в поле порядка поставить 127, получится 0x3FB75A86. Библиотека GRISU2, полностью целочисленная и не зависящая от тонкостей сопроцессора, говорит, что 0x3FB75A86 ↔ 1,43245 — это кратчайшее десятичное число, преобразующееся в данный float . Однако всё-таки единица младшего разряда равняется 1,19·10 ⁻⁷ , и 0x3FB75A86 = 1,432450056 ≈ 1,4324501. Следующее дробное 0x3FB75A87 ↔ 1,4324502 без всяких тонкостей. Отсюда неинтуитивное округление 1,43245008 до 1,4324500.

Примечания

↑ (неопр.) . Дата обращения: 25 августа 2019. 6 февраля 2009 года.
↑
↑
↑ (неопр.) . Дата обращения: 1 февраля 2017. 13 января 2017 года.
(неопр.) . Дата обращения: 9 апреля 2023. 13 апреля 2023 года.
(неопр.) . Дата обращения: 25 августа 2019. 25 августа 2019 года.
(неопр.) . Дата обращения: 9 апреля 2023. 9 апреля 2023 года.
(неопр.) . Дата обращения: 17 августа 2022. 17 августа 2022 года. . от 9 апреля 2023 на Wayback Machine
↑ (неопр.) . Дата обращения: 17 августа 2022. 11 октября 2022 года.
(неопр.) . Дата обращения: 4 октября 2019. 10 апреля 2017 года.
(неопр.) . Дата обращения: 4 октября 2019. 16 октября 2019 года.
(неопр.) . Дата обращения: 6 октября 2019. 12 августа 2019 года.
(неопр.) . Дата обращения: 9 апреля 2023. 9 апреля 2023 года.
↑ (неопр.) . Дата обращения: 12 июня 2022. 17 апреля 2022 года.
↑ (неопр.) . Дата обращения: 4 июля 2022. 10 июля 2020 года.
↑ (неопр.) . Дата обращения: 8 апреля 2023. 8 апреля 2023 года.

Ссылки

C. Lomont, , Technical Report, 2003.
by Matthew Robertson
, further investigations into accuracy and generalizability of the algorithm by Christian Plesner Hansen