0x5f375a86

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

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


Алгоритм[ | ]

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

  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-байтового дробного такова:

Знак
Порядок Мантисса
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.

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

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

Целочисленное представление числа можно приблизить как

Эта операция сложна, но почти «бесплатна»: в зависимости от архитектуры процессора и соглашений вызова, нужно или ничего не делать, или переместить число из дробного регистра в целочисленный.

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

Для будет , и

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

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

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

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

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

Метод Ньютона даёт[5] , , и .

Неизвестно, откуда взялась константа 0x5F3759DF = 1,4324301·263. Перебором Крис Ломонт и Мэттью Робертсон выяснили[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 принёс Майкл Абраш, специалист по графике, или Терье Матисен, специалист по ассемблеру[6]. Изучение вопроса показало, что имел более глубокие корни как в аппаратной, так и в программной сферах компьютерной графики. Исправления и изменения производились как Silicon Graphics, так и 3dfx Interactive, при этом самая ранняя известная версия написана Гэри Таролли для SGI Indigo. Возможно, алгоритм придумали Грег Уолш и Клив Моулер, коллеги Гэри по Ardent Computer[7].

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

Мотивация[ | ]

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

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

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

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

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

Примечания[ | ]

Ссылки[ | ]