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

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

Швидкий обернений квадратний корінь (іноді згадуваний як Fast InvSqrt() або за шістнадцятковою сталою 0x5f3759df) — це метод обчислення , оберненого квадратного кореня для 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;
}

Для визначення оберненого квадратного кореня визначається наближення для , тоді за допомогою чисельного методу це наближення переглядається, щоб отримати прийнятну похибку у кінцевому результаті. Звичайні програмні методи на початку 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 як ціле, як спосіб приблизного обчислення log2(x)
  2. Використовує це наближення для обчислення наближення log2(1/x)
  3. Знов інтерпретує як число з рухомою комою, як спосіб для обчислення наближення 1/x
  4. Уточнює наближення використовуючи метод Ньютона.

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

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

де показник ex є цілим, mx ∈ [0, 1), і 1.b1b2b3... це двійкове представлення мантиси (1 + mx). Варто зазначити, що оскільки єдиний біт перед комою у мантисі завжди 1, то немає потреби його зберігати. З цієї форми маємо три беззнакові цілі числа:

  • Sx, знаковий біт, це 0 якщо x > 0, і 1 якщо x < 0 (1 біт)
  • Ex = ex + B — це зміщена експонента, де B = 127 — зсув[note 2] (8 бітів)
  • Mx = mx × L, де L = 223[note 3] (23 bits)

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

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

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

  • S = 0
  • E = −3 + 127 = 124 = 011111002
  • M = 0.25 × 223 = 2097152 = 010000000000000000000002

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

Float w significand 2.svg

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

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

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

тоді ми маємо

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

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

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

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

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

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

Перше наближення результату[ред. | ред. код]

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

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

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

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

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

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

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

Метод Ньютона[ред. | ред. код]

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

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

 — це обернений квадратний корінь, або, як функція від y,
.
As представляє загальне вираження методу Ньютона із як перше наближення,
де і .
Тому y = y*(1.5f - xhalf*y*y); є тим самим, що

Виноски[ред. | ред. код]

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

Примітки[ред. | ред. код]

  1. а б Sommefeldt, Rys (November 29, 2006). Origin of Quake3's Fast InvSqrt(). Beyond3D. Процитовано 2009-02-12. 
  2. Discussion on CSDN. 
  3. Eberly, 2001, p. 504.
  4. Lomont, 2003, p. 1.
  5. McEniry, 2007, p. 1.
  6. Lomont, 2003, p. 3.
  7. McEniry, 2007, p. 2, 16.
  8. Eberly, 2002, p. 2.
  9. Ruskin, Elan (2009-10-16). Timing square root. Some Assembly Required. Процитовано 2015-05-07. 
  10. McEniry, 2007, p. 3.
  11. Hennessey & Patterson 1998, p. 305.

Документи[ред. | ред. код]