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

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

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

Алгоритм приймає 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-х отримували перше наближення із таблиці пошуку.Шаблон:Sfn Цей шматок коду виявився швидшим ніж використання таблиці пошуку і приблизно в чотири рази швидший ніж звичайне ділення чисел з рухомою комою.Шаблон:Sfn Хоча деяка втрата точності і відбувалася, але її перекривало значне покращення швидкодії.Шаблон:Sfn Алгоритм був розроблений для специфікації Шаблон:Li 32 бітних чисел з рухомою комою, але подальші дослідження Кріса Ломонта і Чарльза Макінері показали, що його можна реалізувати і для інших специфікацій.

Переваги у швидкості пропоновані швидким оберненим квадратним коренем з'явились завдяки трактуванню довгого слова, що містить число з рухомою комою як цілого і віднімання його від специфічної сталої, 0x5f3759df. Ціль цієї сталої не одразу очевидна для читача коду, отже, як і багато інших сталих знайдених у коді, її називають магічним числом.Шаблон:SfnШаблон:SfnШаблон:Sfn Це цілочисельне віднімання і бітовий зсув дають довге слово, яке знов трактується як число з рухомою комою і є грубим наближенням оберненого квадратного кореня вхідного числа. Одна ітерація методу Ньютона виконується для отримання більшої точності, і код завершується. Алгоритм генерує прийнятно точні результати використовуючи унікальне перше наближення для методу Ньютона; однак, він набагато повільніший ніж використання SSE інструкції rsqrtss на x86 процесорах також випущеної у 1999.

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

Як приклад, розглянемо число Шаблон:Math, для якого ми хочемо обчислити Шаблон:Math. Перші кроки алгоритму проілюстровані нижче:

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

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

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

Алгоритм обчислює Шаблон:Math виконуючи такі кроки:

1. Інтерпретує аргумент Шаблон:Math як ціле, як спосіб приблизного обчислення Шаблон:Math
2. Використовує це наближення для обчислення наближення Шаблон:Math
3. Знов інтерпретує як число з рухомою комою, як спосіб для обчислення наближення Шаблон:Math
4. Уточнює наближення використовуючи метод Ньютона.

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

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

$\begin{array}{cc}x& =\pm 1.b_1{b}_2{b}_3\dots \times 2^{e_x}\\ & =\pm 2^{e_x}(1+m_x)\end{array}$

де показник Шаблон:Math є цілим, Шаблон:Math, і Шаблон:Math це двійкове представлення мантиси Шаблон:Math. Варто зазначити, що оскільки єдиний біт перед комою у мантисі завжди 1, то немає потреби його зберігати. З цієї форми маємо три беззнакові цілі числа:

• Шаблон:Math, знаковий біт, це 0 якщо Шаблон:Math, і 1 якщо Шаблон:Math (1 біт)
• Шаблон:Math — це зміщена експонента, де Шаблон:Math — зсув (8 бітів)
• Шаблон:Math, де Шаблон:Math (23 bits)

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

Як приклад розглянемо число Шаблон:Math. Нормалізація Шаблон:Math дає:

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

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

• Шаблон:Math
• Шаблон:Math
• Шаблон:Math

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

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

Якби комусь довелось порахувати Шаблон:Math без комп'ютера чи калькулятора, то йому б стала в пригоді таблиця логарифмів разом із тотожністю Шаблон:Math, яка дійсна для кожної основи Шаблон:Math. Швидкий обернений квадратний корінь базується на цій тотожності і на факті, що інтерпретація float32 у ціле число дає грубе наближення цього логарифма. Ось як:

Якщо Шаблон:Math це додатне нормальне число:

$x=2^{e_x}(1+m_x)$

тоді ми маємо

${\log }_2(x)=e_x+{\log }_2(1+m_x)$

але оскільки Шаблон:Math, логарифм праворуч можна приблизно порахувати через Шаблон:Sfn

${\log }_2(1+m_x)\approx m_x+\sigma$

де Шаблон:Math — це вільний параметр використовуваний для налаштування наближення. Наприклад, Шаблон:Math дає точний результат на обох кінцях інтервалу, тоді як Шаблон:Math дає оптимальне наближення (найкраще у сенсі рівномірної норми похибки).

Інтерпретування числа з рухомою комою IEEE 754 x як цілого I_x (як-от C: float x = ...; int32_t i = * (int32_t *) &x;) дає масштабоване і зсунуте наближення логарифму з основою 2.

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

${\log }_2(x)\approx e_x+m_x+\sigma .$

З іншого боку, інтерпретування бітового представлення Шаблон:Math як цілого дає

$\begin{array}{cc}{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{array}$

Тоді виявляється, що Шаблон:Math є масштабованим і зсунутим кусково-лінійним наближенням Шаблон:Math, як показано на зображенні праворуч. Інакше кажучі, Шаблон:Math наближується за допомогою

${\log }_2(x)\approx \frac{I_x}{L}-(B-\sigma ).$

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

Обчислення Шаблон:Math базується на тотожності

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

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

$\frac{I_y}{L}-(B-\sigma )\approx -\frac{1}{2}(\frac{I_x}{L}-(B-\sigma ))$

З цього, наближення для Шаблон:Math таке:

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

з якого можна зробити висновок, що Шаблон:Math. Другий доданок, Шаблон:Math, обрахований через бітовий зсув Шаблон:Math на одну позицію праворуч.

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

Шаблон:Main Після використання цих цілочисельних операцій, алгоритм знов розглядає довге слово як число з рухомою комою (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^{\prime }(y_n)}$ представляє загальне вираження методу Ньютона із $y_n$ як перше наближення,

$y_{n+1}=\frac{y_n(3-x{y}_n^2)}{2},$ де $f(y)=\frac{1}{y^2}-x$ і $f^{\prime }(y)=\frac{-2}{y^3}$.

Тому y = y*(1.5f - xhalf*y*y); є тим самим, що $y_{n+1}=y_n(1.5-\frac{x{y}_n^2}{2})=\frac{y_n(3-x{y}_n^2)}{2}$

Примітки

Документи

• Шаблон:Cite journal
• Шаблон:Cite book
• Шаблон:Cite book
• Шаблон:Cite book
• Шаблон:Cite journal
• Шаблон:Cite web
• Шаблон:Cite web
• Шаблон:Cite conference
• Шаблон:Cite web