Швидкий обернений квадратний корінь (іноді згадуваний як Fast InvSqrt() або за шістнадцятковою сталою 0x5f3759df) — це метод обчислення , оберненого квадратного кореня для 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; }
Для визначення оберненого квадратного кореня визначається наближення для , тоді за допомогою чисельного методу це наближення переглядається, щоб отримати прийнятну похибку у кінцевому результаті. Звичайні програмні методи на початку 1990-х отримували перше наближення із таблиці пошуку. Цей шматок коду виявився швидшим ніж використання таблиці пошуку і приблизно в чотири рази швидший ніж звичайне ділення чисел з рухомою комою. Хоча деяка втрата точності і відбувалася, але її перекривало значне покращення швидкодії. Алгоритм був розроблений для специфікації 32 бітних чисел з рухомою комою, але подальші дослідження Кріса Ломонта і Чарльза Макінері показали, що його можна реалізувати і для інших специфікацій.
Переваги у швидкості пропоновані швидким оберненим квадратним коренем з'явились завдяки трактуванню довгого словаВикористання типа long зменшує переносність цього коду на сучасні системи. Для того, щоб код виконався правильно, sizeof(long) повинен бути 4 байти, інакше можна отримати від'ємний результат. На багатьох сучасних 64-бітних системах, sizeof(long) становить 8 байтів., що містить число з рухомою комою як цілого і віднімання його від специфічної сталої, 0x5f3759df. Ціль цієї сталої не одразу очевидна для читача коду, отже, як і багато інших сталих знайдених у коді, її називають магічним числом. Це цілочисельне віднімання і бітовий зсув дають довге слово, яке знов трактується як число з рухомою комою і є грубим наближенням оберненого квадратного кореня вхідного числа. Одна ітерація методу Ньютона виконується для отримання більшої точності, і код завершується. Алгоритм генерує прийнятно точні результати використовуючи унікальне перше наближення для методу Ньютона; однак, він набагато повільніший ніж використання SSE інструкції rsqrtss на x86 процесорах також випущеної у 1999.
Робочий приклад
Як приклад, розглянемо число , для якого ми хочемо обчислити . Перші кроки алгоритму проілюстровані нижче:
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
Інтерпретування останнього бітового представлення як числа з рухомою комою дає наближення , яке має похибку близько 3.4%. Після однієї ітерації метода Ньютона, кінцевим результатом є , і помилка становить лише 0.17%.
Перебіг алгоритму
Алгоритм обчислює виконуючи такі кроки:
Інтерпретує аргумент як ціле, як спосіб приблизного обчислення Використовує це наближення для обчислення наближення Знов інтерпретує як число з рухомою комою, як спосіб для обчислення наближення Уточнює наближення використовуючи метод Ньютона.
Представлення чисел з рухомою комою
Оскільки алгоритм сильно покладається на представлення чисел одинарної точності з рухомою комою на бітовому рівні, короткий огляд цього представлення наведений тут. Для того, щоб закодувати ненульове дійсне число як число із рухомою комою одинарної точності, перший крок полягає в записуванні як нормалізованого двійкового числа:
де показник є цілим, , і це двійкове представлення мантиси . Варто зазначити, що оскільки єдиний біт перед комою у мантисі завжди 1, то немає потреби його зберігати. З цієї форми маємо три беззнакові цілі числа:
, знаковий біт, це 0 якщо , і 1 якщо (1 біт) — це зміщена експонента, де — зсув має бути в діапазоні для , щоб бути представна як нормальне число. (8 бітів) , де Єдиними числами представними точно як числа з рухомою комою це ті у яких є цілим. Інші числа можна представити лише приблизно, округлюючи їх до найближчого цілого. (23 bits)
Ці поля пакуються зліва направо у 32 бітовий контейнер.
Як приклад розглянемо число . Нормалізація дає:
і отже, три беззнакові цілочисельні поля такі:
ці поля пакуються як показано нижче:
Інтерпретування цілим як приблизний логарифм
Якби комусь довелось порахувати без комп'ютера чи калькулятора, то йому б стала в пригоді таблиця логарифмів разом із тотожністю , яка дійсна для кожної основи . Швидкий обернений квадратний корінь базується на цій тотожності і на факті, що інтерпретація float32 у ціле число дає грубе наближення цього логарифма. Ось як:
Якщо це додатне нормальне число:
тоді ми маємо
але оскільки , логарифм праворуч можна приблизно порахувати через
де — це вільний параметр використовуваний для налаштування наближення. Наприклад, дає точний результат на обох кінцях інтервалу, тоді як дає оптимальне наближення (найкраще у сенсі рівномірної норми похибки).
Отже, ми маємо наближення
З іншого боку, інтерпретування бітового представлення як цілого дає оскільки .
Тоді виявляється, що є масштабованим і зсунутим кусково-лінійним наближенням , як показано на зображенні праворуч. Інакше кажучі, наближується за допомогою
Перше наближення результату Обчислення базується на тотожності
Використовуючи наближення логарифму наведене вище, застосоване до обох і , рівняння дає:
З цього, наближення для таке:
що записано в коді як
i = 0x5f3759df - ( i >> 1 );
Перший доданок вище це магічне число
з якого можна зробити висновок, що . Другий доданок, , обрахований через бітовий зсув на одну позицію праворуч.Hennessey & Patterson 1998, p. 305.
Метод Ньютона
Після використання цих цілочисельних операцій, алгоритм знов розглядає довге слово як число з рухомою комою (y = *(float*)&i;) і виконує операцію множення із рухомою комою (y = y*(1.5f - xhalf*y*y);). Ця операція представляє одну ітерацію методу Ньютона. Тут ми маємо:
— це обернений квадратний корінь, або, як функція від y,
.
As представляє загальне вираження методу Ньютона із як перше наближення,
де і .
Тому y = y*(1.5f - xhalf*y*y); є тим самим, що
Виноски
Примітки
Документи