Швидкий обернений квадратний корінь

Для обчислення освітлення і віддзеркалення (показано у шутері від першої особи *OpenArena*) використовуються швидкий обернений квадратний корінь для обчислення кутів падіння і відбиття.

Швидкий обернений квадратний корінь (іноді згадуваний як Fast InvSqrt() або за шістнадцятковою сталою 0x5f3759df) — це метод обчислення $f(x)={\frac {1}{\sqrt {x}}}$ , оберненого квадратного кореня для 32-бітного числа у форматі чисел з рухомою комою IEEE 754. Алгоритм ймовірно розробили у Silicon Graphics на початку 1990-х, і реалізація з'явилась 1999 року в сирцевому коді Quake III Arena, але метод не з'являвся на публічних форумах як-от Usenet до 2002 чи 2003.^[1] (Існує обговорення на китайському форумі розробників CSDN у 2000.^[2]) На той час, основна перевага алгоритму полягала у використанні замість обчислювально дорогих операцій над числами з рухомою комою операцій над цілими числами. Обернений квадратний корінь використовують для обчислення кутів падіння і відбивання для освітлення і шейдинга в комп'ютерній графіці.

Алгоритм приймає 32-бітне число з рухомою комою і зберігає його половинне значення для подальшого використання. Тоді, трактуючи числа з рухомою комою як цілі, виконується логічний зсув вправо на один біт і результат віднімається від магічного числа 0x5f3759df. Це буде першим наближенням до оберненого квадратного кореня вхідного числа. Знов трактуючи біти як число з рухомою комою проводиться одна ітерація методу Ньютона, щоб результат був точнішим. Так обчислення наближеного значення оберненого квадратного кореня для числа з рухомою комою відбувається приблизно вчетверо швидше ніж із використанням ділення чисел з рухомою комою.

Огляд коду

Наступний код є реалізацією оберненого квадратного кореня з 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;                       // злий хак із рухомою комою на бітовому рівні
	i  = 0x5f3759df - ( i >> 1 );               // що за чортівня? 
	y  = * ( float * ) &i;
	y  = y * ( threehalfs - ( x2 * y * y ) );   // 1-ша ітерація
//	y  = y * ( threehalfs - ( x2 * y * y ) );   // 2-га ітерація, це можна видалити

	return y;
}

Для визначення оберненого квадратного кореня визначається наближення для $x^{-1/2}$ , тоді за допомогою чисельного методу це наближення переглядається, щоб отримати прийнятну похибку у кінцевому результаті. Звичайні програмні методи на початку 1990-х отримували перше наближення із таблиці пошуку.^[3] Цей шматок коду виявився швидшим ніж використання таблиці пошуку і приблизно в чотири рази швидший ніж звичайне ділення чисел з рухомою комою.^[4] Хоча деяка втрата точності і відбувалася, але її перекривало значне покращення швидкодії.^[5] Алгоритм був розроблений для специфікації IEEE 754-1985(інші мови) 32 бітних чисел з рухомою комою, але подальші дослідження Кріса Ломонта і Чарльза Макінері показали, що його можна реалізувати і для інших специфікацій.

Переваги у швидкості пропоновані швидким оберненим квадратним коренем з'явились завдяки трактуванню довгого слова^{[note 1]}, що містить число з рухомою комою як цілого і віднімання його від специфічної сталої, 0x5f3759df. Ціль цієї сталої не одразу очевидна для читача коду, отже, як і багато інших сталих знайдених у коді, її називають магічним числом.^[1]^[6]^[7]^[8] Це цілочисельне віднімання і бітовий зсув дають довге слово, яке знов трактується як число з рухомою комою і є грубим наближенням оберненого квадратного кореня вхідного числа. Одна ітерація методу Ньютона виконується для отримання більшої точності, і код завершується. Алгоритм генерує прийнятно точні результати використовуючи унікальне перше наближення для методу Ньютона; однак, він набагато повільніший ніж використання SSE інструкції rsqrtss на x86 процесорах також випущеної у 1999.^[9]

Робочий приклад

Як приклад, розглянемо число $x = 0.15625$ , для якого ми хочемо обчислити $1/ \sqrt x \approx 2.52982$ . Перші кроки алгоритму проілюстровані нижче:

0011_1110_0010_0000_0000_0000_0000_0000  Вигляд x та i на бітовому рівні
0001_1111_0001_0000_0000_0000_0000_0000  Зсув вправо на одну позицію: (i >> 1)
0101_1111_0011_0111_0101_1001_1101_1111  Магічне число 0x5f3759df
0100_0000_0010_0111_0101_1001_1101_1111  Результат 0x5f3759df — (i >> 1)

Використовуючи IEEE 32 бітове представлення:

0_01111100_01000000000000000000000  1.25 * 2^-3
0_00111110_00100000000000000000000  1.125 * 2^-65
0_10111110_01101110101100111011111  1.432430... * 2^+63
0_10000000_01001110101100111011111  1.307430... * 2^+1

Інтерпретування останнього бітового представлення як числа з рухомою комою дає наближення $y = 2.61486$ , яке має похибку близько 3.4%. Після однієї ітерації метода Ньютона, кінцевим результатом є $y = 2.52549$ , і помилка становить лише 0.17%.

Перебіг алгоритму

Алгоритм обчислює $1/ \sqrt x$ виконуючи такі кроки:

Інтерпретує аргумент $x$ як ціле, як спосіб приблизного обчислення $log 2 (x)$
Використовує це наближення для обчислення наближення $log 2 (1/ \sqrt x)$
Знов інтерпретує як число з рухомою комою, як спосіб для обчислення наближення $1/ \sqrt x$
Уточнює наближення використовуючи метод Ньютона.

Представлення чисел з рухомою комою

Докладніше: Число одинарної точності

Оскільки алгоритм сильно покладається на представлення чисел одинарної точності з рухомою комою на бітовому рівні, короткий огляд цього представлення наведений тут. Для того, щоб закодувати ненульове дійсне число $x$ як число із рухомою комою одинарної точності, перший крок полягає в записуванні $x$ як нормалізованого двійкового числа:

{\begin{aligned}x&=\pm 1.b_{1}b_{2}b_{3}\ldots \times 2^{e_{x}}\\&=\pm 2^{e_{x}}(1+m_{x})\end{aligned}}

де показник $e x$ є цілим, $m x \in [0, 1)$ , і $1.b 1 b 2 b 3 ...$ це двійкове представлення мантиси $(1 + m x)$ . Варто зазначити, що оскільки єдиний біт перед комою у мантисі завжди 1, то немає потреби його зберігати. З цієї форми маємо три беззнакові цілі числа:

$S x$ , знаковий біт, це 0 якщо $x > 0$ , і 1 якщо $x < 0$ (1 біт)
$E x = e x + B$ — це зміщена експонента, де $B = 127$ — зсув^{[note 2]} (8 бітів)
$M x = m x \times L$ , де $L = 2 23$ ^{[note 3]} (23 bits)

Ці поля пакуються зліва направо у 32 бітовий контейнер.

Як приклад розглянемо число $x = 0.15625 = 0.00101 2$ . Нормалізація $x$ дає:

x=+2^{-3}(1+0.25)

і отже, три беззнакові цілочисельні поля такі:

$S = 0$
$E = -3 + 127 = 124 = 01111100 2$
$M = 0.25 \times 2 23 = 2097152 = 01000000000000000000000 2$

ці поля пакуються як показано нижче:

Інтерпретування цілим як приблизний логарифм

Якби комусь довелось порахувати $1/ \sqrt x$ без комп'ютера чи калькулятора, то йому б стала в пригоді таблиця логарифмів разом із тотожністю $log b (1/ \sqrt x) = -½ log b (x)$ , яка дійсна для кожної основи $b$ . Швидкий обернений квадратний корінь базується на цій тотожності і на факті, що інтерпретація float32 у ціле число дає грубе наближення цього логарифма. Ось як:

Якщо $x$ це додатне нормальне число:

x=2^{e_{x}}(1+m_{x})

тоді ми маємо

\log _{2}(x)=e_{x}+\log _{2}(1+m_{x})

але оскільки $m x \in [0, 1)$ , логарифм праворуч можна приблизно порахувати через ^[10]

\log _{2}(1+m_{x})\approx m_{x}+\sigma

де $σ$ — це вільний параметр використовуваний для налаштування наближення. Наприклад, $σ = 0$ дає точний результат на обох кінцях інтервалу, тоді як $σ \approx 0.0430357$ дає оптимальне наближення (найкраще у сенсі рівномірної норми похибки).

Отже, ми маємо наближення

\log _{2}(x)\approx e_{x}+m_{x}+\sigma .

З іншого боку, інтерпретування бітового представлення $x$ як цілого дає^{[note 4]}

{\begin{aligned}I_{x}&=E_{x}L+M_{x}\\&=L(e_{x}+B+m_{x})\\&=L(e_{x}+m_{x}+\sigma +B-\sigma )\\&\approx L\log _{2}(x)+L(B-\sigma ).\end{aligned}}

Тоді виявляється, що $I x$ є масштабованим і зсунутим кусково-лінійним наближенням $log 2 (x)$ , як показано на зображенні праворуч. Інакше кажучі, $log 2 (x)$ наближується за допомогою

\log _{2}(x)\approx {\frac {I_{x}}{L}}-(B-\sigma ).

Перше наближення результату

Обчислення $y = 1/ \sqrt x$ базується на тотожності

\log _{2}(y)=-{\frac {1}{2}}\log _{2}(x)

Використовуючи наближення логарифму наведене вище, застосоване до обох $x$ і $y$ , рівняння дає:

{\frac {I_{y}}{L}}-(B-\sigma )\approx -{\frac {1}{2}}{\biggl (}{\frac {I_{x}}{L}}-(B-\sigma ){\biggr )}

З цього, наближення для $I y$ таке:

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

що записано в коді як

i  = 0x5f3759df - ( i >> 1 );

Перший доданок вище це магічне число

{\frac {3}{2}}L(B-\sigma )={\text{0x5f3759df}}

з якого можна зробити висновок, що $σ \approx 0.0450466$ . Другий доданок, $½ I x$ , обрахований через бітовий зсув $I x$ на одну позицію праворуч.^[11]

Метод Ньютона

Докладніше: Метод Ньютона

Після використання цих цілочисельних операцій, алгоритм знов розглядає довге слово як число з рухомою комою (y = *(float*)&i;) і виконує операцію множення із рухомою комою (y = y*(1.5f - xhalf*y*y);). Ця операція представляє одну ітерацію методу Ньютона. Тут ми маємо:

y={\frac {1}{\sqrt {x}}}

— це обернений квадратний корінь, або, як функція від y,

f(y)={\frac {1}{y^{2}}}-x=0

.

As

y_{n+1}=y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}

представляє загальне вираження методу Ньютона із

\,y_{n}

як перше наближення,

y_{n+1}={\frac {y_{n}(3-xy_{n}^{2})}{2}},

де

f(y)={\frac {1}{y^{2}}}-x

і

f'(y)={\frac {-2}{y^{3}}}

.

Тому y = y*(1.5f - xhalf*y*y); є тим самим, що

\,y_{n+1}=y_{n}\left(1.5-{\frac {xy_{n}^{2}}{2}}\right)={\frac {y_{n}(3-xy_{n}^{2})}{2}}

Виноски

↑ Використання типа long зменшує переносність цього коду на сучасні системи. Для того, щоб код виконався правильно, sizeof(long) повинен бути 4 байти, інакше можна отримати від'ємний результат. На багатьох сучасних 64-бітних системах, sizeof(long) становить 8 байтів.
↑ $E x$ має бути в діапазоні $[1, 254]$ для $x$ , щоб бути представна як нормальне число.
↑ Єдиними числами представними точно як числа з рухомою комою це ті у яких $M x$ є цілим. Інші числа можна представити лише приблизно, округлюючи їх до найближчого цілого.
↑ $S x = 0$ оскільки $x > 0$ .

Примітки

↑ ^а ^б Sommefeldt, Rys (29 листопада 2006). Origin of Quake3's Fast InvSqrt(). Beyond3D. Архів оригіналу за 9 лютого 2009. Процитовано 12 лютого 2009.
↑ Discussion on CSDN. Архів оригіналу за 2 липня 2015. Процитовано 8 травня 2016.
↑ Eberly, 2001, с. 504.
↑ Lomont, 2003, с. 1.
↑ McEniry, 2007, с. 1.
↑ Lomont, 2003, с. 3.
↑ McEniry, 2007, с. 2, 16.
↑ Eberly, 2002, с. 2.
↑ Ruskin, Elan (16 жовтня 2009). Timing square root. Some Assembly Required. Архів оригіналу за 18 травня 2015. Процитовано 7 травня 2015. [Архівовано 2015-05-18 у Wayback Machine.]
↑ McEniry, 2007, с. 3.
↑ Hennessey & Patterson 1998, p. 305.

Документи

Blinn, Jim (July 1997). Floating Point Tricks. Computer Graphics & Applications, IEEE. 17 (4): 80. doi:10.1109/38.595279.
Blinn, Jim (2003). Jim Blinn's Corner: Notation, notation notation. Morgan Kaufmann. ISBN 1-55860-860-5.
Eberly, David (2001). 3D Game Engine Design. Morgan Kaufmann. ISBN 978-1-55860-593-0.
Hennessey, John; Patterson, David A. (1998). Computer Organization and Design (вид. 2nd). San Francisco, CA: Morgan Kaufmann Publishers. ISBN 978-1-55860-491-9.
Kushner, David (August 2002). The wizardry of Id. IEEE Spectrum. 39 (8): 42—47. doi:10.1109/MSPEC.2002.1021943.
Lomont, Chris (February 2003). Fast Inverse Square Root (PDF). Архів оригіналу (PDF) за 6 лютого 2009. Процитовано 13 лютого 2009.
McEniry, Charles (August 2007). The Mathematics Behind the Fast Inverse Square Root Function Code (PDF). Архів оригіналу (PDF) за 11 травня 2015. Процитовано 13 лютого 2009.
Middendorf, Lars; Mühlbauer, Felix; Umlauf, George; Bodba, Christophe (1 червня 2007). Embedded Vertex Shader in FPGA. У Rettberg, Achin (ред.). Embedded System Design: Topics, Techniques and Trends. IFIP TC10 Working Conference:International Embedded Systems Symposium (IESS). et al. Irvine, California: Springer. ISBN 978-0-387-72257-3.
Striegel, Jason (4 грудня 2008). Quake's fast inverse square root. Hackszine. O'Reilly Media. Архів оригіналу за 15 лютого 2009. Процитовано 7 січня 2013.

[6] Використання типа long зменшує переносність цього коду на сучасні системи. Для того, щоб код виконався правильно, sizeof(long) повинен бути 4 байти, інакше можна отримати від'ємний результат. На багатьох сучасних 64-бітних системах, sizeof(long) становить 8 байтів.

[11] $E x$ має бути в діапазоні $[1, 254]$ для $x$ , щоб бути представна як нормальне число.

[12] Єдиними числами представними точно як числа з рухомою комою це ті у яких $M x$ є цілим. Інші числа можна представити лише приблизно, округлюючи їх до найближчого цілого.

[14] $S x = 0$ оскільки $x > 0$ .

[Beyond3D-1] а ^б Sommefeldt, Rys (29 листопада 2006). Origin of Quake3's Fast InvSqrt(). Beyond3D. Архів оригіналу за 9 лютого 2009. Процитовано 12 лютого 2009.

[csdn-2] Discussion on CSDN. Архів оригіналу за 2 липня 2015. Процитовано 8 травня 2016.

[FOOTNOTEEberly2001504-3] Eberly, 2001, с. 504.

[FOOTNOTELomont20031-4] Lomont, 2003, с. 1.

[FOOTNOTEMcEniry20071-5] McEniry, 2007, с. 1.

[FOOTNOTELomont20033-7] Lomont, 2003, с. 3.

[FOOTNOTEMcEniry20072,_16-8] McEniry, 2007, с. 2, 16.

[FOOTNOTEEberly20022-9] Eberly, 2002, с. 2.

[ruskin-10] Ruskin, Elan (16 жовтня 2009). Timing square root. Some Assembly Required. Архів оригіналу за 18 травня 2015. Процитовано 7 травня 2015. [Архівовано 2015-05-18 у Wayback Machine.]

[FOOTNOTEMcEniry20073-13] McEniry, 2007, с. 3.

[15] Hennessey & Patterson 1998, p. 305.

[1]

[2]

[3]

[4]

[5]

[note 1]

[6]

[7]

[8]

[9]

[note 2]

[note 3]

[10]

[note 4]

[11]