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


112

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


10
Pisano to o milionach razy. Zobacz: google.com/search?q=0x5f3759df
Greg Hewgill

15
W każdym razie dzięki. To było o wiele ciekawsze pytanie niż „jak zrobić ujemną liczbę dodatnią w C #?”
MusiGenesis


7
O cholera, to tylko hack oparty na metodzie Newtona, to nie jest jakiś święty Graal algorytmów, przestań o tym mówić, proszę: P
ldog

Odpowiedzi:


75

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.


4
Co ma znaczyć „matematycznie czystszy”?
Tara

1
Wyobrażałbym sobie, gdzie pierwsze przypuszczenie można wyprowadzić z uzasadnionych stałych, zamiast być pozornie arbitralne. Chociaż jeśli potrzebujesz opisu technicznego, możesz go sprawdzić. Nie jestem matematykiem, a dyskusja semantyczna na temat terminologii matematycznej nie należy do SO.
Rushyo

7
To jest dokładnie powód I obudowane to słowo w cudzysłowie przestraszyć, aby uniknąć tego rodzaju bzdury. Zakładam, że czytelnik jest zaznajomiony z potocznym angielskim pismem, jak sądzę. Można by pomyśleć, że wystarczy zdrowy rozsądek. Nie użyłem niejasnego terminu, ponieważ pomyślałem „wiesz co, naprawdę chcę, aby ktoś zapytał mnie o to, komu nie chce się znaleźć oryginalnego źródła, co zajęłoby dwie sekundy w Google”.
Rushyo,

2
Cóż, właściwie nie odpowiedziałeś na pytanie.
BJovke

1
Dla tych, którzy chcieli wiedzieć, gdzie go znalazł: outside3d.com/content/articles/8
mr5

52

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.


4
Jest jednak dużo szybszy niż std :: sqrt ().
Tara

2
Czy masz źródło? Chcę przetestować środowiska wykonawcze, ale nie mam zestawu deweloperskiego dla konsoli Xbox 360.
DucRP

31

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ń.


2
Rzutowanie wskaźnika na długość jest przybliżeniem log_2 (liczba zmiennoprzecinkowa). Odrzucenie go z powrotem jest przybliżeniem długości 2 ^. Oznacza to, że stosunek można uczynić w przybliżeniu liniowym.
wizzwizz4

22

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ęść.


19

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


6
To zasługuje na więcej uwagi. Wszystko ma sens po uświadomieniu sobie, że to tylko pierwiastek kwadratowy z 2 ^ 127 ...
u8y7541
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.