Швидкий обернений квадратний корінь (іноді згадуваний як 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]) На той час, основна перевага алгоритму полягала у використанні замість обчислювально дорогих операцій над числами з рухомою комою операцій над цілими числами. Обернений квадратний корінь використовують для обчислення кутів падіння і відбивання для освітлення і шейдинга в комп'ютерній графіці.

Алгоритм приймає 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 дає оптимальне наближення (найкраще у сенсі рівномірної норми похибки).

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

${\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.

• 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.