Fast Quake Approximation of Inverse Square Root

readme.md

Немного о Quake III

Я хочу написать небольшой отзыв из своей головы о этом небезызвестном примере оптимизации процессов. Q_rsqrt Функция, которая вычисляет обратный квадратный корень (схематично: 1/√x). В компьютерной графике Возможно часто в играх приходится вычислять нормализацию векторов. Для нормализации вектора, как и в нормальном уравнении прямой есть подобный способ - нахождение нормирующего элемента.

Будем считать, что существует вектор $\overrightarrow{v}$, тогда для его нормализации следует рассчитать такую вот структуру:

$$n = \frac {1} { \sqrt{v^2_x, x^2_y} }$$

Такие операции, поскольку речь идет о реальных условиях, где важна любая дробная часть числа, (как я понимаю), считающему устройству (CPU) даются не просто.

Почему не просто????

Ну, предположим, если брать в пример Архитектуру процессоров Intel: IA32 (или ее более известное имя: x86), Вычисление (всем известных) float значений, происходит приблизительно так

Опрелеление понятия чисел. Восприятие чисел процессором

Если я помню правильно, то давным давно я много читал и смотрел, есть стандарт IEEE 754, который подробно объясняет структуру сегментов дробного числа (или числа с плавающей точкой, как удобно). Внутри стандарта выделю пару моментов:

• Сегмент знака 1 бит
• Экспонента.
• Мантисса.

Сложно быстро обойти этот разговор, но на счет второго и третьего пункта прийдется мне "вставить свои 5 копеек".

Экспонента это некоторое количество сегментов, которое показывает в какую степень необходимо возвести число, чтобы получить нужное значение.

Мягко говоря это показатель степени, в который нужно возвести основание (2), чтобы получить нужное значение.

$$ 2^{e - f} $$

Что ж такое f... Буквой f я обозначил двоичное смещение, которое "поворачивает" отрицательное значение к положительному Раз уж речь про дроби, то про вид $10^+n$ и речи быть не может, согласны же...?

Мантисса - как я понимаю, это вся значимая дробная часть числа

Поскольку для IA32 определяют разные точности дробного числа, дробная часть числа (или "significant") будет тоже разная

• float 23
• double 52 (больше не знаю)

// [0.75]
// Переводится в двоичную систему:
t = 0.5 + 0.25 = 1 * 2^-1 + 1 * 2^-2 = 0.11; // в двоичном виде.
// Нормализуется:
0.11 = 1.1 * 2^-1;
// Ищется дробная часть: 10000000000000000000000 (23 знака)
// Экспонента (биты) -1 со смещением: 126 -> 01111110 (8 знаков) 
// Знак: 0 (положительно)

// Полное представление (float)0.75:
0 01111110 10000000000000000000000;

А теперь вообразите, как считать какие-нибудь нано-отклонения (например: -123,0177352458) по такому же стандарту. Далее должен быть нормальный алгоритм вычисления каких-нибудь выражений, например тот же обратный квадратный корень суммы квадратов координат, но этому и так будет посвящен ни один раздел здесь.

Что происходит в Q_rsqrt?

Долго я изучал этот пример, и несколько компиляторов ведут с этим примером себя совершенно по разному. Разберу все по частям

Преобразование float в long:

i = *(long*)&y;

В памяти компьютера числа с плавающей точкой (float) хранятся в виде набора битов, как и целые числа (long). Эта строчка становится битовым представлением float и интерпретирует его как целое число i. Это работает, потому что и float и long обычно имеют одинаковый размер 32 бита НО в 32-разрядной Операционной системе... или 32-разрядной архитектуре (Раз уж речь про x86, то именно x86, а не x86-64)

Это не приведение типов, а реинтерпретация битов. (Возможно известно всем как Strict Aliasing). С таким я лично встречаюсь очень редко, но все же помню об опасности таких трюков.

Собственно, раскрываю свои предположения... Что же все-таки такое 0x5f3759df?! и почему мы сдвигаем i вправо на 1 бит (i >> 1) и вычитаем результат?

0x5f3759df - Это константа - подобранное значение, которое дает очень хорошее приближение к обратному квадратному корню. По другому не могу ничего предположить, потому что задача буквально одна, и решений ее тоже немного.

Она базируется на анализе формата представления чисел с плавающей запятой IEEE 754, о котором я писал выше...

i >> 1 - Сдвиг вправо на 1 бит.

Побайтовые сдвиги (i >> 1) еще можно прочесть как i /= 2.

Почему это работает?

Если коротко, то эта операция с битовым представлением float приблежает преобразование, которое, в грубом приближении, соответсвует

$$-\log_2{\sqrt{x}}$$

где x - исходное число. Эта операция основана на том, как в IEEE 754 хранится значение экспоненты числа.

Обратное преобразование:

y  = * ( float * ) &i;

Здесь же все буквально наоборот. (с маленькой оговоркой о разных внутренностях типа данных) Берется битовое представление целого числа i и реинтерпритируется как битовое представление дроби (float)

В целом, выглядит умопомрачающе. Но, как можно заметить с оговорками, поэтому...

Подводные камни?

Как раньше уже было написано, Strict Aliasing опасный прием в нетипобезопасной среде. Что я точно могу утверждать, это очень важно вспомнить замечания в стандартах C/++

Если у вас есть указатель на один тип данных (например, float* ), вы не имеете права обращаться к памяти, на которую он указывает, как если бы она содержала данные другого типа (например, long).

Это очень важное замечание, потому что в коде Quake как раз это нарушается, и соответственно ведет к неопределенному поведению. (нельзя обращаться к адресу дробного числа (float) указателем long*)

Всмысле нельзя? Все будет нормально.

Думаю, что все-таки не будет нормально. Потому что Разные архитектуры могут иметь разные требования к выравниванию данных. При реинтерпретации через указатель можно прочитать данные, которые не выровнены по правильным границам, что приведет к исключениям (?) или неверному ответу.

Поэтому Quake III может завестись нормально на 32-разрядной машине с Intel i386 архитектурой на борту, а вот в 64-разрядной системе (без слоя System WOW64), полагаю уже нет. (даже не смотря на то, что много кто вспомнит о том, что x86-64 это просто расширение x86 архитектуры) Мало того, достаточно вспомнить чем опасно UB.

Вставлю "свои 5 копеек" еще раз, потому что тоже считаю очень важным момент, что при сборке используя C/++, проект проходит минимум 2 большие пытки, такие как "Компиляция" и сама "Сборка" (Линковка). И все неопределенное поведение скорее всего решается уже на этапе линковки, потому что флаги оптимизации и проверок тоже очень сильно меняют результат.

У CLang, например есть три уровня представления кода:

• LLVM IR
• Machine IR
• Сам Ассемблер

И по ходу дела, компилятор и линковщик, с определенными флагами ходят проверяют и переводят проект в бинарный файл. Современные компиляторы могут (часто будут) вырезать такой код в целях оптимизации.

Вывод

Озадачился я вопросом о этом примере очень давно, но только спустя несколько месяцев добрались руки до детального анализа и такого щитпоста. Эксперементы с реинтерпритацией битов были сделаны на Windows 11, Open SUSE 15.6, Windows 2000 В инструментах

• Руки
• VirtualBox
• CLang / LLVM
• GCC
• Mingw

Было очень интересно напрячь голову и с горем пополам но разобраться, чем же так примечателен этот пример оптимизации.

quake_approximation_of_inv_sqrt.c

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));
  // y = y = y * (threehalfs - (x2 * y * y));
  return y;
}