Áttekintés

A Quake egy 3 dimenziós lövöldözős játék, ahol a képpontok meghatározásához nagyon sokszor kell egy szám négyzetgyökének inverzét kiszámolni. A Quake III Arena forráskódja a következőt tartalmazza:

float Q_rsqrt(float number)
{
  long i;
  float x2, y;
  const float threehalfs = 1.5F;

  x2 = number * 0.5F;
  y  = number;
  i  = * ( long * ) &y;                       // evil floating point bit level hacking
  i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
  y  = * ( float * ) &i;
  y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
  // y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

  return y;
}

A megvalósítás több szempontból is igénytelen, viszont egyúttal a szoftverfejlesztés történelmének egyik legkülönlegesebb megoldását tartalmazza. A célt elérte: olyan megoldást talált, ami rendkívül gyors, és meglehetősen jó közelítéssel határozza meg az inverz négyzetgyököt.

A leírás a https://en.wikipedia.org/wiki/Fast_inverse_square_root alapján készült.

A kód felbontása

A kód maga C nyelven íródott. Írjuk át C++-ra, kicsit javítsunk fel, és bontsuk két részre.

Első rész:

float y = number;
long i = *(long *)&y;
i = 0x5f3759df - (i >> 1);
y = *(float *)&i;

Itt a következő történik:

• A bemenő, lebegőpontos (float) formátumú számot átértelmezi 32 bites előjeles egész számmá (a C-ben és C++-ban long). Tehát nem átalakítja, hanem bitről bitre átértelmezi.
• Az egész számon elvégez egy olyan műveletet, amelyet ha visszaértelmezünk lebegőpontos számmá, akkor az inverz négyzetgyök közelítő értékét kapjuk.
• Végül visszaértelmezi. Ezáltal az y az inverz gyök meglehetősen jó közelítését tartalmazza.

Második rész:

const float threehalfs = 1.5F;
const float x2 = number * 0.5F;
y = y * (threehalfs - (x2 * y * y));
y = y * (threehalfs - (x2 * y * y));

Ez az ún. Newton-módszer, amellyel bármilyen kezdőértékből kiindulva pár lépés alatt meglehetősen jó közelítést kapunk. Mivel az előző lépésben már kaptunk egy nagyon jó közelítést, igazából elég egy, maximum kettő lépés a pontosításhoz. Ennél sokkal többet már csak amiatt sem érdemes végrehajtani, mert a lebegőpontos számábrázolás pontosságának a határait is feszegetjük.

Fontos kiemelni, hogy a függvény nem tartalmaz lassú műveletet: még általános osztást se, gyökvonást meg pláne nem. A bitmozgatás, a kivonás és a szorzás gyors műveletek.

Az alábbiakban megnézzük, hogy melyik kódrészlet mit jelent.

A lebegőpontos számábrázolás

Mielőtt belefognánk az elméleti háttér és a kód tulajdonképpen áttekintésébe, meg kell ismerkednünk a 32 bites lebegőpontos számábrázolással. Mivel több szabvány létezik, ez az algoritmus egy nemdefiniált működést (https://en.wikipedia.org/wiki/Undefined_behavior) használ ki.

A lebegőpontos szám felépítése az alábbi:

$x=\pm 2^{e}\cdot(1+m)$

ahol $m=[0...1)$

Az egyes bitek jelentése:

• Balról az első bit (S): előjel; 0 ha pozitív és 1, ha negatív. Ebben a példában csak pozitív számokkal dolgozunk, emiatt ez mindig 0 lesz.
• 8 bit (balról a második-kilencedik; E): kitevő. Itt $E = e + B$, ahol $B=127$ a kitevő eltolás (exponent bias).
• 23 bit (balról a tizedik-harminckettedik; M): alap. Itt $M = m\cdot L$, ahol $L=2^{23}$.

Például a 25 átalakítása:

• A legnagyobb kettő hatvány ami nem nagyobb mint 25 a 16, azaz 2 a negyediken.
• Az alap 25 / 16 = 1,5625. Ebből az m a 0,5625 lesz.

Tehát:

$x = 25 = 2^4\cdot(1 + 0.5625)$

A bitek:

• S = 0
• E = 4 + 127 = 131 = 10000011 bináris formában
• M = $0.5625\cdot 2^{23} = 4718592$, ami 10010000000000000000000 bináris formában.

Összesítve tehát: 01000001110010000000000000000000.

Az első közelítés kiszámolása

Logaritmus közelítése

A későbbiekben szükség lesz a kettes alapú logaritmus közelítésére, így először nézzük át a használt logaritmus azonosságokat.

$\log_2{\frac{1}{\sqrt{x}}} = \log_2{x^{-\frac{1}{2}}} = -\frac{1}{2}\log_2 x$

Mivel az x lebegőpontos formában az alábbi:

$x = 2^e (1 + m)$

ennek kiszámolja a logaritmusát, ezt kapjuk:

$\log_2 x = \log_2 (2^e (1 + m)) = \log_2 2^e + \log_2(1+m) = e + log_2(1+m)$

Itt jön egy matematikai hack: az 1+m kettes alapú logaritmusát m-mel, egészen pontosan egy m-hez közeli értékkel helyettesítjük. Lássunk pár példát:

• $\log_2(1+0) = 0$
• $\log_2(1+0.1) \approx 0.1375$
• $\log_2(1+0.2) \approx 0.26303$
• $\log_2(1+0.5) \approx 0.585$
• $\log_2(1+0.8) \approx 0.848$

A 0-nál az érték pontos, afelett valamelyest alábecsüljük. A hibák a 3 esetben kerekítve 0.0375, 0.063, 0.085, 0.048. Úgy tűnik, hogy érdemes egy 0.05 körüli értékkel megnövelni, hogy még pontosabb eredményt kapjunk. Ezt alkalmazzák is, és a jele $\sigma$. Az optimális értéke:

$\sigma = \frac{1}{2} - \frac{1 + \ln(\ln 2)}{2\ln 2} \approx 0.0430357$

A Quake algoritmusában egy ennél valamelyest magasabb értéket használnak: $\sigma \approx 0.0450466$.

Összerakva, tehát az x kettes alapú logaritmusának a közelítése

$\log_2 x \approx e + m + \sigma$

A lebegőpontos szám átértelmezése egésszé

Használjuk itt és a továbbiakban az alábbi jelöléseket:

• x: az eredeti szám, aminek az inverz négyzetgyökét keressük
• y: az inverz négyzetgyök közelítése, azaz az eredmény
• I: az x ill y átértelmezése egésszé
• E: a lebegőpontos kitevő, azaz I-nek a 23-30. bitjei
• M: a lebegőpontos alap, azaz az I-nek a 0-22. bitjei
• L: $2^{23}$
• B: 127
• $\sigma \approx 0.0450466$

Azonosságok, összefoglalva:

• $E = e + B$
• $M = m\cdot L$
• $\log_2 x \approx e + m + \sigma$

Ezeket felhasználva:

$I = EL + M$

Az EL az exponens biteket "feltolja" a 23-30. pozíciókba azzal, hogy megszoroztuk L-lel, és utána hozzáadjuk az M-et. A 25 példájában ez a következőt jelenti:

• E = 10000011
• L = 100000000000000000000000
• M = 10010000000000000000000

Ez esetben I = 10000011 10010000000000000000000 lesz, amit pont keresünk.

Folytassuk az átalakítást:

$I = EL + M = EL + mL = L(E + m) = L(e + B + m) = L(e + B + m + \sigma - \sigma) = L(e + m + \sigma + B - \sigma) = L(e + m + \sigma) + L(B - \sigma) \approx L\log_2 x + L(B - \sigma)$

Az utolsó lépésben a $\log_2 x \approx e + m + \sigma$ azonosságot használtuk ki. A fentiből fejezzük ki $\log_2 x$-et:

$I = L\log_2 x + L(B - \sigma)$

$I - L(B - \sigma) = L\log_2 x$

$\log_2 x = \frac{I}{L} - (B - \sigma)$

Az első közelítés

Most jelölje $I_x$ az $x$ érték egész alakját, míg $I_y$ az $y$ érték egész alakját. Induljuk ki az alábbiakból:

$\log_2 x \approx \frac{I_x}{L} - (B - \sigma)$
$\log_2 y \approx \frac{I_y}{L} - (B - \sigma)$
$\log_2 y = -\frac{1}{2}\log_2 x$

Behelyettesítve:

$\frac{I_y}{L} - (B - \sigma) \approx -\frac{1}{2}(\frac{I_x}{L} - (B - \sigma)) = -\frac{1}{2}\frac{I_x}{L} + \frac{1}{2}(B - \sigma)$

Ebből fejezzük ki $I_y$-t:

$\frac{I_y}{L} \approx -\frac{1}{2}\frac{I_x}{L} + \frac{1}{2}(B - \sigma) + (B - \sigma) = -\frac{1}{2}\frac{I_x}{L} + \frac{3}{2}(B - \sigma)$
$I_y \approx -\frac{1}{2}I_x + \frac{3}{2}L(B - \sigma)$

Tehát:

$I_y \approx \frac{3}{2}L(B - \sigma) - \frac{1}{2}I_x$

Itt az első tag konstans:

$\frac{3}{2}L(B - \sigma) = \frac{3}{2}2^{23}(127 - 0.0450466) \approx 1597463007$ = 0x5F3759DF

A programban a második tagot, azaz az osztást bitmozgatással érik el. Ez tehát megmagyarázza a függvény legkritikusabb részét:

i = 0x5f3759df - (i >> 1);

Hogy mennyire erős ez az algoritmus, azt jól illusztrálja a 25 példája, melynek négyzetgyöke 5, inverz gyöke 0.2, a közelített érték pedig 0.206398.

A Newton-módszer

A New-ton módszer (https://hu.wikipedia.org/wiki/Newton-módszer) az alábbi képletből indul ki (átírva x-et y-ra):

$y_{n+1} = y_n - \frac{f(y_n)}{f'(y_n)}$

Az alábbiból indulunk ki:

$y = \frac{1}{\sqrt{x}}$
$y^2 = \frac{1}{x}$
$\frac{1}{y^2} = x$
$\frac{1}{y^2} - x = 0$

Tehát:

$f(y) = \frac{1}{y^2} - x = 0$
$f'(y) = (y^{-2} - x)' = -2y^{-3} = -\frac{2}{y^3}$

Behelyettesítve:

$y_{n+1} = y_n - \frac{y_n^{-2} - x}{-\frac{2}{y_n^3}} = y_n + y_n^3\frac{y_n^{-2} - x}{2} = y_n + \frac{y_n - xy_n^3}{2} = y_n\left(1 + \frac{1 - xy_n^2}{2}\right) = y_n\left(\frac{3 - xy_n^2}{2}\right) = y_n\left(\frac{3}{2} - \frac{x}{2}y_n^2\right)$

Ez utóbbi pedig pont megegyezik ezzel:

y = y * (threehalfs - (x2 * y * y));

Egy Newton-iterációval az eredmény 0.19969, kettővel pedig 0.199999, ami már szinte pontos.