Rădăcină pătrată inversă rapidă

Versiunea actuală a paginii nu a fost încă examinată de colaboratori experimentați și poate diferi semnificativ de versiunea revizuită la 17 ianuarie 2022; verificările necesită 76 de modificări .

Rădăcină pătrată inversă rapidă (de asemenea rapid InvSqrt() sau 0x5F3759DF după constanta „magică” folosită , în zecimală 1 597 463 007) este un algoritm rapid aproximativ pentru calcularea rădăcinii pătrate inverse pentru numerele pozitive cu virgulă mobilă de 32 de biți . Algoritmul folosește operații cu numere întregi „scădere” și „ deplasare de biți ”, precum și „scădere” și „înmulțire” fracțională - fără operații lente „împărțire” și „rădăcină pătrată”. În ciuda faptului că este „hacky” la nivel de biți, aproximarea este monotonă și continuă : argumentele apropiate dau rezultate apropiate. Precizia (mai puțin de 0,2% în jos și niciodată în sus) [1] [2] nu este suficientă pentru calcule numerice reale, dar este destul de suficientă pentru grafica tridimensională . $y={\frac {1}{{\sqrt {x))))$

Algoritm

Algoritmul ia ca intrare un număr în virgulă mobilă de 32 de biți ( precizie unică în format IEEE 754 ) și efectuează următoarele operații pe acesta:

Tratând un număr fracționar de 32 de biți ca un întreg, efectuați operația y 0 = 5F3759DF 16 − (x >> 1) , unde >> este o deplasare de biți la dreapta. Rezultatul este din nou tratat ca un număr fracționar de 32 de biți.
Pentru clarificare, poate fi efectuată o iterație a metodei lui Newton : y 1 \u003d y 0 (1,5 - 0,5xy 0 ²) .

Algoritmul original:

float Q_rsqrt ( număr float ) { i lung ; float x2 , y ; const float threehalfs = 1.5F ; x2 = număr * 0,5F ; y = număr ; i = * ( lung * ) & y ; // hacking evil flating point la nivel de biți i = 0x5f3759df - ( i >> 1 ); // ce dracu? y = * ( float * ) & i ; y = y * ( trei jumătăți - ( x2 * y * y ) ); // Prima iterație // y = y * ( trei jumătăți - ( x2 * y * y ) ); // A doua iterație, aceasta poate fi eliminată returnează y ; }

Corect după standardele implementării moderne C, ținând cont de posibilele optimizări și cross-platform :

float Q_rsqrt ( număr float ) { const float x2 = număr * 0.5F ; const float threehalfs = 1.5F ; unire { plutire f ; uint32_t i ; } conv = { număr }; // membru „f” setat la valoarea „număr”. conv . i = 0x5f3759df - ( conv . i >> 1 ); conv . f *= trei jumătăți - x2 * conv . f * conv . f ; returnare conv . f ; }

Implementarea de la Quake III: Arena [3] consideră că lungimea este egală cu , și folosește pointeri pentru conversie (optimizarea „dacă , nu s-a schimbat nimic” poate funcționa eronat; pe GCC, la compilare pentru „lansare”, este un avertisment declanșat). Și conține, de asemenea, un cuvânt obscen - John Carmack , așezând jocul în domeniul public, nu a înțeles ce se face acolo. floatlongfloatlong

În C++20 , puteți utiliza noua funcție . bit_cast

#include <bit> #include <limite> #include <cstdint> constexpr float Q_rsqrt ( număr float ) noexcept { static_assert ( std :: numeric_limits < float >:: is_iec559 ); // Verificați dacă mașina țintă este compatibilă float const y = std :: bit_cast < float > ( 0x5f3759df - ( std :: bit_cast < std :: uint32_t > ( număr ) >> 1 )); returnează y * ( 1,5f - ( număr * 0,5f * y * y )); }

Istorie

Algoritmul a fost probabil dezvoltat la Silicon Graphics în anii 1990, iar o implementare a apărut în 1999 în codul sursă pentru jocul de calculator Quake III Arena , dar metoda nu a apărut pe forumuri publice precum Usenet până în 2002-2003. Algoritmul generează rezultate rezonabil de precise folosind prima aproximare unică a metodei lui Newton . La acea vreme, principalul avantaj al algoritmului era respingerea calculelor costisitoare în virgulă mobilă în favoarea operațiilor cu numere întregi. Rădăcinile pătrate inverse sunt folosite pentru a calcula unghiurile de incidență și reflexie pentru iluminare și umbrire în grafica computerizată.

Algoritmul a fost atribuit inițial lui John Carmack , dar el a sugerat că Michael Abrash , specialist în grafică, sau Terje Mathisen, specialist în limbaj de asamblare, l- au adus la id Software . Studiul problemei a arătat că codul are rădăcini mai profunde atât în domeniul hardware, cât și în cel software al graficii computerizate. Corecțiile și modificările au fost făcute atât de Silicon Graphics, cât și de 3dfx Interactive , cea mai veche versiune cunoscută fiind scrisă de Gary Tarolli pentru SGI Indigo . Poate că algoritmul a fost inventat de Greg Walsh și Clive Moler , colegii lui Gary de la Ardent Computer [5] .

Jim Blinn, specialist în grafică 3D, a propus o metodă similară tabelară a rădăcinii pătrate inverse [6] care numără până la 4 cifre (0,01%) și a fost găsită la dezasamblarea jocului Interstate '76 (1997) [7] .

Odată cu lansarea 3DNow! Procesoarele AMD au introdus instrucțiunea de asamblare PFRSQRT [8] pentru calcularea rapidă aproximativă a rădăcinii pătrate inverse. Versiunea pentru dublu este lipsită de sens - acuratețea calculelor nu va crește [2] - prin urmare nu a fost adăugată. În 2000, funcția RSQRTSS [9] a fost adăugată la SSE2 , care este mai precisă decât acest algoritm (0,04% față de 0,2%).

Analiză și eroare

Reprezentarea pe biți a unui număr fracționar de 4 octeți în format IEEE 754 arată astfel:

Semn
	Ordin								mantisa
0	0	unu	unu	unu	unu	unu	0	0	0	unu	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	$=(1+2^{-2})\cdot 2^{-3}=0{,}15625$
31				24				23				16				cincisprezece				opt				7				0

Avem de-a face doar cu numere pozitive (bitul de semn este zero), nu denormalizate , nu ∞ și nu NaN . Astfel de numere sunt scrise în formă standard ca 1,mmmm 2 ·2 e . Partea 1,mmmm se numește mantisa , e este ordinul . Unitatea principală nu este stocată ( unitatea implicită ), așa că să numim valoarea 0,mmmm partea explicită a mantisei. În plus, numerele fracționale ale mașinii au o ordine deplasată : 2 0 se scrie ca 011.1111.1 2 .

Pe numerele pozitive , bijecția „fracțional ↔ întreg” (notat mai jos ca ) este continuă ca funcție liniară pe bucăți și monotonă . Din aceasta putem afirma imediat că rădăcina inversă rapidă, ca o combinație de funcții continue, este continuă. Și prima sa parte - deplasare-scădere - este, de asemenea, monotonă și liniară pe bucăți. Bijecția este complicată, dar aproape „gratuită”: în funcție de arhitectura procesorului și de convențiile de apelare , fie trebuie să nu faci nimic, fie să muți numărul din registrul fracționar în cel întreg. $I_{x}$

De exemplu, reprezentarea binară a întregului hexazecimal 0x5F3759DF este 0|101.1111.0|011.0111.0101.1001.1101.1111 2 (punctele sunt granițe de ronțăit, liniile verticale sunt limite de câmp fracțional de computer). Ordinul lui 101 1111 0 2 este 190 10 , după scăderea decalajului 127 10 obținem exponentul 63 10 . Partea explicită a mantisei 01 101 110 101 100 111 011 111 2 după adăugarea unității conducătoare implicite devine 1.011 011 101 011 001 110 111 11 2 = 1.432 8… 4301 4 . Ținând cont de precizia reală a fracționalului computerizat 0x5F3759DF ↔ 1,4324301 10 2 63 .

Notăm partea explicită a mantisei a numărului , - ordinul neschimbat, - lungimea biților a mantisei, - offset al ordinului. Numărul , scris într-o rețea de biți liniar-logaritmică a fracțiilor computerizate, poate fi [10] [3] aproximat printr-o rețea logaritmică ca , unde este un parametru utilizat pentru a ajusta acuratețea aproximării. Acest parametru variază de la 0 (exact la și ) la 0,086 (exact la un punct, ) $m_{x}\in [0,1)$ $X$ $e_{x}\in \mathbb {Z}$ $L=2^{23}$ $B=127$ $x\equiv 2^{e_{x}}(1+m_{x})$ $\log _{2}x\equiv e_{x}+\log _{2}(1+m_{x})\aprox e_{x}+m_{x}+\sigma$ $\sigma$ $m_{x}=0$ $unu$ $m_{x}=0{,}443$

Folosind această aproximare, reprezentarea întregului număr poate fi aproximată ca $X$

I_{x}\equiv L(e_{x}+B+m_{x})\aprox L\log _{2}x+L(B-\sigma)

În consecință, . $\log _{2}x\aprox {\frac {I_{x}}{L}}-(B-\sigma )$

Să facem același lucru [3] pentru (respectiv ) și să obținem $y={\tfrac {1}{\sqrt {x}))$ $\log _{2}y=-{\tfrac {1}{2}}\log _{2}x$

I_{y}\approx {\tfrac {3}{2}}L(B-\sigma )-{\tfrac {1}{2}}I_{x}

y\aprox I^{-1}\left[{\tfrac {3}{2}}L(B-\sigma )-{\tfrac {1}{2}}I_{x}\right]

Constanta magică , ținând cont de granițe , în aritmetică fracțională va avea o ordine imparțială și mantisă ), iar în notație binară - 0|101.1111.0|01 1 ... (1 este o unitate implicită; 0,5 a venit din ordinea ; o unitate mică corespunde intervalului [1,375; 1,5) și, prin urmare, este foarte probabilă, dar nu este garantată de estimările noastre.) ${\tfrac {3}{2}}L(B-\sigma )$ $\sigma$ $\left\lfloor {\tfrac {3L}{2L}}(B-\sigma )\right\rfloor -B=63$ $c=1+0{,}5-{\tfrac {3}{2}}\sigma \in (1{,}37;1{,}5$

Este posibil să se calculeze cu ce este egală prima aproximare liniară pe bucăți [11] (în sursă, nu mantisa în sine, dar este folosită partea sa explicită ): $t=c-1$

Pentru : ; $x\in [0{,}5;\;c-0{,}5)$ $y_{01}=-x+t+{\tfrac {3}{2}}=-x+c+{\tfrac {1}{2}}$
Pentru : ; $x\in [c-0{,}5;\;1)$ $y_{02}=-{\tfrac {1}{2}}x+{\tfrac {1}{2}}t+{\tfrac {5}{4}}=-{\tfrac {1}{} 2}}x+{\tfrac {1}{2}}c+{\tfrac {3}{4}}$
Pentru : . $x\in [1;\;2)$ $y_{03}=-{\tfrac {1}{4}}x+{\tfrac {1}{2}}t+1=-{\tfrac {1}{4}}x+{\tfrac { 1}{2}}c+{\tfrac {1}{2}}$

La mai mare sau mai mică , rezultatul se modifică proporțional: atunci când este de patru ori , rezultatul este exact înjumătățit. $X$ $X$

Metoda lui Newton dă [11] , , și . Funcția este descrescătoare și convexă în jos; pe astfel de funcții, metoda lui Newton se apropie de valoarea adevărată din stânga - prin urmare, algoritmul subestimează întotdeauna răspunsul. $f(y)={\frac {1}{y^{2)}}-x$ $f'(y)=-{\frac {2}{y^{3))}$ $y_{n+1}=y_{n}-{\frac {f(y_{n})}{f'(y_{n))}}={\frac {y_{n}(3- xy_{n}^{2})}{2}}=y_{n}(1{,}5-0{,}5\,xy_{n}^{2})$ $f(y)$

Nu se știe de unde provine constanta 0x5F3759DF ↔ [a] 1,4324301 2 63 . Prin forța brută, Chris Lomont și Matthew Robertson au descoperit [1] [2] că cea mai bună constantă [b] în termeni de eroare relativă marginală pentru float este 0x5F375A86 ↔ 1.4324500 2 63 , pentru dublu este 0x5FE6EB50C7B537A9. Adevărat, pentru dublu algoritmul este lipsit de sens (nu oferă un câștig în precizie în comparație cu float) [2] . Constanta Lomont a fost obținută și analitic ( c = 1,432450084790142642179 ) [b] , dar calculele sunt destul de complicate [11] [2] .

După un pas al metodei lui Newton, rezultatul este destul de precis ( +0 % -0,18 % ) [1] [2] , ceea ce este mai mult decât potrivit pentru scopuri de grafică pe computer ( 1 ⁄ 256 ≈ 0,39 % ). O astfel de eroare este păstrată pe întregul interval de numere fracționale normalizate. Doi pași dau o precizie de 5 cifre [1] , după patru pași se ajunge la eroarea de dublu . Dacă doriți, puteți reechilibra eroarea înmulțind coeficienții 1,5 și 0,5 cu 1,0009, astfel încât metoda să dea simetric ±0,09% - asta au făcut în jocul Interstate '76 și metoda Blinn, care repetă și metoda lui Newton [7] ] .

Metoda lui Newton nu garantează monotonitatea, dar enumerarea computerizată arată că monotonitatea există.

Cod sursă (C++) #include <iostream> unire FloatInt { float asFloat ; int32_t asInt ; }; int floatToInt ( float x ) { FloatInt r ; r . asFloat = x ; întoarce r . asInt ; } float intToFloat ( int x ) { FloatInt r ; r . asInt = x ; întoarce r . asFloat ; } float Q_rsqrt ( număr float ) { i lung ; float x2 , y ; const float threehalfs = 1.5F ; x2 = număr * 0,5F ; y = număr ; i = * ( lung * ) & y ; // evil flating point bit level hacking i = 0x5f3759df - ( i >> 1 ); // ce dracu? y = * ( float * ) & i ; // nu stiu ce naiba! y = y * ( trei jumătăți - ( x2 * y * y ) ); // Prima iterație returnează y ; } int main () { int iStart = floatToInt ( 1.0 ); int iEnd = floatToInt ( 4.0 ); std :: cout << "Numerele de parcurs: " << iEnd - iStart << std :: endl ; int nProbleme = 0 ; float oldResult = std :: numeric_limits < float >:: infinity (); pentru ( int i = iStart ; i <= iEnd ; ++ i ) { float x = intToFloat ( i ); rezultat float = Q_rsqrt ( x ); if ( rezultat > oldResult ) { std :: cout << "Am găsit o problemă pe " << x << std :: endl ; ++ nProbleme ; } } std :: cout << "Total probleme: " << nProblems << std :: endl ; returnează 0 ; }

Există algoritmi similari pentru alte puteri, cum ar fi rădăcina pătrată sau cubă [3] .

Motivație

Suprapunerea „directă” a luminii pe un model tridimensional , chiar și unul high-poli, chiar și luând în considerare legea lui Lambert și alte formule de reflexie și împrăștiere, va da imediat un aspect poligonal - privitorul va vedea diferența de iluminare de-a lungul muchiile poliedrului [12] . Uneori, acest lucru este necesar - dacă obiectul este cu adevărat unghiular. Și pentru obiectele curbilinie, ele fac acest lucru: într-un program tridimensional, ele indică dacă o muchie ascuțită sau una netedă [12] . În funcție de aceasta, chiar și atunci când exportați modelul în joc, colțurile triunghiurilor calculează normala lungimii unității la suprafața curbată. În timpul animației și rotației, jocul transformă aceste valori normale împreună cu restul datelor 3D; la aplicarea luminii, se interpolează pe întregul triunghi și se normalizează (o aduce la o unitate de lungime).

Pentru a normaliza un vector, trebuie să împărțim toate cele trei componente ale acestuia la lungime. Sau, mai bine, înmulțiți-le cu reciproca lungimii: . Milioane din aceste rădăcini trebuie calculate pe secundă. Înainte ca hardware-ul dedicat transformării și procesării luminii să fie creat , software-ul de calcul ar putea fi lent. În special, la începutul anilor 1990, când a fost dezvoltat codul, majoritatea calculelor în virgulă mobilă au rămas în urma operațiilor cu numere întregi în performanță. $(x',y',z')=(x,y,z){\frac {1}{\sqrt {x^{2}+y^{2}+z^{2}}} }$

Quake III Arena folosește un algoritm rapid de rădăcină pătrată inversă pentru a accelera procesarea grafică de către CPU, dar algoritmul a fost implementat de atunci în unele vertex shadere hardware specializate folosind matrice programabile dedicate (FPGA).

Chiar și pe computerele anilor 2010, în funcție de sarcina coprocesorului fracționat , viteza poate fi de trei până la patru ori mai mare decât utilizarea funcțiilor standard [11] .

Comentarii

↑ Aici săgeata ↔ înseamnă bijecția reprezentării binare a unui număr întreg și reprezentarea binară a unui număr în virgulă mobilă în formatul IEEE 754 explicat mai sus .
↑ 1 2 Dacă puneți 127 în câmpul de comandă, obțineți 0x3FB75A86. Biblioteca GRISU2, care este în întregime întreagă și independentă de subtilitățile coprocesorului, spune că 0x3FB75A86 ↔ 1,43245 este cel mai scurt număr zecimal care se convertește într-un anumit float . Cu toate acestea, unitatea cea mai puțin semnificativă este încă 1,19 10 −7 și 0x3FB75A86 = 1,432450056 ≈ 1,4324501. Următorul fracționar este 0x3FB75A87 ↔ 1,4324502 fără subtilități. De aici și rotunjirea neintuitivă a 1,43245008 la 1,4324500.

Note

↑ 1 2 3 4 Copie arhivată . Preluat la 25 august 2019. Arhivat din original pe 6 februarie 2009. (nedefinit)
↑ 1 2 3 4 5 6 https://web.archive.org/web/20140202234227/http://shelfflag.com/rsqrt.pdf
↑ 1 2 3 4 Hummus și magneți . Preluat la 1 februarie 2017. Arhivat din original la 13 ianuarie 2017. (nedefinit)
↑ Beyond3D - Originea lui Quake3 Fast InvSqrt() . Consultat la 4 octombrie 2019. Arhivat din original la 10 aprilie 2017. (nedefinit)
↑ Beyond3D - Originea lui Quake3's Fast InvSqrt() - Partea a doua . Preluat la 25 august 2019. Arhivat din original la 25 august 2019. (nedefinit)
↑ Trucuri cu virgulă mobilă | Reviste și reviste IEEE | IEEE Explore
↑ 1 2 Rădăcină pătrată reciprocă rapidă... în 1997?! — Blogul lui Shane Peelar
↑ PFRSQRT - Calculați valoarea aproximativă a reciprocei rădăcinii pătrate a unei valori reale scurte - Club155.ru . Preluat la 4 octombrie 2019. Arhivat din original la 16 octombrie 2019. (nedefinit)
↑ RSQRTSS - Calculează reciprocul rădăcinii pătrate a valorii scalare cu precizie unică în virgulă mobilă . Preluat la 6 octombrie 2019. Arhivat din original la 12 august 2019. (nedefinit)
↑ https://web.archive.org/web/20150511044204/http://www.daxia.com/bibis/upload/406Fast_Inverse_Square_Root.pdf
↑ 1 2 3 4 Calculul lui Schwidke al rădăcinii pătrate înfășurate cu constantele magice - pidkhid analitic . Preluat la 12 iunie 2022. Arhivat din original la 17 aprilie 2022. (nedefinit)
↑ 1 2 3 Aceasta este norma: ce sunt hărțile normale și cum funcționează / Sudo Null IT News

Link -uri

C. Lomont, Rădăcină pătrată inversă rapidă , Raport tehnic, 2003.
O scurtă istorie a InvSqrt de Matthew Robertson
0x5f3759df , investigații suplimentare privind acuratețea și generalizarea algoritmului de către Christian Plesner Hansen