Quake-III Arena (雷神之錘3)是90年代的經(jīng)典游戲之一。該系列的游戲不但畫面和內(nèi)容不錯(cuò)，而且即使計(jì)算機(jī)配置低，也能極其流暢地運(yùn)行。這要?dú)w功于它3D引擎的開發(fā)者約翰-卡馬克(John Carmack)。事實(shí)上早在90年代初DOS時(shí)代，只要能在PC上搞個(gè)小動(dòng)畫都能讓人驚嘆一番的時(shí)候，John Carmack就推出了石破天驚的Castle Wolfstein, 然后再接再勵(lì)，doom, doomII, Quake...每次都把3-D技術(shù)推到極致。他的3D引擎代碼資極度高效，幾乎是在壓榨PC機(jī)的每條運(yùn)算指令。當(dāng)初MS的Direct3D也得聽取他的意見，修改了不少API。
最近，QUAKE的開發(fā)商ID SOFTWARE 遵守GPL協(xié)議，公開了QUAKE-III的原代碼，讓世人有幸目睹Carmack傳奇的3D引擎的原碼。
這是QUAKE-III原代碼的下載地址：
http://www.fileshack.com/file.x?fid=7547
(下面是官方的下載網(wǎng)址，搜索 “quake3-1.32b-source.zip” 可以找到一大堆中文網(wǎng)頁的
ftp://ftp.idsoftware.com/idstuff/source/quake3-1.32b-source.zip)
我們知道，越底層的函數(shù)，調(diào)用越頻繁。3D引擎歸根到底還是數(shù)學(xué)運(yùn)算。那么找到最底層的數(shù)學(xué)運(yùn)算函數(shù)(在game/code/q_math.c)， 必然是精心編寫的。里面有很多有趣的函數(shù)，很多都令人驚奇，估計(jì)我們幾年時(shí)間都學(xué)不完。
在game/code/q_math.c里發(fā)現(xiàn)了這樣一段代碼。它的作用是將一個(gè)數(shù)開平方并取倒，經(jīng)測(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;
}
函數(shù)返回1/sqrt(x)，這個(gè)函數(shù)在圖像處理中比sqrt(x)更有用。
注意到這個(gè)函數(shù)只用了一次疊代!(其實(shí)就是根本沒用疊代，直接運(yùn)算)。編譯，實(shí)驗(yàn)，這個(gè)函數(shù)不僅工作的很好，而且比標(biāo)準(zhǔn)的sqrt()函數(shù)快4倍!要知道，編譯器自帶的函數(shù)，可是經(jīng)過嚴(yán)格仔細(xì)的匯編優(yōu)化的啊!
這個(gè)簡潔的函數(shù)，最核心，也是最讓人費(fèi)解的，就是標(biāo)注了“what the fuck?”的一句
i = 0x5f3759df - ( i >> 1 );
再加上y = y ( threehalfs - ( x2 y y ) );
兩句話就完成了開方運(yùn)算!而且注意到，核心那句是定點(diǎn)移位運(yùn)算，速度極快!特別在很多沒有乘法指令的RISC結(jié)構(gòu)CPU上，這樣做是極其高效的。
算法的原理其實(shí)不復(fù)雜,就是牛頓迭代法,用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，現(xiàn)在我們選a=5,選一個(gè)猜測(cè)值比如2，
那么我們可以這么算
5/2 = 2.5; (2.5+2)/2 = 2.25; 5/2.25 = xxx; (2.25+xxx)/2 = xxxx ...
這樣反復(fù)迭代下去，結(jié)果必定收斂于sqrt(5)，沒錯(cuò)，一般的求平方根都是這么算的
但是卡馬克(quake3作者)真正牛B的地方是他選擇了一個(gè)神秘的常數(shù)0x5f3759df 來計(jì)算那個(gè)猜測(cè)值
就是我們加注釋的那一行,那一行算出的值非常接近1/sqrt(n),這樣我們只需要2次牛 頓迭代就可以達(dá)到我們所需要的精度.
好吧 如果這個(gè)還不算NB,接著看:
普渡大學(xué)的數(shù)學(xué)家Chris Lomont看了以后覺得有趣，決定要研究一下卡馬克弄出來的
這個(gè)猜測(cè)值有什么奧秘。Lomont也是個(gè)牛人，在精心研究之后從理論上也推導(dǎo)出一個(gè)
最佳猜測(cè)值，和卡馬克的數(shù)字非常接近, 0x5f37642f?？R克真牛，他是外星人嗎?
傳奇并沒有在這里結(jié)束。Lomont計(jì)算出結(jié)果以后非常滿意，于是拿自己計(jì)算出的起始
值和卡馬克的神秘?cái)?shù)字做比賽，看看誰的數(shù)字能夠更快更精確的求得平方根。結(jié)果是
卡馬克贏了... 誰也不知道卡馬克是怎么找到這個(gè)數(shù)字的。
最后Lomont怒了，采用暴力方法一個(gè)數(shù)字一個(gè)數(shù)字試過來，終于找到一個(gè)比卡馬克數(shù)
字要好上那么一丁點(diǎn)的數(shù)字，雖然實(shí)際上這兩個(gè)數(shù)字所產(chǎn)生的結(jié)果非常近似，這個(gè)暴
力得出的數(shù)字是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()函數(shù)：
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機(jī)、51、AVR、430、ARM、上面編譯并實(shí)驗(yàn)，驚訝一下它的工作效率。
。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。
維基百科參考：
論文：http://www.daxia.com/bibis/upload/406Fast_Inverse_Square_Root.pdf
以上為R的存在說明;
。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。
以下是R的計(jì)算
http://www.guokr.com/post/91389/
基礎(chǔ)知識(shí)1：i>>1
操作i>>1表示將二進(jìn)制數(shù)i向右移動(dòng)一位，也就是將最后一位刪掉并在最前一位添加0
注意到我們將一個(gè)n位的十進(jìn)制數(shù)M刪掉最后一位之后就變成了n-1位，可以看做將這個(gè)十進(jìn)制數(shù)除以10之后向下取整floor(M/10)
同樣的，講一個(gè)二進(jìn)制數(shù)刪掉最后一位之后相當(dāng)于將這個(gè)數(shù)除以2并向下取整floor(M/2)
這樣看上去i>>1就像是floor(i/2)，因?yàn)楹瘮?shù)f(x)=1/sqrt(x)的一階導(dǎo)數(shù)是-1/2(x)^-3/2，正好有個(gè)-1/2在前面，不禁讓人感覺 0x5f3759df - ( i >> 1 )是函數(shù)1/sqrt(x)在某一個(gè)點(diǎn)的一階泰勒展開。在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的解釋是說這相當(dāng)于做了線性近似，但是這個(gè)解釋是不對(duì)的。我們會(huì)看到這個(gè)估計(jì)值是分段線性的，并且這個(gè)近似函數(shù)在各種情況下也并不相同?！?/p>
為什么這么說呢?這里需要用到基礎(chǔ)知識(shí)2：浮點(diǎn)數(shù)存儲(chǔ)方式
各種類型浮點(diǎn)數(shù)的存儲(chǔ)方式可以通過查看IEEE745(完全不知道是什么東西)了解
這里用到的是32位單精度浮點(diǎn)數(shù)，并且總是正數(shù)，表示方式為：
0|E|M
其中0代表正數(shù)
E是指數(shù)，E=0相當(dāng)于2^-127
M表示一個(gè)絕對(duì)值小于1的數(shù)，但是注意到這里省略了一位。也就是說，當(dāng)轉(zhuǎn)化位十進(jìn)制的時(shí)候應(yīng)該表示為(1+M)
那么換算之后的數(shù)就應(yīng)該是：(1+M)2^(E-127)
這樣我們發(fā)現(xiàn)其實(shí)i>>1并不完全是floor(i/2)而是將一個(gè)數(shù)變成(floor(M/2)+1)*2^floor((E-127)/2)
而且根據(jù)E的奇偶性尾數(shù)可能還需要加上1/2
不妨令e=E-127，注意到1/sqrt(x)是讓原數(shù)的指數(shù)變?yōu)?e/2，這么說來卡馬克可能僅僅是希望產(chǎn)生一個(gè)e/2而用上了位移，接下來就是要找到一個(gè)相減之后產(chǎn)生-e/2并且讓尾數(shù)盡量和(1+M)^-1/2相近
由于這個(gè)數(shù)必然為正數(shù)假設(shè)這個(gè)數(shù)值為：
0|R1|R2
接下來便是要分情況討論：
假設(shè)E為偶數(shù)，這時(shí)候指數(shù)的右移并不會(huì)影響到尾數(shù)的數(shù)值：
這時(shí)候e是奇數(shù)，令e=2d+1
那么相減之后指數(shù)部分變?yōu)椋?/p>
注意這里的相減其實(shí)是將浮點(diǎn)數(shù)轉(zhuǎn)化為整型(也就是正則化)之后再相減，而不是普通的浮點(diǎn)數(shù)加減。
由于初始值一定要是正數(shù)，所以我們需要上式一定為正，因?yàn)?=<e<=256，所以r1>=128
因?yàn)槲覀冇懻摰腅為偶數(shù)，也就是末位數(shù)一定是0，所以不用考慮他右移后對(duì)M的影響，所以兩數(shù)相減之后的結(jié)果是：
注意這里用M/2而不是floor(M/2)因?yàn)檫@一點(diǎn)點(diǎn)的誤差相較于其他誤差來說太小了
當(dāng)然，還存在一種情況，那就是R2<m p="" 2，其實(shí)二進(jìn)制的加減和十進(jìn)制差不多，如果尾數(shù)小了，那么就要像更高位數(shù)“借位”，也就是說這種情況下相減之后的結(jié)果是：<="">
我們定義：
那么我們可以將這兩種情況合并為一個(gè)函數(shù)：
這個(gè)函數(shù)就是我們對(duì)函數(shù)1/sqrt(x)的近似了，那么我們的目標(biāo)就是讓這個(gè)函數(shù)的相對(duì)誤差|(y-y0)/y|盡量小：
這樣我們得到一個(gè)誤差方程：
之后：
注意這里的1/8其實(shí)是湊出來的，具體湊法可以先假設(shè)一個(gè)小于一的正數(shù)a，由于0<r2<1，0<m<1，我們可以通過展開得到r1關(guān)于a的一個(gè)區(qū)間。讓a盡量小，使得這個(gè)區(qū)間范圍小于一。根據(jù)r1一定是整數(shù)的特性，我們可以確定使得誤差最小的r1。這里得出r1=190，帶入十六進(jìn)制里面并右移(注意開頭有一個(gè)表示符號(hào)的0)就是0x5f，正好是黑魔法常數(shù)的頭幾位。< p="">
那當(dāng)E為奇數(shù)怎么辦呢?其實(shí)是一樣的辦法，如果E為奇數(shù)，那么M/2就需要加上1/2(尾數(shù)的第一位相當(dāng)于1/2)，根據(jù)同樣的方法，我們可以得到另外一個(gè)相對(duì)誤差函數(shù)：
其中：
有興趣可以算一算這種情況下R1應(yīng)該為多少，作者十分偷懶地說由于需要讓常數(shù)同時(shí)應(yīng)用于兩種情況，所以就讓R1=190了。
之后就是確定R2的值了，各種分段討論，過于糾結(jié)我們就不看了(反正最后也沒算對(duì)，攤手)，確定下來大約在0.45左右，再通過軟件算得最優(yōu)解。
