Nietypowy szybki odwrócony pierwiastek kwadratowy Johna Carmacka (Quake III)

Question

John Carmack pełni specjalną funkcję w kodzie źródłowym Quake III, który oblicza odwrotność pierwiastka kwadratowego float, 4x szybciej niż zwykły (float)(1.0/sqrt(x)), w tym dziwną stałą 0x5f3759df. Zobacz poniższy kod. Czy ktoś może wyjaśnić linijkę po linii, co dokładnie się tutaj dzieje i dlaczego działa to znacznie szybciej niż zwykła implementacja?

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. Carmack tego nie napisał. Terje Mathisen i Gary Tarolli biorą na siebie częściowy (i bardzo skromny) kredyt, jak również kredytują inne źródła.

To, jak powstała mityczna stała, jest czymś tajemniczym.

Cytując Gary Tarolli:

co faktycznie robi pływającą punkt obliczeń w całkowitej - zajęło dużo czasu, aby dowiedzieć się, w jaki sposób i dlaczego to działa, a ja nie pamiętam szczegóły już.

Nieco lepiej na stałym poziomie, opracowana przez matematyka ekspertów (Chris Lomont) stara się wypracować jak oryginalny algorytm pracował to:

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

Pomimo tego, jego początkowa próba matematycznie „superior "Wersja sqrt id (która osiągnęła prawie taką samą wartość 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 jest tak znakomity iirc.

Answer 2

Według to this nice article napisany jakiś czas temu ...

Magia kodu, nawet jeśli nie może obserwować go, wyróżnia się jako i = 0x5f3759df - (I 1 >>) ; linia. Uproszczony, Newton-Raphson to przybliżenie , które rozpoczyna się od zgadywania, a udoskonala je w iteracji. Biorąc wykorzystać charakter 32-bitowych procesorów x86 , I, liczbą całkowitą, jest początkowo ustawiona na wartość liczby zmiennoprzecinkową chcesz wziąć plac odwrotny, za pomocą obsady całkowitą. i jest ustawiony na 0x5f3759df, minus sam przesunął jeden bit o w prawo. Prawe przesunięcie upuszcza najmniej znaczący bit i, , zasadniczo zmniejszając go o połowę.

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

Answer 3

Oczywiście w dzisiejszych czasach okazuje się, że jest znacznie wolniejszy niż użycie sqrt'a FPU (zwłaszcza na 360/PS3), ponieważ zamiana pomiędzy rejestrami float i int powoduje obciążenie sklepu, a jednostka zmiennoprzecinkowa potrafi wykonać odwrotność pierwiastka kwadratowego w sprzęcie.

Pokazuje tylko, w jaki sposób zmiany muszą ewoluować wraz ze zmianami sprzętu.

Answer 4

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

Każda funkcja matematyczna, z pewnymi wyjątkami, mogą być reprezentowane przez sumę wielomianu:

y = f(x) 

może być dokładnie przekształcony:

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

Gdzie A0, A1, A2 ,. .. są stałych. Problem polega na tym, że dla wielu funkcji, takich jak pierwiastek kwadratowy, dla dokładnej wartości ta suma ma nieskończoną liczbę członków, to nie kończy się na jakiejś x^n. Ale jeśli zatrzymamy się na jakimś x^n nadal będziemy mieli wynik do pewnej precyzji.

Tak więc, jeśli mamy:

y = 1/sqrt(x) 

W tym konkretnym przypadku postanowili odrzucić wszystkie wielomianowych członków powyższe sekundę, prawdopodobnie ze względu na szybkość obliczeń:

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

a zadanie ma teraz przyszedł w dół, aby obliczyć a0 i a1, aby y miały najmniejszą różnicę w stosunku do dokładnej wartości. Oni obliczyli, że najbardziej odpowiednie są następujące wartości:

a0 = 0x5f375a86 
a1 = -0.5 

Więc kiedy można umieścić to w równaniu otrzymasz:

y = 0x5f375a86 - 0.5*x 

który jest taki sam jak w wierszu widać w kodzie:

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

Edytuj: faktycznie tutaj y = 0x5f375a86 - 0.5*x to nie to samo, co i = 0x5f375a86 - (i >> 1);, ponieważ zmiana liczby zmiennoprzecinkowej jako liczby całkowitej nie tylko dzieli się przez dwa, ale również dzieli wykładnik o dwie 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 wystarcza do tego celu. Więc dodatkowo zrobił tylko jeden krok iteracji Newtona do poprawy dokładności wyników:

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

Mogli zrobić trochę więcej iteracji w pętli, każdy poprawiając wynik, aż wymagana dokładność jest spełniony. Dokładnie 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 prędkości. CPU/FPU wykonuje tyle iteracji, ile potrzeba, aby osiągnąć dokładność liczby zmiennoprzecinkowej, w której zapisany jest wynik, i ma bardziej ogólny algorytm, który działa we wszystkich przypadkach.

---

Tak w skrócie, to co zrobili to:

Zastosowanie (prawie) taki sam algorytm jak CPU/FPU, wykorzystać poprawę warunków początkowych dla szczególnego przypadku 1/sqrt (x) i nie obliczaj całej drogi do precyzji CPU/FPU, ale zatrzymaj się wcześniej, zyskując w ten sposób prędkość obliczeń.

Answer 5

Byłem ciekawy, jak stała była zmienna, więc napisałem po prostu ten kod i wyszukał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 stałe jest „całkowita przybliżeniem pierwiastka kwadratowego 2^127 lepiej znany przez szesnastkowym postaci jego zmiennoprzecinkowej reprezentacji 0x5f3759df” https://mrob.com/pub/math/numbers-18.html

W tym samym miejscu it wyjaśnia całą sprawę. https://mrob.com/pub/math/numbers-16.html#le009_16