Hurtig omvendt kvadratrod

Den aktuelle version af siden er endnu ikke blevet gennemgået af erfarne bidragydere og kan afvige væsentligt fra den version , der blev gennemgået den 17. januar 2022; kontroller kræver 76 redigeringer .

Hurtig tilbagevendende kvadratrod (også hurtig InvSqrt() eller 0x5F3759DF med den "magiske" konstant brugt , i decimal 1 597 463 007) er en hurtig tilnærmet algoritme til at beregne den omvendte kvadratrod for positive 32-bit flydende kommatal . Algoritmen bruger heltalsoperationer "subtrahere" og " bitforskydning ", samt brøkdele "fratrække" og "multiplicere" - uden langsomme operationer "divide" og "kvadratrod". På trods af at den er "hacky" på bitniveau, er tilnærmelsen monoton og kontinuerlig : tætte argumenter giver tætte resultater. Nøjagtighed (mindre end 0,2 % ned og aldrig op) [1] [2] er ikke nok til rigtige numeriske beregninger, men det er ganske nok til tredimensionel grafik . $y={\frac {1}{{\sqrt {x))))$

Algoritme

Algoritmen tager et 32-bit flydende kommanummer ( enkelt præcision i IEEE 754 -format ) som input og udfører følgende operationer på det:

Behandl et 32-bit brøktal som et heltal, udfør operationen y 0 = 5F3759DF 16 − (x >> 1) , hvor >> er en bitforskydning til højre. Resultatet behandles igen som et 32-bit brøktal.
Til afklaring kan en iteration af Newtons metode udføres : y 1 \u003d y 0 (1,5 - 0,5 xy 0 ²) .

Original algoritme:

float Q_rsqrt ( float nummer ) { lang i ; flyde x2 , y ; const float trehalvdele = 1,5F ; x2 = tal * 0,5F ; y = tal ; i = * ( lang * ) & y ; // ond floating point bit niveau hacking i = 0x5f3759df - ( i >> 1 ); // hvad fanden? y = * ( flydende * ) & i ; y = y * ( tre halvdele - ( x2 * y * y ) ); // 1. iteration // y = y * (tre halvdele - (x2 * y * y)); // 2. iteration, denne kan fjernes returnere y ; }

Korrigeres efter standarderne for moderne C-implementering under hensyntagen til mulige optimeringer og tværplatform :

float Q_rsqrt ( float nummer ) { const float x2 = tal * 0,5F ; const float trehalvdele = 1,5F ; union { flyde f ; uint32_t i ; } konv = { tal }; // medlem 'f' sat til værdien af 'nummer'. konv . i = 0x5f3759df - ( konv . i >> 1 ); konv . f *= trehalvdele - x2 * konv . f * konv . f ; returnere konv . f ; }

Implementeringen fra Quake III: Arena [3] vurderer, at længden er lig med , og bruger pointere til konvertering (optimeringen "hvis , ingen har ændret sig" kan fungere fejlagtigt; på GCC, når der kompileres for at "frigive", er en advarsel udløst). Og det indeholder også et uanstændigt ord - John Carmack , der lagde spillet ud i det offentlige domæne, forstod ikke, hvad der blev gjort der. floatlongfloatlong

I C++20 kan du bruge den nye funktion . bit_cast

#inkluder <bit> #inkluder <grænser> #include <cstdint> constexpr float Q_rsqrt ( floatnummer ) noexcept _ { static_assert ( std :: numeric_limits < float >:: is_iec559 ); // Tjek om målmaskinen er kompatibel float const y = std :: bit_cast < float > ( 0x5f3759df - ( std :: bit_cast < std :: uint32_t > ( nummer ) >> 1 )); returner y * ( 1,5f - ( tal * 0,5f * y * y )); }

Historie

Algoritmen er formentlig udviklet hos Silicon Graphics i 1990'erne, og en implementering dukkede op i 1999 i kildekoden til computerspillet Quake III Arena , men metoden dukkede først op på offentlige fora som Usenet i 2002-2003. Algoritmen genererer rimeligt nøjagtige resultater ved hjælp af den unikke første tilnærmelse af Newtons metode . På det tidspunkt var den største fordel ved algoritmen afvisningen af dyre floating-point-beregninger til fordel for heltalsoperationer. Inverse kvadratrødder bruges til at beregne indfaldsvinkler og refleksion til belysning og skygge i computergrafik.

Algoritmen blev oprindeligt tilskrevet John Carmack , men han foreslog, at Michael Abrash , en grafikspecialist, eller Terje Mathisen, en assemblersprogspecialist, bragte den til id Software . Undersøgelsen af problemet viste, at koden havde dybere rødder i både hardware- og softwareområderne inden for computergrafik. Rettelser og ændringer er foretaget af både Silicon Graphics og 3dfx Interactive , med den tidligst kendte version skrevet af Gary Tarolli for SGI Indigo . Måske er algoritmen opfundet af Greg Walsh og Clive Moler , Garys kolleger hos Ardent Computer [5] .

Jim Blinn, en specialist i 3D-grafik, har foreslået en lignende tabelformet invers kvadratrodsmetode [6] , der tæller op til 4 cifre (0,01%) og blev fundet ved adskillelse af spillet Interstate '76 (1997) [7] .

Med udgivelsen af 3DNow! AMD - processorer introducerede PFRSQRT [8] assembler -instruktionen til hurtig omtrentlig beregning af den inverse kvadratrod. Versionen for dobbelt er meningsløs - nøjagtigheden af beregninger vil ikke øges [2] - derfor blev den ikke tilføjet. I 2000 blev RSQRTSS-funktionen [9] tilføjet til SSE2 , hvilket er mere nøjagtigt end denne algoritme (0,04% mod 0,2%).

Analyse og fejl

Bitrepræsentationen af et 4-byte brøktal i IEEE 754 -format ser sådan ud:

Skilt
	Bestille								Mantissa
0	0	en	en	en	en	en	0	0	0	en	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				femten				otte				7				0

Vi beskæftiger os kun med positive tal (tegnbit er nul), ikke denormaliseret , ikke ∞ og ikke NaN . Sådanne tal skrives i standardform som 1,mmmm 2 ·2 e . Delen 1,mmmm kaldes mantissen , e er rækkefølgen . Hovedenheden er ikke lagret ( implicit enhed ), så lad os kalde værdien 0, mmmm den eksplicitte del af mantissen. Derudover har maskinbrøktal en forskudt rækkefølge : 2 0 skrives som 011.1111.1 2 .

På positive tal er "brøkdelen ↔ heltal" -bijektionen (betegnet nedenfor som ) kontinuerlig som en stykkevis lineær funktion og monoton . Herfra kan vi umiddelbart konstatere, at den hurtige inverse rod, som en kombination af kontinuerte funktioner, er kontinuert. Og dens første del - skift-subtraktion - er også monoton og stykkevis lineær. Bijektionen er kompliceret, men næsten "gratis": afhængigt af processorarkitekturen og kaldekonventioner skal du enten ikke gøre noget eller flytte tallet fra brøkregisteret til det heltal. ${\displaystyle I_{x))$

For eksempel er den binære repræsentation af det hexadecimale heltal 0x5F3759DF 0|101.1111.0|011.0111.0101.1001.1101.1111 2 (punkter er nibble-grænser, lodrette linjer er computerbrøkfeltgrænser). Rækkefølgen af 101 1111 0 2 er 190 10 , efter at have trukket forskydningen 127 10 fra får vi eksponenten 63 10 . Den eksplicitte del af mantissen 01 101 110 101 100 111 011 111 2 efter tilføjelse af den implicitte ledende enhed bliver 1,011 011 101 011 001 110 111 11 2 = 1,432 4... 3 . Under hensyntagen til den reelle præcision af computerbrøk 0x5F3759DF ↔ 1.4324301 10 2 63 .

Vi betegner den eksplicitte del af mantissen af tallet , - uforskudt rækkefølge, - bitlængde af mantissen, - forskydning af rækkefølgen. Tallet , skrevet i et lineært-logaritmisk bitgitter af computerbrøker, kan [10] [3] tilnærmes af et logaritmisk gitter som , hvor er en parameter, der bruges til at justere nøjagtigheden af tilnærmelsen. Denne parameter går fra 0 (nøjagtig ved og ) til 0,086 (nøjagtig på et 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$ $en$ $m_{x}=0{,}443$

Ved at bruge denne tilnærmelse kan heltalsrepræsentationen af et tal tilnærmes som $x$

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

Følgelig . $\log _{2}x\approx {\frac {I_{x}}{L}}-(B-\sigma )$

Lad os gøre det samme [3] for (henholdsvis ) og opnå ${\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 magiske konstant , under hensyntagen til grænserne , vil i brøkregning have en upartisk rækkefølge og mantisse ), og i binær notation - 0|101.1111.0|01 1 ... (1 er en implicit enhed; 0,5 kom fra rækkefølgen ; en lille enhed svarer til området [1,375; 1,5) og er derfor højst sandsynlig, men ikke garanteret af vores estimater.) ${\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 er muligt at beregne, hvad den første stykkevise lineære tilnærmelse [11] er lig (i kilden, ikke selve mantissen, men dens eksplicitte del er brugt ): $t=c-1$

For : ; $x\in [0{,}5;\;c-0{,}5)$ $y_{01}=-x+t+{\tfrac {3}{2}}=-x+c+{\tfrac {1}{2}}$
For : ; $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}}$
For :. _ $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 ændres resultatet proportionalt: når det firdobles, halveres resultatet nøjagtigt. $x$ $x$

Newtons metode giver [11] , , og . Funktionen er aftagende og konveks nedad; på sådanne funktioner nærmer Newtons metode sig den sande værdi fra venstre - derfor undervurderer algoritmen altid 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 vides ikke, hvor konstanten 0x5F3759DF ↔ [a] 1.4324301 2 63 kom fra . Ved brute force fandt Chris Lomont og Matthew Robertson ud af [1] [2] at den bedste konstant [b] i form af marginal relativ fejl for float er 0x5F375A86 ↔ 1,4324500 2 63 , for dobbelt er den 0x5FE6EB50C7B537A9. Sandt nok, for dobbelt er algoritmen meningsløs (giver ikke en gevinst i nøjagtighed sammenlignet med float) [2] . Lomont-konstanten blev også opnået analytisk ( c = 1,432450084790142642179 ) [b] , men beregningerne er ret komplicerede [11] [2] .

Efter et trin i Newtons metode er resultatet ret nøjagtigt ( +0 % -0,18 % ) [1] [2] , hvilket er mere end velegnet til computergrafikformål ( 1 ⁄ 256 ≈ 0,39 % ). En sådan fejl bevares over hele området af normaliserede brøktal. To trin giver en nøjagtighed på 5 cifre [1] , efter fire trin er fejlen dobbelt nået . Hvis det ønskes, kan man rebalancere fejlen ved at gange koefficienterne 1,5 og 0,5 med 1,0009, så metoden giver symmetrisk ±0,09% - det gjorde de i spillet Interstate '76 og Blinn-metoden, som også gentager Newtons metode [7 ] .

Newtons metode garanterer ikke monotonitet, men computeroptælling viser, at monotoni eksisterer.

Kildekode (C++) #include <iostream> union FloatInt { flyde somFloat ; int32_t asInt ; }; int floatToInt ( float x ) { FloatInt r ; r . asFloat = x ; returnere r . asInt ; } float intToFloat ( int x ) { FloatInt r ; r . asInt = x ; returnere r . asFloat ; } float Q_rsqrt ( float nummer ) { lang i ; flyde x2 , y ; const float trehalvdele = 1,5F ; x2 = tal * 0,5F ; y = tal ; i = * ( lang * ) & y ; // ond floating point bit niveau hacking i = 0x5f3759df - ( i >> 1 ); // hvad fanden? y = * ( flydende * ) & i ; // jeg ved ikke hvad fanden! y = y * ( tre halvdele - ( x2 * y * y ) ); // 1. iteration returnere y ; } int main () { int iStart = floatToInt ( 1.0 ); int iEnd = floatToInt ( 4.0 ); std :: cout << "Numbers to go: " << iEnd - iStart << std :: endl ; int nProblemer = 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 > oldResult ) { std :: cout << "Fundet et problem på " << x << std :: endl ; ++ nProblemer ; } } std :: cout << "Totale problemer: " << nProblemer << std :: endl ; returnere 0 ; }

Der er lignende algoritmer for andre potenser, såsom kvadrat- eller terningrod [3] .

Motivation

"Direkte" overlejring af belysning på en tredimensionel model , selv en højpolymodel, selv under hensyntagen til Lamberts lov og andre refleksions- og spredningsformler, vil straks give et polygonalt udseende - seeren vil se forskellen i belysning langs med kanter af polyhedron [12] . Nogle gange er dette nødvendigt - hvis objektet er virkelig kantet. Og for krumlinjede objekter gør de dette: I et tredimensionelt program angiver de, om det er en skarp kant eller en glat [12] . Afhængigt af dette, selv når modellen eksporteres til spillet, beregner hjørnerne af trekanter normalen af enhedslængde til den buede overflade. Under animation og rotation transformerer spillet disse normaler sammen med resten af 3D-dataene; ved anvendelse af belysning interpolerer den over hele trekanten og normaliserer (bringer den til en enhedslængde).

For at normalisere en vektor skal vi dividere alle tre dens komponenter med længden. Eller, bedre, gange dem med den gensidige af længden: . Millioner af disse rødder skal beregnes pr. sekund. Før dedikeret transformations- og lysbehandlingshardware blev skabt , kunne computersoftware være langsom. Især i begyndelsen af 1990'erne, da koden blev udviklet, haltede de fleste flydende kommaberegninger efter heltalsoperationer i ydeevne. $(x',y',z')=(x,y,z){\frac {1}{\sqrt {x^{2}+y^{2}+z^{2}}} }$

Quake III Arena bruger en hurtig invers kvadratrodsalgoritme til at fremskynde grafikbehandlingen af CPU'en, men algoritmen er siden blevet implementeret i nogle specialiserede hardwarevertex shaders ved hjælp af dedikerede programmerbare arrays (FPGA'er).

Selv på computere fra 2010'erne, afhængigt af belastningen af den fraktionelle coprocessor , kan hastigheden være tre til fire gange højere end ved brug af standardfunktioner [11] .

Kommentarer

↑ Her betyder pilen ↔ bijektionen af den binære repræsentation af et heltal og den binære repræsentation af et flydende kommatal i IEEE 754 -formatet forklaret ovenfor .
↑ 1 2 Hvis du sætter 127 i ordrefeltet, får du 0x3FB75A86. GRISU2-biblioteket, som er fuldstændigt heltal og uafhængigt af coprocessor-finesser, siger, at 0x3FB75A86 ↔ 1,43245 er det korteste decimaltal, der konverteres til en given float . Den mindst signifikante enhed er dog stadig 1,19 10 −7 , og 0x3FB75A86 = 1,432450056 ≈ 1,4324501. Den næste brøk er 0x3FB75A87 ↔ 1.4324502 uden nogen finesser. Derfor den ikke-intuitive afrunding af 1,43245008 til 1,4324500.

Noter

↑ 1 2 3 4 Arkiveret kopi . Hentet 25. august 2019. Arkiveret fra originalen 6. februar 2009. (ubestemt)
↑ 1 2 3 4 5 6 https://web.archive.org/web/20140202234227/http://shelfflag.com/rsqrt.pdf
↑ 1 2 3 4 Hummus og magneter . Hentet 1. februar 2017. Arkiveret fra originalen 13. januar 2017. (ubestemt)
↑ Beyond3D - Origin of Quake3's Fast InvSqrt() . Hentet 4. oktober 2019. Arkiveret fra originalen 10. april 2017. (ubestemt)
↑ Beyond3D - Origin of Quake3's Fast InvSqrt() - Part Two . Hentet 25. august 2019. Arkiveret fra originalen 25. august 2019. (ubestemt)
↑ Flydende komma-tricks | IEEE tidsskrifter og magasiner | IEEE Xplore
↑ 1 2 Hurtig gensidig kvadratrod... i 1997?! — Shane Peelars blog
↑ PFRSQRT - Beregn den omtrentlige værdi af den reciproke af kvadratroden af en kort reel værdi - Club155.ru . Hentet 4. oktober 2019. Arkiveret fra originalen 16. oktober 2019. (ubestemt)
↑ RSQRTSS - Beregn reciprok af kvadratrod af skalar Single-Precision Floating-Point-værdi . Hentet 6. oktober 2019. Arkiveret fra originalen 12. august 2019. (ubestemt)
↑ https://web.archive.org/web/20150511044204/http://www.daxia.com/bibis/upload/406Fast_Inverse_Square_Root.pdf
↑ 1 2 3 4 Schwidkes beregning af den omviklede kvadratrod med de magiske konstanter - analytisk pidkhid . Hentet 12. juni 2022. Arkiveret fra originalen 17. april 2022. (ubestemt)
↑ 1 2 3 Dette er normen: hvad er normale kort, og hvordan fungerer de / Sudo Null IT News