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

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

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

Алгоритм приймає 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;
}

Для визначення оберненого квадратного кореня визначається наближення для ${\displaystyle 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/√x ≈ 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/√x виконуючи такі кроки:

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

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

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

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

${\displaystyle {\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 ∈ [0, 1), і 1.b₁b₂b₃... це двійкове представлення мантиси (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 × L, де L = 2^{23[note 3]} (23 bits)

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

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

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

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

• S = 0
• E = −3 + 127 = 124 = 01111100₂
• M = 0.25 × 2²³ = 2097152 = 01000000000000000000000₂

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

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

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

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

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

тоді ми маємо

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

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

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

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

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

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

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

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

${\displaystyle {\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₂(x), як показано на зображенні праворуч. Інакше кажучі, log₂(x) наближується за допомогою

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

${\displaystyle y={\frac {1}{\sqrt {x}}}}$ — це обернений квадратний корінь, або, як функція від y,

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

As ${\displaystyle y_{n+1}=y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}}$ представляє загальне вираження методу Ньютона із ${\displaystyle \,y_{n}}$ як перше наближення,

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

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

Виноски

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