關於平方根倒數速演算法（雷神之錘3，牛B）

關於平方根倒數速演算法（雷神之錘3，牛B）

Quake-III Arena (雷神之錘3)是90年代的經典遊戲之一。該系列的遊戲不但畫面和內容不錯，而且即使計算機配置低，也能極其流暢地執行。這要歸功於它3D引擎的開發者約翰-卡馬克（John Carmack）。事實上早在90年代初DOS時代，只要能在PC上搞個小動畫都能讓人驚歎一番的時候，John Carmack就推出了石破天驚的Castle Wolfstein, 然後再接再勵，doom, doomII, Quake...每次都把3-D技術推到極致。他的3D引擎程式碼資極度高效，幾乎是在壓榨PC機的每條運算指令。當初MS的Direct3D也得聽取他的意見，修改了不少API。

最近，QUAKE的開發商ID SOFTWARE 遵守GPL協議，公開了QUAKE-III的原始碼，讓世人有幸目睹Carmack傳奇的3D引擎的原碼。
這是QUAKE-III原始碼的下載地址：
http://www.fileshack.com/file.x?fid=7547

(下面是官方的下載網址，搜尋 “quake3-1.32b-source.zip” 可以找到一大堆中文網頁的
ftp://ftp.idsoftware.com/idstuff/source/quake3-1.32b-source.zip)

我們知道，越底層的函式，呼叫越頻繁。3D引擎歸根到底還是數學運算。那麼找到最底層的數學運算函式（在game/code/q_math.c）， 必然是精心編寫的。裡面有很多有趣的函式，很多都令人驚奇，估計我們幾年時間都學不完。

在game/code/q_math.c裡發現了這樣一段程式碼。它的作用是將一個數開平方並取倒，經測試這段程式碼比(float)(1.0/sqrt(x))快4倍：
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

#ifndef Q3_VM
#ifdef __linux__
assert( !isnan(y) ); // bk010122 - FPE?
#endif
#endif
return y;
}

函式返回1/sqrt(x)，這個函式在影像處理中比sqrt(x)更有用。
注意到這個函式只用了一次疊代！（其實就是根本沒用疊代，直接運算）。編譯，實驗，這個函式不僅工作的很好，而且比標準的sqrt()函式快4倍！要知道，編譯器自帶的函式，可是經過嚴格仔細的彙編優化的啊！

這個簡潔的函式，最核心，也是最讓人費解的，就是標註了“what the fuck?”的一句
i = 0x5f3759df - ( i >> 1 );

再加上y = y * ( threehalfs - ( x2 * y * y ) );
兩句話就完成了開方運算！而且注意到，核心那句是定點移位運算，速度極快！特別在很多沒有乘法指令的RISC結構CPU上，這樣做是極其高效的。

演算法的原理其實不復雜,就是牛頓迭代法,用x-f(x)/f'(x)來不斷的逼近f(x)=a的根。

簡單來說比如求平方根,f(x)=x^2=a ,f'(x)= 2*x,f(x)/f'(x)=x/2,把f(x)代入

x-f(x)/f'(x)後有(x+a/x)/2，現在我們選a=5,選一個猜測值比如2，
那麼我們可以這麼算
5/2 = 2.5; (2.5+2)/2 = 2.25; 5/2.25 = xxx; (2.25+xxx)/2 = xxxx ...
這樣反覆迭代下去，結果必定收斂於sqrt(5)，沒錯，一般的求平方根都是這麼算的
但是卡馬克(quake3作者)真正牛B的地方是他選擇了一個神祕的常數0x5f3759df 來計算那個猜測值
就是我們加註釋的那一行,那一行算出的值非常接近1/sqrt(n),這樣我們只需要2次牛 頓迭代就可以達到我們所需要的精度.
好吧 如果這個還不算NB,接著看:

普渡大學的數學家Chris Lomont看了以後覺得有趣，決定要研究一下卡馬克弄出來的
這個猜測值有什麼奧祕。Lomont也是個牛人，在精心研究之後從理論上也推匯出一個
最佳猜測值，和卡馬克的數字非常接近, 0x5f37642f。卡馬克真牛，他是外星人嗎？

傳奇並沒有在這裡結束。Lomont計算出結果以後非常滿意，於是拿自己計算出的起始
值和卡馬克的神祕數字做比賽，看看誰的數字能夠更快更精確的求得平方根。結果是
卡馬克贏了... 誰也不知道卡馬克是怎麼找到這個數字的。

最後Lomont怒了，採用暴力方法一個數字一個數字試過來，終於找到一個比卡馬克數
字要好上那麼一丁點的數字，雖然實際上這兩個數字所產生的結果非常近似，這個暴
力得出的數字是0x5f375a86。

Lomont為此寫下一篇論文，"Fast Inverse Square Root"。

論文下載地址：
http://www.math.purdue.edu/~clomont/Math/Papers/2003/InvSqrt.pdf
http://www.matrix67.com/data/InvSqrt.pdf

參考：<IEEE Standard 754 for Binary Floating-Point Arithmetic><FAST INVERSE SQUARE ROOT>

最後，給出最精簡的1/sqrt()函式：
float InvSqrt(float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x; // get bits for floating VALUE
i = 0x5f375a86- (i>>1); // gives initial guess y0
x = *(float*)&i; // convert bits BACK to float
x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
return x;
} 
大家可以嘗試在PC機、51、AVR、430、ARM、上面編譯並實驗，驚訝一下它的工作效率。

。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。

維基百科參考：

http://zh.wikipedia.org/wiki/%E5%B9%B3%E6%96%B9%E6%A0%B9%E5%80%92%E6%95%B0%E9%80%9F%E7%AE%97%E6%B3%95

論文：http://www.daxia.com/bibis/upload/406Fast_Inverse_Square_Root.pdf

以上為R的存在說明；

。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。

以下是R的計算

http://www.guokr.com/post/91389/

基礎知識1：i>>1
操作i>>1表示將二進位制數i向右移動一位，也就是將最後一位刪掉並在最前一位新增0
注意到我們將一個n位的十進位制數M刪掉最後一位之後就變成了n-1位，可以看做將這個十進位制數除以10之後向下取整floor(M/10)
同樣的，講一個二進位制數刪掉最後一位之後相當於將這個數除以2並向下取整floor(M/2)

這樣看上去i>>1就像是floor(i/2)，因為函式f(x)=1/sqrt(x)的一階導數是-1/2*(x)^-3/2，正好有個-1/2在前面，不禁讓人感覺 0x5f3759df - ( i >> 1 )是函式1/sqrt(x)在某一個點的一階泰勒展開。在Fast Inverse Square Root 裡面有這樣一段：

Eberly’s explanation was that this produced linear approximation, but is incorrect; we’ll see the guess is piecewise linear, and the function being approximated is not the same in all cases.

“Eberly的解釋是說這相當於做了線性近似，但是這個解釋是不對的。我們會看到這個估計值是分段線性的，並且這個近似函式在各種情況下也並不相同。”

為什麼這麼說呢？這裡需要用到基礎知識2：浮點數儲存方式
各種型別浮點數的儲存方式可以通過檢視IEEE745（完全不知道是什麼東西）瞭解
這裡用到的是32位單精度浮點數，並且總是正數，表示方式為：
0|E|M
其中0代表正數
E是指數，E=0相當於2^-127
M表示一個絕對值小於1的數，但是注意到這裡省略了一位。也就是說，當轉化位十進位制的時候應該表示為（1+M）
那麼換算之後的數就應該是：（1+M）2^（E-127）
這樣我們發現其實i>>1並不完全是floor(i/2)而是將一個數變成(floor(M/2)+1)*2^floor((E-127)/2)
而且根據E的奇偶性尾數可能還需要加上1/2

不妨令e=E-127，注意到1/sqrt(x)是讓原數的指數變為-e/2，這麼說來卡馬克可能僅僅是希望產生一個e/2而用上了位移，接下來就是要找到一個相減之後產生-e/2並且讓尾數儘量和(1+M)^-1/2相近

由於這個數必然為正數假設這個數值為：
0|R1|R2
接下來便是要分情況討論：
假設E為偶數，這時候指數的右移並不會影響到尾數的數值：
這時候e是奇數，令e=2d+1
那麼相減之後指數部分變為：

注意這裡的相減其實是將浮點數轉化為整型（也就是正則化）之後再相減，而不是普通的浮點數加減。
由於初始值一定要是正數，所以我們需要上式一定為正，因為0=<E<=256，所以R1>=128
因為我們討論的E為偶數，也就是末位數一定是0，所以不用考慮他右移後對M的影響，所以兩數相減之後的結果是：

注意這裡用M/2而不是floor(M/2)因為這一點點的誤差相較於其他誤差來說太小了
當然，還存在一種情況，那就是R2<M/2，其實二進位制的加減和十進位制差不多，如果尾數小了，那麼就要像更高位數“借位”，也就是說這種情況下相減之後的結果是：

我們定義：

那麼我們可以將這兩種情況合併為一個函式：

這個函式就是我們對函式1/sqrt(x)的近似了，那麼我們的目標就是讓這個函式的相對誤差|(y-y0)/y|儘量小：

這樣我們得到一個誤差方程：

之後：

注意這裡的1/8其實是湊出來的，具體湊法可以先假設一個小於一的正數a，由於0<R2<1，0<M<1，我們可以通過展開得到R1關於a的一個區間。讓a儘量小，使得這個區間範圍小於一。根據R1一定是整數的特性，我們可以確定使得誤差最小的R1。這裡得出R1=190，帶入十六進位制裡面並右移（注意開頭有一個表示符號的0）就是0x5f，正好是黑魔法常數的頭幾位。

那當E為奇數怎麼辦呢？其實是一樣的辦法，如果E為奇數，那麼M/2就需要加上1/2（尾數的第一位相當於1/2），根據同樣的方法，我們可以得到另外一個相對誤差函式：

其中：

有興趣可以算一算這種情況下R1應該為多少，作者十分偷懶地說由於需要讓常數同時應用於兩種情況，所以就讓R1=190了。

之後就是確定R2的值了，各種分段討論，過於糾結我們就不看了（反正最後也沒算對，攤手），確定下來大約在0.45左右，再通過軟體算得最優解。