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

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

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

Алгоритм[править | править код]

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

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

Оригинальный алгоритм[3]:

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 при компиляции в «выпуск» срабатывает предупреждение). По комментариям видно, что Джон Кармак, выкладывая игру в открытый доступ, не понял, что там делается.

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

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

История[править | править код]

Алгоритм был, вероятно, разработан в Silicon Graphics в 1990-х, наиболее известная реализация появилась в 1999 году в исходном коде компьютерной игры Quake III Arena, но данный метод не появлялся на общедоступных форумах, таких как Usenet, до 2002—2003-х годов. Алгоритм генерирует достаточно точные результаты, используя уникальное первое приближение метода Ньютона. В то время основным преимуществом алгоритма был отказ от дорогих вычислительных операций с плавающей запятой в пользу целочисленных операций. Обратные квадратные корни используются для расчета углов падения и отражения для освещения и затенения в компьютерной графике.

Алгоритм изначально приписывался Джону Кармаку, но тот предположил, что его в id Software принёс Майкл Абраш, специалист по графике, или Терье Матисен, специалист по ассемблеру[4]. Изучение вопроса показало, что код имел более глубокие корни как в аппаратной, так и в программной сферах компьютерной графики. Исправления и изменения производились как Silicon Graphics, так и 3dfx Interactive, при этом самая ранняя известная версия написана Гэри Таролли для SGI Indigo. Возможно, алгоритм придумали Грег Уолш и Клив Моулер, коллеги Гэри по Ardent Computer[5].

Джим Блинн, специалист по 3D-графике, предложил похожий табличный метод вычисления обратного квадратного корня[6], который считает до 4 знаков (0,01 %) и найден при дизассемблировании игры Interstate ’76 (1997)[7].

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

Анализ и погрешность[править | править код]

Битовое представление 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
31 24 23 16 15 8 7 0
. Приведены крайние случаи — σ = 0 и 0,086.

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

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

Например, двоичное представление 16-ричного целого числа 0x5F3759DF есть 0|101.1111.0|011.0111.0101.1001.1101.11112 (точки — границы полубайтов, вертикальные линии — границы полей компьютерного дробного). Порядок 101 1111 02 равен 19010, после вычитания смещения 12710 получаем показатель степени 6310. Явная часть мантиссы 01 101 110 101 100 111 011 1112 после добавления неявной ведущей единицы превращается в 1,011 011 101 011 001 110 111 112 = 1,432 430 148…10. С учётом реальной точности компьютерных дробных 0x5F3759DF ↔ 1,432430110·263.

Обозначим явную часть мантиссы числа ,  — несмещённый порядок,  — разрядность мантиссы,  — смещение порядка. Число , записанное в линейно-логарифмической разрядной сетке компьютерных дробных, можно[10][3] приблизить логарифмической сеткой как , где  — параметр, используемый для настройки точности приближения. Этот параметр варьируется от 0 (формула точна при и ) до 0,086 (точна в одной точке, )

Воспользовавшись этим приближением, целочисленное представление числа можно приблизить как

Соответственно, .

Проделаем это же[3] для (соответственно ), и получим

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

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

Можно вычислить, чему равняется первое кусочно-линейное приближение[11] (в источнике используется не сама мантисса, а её явная часть ):

  • Для : ;
  • Для : ;
  • Для : .

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

Метод Ньютона даёт[11] , , и . Функция убывает и выпукла вниз, на таких функциях метод Ньютона подбирается к истинному значению слева — потому алгоритм всегда занижает ответ.

Неизвестно, откуда взялась константа 0x5F3759DF ↔[a] 1,4324301·263. Перебором Крис Ломонт и Мэттью Робертсон выяснили[1][2], что наилучшая по предельной относительной погрешности константа[b] для float — 0x5F375A86 ↔ 1,4324500·263, для double — 0x5FE6EB50C7B537A9. Правда, для double алгоритм бессмысленный (не даёт выигрыша в точности по сравнению с float)[2]. Константу Ломонта удалось получить и аналитически (c = 1,432450084790142642179)[b], но расчёты довольно сложны[11][2].

После одного шага метода Ньютона результат получается довольно точный (+0 % −0,18 %)[1][2], что для целей компьютерной графики более чем подходит (1256 ≈ 0,39 %). Такая погрешность сохраняется на всём диапазоне нормированных дробных чисел. Два шага дают точность в 5 цифр[1], после четырёх достигается погрешность double. При желании можно перебалансировать погрешность, умножив коэффициенты 1,5 и 0,5 на 1,0009, чтобы метод давал симметрично ±0,09 % — так поступили в игре Interstate ’76 и методе Блинна, которые также делают итерацию метода Ньютона[7].

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

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

Мотивация[править | править код]

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

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

Чтобы нормализовать вектор, надо разделить все три его компонента на длину. Или, что лучше, умножить их на величину, обратную длине: . За секунду должны вычисляться миллионы этих корней. До того как было создано специальное аппаратное обеспечение для обработки трансформаций и освещения, программное обеспечение вычислений могло быть медленным. В частности, в начале 1990-х, когда код был разработан, большинство вычислений с плавающей запятой отставало по производительности от операций с целыми числами.

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

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

Комментарии[править | править код]

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

Примечания[править | править код]

  1. 1 2 3 4 Архивированная копия. Дата обращения: 25 августа 2019. Архивировано 6 февраля 2009 года.
  2. 1 2 3 4 5 6 https://web.archive.org/web/20140202234227/http://shelfflag.com/rsqrt.pdf
  3. 1 2 3 4 Hummus and Magnets. Дата обращения: 1 февраля 2017. Архивировано 13 января 2017 года.
  4. Beyond3D — Origin of Quake3’s Fast InvSqrt(). Дата обращения: 4 октября 2019. Архивировано 10 апреля 2017 года.
  5. Beyond3D — Origin of Quake3’s Fast InvSqrt() — Part Two. Дата обращения: 25 августа 2019. Архивировано 25 августа 2019 года.
  6. Floating-point tricks | IEEE Journals & Magazine | IEEE Xplore
  7. 1 2 Fast reciprocal square root… in 1997?! — Shane Peelar’s Blog
  8. PFRSQRT — Вычислить приблизительное значение обратной величины квадратного корня от короткого вещественного значения — Club155.ru. Дата обращения: 4 октября 2019. Архивировано 16 октября 2019 года.
  9. RSQRTSS — Compute Reciprocal of Square Root of Scalar Single-Precision Floating-Point Value. Дата обращения: 6 октября 2019. Архивировано 12 августа 2019 года.
  10. https://web.archive.org/web/20150511044204/http://www.daxia.com/bibis/upload/406Fast_Inverse_Square_Root.pdf
  11. 1 2 3 4 Швидке обчислення оберненого квадратного кореня з використанням магічної константи — аналітичний підхід. Дата обращения: 12 июня 2022. Архивировано 17 апреля 2022 года.
  12. 1 2 3 Это норма: что такое карты нормалей и как они работают / Хабр

Ссылки[править | править код]