Rychlá inverzní odmocnina

Aktuální verze stránky ještě nebyla zkontrolována zkušenými přispěvateli a může se výrazně lišit od verze recenzované 17. ledna 2022; kontroly vyžadují 76 úprav .

Rychlá obrácená druhá odmocnina (také rychlá InvSqrt() nebo 0x5F3759DF podle použité "magické" konstanty , v desítkové soustavě 1 597 463 007) je rychlý přibližný algoritmus pro výpočet inverzní odmocniny pro kladná 32bitová čísla s plovoucí desetinnou čárkou . Algoritmus používá celočíselné operace "odčítání" a " bitový posun ", stejně jako zlomkové "odečítání" a "násobení" - bez pomalých operací "rozdělení" a "druhá odmocnina". Navzdory tomu, že je na bitové úrovni „hacky“, je aproximace monotónní a spojitá : blízké argumenty dávají těsné výsledky. Přesnost (méně než 0,2 % dolů a nikdy nahoru) [1] [2] nestačí pro skutečné numerické výpočty, ale pro trojrozměrnou grafiku je dostačující . $y={\frac {1}{{\sqrt {x))))$

Algoritmus

Algoritmus používá jako vstup 32bitové číslo s plovoucí desetinnou čárkou ( s jednoduchou přesností ve formátu IEEE 754 ) a provádí s ním následující operace:

Považujte 32bitové zlomkové číslo za celé číslo a proveďte operaci y 0 = 5F3759DF 16 − (x >> 1) , kde >> je bitový posun doprava. Výsledek je opět považován za 32bitové zlomkové číslo.
Pro upřesnění lze provést jednu iteraci Newtonovy metody : y 1 \u003d y 0 (1,5 - 0,5xy 0 ²) .

Původní algoritmus:

float Q_rsqrt ( float number ) { dlouhé i ; float x2 , y ; const float tři poloviny = 1,5F ; x2 = číslo * 0,5F ; y = číslo ; i = * ( dlouhé * ) & y ; // zlé hackování bitové úrovně s plovoucí desetinnou čárkou i = 0x5f3759df - ( i >> 1 ); // co to sakra? y = * ( float * ) & i ; y = y * ( tři poloviny - ( x2 * y * y ) ); // 1. iterace // y = y * ( tři poloviny - ( x2 * y * y)); // 2. iterace, toto lze odstranit vrátit y ; }

Opravte podle standardů moderní implementace C, s ohledem na možné optimalizace a multiplatformní :

float Q_rsqrt ( float number ) { const float x2 = číslo * 0,5F ; const float tři poloviny = 1,5F ; unie { plovák f ; uint32_t i ; } konv = { číslo }; // člen 'f' nastaven na hodnotu 'číslo'. konv . i = 0x5f3759df - ( konv . i >> 1 ); konv . f *= tři poloviny - x2 * konv . f * konv . f ; zpětná konv . f ; }

Implementace z Quake III: Arena [3] uvažuje, že délka je rovna , a pro převod používá ukazatele (optimalizace „pokud se , nic se nezměnilo“ může fungovat chybně; na GCC se při kompilaci do „release“ zobrazí varování spuštěno). A také obsahuje obscénní slovo - John Carmack , který hru veřejně vystavil, nechápal, co se tam dělá. floatlongfloatlong

V C++20 můžete použít novou funkci . bit_cast

#include <bit> #include <limity> #include <cstdint> constexpr float Q_rsqrt ( float number ) noexcept { static_assert ( std :: numeric_limits < float >:: is_iec559 ); // Zkontrolujte, zda je cílový počítač kompatibilní float const y = std :: bit_cast < float > ( 0x5f3759df - ( std :: bit_cast < std :: uint32_t > ( číslo ) >> 1 )); návrat y * ( 1,5f - ( číslo * 0,5f * y * y )); }

Historie

Algoritmus byl pravděpodobně vyvinut v Silicon Graphics v 90. letech 20. století a implementace se objevila v roce 1999 ve zdrojovém kódu počítačové hry Quake III Arena , ale metoda se na veřejných fórech jako Usenet neobjevila až do roku 2002-2003. Algoritmus generuje přiměřeně přesné výsledky pomocí jedinečné první aproximace Newtonovy metody . V té době bylo hlavní výhodou algoritmu odmítnutí drahých výpočtů s plovoucí desetinnou čárkou ve prospěch celočíselných operací. Inverzní druhé odmocniny se používají k výpočtu úhlů dopadu a odrazu pro osvětlení a stínování v počítačové grafice.

Algoritmus byl původně připisován Johnu Carmackovi , ale ten navrhl, aby jej do id Software přinesli Michael Abrash , specialista na grafiku, nebo Terje Mathisen, specialista na jazyk symbolických instrukcí . Studie problému ukázala, že kód měl hlubší kořeny v hardwarové i softwarové oblasti počítačové grafiky. Opravy a změny provedly jak Silicon Graphics, tak 3dfx Interactive , přičemž nejstarší známou verzi napsal Gary Tarolli pro SGI Indigo . Možná tento algoritmus vynalezli Greg Walsh a Clive Moler , Garyho kolegové z Ardent Computer [5] .

Jim Blinn, specialista na 3D grafiku, navrhl podobnou tabulkovou inverzní odmocninu [6] , která počítá až 4 číslice (0,01 %) a byla nalezena při rozebírání hry Interstate '76 (1997) [7] .

S vydáním 3DNow! Procesory AMD zavedly instrukci assembleru PFRSQRT [8] pro rychlý přibližný výpočet inverzní odmocniny. Verze pro double je nesmyslná - přesnost výpočtů se nezvýší [2] - proto nebyla přidána. V roce 2000 byla do SSE2 přidána funkce RSQRTSS [9] , která je přesnější než tento algoritmus (0,04 % oproti 0,2 %).

Analýza a chyba

Bitová reprezentace 4bajtového zlomkového čísla ve formátu IEEE 754 vypadá takto:

Podepsat
	Objednat								Mantisa
0	0	jeden	jeden	jeden	jeden	jeden	0	0	0	jeden	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				patnáct				osm				7				0

Zabýváme se pouze kladnými čísly (znaménkový bit je nula), nedenormalizovanými , ne ∞ a ne NaN . Taková čísla se zapisují ve standardním tvaru jako 1,mmmm 2 ·2 e . Část 1,mmmm se nazývá mantisa , e je řád . Hlavní jednotka není uložena ( implicitní jednotka ), proto nazvěme hodnotu 0,mmmm explicitní částí mantisy. Strojová zlomková čísla mají navíc posunuté pořadí : 2 0 se zapisuje jako 011.1111.1 2 .

Na kladných číslech je bijekce "zlomek ↔ celé číslo" (níže označená jako ) spojitá jako po částech lineární funkce a monotónní . Z toho můžeme ihned konstatovat, že rychlá inverzní odmocnina je jako kombinace spojitých funkcí spojitá. A jeho první část - posun-odčítání - je také monotónní a po částech lineární. Bijekce je komplikovaná, ale téměř „zdarma“: v závislosti na architektuře procesoru a konvencích volání nemusíte buď nic dělat, nebo přesunout číslo z zlomkového registru na celočíselné. ${\displaystyle I_{x))$

Například binární reprezentace hexadecimálního celého čísla 0x5F3759DF je 0|101.1111.0|011.0111.0101.1001.1101.1111 2 (body jsou nibble hranice, svislé čáry jsou počítačové zlomkové hranice pole). Pořadí 101 1111 0 2 je 190 10 , po odečtení offsetu 127 10 dostaneme exponent 63 10 . Explicitní část mantisy 01 101 110 101 100 111 011 111 2 po přidání implicitní vedoucí jednotky se stane 1,011 011 101 011 001 110 111 11 2 = 1,432 111 430 1 S přihlédnutím ke skutečné přesnosti počítačového zlomku 0x5F3759DF ↔ 1,4324301 10 2 63 .

Označujeme explicitní část mantisy čísla , - neposunutý řád, - bitová délka mantisy, - offset řádu. Číslo , zapsané v lineárně-logaritmické bitové mřížce počítačových zlomků, lze [10] [3] aproximovat logaritmickou mřížkou jako , kde je parametr používaný k úpravě přesnosti aproximace. Tento parametr se pohybuje od 0 (přesný v a ) do 0,086 (přesný v jednom bodě, ) $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})\přibližně e_{x}+m_{x}+\sigma$ $\sigma$ $m_{x}=0$ $jeden$ $m_{x}=0{,}443$

Pomocí této aproximace lze celočíselnou reprezentaci čísla aproximovat jako $X$

I_{x}\equiv L(e_{x}+B+m_{x})\přibližně L\log _{2}x+L(B-\sigma )

V souladu s tím, . $\log _{2}x\cca {\frac {I_{x}}{L}}-(B-\sigma )$

Udělejme totéž [3] pro (respektive ) a získáme ${\displaystyle 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\approx I^{-1}\left[{\tfrac {3}{2}}L(B-\sigma )-{\tfrac {1}{2}}I_{x}\right]

Magická konstanta s přihlédnutím k hranicím bude mít ve zlomkové aritmetice nezaujaté pořadí a mantisu a v binárním zápisu - 0|101.1111.0|01 1 ... (1 je implicitní jednotka; 0,5 pochází z řádu ; malá jednotka odpovídá rozsahu [1,375; 1,5), a je tedy vysoce pravděpodobná, ale není zaručena našimi odhady.) ${\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$

Je možné spočítat, čemu se rovná první dílčí lineární aproximace [11] (ve zdroji není použita mantisa samotná, ale její explicitní část ): $t=c-1$

Pro : ; $x\in [0{,}5;\;c-0{,}5)$ $y_{01}=-x+t+{\tfrac {3}{2}}=-x+c+{\tfrac {1}{2}}$
Pro : ; $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}}$
Pro : . $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}}$

Na větším nebo menším se výsledek mění proporcionálně: při čtyřnásobku je výsledek přesně poloviční. $X$ $X$

Newtonova metoda dává [11] , , a . Funkce je klesající a konvexní směrem dolů, na takových funkcích se Newtonova metoda blíží skutečné hodnotě zleva - proto algoritmus vždy podcení odpověď. $f(y)={\frac {1}{y^{2))}-x$ ${\displaystyle 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)$

Není známo, odkud se vzala konstanta 0x5F3759DF ↔ [a] 1,4324301 2 63 . Hrubou silou Chris Lomont a Matthew Robertson zjistili [1] [2] , že nejlepší konstanta [b] z hlediska mezní relativní chyby pro float je 0x5F375A86 ↔ 1,4324500 2 63 , pro double je to 0x5FE6EB50C7B537A9. Pravda, pro dvojnásobek je algoritmus nesmyslný (nepřináší zisk v přesnosti ve srovnání s plovoucím) [2] . Lomontova konstanta byla získána také analyticky ( c = 1,432450084790142642179 ) [b] , ale výpočty jsou poměrně komplikované [11] [2] .

Po jednom kroku Newtonovy metody je výsledek celkem přesný ( +0 % -0,18 % ) [1] [2] , což je více než vhodné pro účely počítačové grafiky ( 1 ⁄ 256 ≈ 0,39 % ). Taková chyba je zachována v celém rozsahu normalizovaných zlomkových čísel. Dva kroky dávají přesnost 5 číslic [1] , po čtyřech krocích je dosaženo chyby dvojnásobku . Pokud chcete, můžete chybu vyvážit vynásobením koeficientů 1,5 a 0,5 1,0009 tak, aby metoda dávala symetricky ±0,09 % - to udělali ve hře Interstate '76 a Blinnovu metodu, která také iteruje Newtonovu metodu [7 ] .

Newtonova metoda nezaručuje monotónnost, ale počítačový výčet ukazuje, že monotonie existuje.

Zdrojový kód (C++) #include <iostream> union FloatInt { plovák jakoPloat ; int32_t asInt ; }; int floatToInt ( float x ) { FloatInt r ; r . asFloat = x ; vrátit r . asInt ; } float intToFloat ( int x ) { FloatInt r ; r . asInt = x ; vrátit r . asFloat ; } float Q_rsqrt ( float number ) { dlouhé i ; float x2 , y ; const float tři poloviny = 1,5F ; x2 = číslo * 0,5F ; y = číslo ; i = * ( dlouhé * ) & y ; // hacking na úrovni bitů s plovoucí desetinnou čárkou i = 0x5f3759df - ( i >> 1 ); // co to sakra? y = * ( float * ) & i ; // Nevím, co sakra! y = y * ( tři poloviny - ( x2 * y * y ) ); // 1. iterace vrátit y ; } int main () { int iStart = floatToInt ( 1.0 ); int iEnd = floatToInt ( 4.0 ); std :: cout << "Čísla do konce: " << iEnd - iStart << std :: endl ; int nProblémy = 0 ; float oldResult = std :: numeric_limits < float >:: infinity (); for ( int i = iStart ; i <= iEnd ; ++ i ) { float x = intToFloat ( i ); float vysledek = Q_rsqrt ( x ); if ( vysledek > stary vysledek ) { std :: cout << "Nalezen problém na " << x << std :: endl ; ++ nProblémy ; } } std :: cout << "Totální problémy: " << nProblémy << std :: endl ; návrat 0 ; }

Existují podobné algoritmy pro další mocniny, jako je druhá mocnina nebo krychlová odmocnina [3] .

Motivace

"Přímé" překrytí osvětlení na trojrozměrném modelu , dokonce i vysoce poly, i když vezmeme v úvahu Lambertův zákon a další vzorce odrazu a rozptylu, okamžitě poskytne polygonální vzhled - divák uvidí rozdíl v osvětlení podél okraje mnohostěnu [12] . Někdy je to nutné - pokud je objekt skutečně hranatý. A u křivočarých objektů dělají toto: v trojrozměrném programu označují, zda ostrou hranu nebo hladkou [12] . V závislosti na tom, i při exportu modelu do hry, rohy trojúhelníků vypočítají normálu jednotky délky ke zakřivené ploše. Během animace a rotace hra transformuje tyto normály spolu se zbytkem 3D dat; při aplikaci osvětlení se interpoluje přes celý trojúhelník a normalizuje (uvede jej na jednotku délky).

Abychom vektor normalizovali, musíme všechny tři jeho složky vydělit délkou. Nebo, lépe, vynásobte je převrácenou hodnotou délky: . Miliony těchto kořenů se musí vypočítat za sekundu. Než byl vytvořen specializovaný hardware pro transformaci a zpracování osvětlení , mohl být výpočetní software pomalý. Zejména na počátku 90. let, kdy byl kód vyvinut, většina výpočtů s plovoucí desetinnou čárkou zaostávala za celočíselnými operacemi ve výkonu. $(x',y',z')=(x,y,z){\frac {1}{\sqrt {x^{2}+y^{2}+z^{2}}} }$

Quake III Arena používá rychlý inverzní algoritmus druhé odmocniny pro urychlení grafického zpracování CPU, ale tento algoritmus byl od té doby implementován v některých specializovaných hardwarových vertex shaderech pomocí vyhrazených programovatelných polí (FPGA).

I na počítačích z 10. let 20. století může být rychlost v závislosti na zatížení frakčního koprocesoru třikrát až čtyřikrát vyšší než při použití standardních funkcí [11] .

Komentáře

↑ Šipka ↔ zde znamená bijekci binární reprezentace celého čísla a binární reprezentace čísla s plovoucí desetinnou čárkou ve formátu IEEE 754 vysvětleném výše .
↑ 1 2 Pokud do pole objednávky vložíte 127, dostanete 0x3FB75A86. Knihovna GRISU2, která je zcela celočíselná a nezávislá na jemnosti koprocesoru, říká, že 0x3FB75A86 ↔ 1,43245 je nejkratší desetinné číslo, které se převádí na daný float . Nejméně významnou jednotkou je však stále 1,19 10 −7 a 0x3FB75A86 = 1,432450056 ≈ 1,4324501. Další zlomek je 0x3FB75A87 ↔ 1,4324502 bez jakýchkoliv jemností. Proto neintuitivní zaokrouhlení 1,43245008 na 1,4324500.

Poznámky

↑ 1 2 3 4 Archivovaná kopie . Získáno 25. srpna 2019. Archivováno z originálu 6. února 2009. (neurčitý)
↑ 1 2 3 4 5 6 https://web.archive.org/web/20140202234227/http://shelfflag.com/rsqrt.pdf
↑ 1 2 3 4 Hummus a magnety . Staženo 1. února 2017. Archivováno z originálu 13. ledna 2017. (neurčitý)
↑ Beyond3D – Původ rychlé InvSqrt() od Quake3 . Získáno 4. října 2019. Archivováno z originálu 10. dubna 2017. (neurčitý)
↑ Beyond3D – Původ rychlé InvSqrt() od Quake3 – část druhá . Získáno 25. srpna 2019. Archivováno z originálu dne 25. srpna 2019. (neurčitý)
↑ Triky s plovoucí desetinnou čárkou | IEEE Journals & Magazines | IEEE Xplore
↑ 1 2 Rychlá reciproční druhá odmocnina… v roce 1997?! — Blog Shanea Peelara
↑ PFRSQRT - Vypočítejte přibližnou hodnotu převrácené hodnoty druhé odmocniny krátké reálné hodnoty - Club155.ru . Získáno 4. října 2019. Archivováno z originálu dne 16. října 2019. (neurčitý)
↑ RSQRTSS - Výpočet převrácené hodnoty druhé odmocniny ze skalární s jednoduchou přesností s plovoucí desetinnou čárkou . Získáno 6. října 2019. Archivováno z originálu 12. srpna 2019. (neurčitý)
↑ https://web.archive.org/web/20150511044204/http://www.daxia.com/bibis/upload/406Fast_Inverse_Square_Root.pdf
↑ 1 2 3 4 Schwidkeův výpočet zalomené odmocniny s magickými konstantami - analytický pidkhid . Získáno 12. června 2022. Archivováno z originálu dne 17. dubna 2022. (neurčitý)
↑ 1 2 3 Toto je norma: co jsou normální mapy a jak fungují / Sudo Null IT News

Odkazy

C. Lomont, Rychlá inverzní odmocnina , Technická zpráva, 2003.
Stručná historie InvSqrt od Matthew Robertson
0x5f3759df , další zkoumání přesnosti a zobecnitelnosti algoritmu Christian Plesner Hansen