WikiSort.ru - Не сортированное

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

Алгоритм

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

• вычисляет половину значения числа и сохраняет для дальнейшего использования
• трактуя 32-битное с плавающей запятой как 32-битное целое:
  • производит сдвиг вправо на один бит
  • вычитает число из «магической» константы 5F3759DF₁₆
• (в трактовке результата как 32-битного с плавающей запятой на этом этапе получается первое приближение обратного квадратного корня исходного числа)
• вычисляет одну итерацию метода Ньютона для получения более точного приближения.

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

Структура 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	${\displaystyle =(1+2^{-2})\cdot 2^{-3}=0{,}15625}$
31				24				23				16				15				8				7				0

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

Если ${\displaystyle x\in (2^{k},2^{k+1})}$ (а значит, не точная степень двойки), то машинный порядок числа ${\displaystyle x}$ равняется a = 011.1111.1₂ + k, а порядок числа ${\displaystyle {\frac {1}{\sqrt {x}}}}$ : b = 011.1111.0₂ − [k/2], [x] — целая часть числа.^[2]

Попробуем вывести k через a и b: k = a − 011.1111.1₂ = 111.1110.0₂ − 2b. Отсюда b = (1011.1101.1₂ − a)/2 (при чётном k).

Теперь посмотрим, что будет при нечётном k. В такой ситуации b на единицу больший, и оба случая можно объединить в формулу: b = [(1011.1110.0₂ − a)/2].

Если добавить справа биты мантиссы (обозначенные звёздочками), то конструкции b*** = [(1011.1110.0***0 − a***)/2] и b*** = (101.1111.00*** − [a***/2]) отличаются одним младшим битом мантиссы. Что и даёт: верхний байт 5F, второй в диапазоне 00…3F. Оставшиеся биты мантиссы подбираются так, чтобы минимизировать погрешность.

Оценим магическую константу. Поскольку при ${\displaystyle x'=4x}$ мантисса результата не изменится, достаточно рассмотреть два диапазона: [1, 2) и [2, 4). Легко видеть, что на самом деле сдвиг-вычитание — грубая кусочно-линейная аппроксимация:

• Для ${\displaystyle x\in [1,2)}$ явная часть мантиссы x равняется ${\displaystyle x-1}$ , мантисса результата начинается с 1,1₂ = 1,5 (неявная единица, двоичная запятая, затем единица, пришедшая из порядка), порядок 011.1111.0₂ ≡ 2⁻¹ = 0,5, и ${\displaystyle {\hat {y}}_{1}(x)=0{,}5(1{,}5+(c-1)-0{,}5(x-1))=0{,}5+0{,}5c-0{,}25x}$ (c — мантисса магической константы; по вышеприведённым расчётам ${\displaystyle c<1{,}5}$ ).
• Для ${\displaystyle x\in [2,4)}$  — явная часть мантиссы x равняется ${\displaystyle 0{,}5x-1}$ , мантисса результата начинается с 1,0₂ = 1, порядок тот же, и ${\displaystyle {\hat {y}}_{2}=0{,}5(1+(c-1)-0{,}5(0{,}5x-1))=0{,}25+0{,}5c-0{,}125x}$ .
• Когда вычисленное на предыдущем этапе ${\displaystyle {\hat {y}}}$ окажется меньше 0,5, происходит заём из поля порядка, и ${\displaystyle {\hat {y}}_{3}(x)=0{,}5-{\frac {0{,}5-{\hat {y}}_{2}}{2}}={\frac {3}{8}}+0{,}25c-{\frac {1}{16}}x}$ .

Легко видеть, что ${\displaystyle {\hat {y}}_{1}(2)={\hat {y}}_{2}(2)}$ . ${\displaystyle {\hat {y}}_{3}}$ очевидно «сшивается» с ${\displaystyle {\hat {y}}_{2}}$ . Совпадающий с ${\displaystyle {\hat {y}}_{3}}$ результат получается, если расписать формулы для ${\displaystyle x\in [4,8)}$ . Потому наша кусочно-линейная аппроксимация непрерывна.

Наши расчёты грубые и прикидочные, так что обойдёмся условием: минимизировать относительную погрешность на концах первого отрезка, без учёта метода Ньютона. Поскольку ${\displaystyle {\frac {1}{\sqrt {1}}}=1}$ , ${\displaystyle {\frac {1}{\sqrt {2}}}\approx 0{,}7071}$ и линейная функция на отрезке падает на 0,25, то ${\displaystyle {\hat {y}}(0)=1-{\frac {(1-0{,}7071-0,25)\cdot 1}{1+0{,}7071}}\approx 0{,}9749\equiv 0{,}5+0{,}5c-0{,}25\cdot 1}$ . Таким образом, магическая константа — дробное число с порядком ровно 101.1111.0₂ ≡ 2⁶³ и мантиссой приблизительно 1,45.

Метод Ньютона даёт ${\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}}}$ . Поскольку функция ${\displaystyle f(y)}$ убывает и выпукла вниз, итерация метода Ньютона даёт лучшую точность, если начальное приближение будет меньше реального. Так что реальная константа должна быть меньше наших расчётных 1,45. После одного шага метода Ньютона результат получается довольно точный (+0 % −0,17 %), что для целей компьютерной графики вполне подходит. Такая погрешность сохраняется на всём диапазоне нормированных дробных чисел.

Неизвестно, откуда взялась константа 0x5F3759DF ≈ 1,432430·2⁶³. Перебором Крис Ломонт и Мэттью Робертсон выяснили, что наилучшая константа для float — 0x5F375A86 ≈ 1,432450·2⁶³, для double — 0x5FE6EB50C7B537A9. Правда, для double алгоритм бессмысленный (не даёт выигрыша в точности по сравнению с float). Более тонкие расчёты, учитывающие метод Ньютона, дают результат, совпадающий с эмпирической константой Ломонта^[3].

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

Реализация из Quake III: Arena

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

История

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

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

С выходом в свет в 1998 году набора инструкций 3DNow! в процессорах фирмы AMD появилась ассемблерная инструкция PFRSQRT для быстрого приближенного вычисления инверсного квадратного корня. Также можно получить константу для чисел двойной точности, но это бессмысленно, так как точность вычислений не увеличится, поэтому инструкции для вычисления быстрого обратного квадратного корня для чисел двойной точности не была добавлена.

Мотивация

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

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

${\displaystyle \|{\boldsymbol {v}}\|={\sqrt {{v_{1}}^{2}+{v_{2}}^{2}+{v_{3}}^{2}}}}$  — это евклидова норма вектора, аналог Евклидова расстояния между двумя точками в пространстве.

${\displaystyle {\boldsymbol {\hat {v}}}={\boldsymbol {v}}/\|{\boldsymbol {v}}\|}$  — это нормализованный единичный вектор. Вычисленное ${\displaystyle {v_{1}}^{2}+{v_{2}}^{2}+{v_{3}}^{2}}$ обозначим как ${\displaystyle {\boldsymbol {x}}}$ ,

${\displaystyle {\boldsymbol {\hat {v}}}={\boldsymbol {v}}/{\sqrt {x}}}$ , определяет соотношение единичного вектора и обратного квадратного корня от ${\displaystyle x}$ .

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

Примечания

Ссылки

• Magical Square Root Implementation In Quake III
• C. Lomont, Fast inverse square root, Technical Report, 2003.
• A Brief History of InvSqrt (недоступная ссылка) by Matthew Robertson
• 0x5f3759df, further investigations into accuracy and generalizability of the algorithm by Christian Plesner Hansen

В этой статье не хватает ссылок на источники информации. Информация должна быть проверяема, иначе она может быть поставлена под сомнение и удалена. Вы можете отредактировать эту статью, добавив ссылки на авторитетные источники. Эта отметка установлена 13 мая 2011 года.