Необычный быстрый квадратный корень Джона Кармака (Quake III)

Question

Джон Кармак имеет специальную функцию в исходном коде Quake III, которая вычисляет обратный квадратный корень с плавающей точкой, в 4 раза быстрее, чем обычная (float)(1.0/sqrt(x)), включая странную 0x5f3759df константу. Смотрите код ниже. Может кто-нибудь построчно объяснить, что именно здесь происходит и почему это работает намного быстрее, чем обычная реализация?

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 ) );

  #ifndef Q3_VM
  #ifdef __linux__
    assert( !isnan(y) );
  #endif
  #endif
  return y;
}

Answer 1

FYI. Кармак этого не написал. Терье Матисен и Гари Таролли принимают за это частичный (и очень скромный) кредит, а также зачисляют некоторые другие источники.

Как возникла мифическая константа, остается загадкой.

Цитировать Гэри Таролли:

Который на самом деле делает плавающий вычисление точки в целых числах - потребовалось долго разбираться как и почему это работает, и я не могу вспомнить подробности больше.

Немного лучшая константа, разработанная экспертом-математиком (Крис Ломонт), пытающаяся выяснить, как работал оригинальный алгоритм:

float InvSqrt(float x)
{
    float xhalf = 0.5f * x;
    int i = *(int*)&x;              // get bits for floating value
    i = 0x5f375a86 - (i >> 1);      // gives initial guess y0
    x = *(float*)&i;                // convert bits back to float
    x = x * (1.5f - xhalf * x * x); // Newton step, repeating increases accuracy
    return x;
}

Несмотря на это, его первоначальная попытка математически «превосходящей» версии sqrt id (которая достигла почти той же константы) оказалась ниже той, которая была первоначально разработана Гэри, несмотря на то, что она математически намного «чище». Он не мог объяснить, почему id был таким превосходным, ирк.

Answer 2

Конечно, в наши дни это происходит намного медленнее, чем просто использование sqrt FPU (особенно на 360 / PS3), потому что переключение между регистрами float и int вызывает сохранение нагрузки, в то время как модуль с плавающей запятой сделать взаимный квадратный корень в аппаратных средствах.

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

Answer 3

Грег Хьюгилл и IllidanS4 дали ссылку с превосходным математическим объяснением. Я постараюсь обобщить это здесь для тех, кто не хочет вдаваться в подробности.

Любая математическая функция, за некоторыми исключениями, может быть представлена ​​полиномиальной суммой:

y = f(x)

может быть точно преобразовано в:

y = a0 + a1*x + a2*(x^2) + a3*(x^3) + a4*(x^4) + ...

Где a0, a1, a2, ... являются константами . Проблема состоит в том, что для многих функций, таких как квадратный корень, для точного значения эта сумма имеет бесконечное число членов, она не заканчивается на некотором x ^ n . Но если мы остановимся на некотором x ^ n , мы все равно получим результат с некоторой точностью.

Итак, если у нас есть:

y = 1/sqrt(x)

В этом конкретном случае они решили отбросить все полиномиальные элементы выше секунды, вероятно, из-за скорости вычисления:

y = a0 + a1*x + [...discarded...]

И теперь задача сводится к вычислению a0 и a1, чтобы у y было наименьшее отличие от точного значения. Они подсчитали, что наиболее подходящие значения:

a0 = 0x5f375a86
a1 = -0.5

Итак, когда вы положите это в уравнение, вы получите:

y = 0x5f375a86 - 0.5*x

То же самое, что строка, которую вы видите в коде:

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

Редактировать: на самом деле здесь y = 0x5f375a86 - 0.5*x - это не то же самое, что i = 0x5f375a86 - (i >> 1);, поскольку смещение числа с плавающей точкой как целого не только делит на два, но также делит экспоненту на два и вызывает некоторые другие артефакты, но все же сводится к вычислению некоторые коэффициенты a0, a1, a2 ....

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

x = x * (1.5f - xhalf * x * x)

Они могли бы сделать еще несколько итераций в цикле, каждый из которых улучшал бы результат, пока не была достигнута требуемая точность. Именно так оно и работает в CPU / FPU! Но, похоже, достаточно было только одной итерации, что также было благословением для скорости. CPU / FPU выполняет столько итераций, сколько необходимо для достижения точности числа с плавающей запятой, в котором сохраняется результат, и имеет более общий алгоритм, который работает для всех случаев.

---

Короче говоря, они сделали следующее:

Используйте (почти) тот же алгоритм, что и CPU / FPU, используйте улучшение начальных условий для особого случая 1 / sqrt (x) и не рассчитывайте весь путь к точности, с которой будет идти CPU / FPU но раньше останавливаться, тем самым увеличивая скорость расчета.

Answer 4

Согласно этой прекрасной статье написанной некоторое время назад ...

Волшебство кода, даже если вы не может следовать за ним, выделяется как я = 0x5f3759df - (i >> 1); линия. Упрощенный, Ньютон-Рафсон является приближением это начинается с догадки и уточняет это с итерацией. принятие Преимущество природы 32-битного x86 Процессоры, я, целое число, это изначально установлен на значение число с плавающей запятой, которое вы хотите взять обратный квадрат, используя целочисленное приведение. Затем я настроен на 0x5f3759df, сам минус один сдвинут немного вправо. Правильный сдвиг отбрасывает наименее значимый бит I, по существу, вдвое.

Это действительно хорошее чтение. Это всего лишь маленький кусочек.

Answer 5

Мне было любопытно посмотреть, что такое константа как число с плавающей точкой, поэтому я просто написал этот фрагмент кода и гуглил выдаваемое число.

    long i = 0x5F3759DF;
    float* fp = (float*)&i;
    printf("(2^127)^(1/2) = %f\n", *fp);
    //Output
    //(2^127)^(1/2) = 13211836172961054720.000000

Похоже, константа имеет вид «Целочисленное приближение к квадратному корню из 2 ^ 127, более известное по шестнадцатеричной форме его представления с плавающей точкой, 0x5f3759df» https://mrob.com/pub/math/numbers-18.html

На том же сайте все объясняется. https://mrob.com/pub/math/numbers-16.html#le009_16