[轉(zhuǎn)]Quake-III代碼里神奇的浮點(diǎn)開方函數(shù)
[轉(zhuǎn)]Quake-III代碼里神奇的浮點(diǎn)開方函數(shù)
2010-09-30 19:42:18Quake-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è)小動畫都能讓人驚嘆一番的時(shí)候，John Carmack就推出了石破天驚的Castle Wolfstein, 然后再接再勵，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)測試這段代碼比(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)= 2*x,f(x)/f'(x)=x/2,把f(x)代入
x-f(x)/f'(x)后有(x+a/x)/2，現(xiàn)在我們選a=5,選一個(gè)猜測值比如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è)猜測值
就是我們加注釋的那一行,那一行算出的值非常接近1/sqrt(n),這樣我們只需要2次牛 頓迭代就可以達(dá)到我們所需要的精度.
好吧 如果這個(gè)還不算NB,接著看:
普渡大學(xué)的數(shù)學(xué)家Chris Lomont看了以后覺得有趣，決定要研究一下卡馬克弄出來的
這個(gè)猜測值有什么奧秘。Lomont也是個(gè)牛人，在精心研究之后從理論上也推導(dǎo)出一個(gè)
最佳猜測值，和卡馬克的數(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
參考：<IEEE Standard 754 for Binary Floating-Point Arithmetic><FAST INVERSE SQUARE ROOT>
最后，給出最精簡的1/sqrt()函數(shù)：
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機(jī)、51、AVR、430、ARM、上面編譯并實(shí)驗(yàn)，驚訝一下它的工作效率。
?
前兩天有一則新聞，大意是說 Ryszard Sommefeldt 很久以前看到這麼樣的一段 code (可能出自 Quake III 的 source code)：
float InvSqrt (float x) {
float xhalf = 0.5f*x;
int i = *(int*)&x;
i = 0x5f3759df - (i>>1);
x = *(float*)&i;
x = x*(1.5f - xhalf*x*x);
return x;
}
他一看之下驚為天人，想要拜見這位前輩高人，但是一路追尋下去卻一直找不到人；同時(shí)間也有其他人在找，雖然也沒找到出處，但是 Chris Lomont 寫了一篇論文 (in PDF) 解析這段 code 的演算法 (用的是 Newton’s Method，牛頓法；比較重要的是後半段講到怎麼找出神奇的 0x5f3759df 的)。
PS. 這個(gè) function 之所以重要，是因?yàn)榍?開根號倒數(shù) 這個(gè)動作在 3D 運(yùn)算 (向量運(yùn)算的部份) 裡面常常會用到，如果你用最原始的 sqrt() 然後再倒數(shù)的話，速度比上面的這個(gè)版本大概慢了四倍吧… XD
PS2. 在他們追尋的過程中，有人提到一份叫做 MIT HACKMEM 的文件，這是 1970 年代的 MIT 強(qiáng)者們做的一些筆記 (hack memo)，大部份是 algorithm，有些 code 是 PDP-10 asm 寫的，另外有少數(shù)是 C code (有人整理了一份列表)。
附：牛頓迭代法快速尋找平方根
???????下面這種方法可以很有效地求出根號a的近似值：首先隨便猜一個(gè)近似值x，然后不斷令x等于x和a/x的平均數(shù)，迭代個(gè)六七次后x的值就已經(jīng)相當(dāng)精確了。
???????例如，我想求根號2等于多少。假如我猜測的結(jié)果為4，雖然錯(cuò)的離譜，但你可以看到使用牛頓迭代法后這個(gè)值很快就趨近于根號2了：
( 4 + 2/ 4 ) / 2 = 2.25
( 2.25 + 2/ 2.25 ) / 2 = 1.56944..
( 1.56944..+ 2/1.56944..) / 2 = 1.42189..
( 1.42189..+ 2/1.42189..) / 2 = 1.41423..
....
?
??????
???????這種算法的原理很簡單，我們僅僅是不斷用(x,f(x))的切線來逼近方程x^2-a=0的根。根號a實(shí)際上就是x^2-a=0的一個(gè)正實(shí)根，這個(gè)函數(shù)的導(dǎo)數(shù)是2x。也就是說，函數(shù)上任一點(diǎn)(x,f(x))處的切線斜率是2x。那么，x-f(x)/(2x)就是一個(gè)比x更接近的近似值。代入 f(x)=x^2-a得到x-(x^2-a)/(2x)，也就是(x+a/x)/2。
源地址： http://blog.renren.com/GetEntry.do?id=491777510&owner=245298353?
更多文章、技術(shù)交流、商務(wù)合作、聯(lián)系博主
微信掃碼或搜索：z360901061
微信掃一掃加我為好友
QQ號聯(lián)系： 360901061
您的支持是博主寫作最大的動力，如果您喜歡我的文章，感覺我的文章對您有幫助，請用微信掃描下面二維碼支持博主2元、5元、10元、20元等您想捐的金額吧，狠狠點(diǎn)擊下面給點(diǎn)支持吧，站長非常感激您！手機(jī)微信長按不能支付解決辦法：請將微信支付二維碼保存到相冊，切換到微信，然后點(diǎn)擊微信右上角掃一掃功能，選擇支付二維碼完成支付。
【本文對您有幫助就好】元
