Niezwykłe Fast Inverse Square Root Johna Carmacka (Quake III)

John Carmack ma specjalną funkcję w kodzie źródłowym Quake III, która oblicza odwrotny pierwiastek kwadratowy z liczby zmiennoprzecinkowej, 4x szybciej niż normalnie (float)(1.0/sqrt(x)), włączając dziwną 0x5f3759dfstałą. Zobacz poniższy kod. Czy ktoś może wyjaśnić wiersz po wierszu, co dokładnie się tutaj dzieje i dlaczego działa to znacznie szybciej niż zwykłe wdrożenie?

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

FYI. Carmack tego nie napisał. Terje Mathisen i Gary Tarolli przypisują to częściowo (i bardzo skromnie), a także innym źródłom.

To, jak wyprowadzono stałą mityczną, pozostaje tajemnicą.

Cytując Gary'ego Tarolliego:

Co w rzeczywistości jest obliczeniem zmiennoprzecinkowym na liczbach całkowitych - zajęło dużo czasu, aby dowiedzieć się, jak i dlaczego to działa, i nie pamiętam już szczegółów.

Nieco lepszą stałą, opracowaną przez eksperta matematyka (Chrisa Lomonta), próbującego ustalić, jak działał oryginalny algorytm, jest:

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

Mimo to, jego początkowa próba matematycznie „lepszej” wersji sqrt id (która osiągnęła prawie taką samą stałą), okazała się gorsza od tej, którą początkowo opracował Gary, mimo że matematycznie była dużo „czystsza”. Nie potrafił wyjaśnić, dlaczego id było tak doskonałe iirc.

Oczywiście w dzisiejszych czasach okazuje się, że jest znacznie wolniejszy niż zwykłe użycie sqrt FPU (szczególnie na 360 / PS3), ponieważ zamiana między rejestrami float i int indukuje magazyn obciążenia, podczas gdy jednostka zmiennoprzecinkowa może wykonać odwrotność kwadratu root w sprzęcie.

Pokazuje tylko, jak muszą ewoluować optymalizacje, gdy zmienia się natura sprzętu.

Greg Hewgill i IllidanS4 podali link z doskonałym wyjaśnieniem matematycznym. Spróbuję podsumować to tutaj dla tych, którzy nie chcą wdawać się w szczegóły.

Dowolną funkcję matematyczną, z pewnymi wyjątkami, można przedstawić za pomocą sumy wielomianów:

y = f(x)

można dokładnie przekształcić w:

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

Gdzie a0, a1, a2, ... są stałymi . Problem polega na tym, że dla wielu funkcji, takich jak pierwiastek kwadratowy, dla dokładnej wartości suma ta ma nieskończoną liczbę elementów, nie kończy się na jakimś x ^ n . Ale jeśli zatrzymamy się na jakimś x ^ n , nadal otrzymamy wynik z pewną precyzją.

Więc jeśli mamy:

y = 1/sqrt(x)

W tym konkretnym przypadku postanowili odrzucić wszystkie składowe wielomianu powyżej sekundy, prawdopodobnie z powodu szybkości obliczeń:

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

Zadanie sprowadza się teraz do obliczenia a0 i a1, aby y miało najmniejszą różnicę od dokładnej wartości. Obliczyli, że najbardziej odpowiednie wartości to:

a0 = 0x5f375a86
a1 = -0.5

Więc kiedy umieścisz to w równaniu, otrzymasz:

y = 0x5f375a86 - 0.5*x

Która jest taka sama jak linia, którą widzisz w kodzie:

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

Edycja: właściwie y = 0x5f375a86 - 0.5*xto nie to samo, co i = 0x5f375a86 - (i >> 1);od czasu, gdy przesunięcie liczby zmiennoprzecinkowej jako liczby całkowitej nie tylko dzieli przez dwa, ale również dzieli wykładnik przez dwa i powoduje inne artefakty, ale nadal sprowadza się do obliczenia niektórych współczynników a0, a1, a2 ....

W tym momencie odkryli, że precyzja tego wyniku nie jest wystarczająca do tego celu. Więc dodatkowo wykonali tylko jeden krok iteracji Newtona, aby poprawić dokładność wyników:

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

Mogli wykonać więcej iteracji w pętli, z których każda poprawiała wynik, aż do osiągnięcia wymaganej dokładności. Dokładnie tak to działa w CPU / FPU! Ale wydaje się, że wystarczyła tylko jedna iteracja, co było również błogosławieństwem dla szybkości. CPU / FPU wykonuje tyle iteracji, ile potrzeba, aby osiągnąć dokładność liczby zmiennoprzecinkowej, w której przechowywany jest wynik, i ma bardziej ogólny algorytm, który działa we wszystkich przypadkach.

Krótko mówiąc, zrobili:

Użyj (prawie) tego samego algorytmu, co CPU / FPU, wykorzystaj poprawę warunków początkowych dla specjalnego przypadku 1 / sqrt (x) i nie obliczaj aż do osiągnięcia precyzji CPU / FPU, ale zatrzymaj się wcześniej, w ten sposób przyspieszenie obliczeń.

Według tego fajnego artykułu, napisanego jakiś czas temu ...

Magia kodu, nawet jeśli nie możesz go śledzić, wyróżnia się jako i = 0x5f3759df - (i >> 1); linia. Upraszczając, Newton-Raphson jest przybliżeniem, które zaczyna się od domysłów i poprawia je za pomocą iteracji. Wykorzystując naturę 32-bitowych procesorów x86, i, liczba całkowita, jest początkowo ustawiana na wartość liczby zmiennoprzecinkowej, dla której chcesz wziąć odwrotność kwadratu, używając rzutowania liczb całkowitych. i jest następnie ustawiane na 0x5f3759df, minus samo siebie przesunięte o jeden bit w prawo. Prawe przesunięcie zmniejsza najmniej znaczący bit i, zasadniczo zmniejszając go o połowę.

To naprawdę dobra lektura. To tylko niewielka część.

Byłem ciekawy, jaka jest stała jako zmiennoprzecinkowa, więc po prostu napisałem ten fragment kodu i wygooglowałem liczbę całkowitą, która się pojawiła.

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

Wygląda na to, że stała to „Przybliżenie liczby całkowitej do pierwiastka kwadratowego z 2 ^ 127, lepiej znane w postaci szesnastkowej jej reprezentacji zmiennoprzecinkowej, 0x5f3759df” https://mrob.com/pub/math/numbers-18.html

Na tej samej stronie wyjaśnia całą sprawę. https://mrob.com/pub/math/numbers-16.html#le009_16