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

Ответ 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;
}

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

Ответ 2

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

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

Ответ 3

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

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

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);, поскольку сдвиг float как целого не только делит на два, но также делит показатель на два и вызывает некоторые другие артефакты, но он все же сводится к вычислению некоторых коэффициентов a0, a1, a2....

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

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

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

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

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

Ответ 4

В соответствии с в эту приятную статью, написанную некоторое время назад...

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

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

Ответ 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