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

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

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

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

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

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

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

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

Реализация из Quake III: Arena[3] считает, что float по длине равен long, и использует для преобразования указатели (может ошибочно сработать оптимизация «если изменился float, ни один long не менялся»; на GCC при компиляции в «выпуск» срабатывает предупреждение). А ещё она содержит нецензурное слово — очевидно, Джон Кармак, выкладывая игру в открытый доступ, не понял, что там делается.

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

Структура 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 (Точки — границы полубайтов, вертикальные линии — границы полей компьютерного дробного). Порядок 10 111 1102 равен 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.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Поверхность нормалей широко используются в расчетах освещения и затенения, требующих расчета норм для векторов. Здесь показано поле векторов нормали к поверхности.

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

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

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

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

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

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