Quake-III Arena (雷神之錘3)是90年代的經典游戲之一。該系列的游戲不但畫面和內容不錯，而且即使計算機配置低，也能極其流暢地運行。這要歸功于它3D引擎的開發(fā)者約翰-卡馬克(John Carmack)。事實上早在90年代初DOS時代，只要能在PC上搞個小動畫都能讓人驚嘆一番的時候，John Carmack就推出了石破天驚的Castle Wolfstein, 然后再接再勵，doom, doomII, Quake...每次都把3-D技術推到極致。他的3D引擎代碼資極度高效，幾乎是在壓榨PC機的每條運算指令。當初MS的Direct3D也得聽取他的意見，修改了不少API。
最近，QUAKE的開發(fā)商ID SOFTWARE 遵守GPL協(xié)議，公開了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里發(fā)現了這樣一段代碼。它的作用是將一個數開平方并取倒，經測試這段代碼比(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倍!要知道，編譯器自帶的函數，可是經過嚴格仔細的匯編優(yōu)化的啊!
這個簡潔的函數，最核心，也是最讓人費解的，就是標注了“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)= 2x,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?？R克真牛，他是外星人嗎?
傳奇并沒有在這里結束。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
參考：
最后，給出最精簡的1/sqrt()函數：
float InvSqrt(float x)
{
float xhalf = 0.5fx;
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-xhalfxx); // Newton step, repeating increases accuracy
return x;
}
大家可以嘗試在PC機、51、AVR、430、ARM、上面編譯并實驗，驚訝一下它的工作效率。
。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。
維基百科參考：
論文：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的解釋是說這相當于做了線性近似，但是這個解釋是不對的。我們會看到這個估計值是分段線性的，并且這個近似函數在各種情況下也并不相同?！?/p>
為什么這么說呢?這里需要用到基礎知識2：浮點數存儲方式
各種類型浮點數的存儲方式可以通過查看IEEE745(完全不知道是什么東西)了解
這里用到的是32位單精度浮點數，并且總是正數，表示方式為：
0|E|M
其中0代表正數
E是指數，E=0相當于2^-127
M表示一個絕對值小于1的數，但是注意到這里省略了一位。也就是說，當轉化位十進制的時候應該表示為(1+M)
那么換算之后的數就應該是：(1+M)2^(E-127)
這樣我們發(fā)現其實i>>1并不完全是floor(i/2)而是將一個數變成(floor(M/2)+1)*2^floor((E-127)/2)
而且根據E的奇偶性尾數可能還需要加上1/2
不妨令e=E-127，注意到1/sqrt(x)是讓原數的指數變?yōu)?e/2，這么說來卡馬克可能僅僅是希望產生一個e/2而用上了位移，接下來就是要找到一個相減之后產生-e/2并且讓尾數盡量和(1+M)^-1/2相近
由于這個數必然為正數假設這個數值為：
0|R1|R2
接下來便是要分情況討論：
假設E為偶數，這時候指數的右移并不會影響到尾數的數值：
這時候e是奇數，令e=2d+1
那么相減之后指數部分變?yōu)椋?/p>
注意這里的相減其實是將浮點數轉化為整型(也就是正則化)之后再相減，而不是普通的浮點數加減。
由于初始值一定要是正數，所以我們需要上式一定為正，因為0=<e<=256，所以r1>=128
因為我們討論的E為偶數，也就是末位數一定是0，所以不用考慮他右移后對M的影響，所以兩數相減之后的結果是：
注意這里用M/2而不是floor(M/2)因為這一點點的誤差相較于其他誤差來說太小了
當然，還存在一種情況，那就是R2<m p="" 2，其實二進制的加減和十進制差不多，如果尾數小了，那么就要像更高位數“借位”，也就是說這種情況下相減之后的結果是：<="">
我們定義：
那么我們可以將這兩種情況合并為一個函數：
這個函數就是我們對函數1/sqrt(x)的近似了，那么我們的目標就是讓這個函數的相對誤差|(y-y0)/y|盡量?。?/p>
這樣我們得到一個誤差方程：
之后：
注意這里的1/8其實是湊出來的，具體湊法可以先假設一個小于一的正數a，由于0<r2<1，0<m<1，我們可以通過展開得到r1關于a的一個區(qū)間。讓a盡量小，使得這個區(qū)間范圍小于一。根據r1一定是整數的特性，我們可以確定使得誤差最小的r1。這里得出r1=190，帶入十六進制里面并右移(注意開頭有一個表示符號的0)就是0x5f，正好是黑魔法常數的頭幾位。< p="">
那當E為奇數怎么辦呢?其實是一樣的辦法，如果E為奇數，那么M/2就需要加上1/2(尾數的第一位相當于1/2)，根據同樣的方法，我們可以得到另外一個相對誤差函數：
其中：
有興趣可以算一算這種情況下R1應該為多少，作者十分偷懶地說由于需要讓常數同時應用于兩種情況，所以就讓R1=190了。
之后就是確定R2的值了，各種分段討論，過于糾結我們就不看了(反正最后也沒算對，攤手)，確定下來大約在0.45左右，再通過軟件算得最優(yōu)解。
