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 wywodziła się od unikania 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ą światło na temat możliwych metod.

Uzasadnianie[edytuj | edytuj kod]

Pole wektorów normalnych
Wektory normalne są szeroko wykorzystywane przy obliczaniu oświetlenia i cieniowania, i 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, miliony takich obliczeń muszą być wykonane na sekundę. Zanim wprowadzono moduły sprzętowe do transformacji i oświetlenia, obliczenia w oprogramowaniu mogły być wolne. Zwłaszcza że na początku lat dziewięćdziesiątych XX wieku, kiedy kod był rozwijany, większość zmiennoprzecinkowej mocy obliczeniowej nie nadążała za prędkością obliczeń całkowitoliczbowych[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ę tę długość, to nowy wektor będzie wektorem jednostkowym o takim samym zwrocie.

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

\|\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 relację

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

która łączy wektor jednostkowy z odwrotnością pierwiastka kwadratowego ze składowych odległości.

Quake III Arena wykorzystuje algorytm szybkiej odwrotności pierwiastka kwadratowego aby przyspieszyć obliczenia procesora graficznego, lecz od tego czasu ten algorytm jest już realizowany 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. Powszechne 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 zrekompensowana 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 tej stałej, więc, jak wielu podobnych przypadkach stosowania stałych w kodzie, nazywa się je „liczbami magicznymi”[1][8][9][10]. To całkowite 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 on znacznie wolniejszy i mniej dokładny niż instrukcja SSE rsqrtss na procesorach x86, wydanych także w 1999 roku[11].

Rzutowanie z liczby zmiennoprzecinkowej na całkowitą i odwrotnie[edytuj | edytuj kod]

Przykład liczby zmiennoprzecinkowej.svg

Zrozumienie tego wymaga znajomości w jaki sposób liczba zmiennoprzecinkowa jest przechowywana. Liczba zmiennoprzecinkowa przedstawia liczbę wymierną wyrazoną za pomocą trzech części rozmieszczonych w 32 bitach. Na przykład liczba dziesiętna 0,15625 w zapisie dwójkowym to 0,00101, która następnie jest normalizowana dzięki przemnożeniu przez odpowiednią potęgę liczby 2, tak aby na początku uzyskać bit o wartości 1, stąd otrzymuje się zapis 1,01×2-3, o którym można powiedzieć, że wykładnik to -3 a mantysa to 1,01. Ponieważ wszystkie znormalizowane liczby mają na pierwszej cyfrze mantysy wartość 1, to w tym formacie jest on domniemany i pomijany, a zapisywana jest jedynie część mantysy po przecinku M = ,01. Reprezentacja wewnętrzna ma jeszcze jeden bit, który służy do przechowywania znaku liczby i jest on pierwszym bitem w zapisie. Jego wartość to 0 dla liczb dodatnich i 1 dla ujemnych. Następne 8 bitów zajmuje wykładnik, który jest zapisany ze składową stałą[b], aby zapisać wartości z przedziału od -128 do 127 (jeśli jest odczytywany jako jako 8-bitowa liczba ze znakiem). Wykładnik o wartości -3 jest zapisany jako E = 124 (= -3 + 127). Cyfry znaczące mantysy obejmują pozostałe 23 bity zapisanej liczby. Liczby wyrażone są więc jako \scriptstyle x = (-1)^{Znak}(1 + M)2^{E - B}, gdzie \scriptstyle B=127[12]. Stąd, dla przykładu, \scriptstyle 0.15625_{10} = (-1)^{0}(1 + .01_2)2^{124 - 127}. Wartość wykładnika określa również, czy mantysa zawiera ułamek, czy liczbę całkowitą[13].

bit znaku
0 1 1 1 1 1 1 1 = 127
0 0 0 0 0 0 1 0 = 2
0 0 0 0 0 0 0 1 = 1
0 0 0 0 0 0 0 0 = 0
1 1 1 1 1 1 1 1 = −1
1 1 1 1 1 1 1 0 = −2
1 0 0 0 0 0 0 1 = −127
1 0 0 0 0 0 0 0 = −128

8-bitowe liczby całkowite ze znakiem

Dodatnia liczba całkowita ze znakiem w systemie uzupełnień do dwóch ma pierwszy (najstarszy) bit równy 0, po którym następują bity z jej binarną wartością. Rzutowanie liczby zmiennoprzecinkowej na liczbę całkowitą ze znakiem zwraca wartość \scriptstyle I = E*2^{23} + M, gdzie I to uzyskana liczba całkowita, E to wykładnik w postaci wewnętrznej (z dodaną składową stałą i zinterpretowany jako 8-bitowa liczba bez znaku), a M to mantysa w postaci wewnętrznej (z pominiętym pierwszym bitem równym 1) zinterpretowana jako liczba całkowita a nie ułamek. Ponieważ odwrotność pierwiastka kwadratowego dotyczy tylko liczb dodatnich, bit znaku w liczbie zmiennoprzecinkowej jest zawsze równy 0, i liczba całkowita uzyskana przez rzutowanie jest również dodatnia.

Należy podkreślić, że domniemany bit o wartości 1 jak również pole wykładnika powoduje, że rzutowanie na liczbę całkowitą zwraca wyniki znacznie różniące się od typowej konwersji typu zmiennoprzecinkowego do całkowitego. Dla rozpatrywanego przykładu, E = 124 a M jako liczba całkowita staje się 221 ponieważ jej jedyny „ustawiony” bit jest na pozycji 21. Oznacza to, że jej binarna forma to 01000000000000000000000, a nie ,01000000000000000000000. Stąd \scriptstyle I = 124*2^{23}+ 2^{21}, które równa się 1 040 187 392 + 2 097 152 lub 1 042 284 544. Ta „konwersja” lub raczej reinterpretacja odbywa się natychmiast, bez potrzeby skomplikowanych korekt, których wymagałby niezerowy wykładnik. Następnie wykonywane jest kilka nieskomplikowanych obliczeniowo operacji aby uzyskać początkowe oszacowanie, które zostanie ulepszone przez zastosowanie metody Newtona.

Pierwsza prosta operacja, przesunięcie bitowe o 1 bit w prawo, dzieli liczbę całkowitą przez 2[14], której wynik jest następnie odejmowany od magicznej stałej, a ostateczny rezultat rzutowany ponownie na typ zmiennoprzecinkowy, aby stać się początkowym oszacowaniem. Należy zauważyć, że to przesunięcie dzieli przez dwa wartość pola, które stanie się ponownie mantysą (lecz nie wartość mantysy, ponieważ domniemana a pominięta wartość 1 pozostaje; 1,(x) staje się 1,(x/2) a nie (1,x)/2, wartość pola, które stanie się wykładnikiem (interpretując je jako wartość bez znaku), jak również przesuwa najmniej znaczący bit wykładnika z bitu 23 (o wartości zero w przykładzie) na najstarszy bit w polu mantysy (bit 22). Domniemany ustawiony bit nie bierze udziału w tych przekształceniach, to jest:

Zeeeeeeee(1)mmmmmmmmmmmmmmmmmmmmmmm    1 Znak, 8 bitów wykładnika, (bit domniemany), 23 bity mantysy.
001111100   01000000000000000000000    1 042 284 544 Wartość z przykładu, i
000111110   00100000000000000000000      521 142 272 Przesunięcie w prawo o 1.   i/2
010111110   01101110101100111011111    1 597 463 007 Magiczna stała = 5F3759DF
010000000   01001110101100111011111    1 076 320 735 Oszacowany wynik 5F3759DF - i/2

Uzyskany całkowity wynik jest rzutowany na typ zmiennoprzecinkowy: bit znaku jest niezmieniony i wyzerowany, pole wykładnika zawiera E = 128 a pole mantysy to ,01001110101100111011111 binarnie czyli 2578911/223 lub 0,3074301481247, skąd wynikiem jest

 y = (1 + ,3074301481247)2^{128 - 127} = 2,614

podczas gdy rzeczywista wartość to

\frac{1}{\sqrt{0.15625}} = 2,5298221281347

czyli jest to dobre pierwsze przybliżenie. Zastosowanie jednej iteracji metody Newtona zwraca 2,5255 a drugiej 2,529811. Trzecia iteracja mogłaby wyczerpać precyzję zmiennej.

Jeśli wartość początkowa byłaby podwojona, 0,3125, mantysa byłaby taka sama ale wykładnik byłby większy o jeden, stąd na bit na pozycji 23 byłby ustawiony, a połowienie wartości całkowitej ustawiłoby najstarszy bit mantysy na pozycji 22, i po uwzględnieniu domyślnego bitu 1 wynikiem byłoby 1,807 jako początkowe oszacowanie dla 1,788854; pierwsza iteracja metody Newtona ulepszyłaby go do 1,788564.

Magiczna liczba[edytuj | edytuj kod]

Z(znak) E(wykładnik) M(mantysa)
1 bit b bitów (n-1-b) bitów
n bitów[15]

Wybór stałej 0x5f3759df wywołał wiele spekulacji na temat szybkiej odwrotności pierwiastka kwadratowego. W celu określenia jak programista mógł pierwotnie uzyskać wartość tej stałej jako mechanizm do oszacowania odwrotności pierwiastka kwadratowego, Charles McEniry najpierw zbadał jak wybór dowolnej stałej R może dać pierwsze oszacowanie odwrotności pierwiastka kwadratowego. Powołując się na porównanie typów całkowitego i zmiennoprzecinkowego przedstawione wyżej, zauważył, że \scriptstyle x, zadana liczba zmiennoprzecinkowa, to \scriptstyle x=(1+m_x)2^{e_x}, a \scriptstyle I_x, uzyskana liczba całkowita z liczby zmiennoprzecinkowej, to \scriptstyle I_x=E_xL+M_x[c]. Te tożsamości wprowadzają kilka nowych elementów, które są prostym przekształceniem wartości wykładnika i mantysy.

\textstyle m_x=\frac{M_x}{L} \quad gdzie \quad L=2^{n-1-b}
e_x=E_x-B \quad gdzie \quad B=2^{b-1}-1

Ilustracja z McEniry 2007 wyprowadza:

y=\frac{1}{\sqrt{x}}
\log_2{(y)}=-\frac{1}{2}\log_2{(x)}
\log_2(1+m_y)+e_y=-\frac{1}{2}\log_2{(1+m_x)}-\frac{1}{2}e_x

biorąc logarytm binarny po obu stronach. Logarytm binarny to funkcja odwrotna dla \scriptstyle f(n)=2^n redukująca mnożenia w reprezentacji zmiennoprzecinkowej do dodawania. Zależność między \scriptstyle \log_2{(x)} i \scriptstyle \log_2{(x^{-1/2})} jest liniowa, co pozwala na utworzenie równania wyrażającego x and y0 (argument i pierwsze oszacowanie) jako kombinację liniową wyrazów[12]. McEniry wprowadza wyraz \scriptstyle \sigma, który odpowiada za przybliżenie dla \scriptstyle \log_2{(1+x)} w pośrednim kroku do przybliżenia R[d]. Ponieważ \scriptstyle 0\le x< 1, \scriptstyle \log_2{(1+x)}\approx {x}, można teraz zdefiniować \scriptstyle \log_2{(1+x)}\cong x+\sigma. Ta definicja daje pierwsze przybliżenie logarytmu binarnego. Dla bieżącego celu, \scriptstyle \sigma jest liczbą rzeczywistą ograniczoną do [0, 1/3], dla R równego 0x5f3759df, \scriptstyle \sigma=0.0450461875791687011756[15].

m_y+\sigma+e_y=-\frac{1}{2}m_x-\frac{1}{2}\sigma-\frac{1}{2}e_x

Stosując tożsamości dla \scriptstyle M_x, \scriptstyle E_x, \scriptstyle B i \scriptstyle L:

M_y+(E_y-B)L=-\frac{3}{2}\sigma{L}-\frac{1}{2}M_x-\frac{1}{2}(E_x-B)L

Przekształcając wyrazy:

E_yL+M_y=\frac{3}{2}(B-\sigma)L-\frac{1}{2}(E_xL+M_x)

Wartość całkowita ze ściśle dodatniej liczby zmiennoprzecinkowej x to \scriptstyle I_x=E_xL+M_x. To daje wyrażenie na całkowitą wartość z y (gdzie \scriptstyle \textstyle y=\frac{1}{\sqrt{x}}, to pierwsze przybliżenie dla odwrotności pierwiastka kwadratowego) za pomocą całkowitych komponentów z x. W szczególności

I_y=E_yL+M_y=R-\frac{1}{2}(E_xL+M_x)=R-\frac{1}{2}I_x \quad gdzie \quad R=\frac{3}{2}(B-\sigma)L

Równanie \scriptstyle I_y=R-\frac{1}{2}I_x to linia i = 0x5f3759df - (i>>1); w szybkim InvSqrt(), przybliżenie całkowite dla \scriptstyle y_0 to całkowita wartość z x, przesunięta w prawo i odjęta od R[15]. Dowód tutaj zaprezentowany tylko pokazuje, że stałą R można wykorzystać do przybliżenia wartości całkowitej odwrotności pierwiastka kwadratowego z liczby zmiennoprzecinkowej. Nie znaczy to, że R przyjmuje wartość taką jaką użyto w kodzie.

Względnie proste wyjaśnienie jak przesunięcie bitowe i operacja odejmowania używając oczekiwanej wartości R skutkuje podzieleniem wykładnika E przez minus dwa można znaleźć w pracy Chrisa Lomonta traktującej o tej stałej. Na przykład dla 10000=104, dzielenie wykładnika przez −2 zwróciłoby 10000-1/2=10-2=1/100. Ponieważ wykładnik jest obarczony składową stałą, rzeczywista wartość wykładnika (tutaj e) to e=E-127, dzięki czemu wartość wyniku ze składową stałą to -e/2+127[16]. Odjęcie liczby całkowitej od R („magicznej” wartości 0x5f3759df) zmusza najmniej znaczący bit wykładnika do przeniesienia na mantysę i po powrocie do zapisu zmiennoprzecinkowego, zwraca liczbę zmiennoprzecinkową bardzo bliską odwrotności pierwiastka kwadratowego z argumentu. Określona wartość R została tak dobrana, aby zminimalizować oczekiwany błąd w dzieleniu wykładnika jak również oczekiwany błąd w przesuwaniu mantysy. 0xbe przedstawia liczbę całkowitą, która minimalizuje oczekiwany błąd wynikający z dzielenia wykładnika liczby zmiennoprzecinkowej przez przesunięcie bitowe — w szczególności wartość 0xbe przesunięta w prawo o jeden bit to 0x5f, pierwsze cyfry magicznej stałe R[17].

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.

Metoda Newtona[edytuj | edytuj kod]

 Osobny artykuł: Metoda Newtona.

Po wykonaniu operacji całkowitoliczbowych, algorytm ponownie traktuje długie słowo jako liczbę zmiennoprzecinkową (x = *(float*)&i;) i wykonuje operację mnożenia zmiennoprzecinkowego (x = x*(1.5f - xhalf*x*x);). 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

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

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ą[18]. Dla celów silnika Quake III użyto tylko jednej iteracji. Druga iteracja pozostała w kodzie, lecz w formie komentarza[10].

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[19].

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[20]. 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[20]. Stwierdził, pytając, czy dokładna wartość stałej została wybrana przez jej wyprowadzenie czy metodą prób i błędów[21]. 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[22]. 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[23]. 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”[24].

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ć negatywne. W przypadku wielu współczesnych systemów 64-bitowych, sizeof(long) wynosi 8 bajtów.
  2. Wartości wykładnika zapisane są w kodzie uzupełnień do dwóch, co można interpretować jako kod z przesunięciem. W przypadku zmiennoprzecinkowych liczb 32-bitowych wartość przesunięcia wynosi 127. Dzięki takiemu zabiegowi, wykładnik w zapisie binarnym 000..000, oznacza zero lub liczbę zdenormalizowaną, a 111...111 nieskończoność lub NaN
  3. Liczby zmiennoprzecinkowe są znormalizowane—mantysa jest wyrażona jako \scriptstyle m_x\in[0,1). Więcej szczegółów w David Goldberg. What Every Computer Scientist Should Know About Floating-Point Arithmetic. „ACM Computing Surveys”. 23 (1), s. 5–48, marzec 1991. doi:10.1145/103162.103163 (ang.). .
  4. Lomont 2003 podchodzi do określenia R w inny sposób, dzieląc R na \scriptstyle R_1 i \scriptstyle R_2 dla bitów mantysy i wykładnika z R.

Zobacz też[edytuj | edytuj kod]

Przypisy

Bibliografia[edytuj | edytuj kod]

Linki zewnętrzne[edytuj | edytuj kod]