Szybka odwrotność pierwiastka kwadratowego

Szybka odwrotność pierwiastka kwadratowego, czasami określana jako szybki InvSqrt() lub przez stałą szesnastkową 0x5f3759df – metoda obliczania ${\displaystyle 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 90. XX wieku, a jedna z jego implementacji 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ę.

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 90. 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:

${\displaystyle \|{\boldsymbol {v}}\|={\sqrt {v_{1}^{2}+v_{2}^{2}+v_{3}^{2}}},}$

a wektor unormowany (jednostkowy):

${\displaystyle {\boldsymbol {\hat {v}}}={\boldsymbol {v}}/\|{\boldsymbol {v}}\|.}$

Zastępując wyrażenie ${\displaystyle v_{1}^{2}+v_{2}^{2}+v_{3}^{2}}$ przez ${\displaystyle {\boldsymbol {x}},}$ otrzymujemy związek:

${\displaystyle {\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].

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 ${\displaystyle 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 90. 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 unikatowe 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] .

Rozważmy przykład z liczbą ${\displaystyle x=0{,}15625}$ dla której chcemy wyznaczyć wielkość ${\displaystyle 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 ${\displaystyle y=2{,}61486,}$ które różni się od dokładnego wyniku o 3,4%. Po jednej iteracji metodą Newtona ostateczny wynik wynosi ${\displaystyle y=2{,}52549}$ z uchybem tylko 0,17%.

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

1. interpretuje parametr ${\displaystyle x}$ jako liczbę całkowitą, co może być sposobem na wyznaczenie przybliżenia ${\displaystyle \log _{2}x,}$
2. przekształca wyliczone przybliżenie aby wyznaczyć przybliżoną wartość ${\displaystyle \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.

 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 ${\displaystyle x}$ w postaci zmiennoprzecinkowej o pojedynczej precyzji pierwszym krokiem jest przedstawienie ${\displaystyle x}$ w postaci liczby znormalizowanej:

${\displaystyle {\begin{aligned}x&=\pm 1.b_{1}b_{2}b_{3}\ldots \times 2^{e_{x}}\\&=\pm 2^{e_{x}}(1+m_{x}),\end{aligned}}}$

gdzie wykładnik ${\displaystyle e_{x}}$ jest liczbą całkowitą, ${\displaystyle m_{x}\in [0,1),}$ a ${\displaystyle 1.b_{1}b_{2}b_{3}\ldots }$ binarną wartością mantysy ${\displaystyle (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:

• ${\displaystyle S_{x}}$ – bit znaku, który przyjmuje wartość 0 jeśli ${\displaystyle x>0}$ lub 1 jeśli ${\displaystyle x<0}$ (1 bit)
• ${\displaystyle E_{x}=e_{x}+B}$ – wykładnik ze składową stałą, w którym ${\displaystyle B=127}$^[b] (8 bitów)
• ${\displaystyle M_{x}=m_{x}\times L,}$ gdzie ${\displaystyle 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 ${\displaystyle x=0{,}15625=0{,}00101_{2}.}$ Jej postać znormalizowana to:

${\displaystyle x=+2^{-3}(1+0{,}25),}$

dla której odpowiednie trzy pola bitowe wynoszą

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

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

Chcąc wyznaczyć ${\displaystyle 1/{\sqrt {x}}}$ bez komputera lub kalkulatora, można skorzystać z tablic logarytmów i tożsamości ${\displaystyle \log _{b}(1/{\sqrt {x}})=-{\frac {1}{2}}\log _{b}(x),}$ która zachodzi dla każdej podstawy ${\displaystyle 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 ${\displaystyle x}$ jest normalną liczbą dodatnią:

${\displaystyle x=2^{e_{x}}(1+m_{x}),}$

to zachodzi

${\displaystyle \log _{2}(x)=e_{x}+\log _{2}(1+m_{x}),}$

lecz skoro ${\displaystyle m_{x}\in [0,1),}$ to logarytm po prawej stronie równości można przybliżyć^[12]

${\displaystyle \log _{2}(1+m_{x})\approx m_{x}+\sigma ,}$

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

Stąd otrzymane przybliżenie to

${\displaystyle \log _{2}(x)\approx e_{x}+m_{x}+\sigma .}$

Z drugiej strony interpretacja bitowej reprezentacji ${\displaystyle x}$ jako liczby całkowitej daje^[d]

${\displaystyle {\begin{aligned}I_{x}&=E_{x}L+M_{x}\\&=L(e_{x}+B+m_{x})\\&\approx L\log _{2}(x)+L(B-\sigma ).\end{aligned}}}$

Okazuje się, że ${\displaystyle I_{x}}$ to przeskalowane i przesunięte liniowe przybliżenie ${\displaystyle \log _{2}(x),}$ przedstawione na ilustracji po prawej. Innymi słowy ${\displaystyle \log _{2}(x)}$ jest przybliżony przez

${\displaystyle \log _{2}(x)\approx {\frac {I_{x}}{L}}-(B-\sigma ).}$

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

${\displaystyle \log _{2}(y)=-{\frac {1}{2}}\log _{2}(x).}$

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

${\displaystyle 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

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

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

 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:

${\displaystyle y={\frac {1}{\sqrt {x}}}.}$

Wybieramy funkcję ${\displaystyle f(z),}$ dla której ${\displaystyle f(y)=0}$

${\displaystyle f(z)={\frac {1}{z^{2}}}-x.}$

Ponieważ ogólne wyrażenie metody Newtona z ${\displaystyle y_{n}}$ jako pierwszym przybliżeniem, to

${\displaystyle y_{n+1}=y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}.}$

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

${\displaystyle y_{n+1}={\frac {y_{n}(3-xy_{n}^{2})}{2}},}$

gdzie:

${\displaystyle f(z)={\frac {1}{z^{2}}}-x}$ i ${\displaystyle f'(z)={\frac {-2}{z^{3}}}.}$

Stąd

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

jest tym samym co

${\displaystyle 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 funkcji ${\displaystyle (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].

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.

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 90. 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].

• metody obliczania pierwiastka kwadratowego

• JimJ. Blinn JimJ., Jim Blinn’s Corner: Notation, notation notation, Morgan Kaufmann, 2003, ISBN 1-55860-860-5 .
• DavidD. Eberly DavidD., 3D Game Engine Design, Morgan Kaufmann, 2001, ISBN 978-1-55860-593-0 .
• DavidD. Eberly DavidD., Fast Inverse Square Root [online], Geometric Tools, 2002 [dostęp 2012-08-19] [zarchiwizowane z adresu 2009-02-24]  (ang.).
• JohnJ. Hennessey JohnJ., David A.D.A. Patterson David A.D.A., Computer Organization and Design, wyd. 2nd, San Francisco, CA: Morgan Kaufmann Publishers, 1998, ISBN 978-1-55860-491-9  (ang.).
• ChrisCh. Lomont ChrisCh., Fast Inverse Square Root [online], 2003 [dostęp 2012-08-18]  (ang.).
• CharlesCh. McEniry CharlesCh., The Mathematics Behind the Fast Inverse Square Root Function Code [online], 2007 [dostęp 2009-02-13] [zarchiwizowane z adresu 2015-05-11]  (ang.).
• LarsL. Middendorf LarsL. i inni, Embedded Vertex Shader in FPGA, [w:] AchinA. Rettberg, Embedded System Design: Topics, Techniques and Trends. International Federation for Information Processing (IFIP) TC10 Working Conference:International Embedded Systems Symposium (IESS), Irvine (Kalifornia): Springer, 1 czerwca 2007, ISBN 978-0-387-72257-3 .
• MatthewM. Robertson MatthewM., A Brief History of InvSqrt [online], UNBSJ, 24 kwietnia 2012 [dostęp 2012-08-22] [zarchiwizowane z adresu 2013-03-09]  (ang.).
• ElanE. Ruskin ElanE., Timing square root, [w:] Some Assembly Required [online], 16 października 2009 [dostęp 2010-09-13] [zarchiwizowane z adresu 2014-07-25] .
• RysR. Sommefeldt RysR., Origin of Quake3's Fast InvSqrt() [online], Beyond3D, 29 listopada 2006a [dostęp 2012-08-17]  (ang.).
• RysR. Sommefeldt RysR., Origin of Quake3's Fast InvSqrt() - Part Two [online], Beyond3D, 19 grudnia 2006b [dostęp 2012-08-22]  (ang.).

• Tomer Margolin: Magical Square Root Implementation In Quake III. [w:] CodeMaestro [on-line]. The Coding Experience, 2005-08-27. [dostęp 2012-08-23]. [zarchiwizowane z tego adresu (2012-04-14)]. (ang.).