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

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

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

Процитирую Гэри Таролли:

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

Чуть лучшая константа, , разработанная опытным математиком (Крис Ломонт), пытаясь чтобы выяснить, как работал исходный алгоритм:

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

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

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

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

Greg Hewgill и 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 будет идти, но останавливаться раньше, тем самым увеличивая скорость вычислений.

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

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

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

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

    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