Snabb invers kvadratrot

Den aktuella versionen av sidan har ännu inte granskats av erfarna bidragsgivare och kan skilja sig väsentligt från versionen som granskades den 17 januari 2022; kontroller kräver 76 redigeringar .

Fast Reverse Square Root (även snabb InvSqrt() eller 0x5F3759DF med den "magiska" konstanten som används , i decimal 1 597 463 007) är en snabb approximativ algoritm för att beräkna den inversa kvadratroten för positiva 32-bitars flyttal . Algoritmen använder heltalsoperationer "subtrahera" och " bitskift ", samt bråktal "subtrahera" och "multiplicera" - utan långsamma operationer "divide" och "kvadratrot". Trots att den är "hackig" på bitnivå är approximationen monoton och kontinuerlig $y={\frac {1}{{\sqrt {x))))$ : nära argument ger nära resultat. Noggrannhet (mindre än 0,2 % ner och aldrig upp) [1] [2] räcker inte för riktiga numeriska beräkningar, men det är tillräckligt för tredimensionell grafik .

Algoritm

Algoritmen tar ett 32-bitars flyttal ( enkel precision i IEEE 754 -format ) som indata och utför följande operationer på det:

Behandla ett 32-bitars bråktal som ett heltal, utför operationen y 0 = 5F3759DF 16 − (x >> 1) , där >> är en bitförskjutning åt höger. Resultatet behandlas återigen som ett 32-bitars bråktal.
För förtydligande kan en iteration av Newtons metod utföras : y 1 \u003d y 0 (1,5 - 0,5 xy 0 ²) .

Ursprunglig algoritm:

float Q_rsqrt ( flytnummer ) _ { lång i ; flyta x2 , y ; const float threehalvs = 1,5F ; x2 = antal * 0,5F ; y = nummer ; i = * ( lång * ) & y ; // ond floating point bit level hacking i = 0x5f3759df - ( i >> 1 ); // vad fan? y = * ( flytande * ) & i ; y = y * ( trehalvor - ( x2 * y * y ) ); // Första iterationen // y = y * (trehalvor - (x2 * y * y)); // 2:a iterationen, denna kan tas bort returnera y ; }

Korrigera enligt standarderna för modern C-implementering, med hänsyn till möjliga optimeringar och plattformsoberoende :

float Q_rsqrt ( flytnummer ) _ { const float x2 = antal * 0,5F ; const float threehalvs = 1,5F ; union { flyta f ; uint32_t i ; } konv = { nummer }; // medlem 'f' satt till värdet på 'nummer'. konv . i = 0x5f3759df - ( konv . i >> 1 ); konv . f *= trehalvor - x2 * omv . f * konv . f ; returnera konv . f ; }

Implementeringen från Quake III: Arena [3] anser att längden är lika med , och använder pekare för konvertering (optimeringen "om , ingen har ändrats" kan fungera felaktigt; på GCC, vid kompilering för att "släppa", är en varning upprörd). Och det innehåller också ett obscent ord - John Carmack , som lade ut spelet i det offentliga området, förstod inte vad som gjordes där. floatlongfloatlong

I C++20 kan du använda den nya funktionen . bit_cast

#inkludera <bit> #inkludera <gränser> #include <cstdint> constexpr float Q_rsqrt ( flytnummer ) noexcept _ { static_assert ( std :: numeric_limits < float >:: is_iec559 ); // Kontrollera om målmaskinen är kompatibel float const y = std :: bit_cast < float > ( 0x5f3759df - ( std :: bit_cast < std :: uint32_t > ( nummer ) >> 1 )); return y * ( 1,5f - ( nummer * 0,5f * y * y )); }

Historik

Algoritmen utvecklades troligen på Silicon Graphics på 1990-talet, och en implementering dök upp 1999 i källkoden till dataspelet Quake III Arena , men metoden dök inte upp på offentliga forum som Usenet förrän 2002-2003. Algoritmen genererar någorlunda exakta resultat genom att använda den unika första approximationen av Newtons metod . På den tiden var den största fördelen med algoritmen avvisandet av dyra flyttalsberäkningar till förmån för heltalsoperationer. Omvända kvadratrötter används för att beräkna infallsvinklar och reflektion för belysning och skuggning i datorgrafik.

Algoritmen tillskrevs ursprungligen John Carmack , men han föreslog att Michael Abrash , en grafikspecialist, eller Terje Mathisen, en assemblerspråkspecialist, tog den till id Software . Studien av problemet visade att koden hade djupare rötter i både hårdvaru- och mjukvaruområdena för datorgrafik. Korrigeringar och ändringar har gjorts av både Silicon Graphics och 3dfx Interactive , med den tidigaste kända versionen skriven av Gary Tarolli för SGI Indigo . Kanske uppfanns algoritmen av Greg Walsh och Clive Moler , Garys kollegor på Ardent Computer [5] .

Jim Blinn, en 3D-grafikspecialist, har föreslagit en liknande tabellformad invers kvadratrotmetod [6] som räknar upp till 4 siffror (0,01%) och hittades när spelet Interstate '76 (1997) plockades isär [7] .

Med lanseringen av 3DNow! AMD - processorer introducerade PFRSQRT [8] assemblerinstruktionen för snabb ungefärlig beräkning av den inversa kvadratroten. Versionen för dubbel är meningslös - noggrannheten i beräkningarna kommer inte att öka [2] - därför lades den inte till. År 2000 lades RSQRTSS-funktionen [9] till SSE2 , vilket är mer exakt än denna algoritm (0,04 % mot 0,2 %).

Analys och fel

Bitrepresentationen av ett 4-byte bråktal i IEEE 754 -format ser ut så här:

Tecken
	Ordning								Mantissa
0	0	ett	ett	ett	ett	ett	0	0	0	ett	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				femton				åtta				7				0

Vi behandlar bara positiva tal (teckenbiten är noll), inte denormaliserade , inte ∞ och inte NaN . Sådana nummer skrivs i standardform som 1,mmmm 2 ·2 e . Delen 1,mmmm kallas mantissan , e är ordningen . Huvudenheten lagras inte ( implicit enhet ), så låt oss kalla värdet 0,mmmm den explicita delen av mantissan. Dessutom har maskinbråktal en förskjuten ordning : 2 0 skrivs som 011.1111.1 2 .

På positiva tal är "bråkdelen ↔ heltal" -bijektionen (betecknad nedan som ) kontinuerlig som en bitvis linjär funktion och monoton . Av detta kan vi omedelbart konstatera att den snabba inversa roten, som en kombination av kontinuerliga funktioner, är kontinuerlig. Och dess första del - skift-subtraktion - är också monoton och bitvis linjär. Bijektionen är komplicerad, men nästan "gratis": beroende på processorarkitekturen och anropskonventioner behöver du antingen inte göra någonting eller flytta numret från bråkregistret till heltalsregistret. ${\displaystyle I_{x))$

Till exempel är den binära representationen av det hexadecimala heltal 0x5F3759DF 0|101.1111.0|011.0111.0101.1001.1101.1111 2 (punkter är nibble-gränser, vertikala linjer är datorbråkgränser). Ordningen på 101 1111 0 2 är 190 10 , efter subtrahering av offset 127 10 får vi exponenten 63 10 . Den explicita delen av mantissan 01 101 110 101 100 111 011 111 2 efter tillägg av den implicita ledande enheten blir 1,011 011 101 011 001 110 111 11 2 = 1,432 4... 3 . Med hänsyn till den verkliga precisionen för datorbråk 0x5F3759DF ↔ 1,4324301 10 2 63 .

Vi betecknar den explicita delen av mantissan av numret , - oförskjuten ordning, - bitlängd av mantissan, - förskjutning av ordningen. Antalet , skrivet i ett linjärt-logaritmiskt bitrutnät av datorbråk, kan [10] [3] approximeras av ett logaritmiskt rutnät som , där är en parameter som används för att justera approximationens noggrannhet. Denna parameter sträcker sig från 0 (exakt vid och ) till 0,086 (exakt vid en punkt, ) $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$ $ett$ $m_{x}=0{,}443$

Med hjälp av denna approximation kan heltalsrepresentationen av ett tal approximeras som $x$

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

Följaktligen, . $\log _{2}x\approx {\frac {I_{x}}{L}}-(B-\sigma )$

Låt oss göra samma sak [3] för (respektive ) och erhålla ${\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]

Den magiska konstanten , med hänsyn till gränserna , kommer i bråkaritmetik att ha en opartisk ordning och mantissa ), och i binär notation - 0|101.1111.0|01 1 ... (1 är en implicit enhet; 0,5 kom från ordningen ; en liten enhet motsvarar intervallet [1,375; 1,5) och är därför mycket troligt, men inte garanterat av våra uppskattningar.) ${\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$

Det är möjligt att beräkna vad den första bitvis linjära approximationen [11] är lika med (i källan, inte själva mantissan, utan dess explicita del används ): $t=c-1$

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

På större eller mindre ändras resultatet proportionellt: när det fyrdubblas halveras resultatet exakt. $x$ $x$

Newtons metod ger [11] , , och . Funktionen är avtagande och konvex nedåt, på sådana funktioner närmar sig Newtons metod det sanna värdet från vänster - därför underskattar algoritmen alltid svaret. $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)$

Det är inte känt var konstanten 0x5F3759DF ↔ [a] 1.4324301 2 63 kom ifrån . Med brute force upptäckte Chris Lomont och Matthew Robertson [1] [2] att den bästa konstanten [b] i termer av marginellt relativ fel för float är 0x5F375A86 ↔ 1,4324500 2 63 , för dubbelt är det 0x5FE6EB50C7B537A9. Sant, för dubbel är algoritmen meningslös (ger ingen vinst i noggrannhet jämfört med float) [2] . Lomont-konstanten erhölls också analytiskt ( c = 1,432450084790142642179 ) [b] , men beräkningarna är ganska komplicerade [11] [2] .

Efter ett steg i Newtons metod är resultatet ganska exakt ( +0 % -0,18 % ) [1] [2] , vilket är mer än lämpligt för datorgrafikändamål ( 1 ⁄ 256 ≈ 0,39 % ). Ett sådant fel bevaras över hela området av normaliserade bråktal. Två steg ger en noggrannhet på 5 siffror [1] , efter fyra steg uppnås dubbelfelet . Om så önskas kan du balansera om felet genom att multiplicera koefficienterna 1,5 och 0,5 med 1,0009 så att metoden ger symmetriskt ±0,09% - det är vad de gjorde i spelet Interstate '76 och Blinn-metoden, som också itererar Newtons metod [7 ] .

Newtons metod garanterar inte monotonitet, men datoruppräkning visar att monotonitet existerar.

Källkod (C++) #include <iostream> union FloatInt { flyta somFlöta ; int32_t asInt ; }; int floatToInt ( float x ) { FloatInt r ; r . asFloat = x ; returnera r . asInt ; } float intToFloat ( int x ) { FloatInt r ; r . asInt = x ; returnera r . asFloat ; } float Q_rsqrt ( flytnummer ) _ { lång i ; flyta x2 , y ; const float threehalvs = 1,5F ; x2 = antal * 0,5F ; y = nummer ; i = * ( lång * ) & y ; // evil flyttal bitnivå hacking i = 0x5f3759df - ( i >> 1 ); // vad fan? y = * ( flytande * ) & i ; // jag vet inte vad fan! y = y * ( trehalvor - ( x2 * y * y ) ); // Första iterationen returnera y ; } int main () { int iStart = floatToInt ( 1.0 ); int iEnd = floatToInt ( 4.0 ); std :: cout << "Number to go: " << iEnd - iStart << std :: endl ; int nProblem = 0 ; float oldResult = std :: numeric_limits < float >:: infinity (); for ( int i = iStart ; i <= iEnd ; ++ i ) { float x = intToFloat ( i ); float resultat = Q_rsqrt ( x ); if ( resultat > gammalt Resultat ) { std :: cout << "Hittade ett problem på " << x << std :: endl ; ++ nProblem ; } } std :: cout << "Totala problem: " << nProblem << std :: endl ; returnera 0 ; }

Det finns liknande algoritmer för andra potenser, såsom kvadrat eller kubrot [3] .

Motivation

"Direkt" överlagring av belysning på en tredimensionell modell , till och med en högpoly-modell, även med hänsyn till Lamberts lag och andra formler för reflektion och spridning, kommer omedelbart att ge ett polygonalt utseende - betraktaren kommer att se skillnaden i belysning längs med polyederns kanter [12] . Ibland är detta nödvändigt - om föremålet är riktigt kantigt. Och för kurvlinjära föremål gör de detta: i ett tredimensionellt program indikerar de om en skarp kant eller en slät kant [12] . Beroende på detta, även när du exporterar modellen till spelet, beräknar trianglarnas hörn normalen för enhetslängd till den krökta ytan. Under animering och rotation omvandlar spelet dessa normaler tillsammans med resten av 3D-data; vid applicering av belysning interpolerar den över hela triangeln och normaliserar (tar den till en enhetslängd).

För att normalisera en vektor måste vi dividera alla tre dess komponenter med längden. Eller, bättre, multiplicera dem med den reciproka av längden: . Miljontals av dessa rötter måste beräknas per sekund. Innan dedikerad hårdvara för transformation och ljusbehandling skapades kunde datormjukvaran vara långsam. I synnerhet i början av 1990-talet, när koden utvecklades, släpade de flesta flyttalsberäkningar efter heltalsoperationer i prestanda. $(x',y',z')=(x,y,z){\frac {1}{\sqrt {x^{2}+y^{2}+z^{2))} }$

Quake III Arena använder en snabb invers kvadratrotsalgoritm för att påskynda grafikbearbetningen av processorn, men algoritmen har sedan dess implementerats i vissa specialiserade hårdvaruvertexshaders med hjälp av dedikerade programmerbara arrayer (FPGA).

Även på 2010-talsdatorer, beroende på belastningen av den fraktionerade samprocessorn , kan hastigheten vara tre till fyra gånger högre än att använda standardfunktioner [11] .

Kommentarer

↑ Här betyder pilen ↔ bijektionen av den binära representationen av ett heltal och den binära representationen av ett flyttalstal i IEEE 754 -formatet som förklaras ovan .
↑ 1 2 Om du sätter 127 i beställningsfältet får du 0x3FB75A86. GRISU2-biblioteket, som är helt och hållet heltal och oberoende av samprocessorfinesser, säger att 0x3FB75A86 ↔ 1,43245 är det kortaste decimaltalet som konverteras till en given float . Den minst signifikanta enheten är dock fortfarande 1,19 10 −7 och 0x3FB75A86 = 1,432450056 ≈ 1,4324501. Nästa bråktal är 0x3FB75A87 ↔ 1.4324502 utan några subtiliteter. Därav den icke-intuitiva avrundningen av 1,43245008 till 1,4324500.

Anteckningar

↑ 1 2 3 4 Arkiverad kopia . Hämtad 25 augusti 2019. Arkiverad från originalet 6 februari 2009. (obestämd)
↑ 1 2 3 4 5 6 https://web.archive.org/web/20140202234227/http://shelfflag.com/rsqrt.pdf
↑ 1 2 3 4 Hummus och magneter . Hämtad 1 februari 2017. Arkiverad från originalet 13 januari 2017. (obestämd)
↑ Beyond3D - Ursprunget till Quake3s snabba InvSqrt() . Hämtad 4 oktober 2019. Arkiverad från originalet 10 april 2017. (obestämd)
↑ Beyond3D - Ursprunget till Quake3s snabba InvSqrt() - Del två . Hämtad 25 augusti 2019. Arkiverad från originalet 25 augusti 2019. (obestämd)
↑ Flyttalstrick | IEEE tidskrifter och tidskrifter | IEEE Xplore
↑ 1 2 Snabb ömsesidig kvadratrot... 1997?! — Shane Peelars blogg
↑ PFRSQRT - Beräkna det ungefärliga värdet av reciproka av kvadratroten av ett kort reellt värde - Club155.ru . Hämtad 4 oktober 2019. Arkiverad från originalet 16 oktober 2019. (obestämd)
↑ RSQRTSS - Beräkna ömsesidigt av kvadratroten av skalärt Single-Precision Flytpunktsvärde . Hämtad 6 oktober 2019. Arkiverad från originalet 12 augusti 2019. (obestämd)
↑ https://web.archive.org/web/20150511044204/http://www.daxia.com/bibis/upload/406Fast_Inverse_Square_Root.pdf
↑ 1 2 3 4 Schwidkes beräkning av den lindade kvadratroten med de magiska konstanterna - analytisk pidkhid . Hämtad 12 juni 2022. Arkiverad från originalet 17 april 2022. (obestämd)
↑ 1 2 3 Detta är normen: vad är normala kartor och hur fungerar de / Sudo Null IT News

Länkar

C. Lomont, Fast invers kvadratrot , Teknisk rapport, 2003.
En kort historia av InvSqrt av Matthew Robertson
0x5f3759df , ytterligare undersökningar av algoritmens noggrannhet och generaliserbarhet av Christian Plesner Hansen