Szybka odwrotność pierwiastka kwadratowego

Z Wikipedii, wolnej encyklopedii
Skocz do: nawigacja, szukaj
Iluminacje i odbicia (ukazane tutaj w grze FPS OpenArena) wykorzystują kod szybkiej odwrotności pierwiastka kwadratowego do obliczenia kątów padania i odbicia

Szybka odwrotność pierwiastka kwadratowego, czasami określana jako szybki InvSqrt() lub przez stałą szesnastkową 0x5f3759df – metoda obliczania x-1/2, odwrotności pierwiastka kwadratowego z 32-bitowej liczby zmiennoprzecinkowej w standardzie IEEE 754. Algorytm został prawdopodobnie opracowany w Silicon Graphics na początku lat dziewięćdziesiątych XX wieku, a implementacja pojawiła się w 1999 roku w kodzie źródłowym gry Quake III: Arena, lecz metoda nie pojawiła się na forach publicznych, takich jak Usenet aż do roku 2002 lub 2003[1]. W ówczesnym czasie podstawowa zaleta algorytmu polegała na unikaniu kosztownych obliczeniowo operacji zmiennoprzecinkowych na korzyść operacji na liczbach całkowitych. Odwrotności pierwiastków kwadratowych są używane do obliczania kątów padania i odbicia dla oświetlenia i cieniowania w grafice komputerowej.

Algorytm na wejściu przyjmuje dodatnią 32-bitową liczbę zmiennoprzecinkową i zachowuje jej połowę do późniejszego wykorzystania. Następnie, interpretując bitową reprezentację liczby zmiennoprzecinkowej jako 32-bitową liczbę całkowitą, dokonuje przesunięcia bitowego w prawo o jeden bit, a wynik odejmuje od „magicznej” wartości 0x5f3759df. Jest to pierwsze przybliżenie odwrotności pierwiastka kwadratowego z argumentu. Interpretując ponownie reprezentację bitową jako liczbę zmiennoprzecinkową stosuje jedną iterację metody Newtona do wyznaczenia dokładniejszego przybliżenia. Ten algorytm oblicza przybliżenie odwrotności pierwiastka kwadratowego z liczby zmiennoprzecinkowej około cztery razy szybciej niż dzielenie zmiennoprzecinkowe.

Pierwotnie algorytm ten został przypisany Johnowi Carmackowi, lecz badania wykazały, że ten kod ma głębsze korzenie, zarówno w sprzęcie, jak i po stronie oprogramowania grafiki komputerowej. Korekty i zmiany wprowadziły tak Silicon Graphics, jak i 3dfx Interactive, w tym samym czasie, co najwcześniej znane zastosowanie implementacji Gary'ego Tarolliego dla SGI Indigo. Nie wiadomo, jak początkowo została wyznaczona stała, lecz badania rzucają pewne światło na tę kwestię.

Uzasadnianie[edytuj | edytuj kod]

Pole wektorów normalnych
Wektory normalne są szeroko wykorzystywane przy obliczaniu oświetlenia i cieniowania; wymaga to wyznaczenia normy wektorów

Odwrotność pierwiastka kwadratowego pozwala wyznaczyć wektor unormowany[2]. Z uwagi na to, że programy do grafiki trójwymiarowej używają tych unormowanych wektorów do określenia oświetlenia i odbicia, w każdej sekundzie muszą być wykonane miliony takich obliczeń. Zanim wprowadzono moduły sprzętowe do transformacji i oświetlenia, obliczenia wykonywane programowo mogły być powolne. Zwłaszcza na początku lat dziewięćdziesiątych XX wieku, kiedy powstawał ten algorytm, prędkość obliczeń zmiennoprzecinkowych pozostawała mocno w tyle za obliczeniami całkowitoliczbowymi[1].

Aby unormować wektor należy wyznaczyć jego długość przez obliczenie jego normy euklidesowej: pierwiastka kwadratowego z sumy kwadratów jego składowych. Jeśli każdą składową wektora podzieli się przez tę długość, powstanie wektor jednostkowy o takim samym kierunku i zwrocie.

Norma euklidesowa wektora, analogiczna do odległości euklidesowej między dwoma punktami w przestrzeni euklidesowej dana jest wzorem:

\|\boldsymbol{v}\| = \sqrt{v_1^2+v_2^2+v_3^2}

a wektor unormowany (jednostkowy):

\boldsymbol{\hat{v}} = \boldsymbol{v} / \|\boldsymbol{v}\|

Zastępując wyrażenie \scriptstyle v_1^2+v_2^2+v_3^2 przez \scriptstyle \boldsymbol{x} otrzymujemy związek:

\boldsymbol{\hat{v}} = \boldsymbol{v} / \sqrt{\boldsymbol{x}}

który łączy wektor jednostkowy z odwrotnością pierwiastka kwadratowego z sumy kwadratów składowych.

Quake III Arena wykorzystuje algorytm szybkiej odwrotności pierwiastka kwadratowego aby przyspieszyć obliczenia procesora graficznego, jednak od tego czasu ten algorytm został już także zrealizowany w niektórych sprzętowych implementacjach Vertex Shader za pomocą FPGA[3].

Przegląd kodu[edytuj | edytuj kod]

Poniższy kod zawiera implementację szybkiej odwrotności pierwiastka kwadratowego z Quake III: Arena, z pominięciem dyrektyw preprocesora C, lecz zawierającym oryginalne komentarze[4]:

float Q_rsqrt( float number )
{
	long i;
	float x2, y;
	const float threehalfs = 1.5F;

	x2 = number * 0.5F;
	y  = number;
	i  = * ( long * ) &y;                       // evil floating point bit level hacking
	i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
	y  = * ( float * ) &i;
	y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//      y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

	return y;
}

W celu wyznaczenia odwrotności pierwiastka kwadratowego w oprogramowaniu zastosowano przybliżenie dla x-1/2, a następnie przy pomocy pewnej metody numerycznej sprawdzano czy uzyskane przybliżenie mieści się w dopuszczalnym zakresie błędu względem rzeczywistego wyniku. Powszechnie stosowane metody numeryczne z początku lat dziewięćdziesiątych XX wieku pierwsze przybliżenie odczytywały z tabeli[5]. Ten fragment kodu okazał się szybszy niż tablicowanie i około cztery razy szybszy niż zwyczajne dzielenie zmiennoprzecinkowe[6]. Występuje pewna utrata dokładności, lecz jest ona rekompensowana znaczącym wzrostem wydajności[7]. Algorytm został zaprojektowany dla 32-bitowej specyfikacji liczb zmiennoprzecinkowych standardu IEEE 754, lecz Chris Lemont, a później Charles McEniry, wykazali w swoich badaniach, że podobne implementacje można uzyskać dla innych reprezentacji zmiennoprzecinkowych.

Korzyści w prędkości oferowane przez szybką odwrotność pierwiastka kwadratowego pochodzą od interpretowania słowa rozmiaru long[a] zawierającego liczbę zmiennoprzecinkową jako liczbę całkowitą a następnie odjęcie jej od określonej stałej, 0x5f3759df. Dla osoby przeglądającej kod nie jest od razu jasne jaki jest cel wprowadzenia tej stałej, tak więc, jak wiele podobnych stałych w kodzie, nazywa się ją „liczbą magiczną”[1][8][9][10]. To całkowitoliczbowe odejmowanie i przesunięcie bitowe zwraca słowo, które zinterpretowane jako liczba zmiennoprzecinkowa jest przybliżeniem odwrotności pierwiastka kwadratowego argumentu funkcji. Po wykonaniu jednej iteracji metodą Newtona wymagana dokładność jest osiągnięta a kod zakończony. Ten algorytm zwraca wyniki o rozsądnej dokładności wykorzystując unikalne pierwsze przybliżenie dla metody Newtona, chociaż jest znacznie wolniejszy i mniej dokładny niż instrukcja SSE rsqrtss na procesorach x86, wydanych także w 1999 roku[11].

Przykład działania[edytuj | edytuj kod]

Rozważmy przykład z liczbą \scriptstyle x = 0,15625 dla której chcemy wyznaczyć wielkość \scriptstyle 1/\sqrt{x} \approx 2,52982. Pierwsze kroki algorytmu wyglądają następująco:

00111110001000000000000000000000  Bitowy wzór x lub i
00011111000100000000000000000000  Przesunięcie w prawo o jedną pozycję: (i >> 1)
01011111001101110101100111011111  Magiczna stała 0x5f3759df
01000000001001110101100111011111  Wynik obliczeń 0x5f3759df - (i >> 1)

Interpretacja ostatniego wzoru bitowego jako liczby zmiennoprzecinkowej daje przybliżenie \scriptstyle y = 2,61486, które różni się od dokładnego wyniku 0 3,4%. Po jednej iteracji metodą Newtona ostateczny wynik wynosi \scriptstyle y = 2.52549 z uchybem tylko 0,17%.

Opis algorytmu[edytuj | edytuj kod]

Algorytm wyznacza \scriptstyle 1/\sqrt{x} w następujących krokach:

  1. interpretuje parametr \scriptstyle x jako liczbę całkowitą, co może być sposobem na wyznaczenie przybliżenia \scriptstyle \log_{2} x
  2. przekształca wyliczone przybliżenie aby wyznaczyć przybliżoną wartość \scriptstyle \log_{2}(1/\sqrt{x})
  3. interpretuje otrzymaną liczbę całkowitą jako zmiennoprzecinkową, co pozwala obliczyć przybliżenie funkcji wykładniczej o podstawie 2
  4. zwiększa dokładność wyniku stosując jedną iterację metody Newtona

Liczby zmiennoprzecinkowe[edytuj | edytuj kod]

 Główny artykuł: Liczba zmiennoprzecinkowa.

Z uwagi na fakt, że algorytm silnie polega na wewnętrznej reprezentacji binarnej formatu zmiennoprzecinkowego niezbędne jest krótkie wprowadzenie. W celu zapisania niezerowej liczby rzeczywistej \scriptstyle x w postaci zmiennoprzecinkowej o pojedynczej precyzji pierwszym krokiem jest przedstawienie \scriptstyle x w postaci liczby znormalizowanej:

\begin{align}
x &= \pm 1.b_1b_2b_3\ldots \times 2^{e_x}\\
  &= \pm 2^{e_x} (1 + m_x)
\end{align}

gdzie wykładnik \scriptstyle e_x jest liczbą całkowitą, \scriptstyle m_x \in [0, 1), a \scriptstyle 1,b_1b_2b_3\ldots binarną wartością mantysy \scriptstyle (1 + m_x). Należy pamiętać, że bit przed przecinkiem jest zawsze równy 1 i nie musi być zapisywany. Z tej postaci wyznaczane są trzy liczby całkowite:

  • \scriptstyle S_x – bit znaku, który przyjmuje wartość 0 jeśli \scriptstyle x > 0 lub 1 jeśli \scriptstyle x < 0 (1 bit)
  • \scriptstyle E_x = e_x + B – wykładnik ze składową stałą, w którym \scriptstyle B = 127[b] (8 bitów)
  • \scriptstyle M_x = m_x \times L, gdzie \scriptstyle L = 2^{23}[c] (23 bity)

Te trzy wartości są następnie umieszczane od lewej do prawej w 32 bitowym pojemniku.

Przykładem niech będzie liczba \scriptstyle x = 0,15625 = 0,00101_2. Jej postać znormalizowana to:

x = +2^{-3}(1 + 0.25)

dla której odpowiednie trzy pola bitowe wynoszą

  • \scriptstyle S = 0
  • \scriptstyle E = -3 + 127 = 124 = 01111100_2
  • \scriptstyle M = 0,25 \times 2^{23} = 01000000000000000000000_2

oraz po upakowaniu w całość dają wartość jak na poniższej ilustracji:

Przykład liczby zmiennoprzecinkowej.svg

Przybliżenie logarytmu[edytuj | edytuj kod]

Chcąc wyznaczyć \scriptstyle 1/\sqrt{x} bez komputera lub kalkulatora można skorzystać z tablic logarytmów i tożsamości \scriptstyle \log_b(1/\sqrt{x}) = -\frac{1}{2} \log_b(x), która zachodzi dla każdej podstawy \scriptstyle b. Algorytm szybkiej odwrotności pierwiastka opiera się na tej zasadzie oraz fakcie, że interpretacja liczby w notacji zmiennoprzecinkowej jako liczby całkowitej daje zgrubne przybliżenie jej logarytmu:

Jeśli \scriptstyle x jest normalną liczbą dodatnią:

x = 2^{e_x} (1 + m_x)

to zachodzi

\log_2(x) = e_x + \log_2(1 + m_x)

lecz skoro \scriptstyle m_x \in [0, 1) to logarytm po prawej stronie równości można przybliżyć[12]

\log_2(1 + m_x) \approx m_x + \sigma

gdzie \scriptstyle \sigma jest współczynnikiem do polepszania przybliżenia. Na przykład \scriptstyle \sigma = 0 ustala dokładne wyniki na obu końcach przedziału, natomiast \scriptstyle \sigma = 0,0430357 zwraca najmniejsze błędy stosując przybliżenie funkcji wielomianem.

Liczba całkowita interpretowana jako zmiennoprzecinkowa (niebieski), w porównaniu do przeskalowanego i przesuniętego logarytmu (szary)

Stąd otrzymane przybliżenie to

\log_2(x) \approx e_x + m_x + \sigma

Z drugiej strony interpretacja bitowej reprezentacji \scriptstyle x jako liczby całkowitej daje[d]

\begin{align}
I_x &= E_x L + M_x\\
    &= L (e_x + B + m_x)\\
    &\approx L \log_2(x) + L (B - \sigma).
\end{align}

Okazuje się, że \scriptstyle I_x to przeskalowane i przesunięte liniowe przybliżenie \scriptstyle \log_2(x), przedstawione na ilustracji po prawej. Innymi słowy \scriptstyle \log_2(x) jest przybliżony przez

\log_2(x) \approx \frac{I_x}{L} - (B - \sigma)

Pierwsze przybliżenie wyniku[edytuj | edytuj kod]

Obliczanie \scriptstyle y=1/\sqrt{x} opiera się na tożsamości

\log_2(y) = - \frac{1}{2}\log_2(x)

Stosując przedstawione wyżej przybliżenie logarytmu do \scriptstyle x i \scriptstyle y w tym równaniu otrzymuje się:

I_y \approx \frac{3}{2} L (B - \sigma) - \frac{1}{2} I_x

który w kodzie źródłowym objawia się jako

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

Pierwszy składnik to magiczna liczba

\frac{3}{2} L (B - \sigma) = \text{0x5f3759df}

z której można wywnioskować, że \scriptstyle \sigma = 0,0450466. Drugi sładnik, \scriptstyle \frac{1}{2} I_x jest obliczony stosując przesunięcie bitowe wobec \scriptstyle I_x o jedną pozycję w prawo[13].

Metoda Newtona[edytuj | edytuj kod]

 Osobny artykuł: Metoda Newtona.

Po wykonaniu operacji całkowitoliczbowych, algorytm ponownie traktuje długie słowo jako liczbę zmiennoprzecinkową (y = *(float*)&i;) i wykonuje operację mnożenia zmiennoprzecinkowego (y = y*(1.5f - x2*y*y);). Operacja zmiennoprzecinkowa przedstawia jedną iterację metody Newtona szukania pierwiastka zadanego równania. W tym przykładzie odwrotność pierwiastka kwadratowego to:

y=\frac{1}{\sqrt{x}}

lub jako funkcja y

f(y)=\frac{1}{y^2}-x=0

Ponieważ ogólne wyrażenie metody Newtona z \scriptstyle \, y_n jako pierwszym przybliżeniem to

y_{n+1} = y_{n} - \frac{f(y_n)}{f'(y_n)}

po podstawieniu rozpatrywanej funkcji wyrażenie przyjmuje postać szczególną

y_{n+1} = \frac{y_{n}(3-xy_n^2)}{2}

gdzie

f(y)=\frac{1}{y^2}-x i f'(y)=\frac{-2}{y^3}

Stąd

y = y*(1.5f - x2*y*y);

jest tym samym co

y_{n+1} = y_{n}(1,5-\frac{xy_n^2}{2}) = \frac{y_{n}(3-xy_n^2)}{2}

Pierwsze przybliżenie jest wygenerowane powyżej przez operacje całkowitoliczbowe i przekazane do ostatnich dwóch linii funkcji. Powtarzane iteracje w tym algorytmie, używając wyniku funcji (\scriptstyle y_{n+1}) jako argumentu w kolejnej iteracji, powoduje, że wyniki są zbieżne do pierwiastka kwadratowego z rosnącą precyzją[14]. Dla celów silnika Quake III użyto tylko jednej iteracji. Druga iteracja pozostała w kodzie, lecz w formie komentarza[10].

Dokładność[edytuj | edytuj kod]

Wykres pokazujący różnice między heurystyczną szybką odwrotnością pierwiastka kwadratowego a funkcją dostarczaną przez libstdc (uwaga na skale logarytmiczne na obu osiach)

Jak jest to przedstawione wyżej, przybliżenie jest zaskakująco dokładne. Wykres po prawej przedstawia funkcję błędu (tj. błędu po ulepszeniu wyniku dzięki wykonaniu jednej iteracji metody Newtona), dla wartości rozpoczynających się od 0,01, w której biblioteka standardowa zwraca 10,0 jako wynik, podczas gdy InvSqrt() zwraca 9,982522, różniąc się o 0,017479 czyli 0,175%. Jest to największy błąd bezwzględny, który maleje dla rosnących argumentów, natomiast błąd względny pozostaje w tych samych granicach dla wszystkich rzędów wielkości.

Historia i badania[edytuj | edytuj kod]

John Carmack, współzałożyciel id Software, jest często kojarzony z kodem, chociaż to nie on go napisał

Kod źródłowy Quake III nie był upubliczniony aż do QuakeCon 2005, lecz kopie kodu szybkiej odwrotności pierwiastka kwadratowego pojawiły się na Usenet i innych forach już w 2002 lub 2003 roku[1]. Początkowe spekulacje wskazywały na Johna Carmacka jako prawdopodobnego autora kodu, lecz on zaprzeczył i zasugerował, że ten kod napisał Terje Mathisen, utalentowany programista asemblerowy, który wcześniej pomagał id Software przy optymalizacji Quake'a. Mathisen napisał implementację podobnego kodu bitowego w późnych latach dziewięćdziesiątych XX wieku, lecz okazało się, że w historii grafiki komputerowej 3D pierwsi autorzy pojawili się znacznie wcześniej. Najwcześniejszą znaną i używaną implementację napisał Gary Tarolli dla SGI Indigo. Rys Sommefeldt stwierdził, że algorytm został opracowany przez Grega Walsha w Ardent Computer w porozumieniu z Clevem Molerem, twórcą programu MATLAB, chociaż ostateczny dowód autorstwa nie istnieje[15].

Nie wiadomo dokładnie w jaki sposób wyznaczono wartość magicznej stałej. Chris Lomont opracował funkcję do minimalizacji błędu przybliżenia przez wybór magicznej liczby R z przedziału. Najpierw wyznaczył optymalną stałą dla etapu liniowej aproksymacji jako 0x5f37642f, blisko 0x5f3759df, lecz ta nowa stała dawała odrobinę mniej dokładne wyniki po wykonaniu jednej iteracji metodą Newtona[16]. Następnie Lomont poszukiwał optymalnej stałej nawet po jednej i dwóch iteracjach metody Newtona, i znalazł 0x5f375a86, która jest bardziej dokładna niż oryginalna dla każdej iteracji[16]. Stwierdził, pytając, czy dokładna wartość stałej została wybrana przez jej wyprowadzenie czy metodą prób i błędów[17]. Lomont wskazał, że „magiczna liczba” dla 64-bitowej wersji liczby zmiennoprzecinkowej IEEE 754 (podwójna precyzja) to 0x5fe6ec85e7de30da, lecz w rzeczywistości okazało się, że to 0x5fe6eb50c7b537a9[18]. Charles McEniry wykonał podobne, lecz bardziej wyrafinowane optymalizacje wobec prawdopodobnych wartości R. Jego początkowa szukająca metoda czołgowa zwróciła tę samą wartość jaką wyznaczył Lomont[19]. Kiedy próbował znaleźć stałą metodą równego podziału, uzyskał specyficzną wartość R, która występuje w funkcji. Utwierdziło to go w przekonaniu, że stała mogła być oryginalnie wyprowadzona „metodą równego podziału do zadanej dokładności”[20].

Uwagi

  1. Użycie typu long zmniejsza przenośność tego kodu na współczesnych maszynach. Aby kod wykonał się prawidłowo sizeof(long) musi zwracać 4 bajty, inaczej wyniki mogą być ujemne. W przypadku wielu współczesnych systemów 64-bitowych, sizeof(long) wynosi 8 bajtów.
  2. \scriptstyle E_x powinien się zmieścić w zakresie \scriptstyle [1, 254] aby \scriptstyle x był zapisywalny jako tak zwana liczba normalna.
  3. Tylko takie liczby rzeczywiste można dokładnie przedstawić w zapisie zmiennoprzecinkowym dla których \scriptstyle M_x jest liczbą całkowitą. Pozostałe są jedynie przybliżone przez zaokrąglenie ich wartości do najbliższej liczby, która daje się przedstawić w tym formacie.
  4. \scriptstyle S_x = 0 ponieważ \scriptstyle x > 0.

Zobacz też[edytuj | edytuj kod]

Przypisy

Bibliografia[edytuj | edytuj kod]

Linki zewnętrzne[edytuj | edytuj kod]