Gyors inverz négyzetgyök

Az oldal jelenlegi verzióját még nem ellenőrizték tapasztalt közreműködők, és jelentősen eltérhet a 2022. január 17-én felülvizsgált verziótól ; az ellenőrzések 76 szerkesztést igényelnek .

A Fast Reverse Square Root (más néven gyors InvSqrt() vagy 0x5F3759DF a használt "mágikus" állandóval , decimálisan 1 597 463 007) egy gyors közelítő algoritmus a pozitív 32 bites lebegőpontos számok fordított négyzetgyökének kiszámítására . Az algoritmus egész számokat használ "kivonás" és " biteltolás ", valamint tört "kivonás" és "szorzás" - lassú "osztás" és "négyzetgyök" műveletek nélkül. Annak ellenére, hogy bitszinten "hacky", a közelítés monoton és folyamatos : a közeli argumentumok szoros eredményeket adnak. A pontosság (kevesebb, mint 0,2% lefelé és soha fel) [1] [2] nem elég valós numerikus számításokhoz, de háromdimenziós grafikákhoz bőven elég . $y={\frac {1}{{\sqrt {x))))$

Algoritmus

Az algoritmus egy 32 bites lebegőpontos számot vesz ( egyetlen pontosságú IEEE 754 formátumban ) bemenetként, és a következő műveleteket hajtja végre rajta:

Ha egy 32 bites tört számot egész számként kezel, hajtsa végre az y 0 = 5F3759DF 16 − (x >> 1) műveletet , ahol a >> egy biteltolódás jobbra. Az eredményt ismét 32 bites törtszámként kezeljük.
Az egyértelműség kedvéért a Newton-módszer egy iterációja végrehajtható : y 1 \u003d y 0 (1,5 - 0,5xy 0 ²) .

Eredeti algoritmus:

float Q_rsqrt ( lebegő szám ) { hosszú i ; úszó x2 , y ; const float threehalfs = 1,5F ; x2 = szám * 0,5F ; y = szám ; i = * ( hosszú * ) & y ; // gonosz lebegőpontos bitszintű hackelés i = 0x5f3759df - ( i >> 1 ); // mi a fasz? y = * ( lebegő * ) & i ; y = y * ( háromfél - ( x2 * y * y ) ); // 1. iteráció // y = y * ( háromfél - ( x2 * y * y ) ); // 2. iteráció, ez eltávolítható vissza y ; }

Korrekt a modern C implementáció szabványai szerint, figyelembe véve a lehetséges optimalizációkat és a platformok közötti átlépést :

float Q_rsqrt ( lebegő szám ) { const float x2 = szám * 0,5F ; const float threehalfs = 1,5F ; szakszervezet { úszó f ; uint32_t i ; } konv = { szám }; // az 'f' tag a 'szám' értékére van állítva. konv . i = 0x5f3759df - ( konv . i >> 1 ); konv . f *= háromfél - x2 * konv . f * konv . f ; visszatérési konv . f ; }

A Quake III: Arena [3] implementációja úgy ítéli meg, hogy a hossza egyenlő , és mutatókat használ a konverzióhoz (az „ha , egyik sem változott” optimalizálás hibásan működhet; GCC-n a „kiadásra” fordításkor figyelmeztetés jelenik meg kiváltva). És egy obszcén szót is tartalmaz - John Carmack , aki nyilvánosságra hozta a játékot, nem értette, mit csinálnak ott. floatlongfloatlong

A C++20 nyelven az új függvény használható . bit_cast

#include <bit> #include <limits> #include <cstdint> constexpr float Q_rsqrt ( float number ) nokivéve { static_assert ( std :: numeric_limits < float >:: is_iec559 ); // Ellenőrizze, hogy a célgép kompatibilis-e float const y = std :: bit_cast < float > ( 0x5f3759df - ( std :: bit_cast < std :: uint32_t > ( szám ) >> 1 )); return y * ( 1.5f - ( szám * 0.5f * y * y )); }

Történelem

Az algoritmust valószínűleg a Silicon Graphics fejlesztette ki a kilencvenes években, és 1999-ben jelent meg egy implementáció a Quake III Arena számítógépes játék forráskódjában , de a módszer csak 2002-2003 között jelent meg nyilvános fórumokon, például a Usenetben . Az algoritmus ésszerűen pontos eredményeket generál a Newton-módszer egyedi első közelítésével . Akkoriban az algoritmus fő előnye a drága lebegőpontos számítások elutasítása volt az egész műveletek javára. Inverz négyzetgyököket használnak a beesési és visszaverődési szögek kiszámítására a számítógépes grafika megvilágításához és árnyékolásához.

Az algoritmust eredetileg John Carmacknek tulajdonították , de ő azt javasolta, hogy Michael Abrash grafikus vagy Terje Mathisen assembly nyelvi specialista vigye el az id Software -hez . A probléma tanulmányozása kimutatta, hogy a kódnak mélyebb gyökerei vannak a számítógépes grafika hardver és szoftver területén egyaránt. Javításokat és változtatásokat hajtott végre a Silicon Graphics és a 3dfx Interactive is, a legkorábbi ismert verziót Gary Tarolli írta az SGI Indigo számára . Az algoritmust talán Greg Walsh és Clive Moler találták ki , Gary munkatársai az Ardent Computernél [5] .

Jim Blinn, egy 3D-s grafika specialistája egy hasonló táblázatos inverz négyzetgyök módszert javasolt [6] , amely 4 számjegyig (0,01%) számol, és az Interstate '76 (1997) [7] szétszedésekor találták meg .

A 3DNow megjelenésével! Az AMD processzorok bevezették a PFRSQRT [8] assembler utasítást az inverz négyzetgyök gyors közelítő kiszámításához. A dupla verzió értelmetlen - a számítások pontossága nem nő [2] - ezért nem került bele. 2000-ben az RSQRTSS függvényt [9] hozzáadták az SSE2 -hez, amely pontosabb ennél az algoritmusnál (0,04% versus 0,2%).

Elemzés és hiba

Egy 4 bájtos tört szám bitreprezentációja IEEE 754 formátumban így néz ki:

Jel
	Rendelés								Mantissa
0	0	egy	egy	egy	egy	egy	0	0	0	egy	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				tizenöt				nyolc				7				0

Csak pozitív számokkal foglalkozunk (az előjelbit nulla), denormalizált nem , nem ∞ és nem NaN . Az ilyen számokat szabványos formában írjuk fel : 1,mmmm 2 ·2 e . Az 1,mmmm részt mantisszának nevezzük , e a sorrend . A fejegység nincs tárolva ( implicit unit ), ezért nevezzük a 0,mmmm értéket a mantissza explicit részének. Ezenkívül a gépi törtszámok eltolt sorrendűek : 2 0 011.1111.1 2 -ként van írva .

Pozitív számokon a "tört ↔ egész" bijekció (a továbbiakban: ) folytonos, mint darabonként lineáris függvény és monoton . Ebből rögtön kijelenthetjük, hogy a gyors inverz gyök, mint folytonos függvények kombinációja folytonos. És az első része - eltolás-kivonás - szintén monoton és darabonként lineáris. A bijekció bonyolult, de szinte „ingyenes”: a processzor architektúrától és a hívási konvencióktól függően vagy semmit sem kell tennie, vagy át kell helyeznie a számot a tört regiszterből az egész számba. $I_{x}$

Például a 0x5F3759DF hexadecimális egész szám bináris megjelenítése 0|101.1111.0|011.0111.0101.1001.1101.1111 2 (a pontok nibble-határok, a függőleges vonalak a számítógépes tört mezők). A 101 1111 0 2 sorrendje 190 10 , a 127 10 eltolás levonása után a 63 10 kitevőt kapjuk . A mantissza 01 101 110 101 100 111 011 111 2 explicit része az implicit vezető egység hozzáadása után 1,011 011 101 011 001 110 111 11 2 = 1,4314 4 … 3014 4 Figyelembe véve a számítógépes törtszám valós pontosságát 0x5F3759DF ↔ 1.4324301 10 2 63 .

Jelöljük a szám mantissza explicit részét , - eltolatlan sorrendet, - a mantissza bithosszát, - a sorrend eltolását. A számítógépes törtek lineáris-logaritmikus bitrácsába írt szám [10] [3] logaritmikus ráccsal közelíthető , mint , ahol a közelítés pontosságának beállítására szolgáló paraméter. Ez a paraméter 0 (pontosan és ) és 0,086 (pontosan egy ponton, ) között mozog. $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})\approx e_{x}+m_{x}+\sigma$ $\sigma$ $m_{x}=0$ $egy$ $m_{x}=0{,}443$

Ezzel a közelítéssel egy szám egész reprezentációja a következőképpen közelíthető meg $x$

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

Ennek megfelelően, . $\log _{2}x\approx {\frac {I_{x}}{L}}-(B-\sigma )$

Tegyük meg ugyanezt [3] -ra (illetve ) és kapjuk meg ${\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]

A mágikus konstansnak , figyelembe véve a határokat , a tört aritmetikában torzítatlan sorrendje és mantisszája lesz , bináris jelölésben pedig - 0|101.1111.0|01 1 ... (1 implicit egység; 0,5 a sorrendből származik ; egy kis egység megfelel az [1,375; 1,5) tartománynak, ezért nagyon valószínű, de becsléseink nem garantálják.) ${\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$

Kiszámolható, hogy az első darabonkénti lineáris közelítés [11] mennyivel egyenlő (a forrásban nem magát a mantisszát, hanem annak explicit részét használjuk ): $t=c-1$

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

Nagyobbnál vagy kisebbnél az eredmény arányosan változik: négyszerezve az eredmény pontosan megfeleződik. $x$ $x$

Newton módszere [11] , , és . A függvény csökkenő és lefelé konvex, ezeken a függvényeken Newton módszere balról közelíti meg a valódi értéket – ezért az algoritmus mindig alulbecsüli a választ. $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)$

Nem ismert, honnan származik a 0x5F3759DF ↔ [a] 1.4324301 2 63 konstans . Brutális erővel Chris Lomont és Matthew Robertson rájöttek [1] [2] , hogy a legjobb állandó [b] a relatív határhibát tekintve lebegésnél 0x5F375A86 ↔ 1,4324500 2 63 , duplánál pedig 0x5FE6EB50C7B537A9. Igaz, duplán az algoritmus értelmetlen (nem ad pontosságnövekedést a floathoz képest) [2] . A Lomont-állandót analitikusan is megkaptuk ( c = 1,432450084790142642179 ) [b] , de a számítások meglehetősen bonyolultak [11] [2] .

Newton módszerének egy lépése után az eredmény meglehetősen pontos ( +0 % -0,18 % ) [1] [2] , ami több mint alkalmas számítógépes grafikai célokra ( 1 ⁄ 256 ≈ 0,39 % ). Ez a hiba a normalizált törtszámok teljes tartományában megmarad. Két lépés 5 számjegy pontosságot ad [1] , négy lépés után elérjük a dupla hibát . Ha szükséges, a hibát újra kiegyenlítheti az 1,5 és 0,5 együtthatók 1,0009-cel való megszorzásával, így a módszer szimmetrikusan ±0,09%-ot ad – ezt tették az Interstate '76 játékban és a Blinn metódusban, amelyek szintén iterálják Newton módszerét [7 ] .

Newton módszere nem garantálja a monotonitást, de a számítógépes felsorolás azt mutatja, hogy létezik monotonitás.

Forráskód (C++) #include <iostream> union FloatInt { float asFloat ; int32_t asInt ; }; int floatToInt ( float x ) { FloatInt r ; r . asFloat = x ; visszatér r . asInt ; } float intToFloat ( int x ) { FloatInt r ; r . asInt = x ; visszatér r . asFloat ; } float Q_rsqrt ( lebegő szám ) { hosszú i ; úszó x2 , y ; const float threehalfs = 1,5F ; x2 = szám * 0,5F ; y = szám ; i = * ( hosszú * ) & y ; // gonosz lebegőpontos bitszintű hackelés i = 0x5f3759df - ( i >> 1 ); // mi a fasz? y = * ( lebegő * ) & i ; // nem tudom mi a fasz! y = y * ( háromfél - ( x2 * y * y ) ); // 1. iteráció vissza y ; } int main () { int iStart = floatToInt ( 1.0 ); int iEnd = floatToInt ( 4.0 ); std :: cout << "Meg kell hagyni a számokat: " << iEnd - iStart << std :: endl ; int nProblémák = 0 ; float oldResult = std :: numeric_limits < float >:: végtelen (); for ( int i = iStart ; i <= iEnd ; ++ i ) { float x = intToFloat ( i ); float result = Q_rsqrt ( x ); if ( eredmény > oldResult ) { std :: cout << "Hibát találtunk a következőn: " << x << std :: endl ; ++ nProblémák ; } } std :: cout << "Összes probléma: " << nProblémák << std :: endl ; return 0 ; }

Hasonló algoritmusok léteznek más hatványokhoz is, például négyzet- vagy kockagyökhöz [3] .

Motiváció

A világítás „közvetlen” átfedése egy háromdimenziós modellen , még egy nagy pólusú modellen is, még a Lambert-törvényt és más tükrözési és szórási képleteket is figyelembe véve , azonnal sokszögletű megjelenést kölcsönöz – a néző látni fogja a világítás különbségét a poliéder élei [12] . Néha ez szükséges - ha az objektum valóban szögletes. A görbe vonalú objektumok esetében pedig ezt teszik: egy háromdimenziós programban jelzik, hogy éles vagy sima élről van szó [12] . Ennek függvényében a modell játékba való exportálásakor is a háromszögek sarkai az egységnyi hosszúságú normált számítják ki az íves felületre. Az animáció és a forgatás során a játék ezeket a normál értékeket a többi 3D adattal együtt átalakítja; világítás alkalmazásakor a teljes háromszögre interpolál és normalizál (egységnyi hosszra hozza).

Egy vektor normalizálásához el kell osztanunk mindhárom komponensét a hosszával. Vagy jobb, ha megszorozzuk őket a hossz reciprokával: . Ezeknek a gyökereknek a millióit kell másodpercenként számolni. A dedikált transzformációs és világítás -feldolgozó hardver létrehozása előtt a számítási szoftverek lassúak voltak. Különösen az 1990-es évek elején, amikor a kódot kifejlesztették, a legtöbb lebegőpontos számítás teljesítményében elmaradt az egész műveletekhez képest. $(x',y',z')=(x,y,z){\frac {1}{\sqrt {x^{2}+y^{2}+z^{2}}} }$

A Quake III Arena gyors inverz négyzetgyök algoritmust használ a CPU grafikus feldolgozásának felgyorsítására, de az algoritmust azóta néhány speciális hardveres vertex shaderben implementálták dedikált programozható tömbök (FPGA) segítségével.

Még a 2010-es évek számítógépein is, a frakcionált koprocesszor terhelésétől függően , a sebesség három-négyszer nagyobb lehet, mint a szabványos funkciók használata [11] .

Megjegyzések

↑ Itt a ↔ nyíl egy egész szám bináris reprezentációjának és egy lebegőpontos szám bináris ábrázolásának bijekcióját jelenti a fentebb ismertetett IEEE 754 formátumban .
↑ 1 2 Ha 127-et ír be a rendelési mezőbe, akkor 0x3FB75A86-ot kap. A GRISU2 könyvtár, amely teljesen egész szám, és független a társprocesszor finomságaitól, azt mondja, hogy a 0x3FB75A86 ↔ 1,43245 a legrövidebb decimális szám, amely egy adott lebegőponttá konvertálódik . A legkisebb szignifikáns egység azonban továbbra is 1,19 10 −7 , és 0x3FB75A86 = 1,432450056 ≈ 1,4324501. A következő tört a 0x3FB75A87 ↔ 1.4324502, minden finomság nélkül. Ezért az 1,43245008 nem intuitív kerekítése 1,4324500-ra.

Jegyzetek

↑ 1 2 3 4 Archív másolat . Letöltve: 2019. augusztus 25. Az eredetiből archiválva : 2009. február 6.. (határozatlan)
↑ 1 2 3 4 5 6 https://web.archive.org/web/20140202234227/http://shelfflag.com/rsqrt.pdf
↑ 1 2 3 4 Hummusz és mágnesek . Letöltve: 2017. február 1. Az eredetiből archiválva : 2017. január 13. (határozatlan)
↑ Beyond3D – A Quake3 Fast InvSqrt() eredete . Letöltve: 2019. október 4. Az eredetiből archiválva : 2017. április 10. (határozatlan)
↑ Beyond3D – A Quake3 Fast InvSqrt() eredete – Második rész . Letöltve: 2019. augusztus 25. Az eredetiből archiválva : 2019. augusztus 25. (határozatlan)
↑ Lebegőpontos trükkök | IEEE folyóiratok és magazinok | IEEE Xplore
↑ 1 2 Gyors reciprok négyzetgyök… 1997-ben?! - Shane Peelar blogja
↑ PFRSQRT - Számítsa ki egy rövid valós érték négyzetgyökének reciprokának közelítő értékét - Club155.ru . Letöltve: 2019. október 4. Az eredetiből archiválva : 2019. október 16. (határozatlan)
↑ RSQRTSS – Számítsa ki a skaláris egypontosságú lebegőpontos érték négyzetgyökének reciprokát . Letöltve: 2019. október 6. Az eredetiből archiválva : 2019. augusztus 12. (határozatlan)
↑ https://web.archive.org/web/20150511044204/http://www.daxia.com/bibis/upload/406Fast_Inverse_Square_Root.pdf
↑ 1 2 3 4 Schwidke számítása a burkolt négyzetgyökről a mágikus állandókkal - analitikus pidkhid . Letöltve: 2022. június 12. Az eredetiből archiválva : 2022. április 17. (határozatlan)
↑ 1 2 3 Ez a norma: mik a normál térképek és hogyan működnek / Sudo Null IT News

Linkek

C. Lomont, Gyors inverz négyzetgyök , Technical Report, 2003.
Az InvSqrt rövid története Matthew Robertsontól
0x5f3759df , Christian Plesner Hansen további vizsgálatai az algoritmus pontosságáról és általánosíthatóságáról