第一篇：卡馬克的求平方根的倒數(shù)的程序(快速平方根倒數(shù)算法)

卡馬克的求平方根的倒數(shù)的程序（快速平方根倒數(shù)算法）

2009-09-18 10:45:01| 分類： 離奇的code |字號(hào)大中小 訂閱

這個(gè)程序，大概在2006年看到的，當(dāng)時(shí)進(jìn)行了分析，主要分析了位運(yùn)算那部分：如何利用浮點(diǎn)數(shù)的位表示法快速計(jì)算一個(gè)近似值。而對(duì)于整個(gè)迭代計(jì)算原理并沒有仔細(xì)看，今天再重新分析一下，并把當(dāng)時(shí)那部分分析寫的更明了些。原貼如下：-------------發(fā)信人: interma(4PZP | 抓緊時(shí)間 | OFG P), 信區(qū): C_Cpp 標(biāo)

題: [zz] 源自Quake3的快速求InvSqrt()函數(shù)

發(fā)信站: 兵馬俑BBS(Mon Dec 4 13:11:36 2006), 本站(202.117.1.8)http://games.solidot.org/article.pl?sid=06/12/04/0210205&from=rss

"人們很早就在Quake3源代碼中發(fā)現(xiàn)了如下的C代碼，它可以快速的求1/sqrt(x)，在3D圖形向量計(jì)算方面應(yīng)用很廣。float InvSqrt(float x){

float xhalf=0.5f*x;

long i=*(long*)&x;

i=0x5f3759df(i>>1);

x=*(float *)&i;

其他地方都比較簡單，這三句的意思，應(yīng)當(dāng)是把x的浮點(diǎn)表示格式直接看出long類型的，然后進(jìn)行0x5f3759df(i>>1)

對(duì)于0x5f3759df如果看成浮點(diǎn)數(shù)的話，其指數(shù)位為(0101)(1111)0->1011111=190

設(shè)原來的xIEEE表示的指數(shù)位為a則新的結(jié)果的指數(shù)位變成了190-a/2；注意如果要轉(zhuǎn)換成

m*2^n格式，需要將指數(shù)a-127；所以得到的新結(jié)果用數(shù)學(xué)表示的指數(shù)位為190-a/2-127= 63-a/2

而x的原來如果從ieee格式變成數(shù)學(xué)表示應(yīng)為a-127

比較63-a/2與a-127

可以看出指數(shù)之間的關(guān)系為63-a/2==（a-127)*(-1/2)

所以經(jīng)過上述三句，x應(yīng)該已經(jīng)是1/sqrt(x)的近似值

然后下面的應(yīng)該是繼續(xù)進(jìn)行了精化，具體依據(jù)還沒搞懂

【 在 interma(4PZP | 抓緊時(shí)間 | OFG P)的大作中提到: 】

☆──────────────────────────────────────☆

xiaocui 于 Tue Dec 5 10:06:34 2006 提到：

厲害。

【 在 phylips(愛立佛)的大作中提到: 】

☆──────────────────────────────────────☆

NemoK 于 Tue Dec 5 11:58:01 2006 提到：

re!

【 在 phylips(愛立佛)的大作中提到: 】

☆──────────────────────────────────────☆

zqfalcon 于 Tue Dec 5 12:01:02 2006 提到：

呵呵大牛啊，佩服

【 在 phylips 的大作中提到: 】

------------------------

首先看牛頓迭代法的原理

若函數(shù)f（x）在點(diǎn)的某一臨域內(nèi)具有直到（n+1）階導(dǎo)數(shù)，則在該鄰域內(nèi)f（x）的n階泰勒公式為：

f(x)= f(x0)+(x－x0)f'(x0)+(x－x0)^2*f''(x0)/2!+… +...fn(x0)(x-x0)n/n!....解非線性方程f(x)=0的牛頓法是把非線性方程線性化的一種近似方法。把f(x)在x0點(diǎn)附近展開成泰勒級(jí)數(shù) f(x)= f(x0)+(x－x0)f'(x0)+(x－x0)^2*f''(x0)/2!+… 取其線性部分，作為非線性方程f(x)= 0的近似方程，即泰勒展開的前兩項(xiàng)，則有f(x0)+f'(x0)(x－x0)=f(x)=0 設(shè)f'(x0)≠0則f(x0)+f'(x0)(x－x0)=0 其解為x1=x0－f(x0)/f'(x0)這樣，得到牛頓法的一個(gè)迭代序列：x(n+1)=x(n)－f(x(n))/f'(x(n))。也就是說x(n+1)=x(n)－f(x(n))/f'(x(n))，是用來計(jì)算f(x)=0的根的迭代公式，如果x(n)與真根足夠靠近的話，那么只需要少數(shù)幾次迭代，就可以得到滿意的解。我們現(xiàn)在要求的是一個(gè)數(shù)的平方根的倒數(shù)，設(shè)這個(gè)數(shù)為a，那么我們要計(jì)算的就是x= a^(-0.5)。首先我們要把這個(gè)計(jì)算需求轉(zhuǎn)換為f(x)=0的形式。x = a^(-0.5)x^(-2)= a x^(-2)-a = 0 這樣就轉(zhuǎn)化成了f(x)= 0的形式，其中x就是a的平方根的倒數(shù)。f(x)= x^(-2)-a 則f'(x)=-2*x^(-3)根據(jù)牛頓迭代公式可以得到

x(n+1)= x(n)a)/(-2*x(n)^(-3))x(n+1)= x(n)+(x(n)^(-2)(a*x(n)^3)/2 x(n+1)= 1.5*x(n)(i>>1);

x=*(float *)&i;其他地方都比較簡單，這三句的意思，應(yīng)當(dāng)是把x的浮點(diǎn)表示格式直接看出long類型的，然后進(jìn)行0x5f3759df(i>>1);對(duì)于0x5f3759df如果看成浮點(diǎn)數(shù)的話，二進(jìn)制位為.10111110.***11011111,其指數(shù)位為10111110=190。設(shè)原來的x的IEEE表示的指數(shù)位為a則新的結(jié)果的指數(shù)位變成了190-a/2；注意如果要轉(zhuǎn)換成m*2^n格式，需要將指數(shù)減去127；所以得到的新結(jié)果用數(shù)學(xué)表示的指數(shù)位為190-a/2-127= 63-a/2，而x的原來的指數(shù)a，如果從ieee格式變成數(shù)學(xué)表示應(yīng)為a-127。

比較63-a/2與a-127，可以看出指數(shù)之間的關(guān)系為63-a/2==（a-127)*(-1/2)所以經(jīng)過上述三句，x應(yīng)該已經(jīng)是1/sqrt(x)的近似值。當(dāng)然上面只是從指數(shù)位進(jìn)行了簡單分析，如果要細(xì)化到尾數(shù)，還需要看具體效果。

--------------以下摘自：http://tbl.javaeye.com/blog/231425 雷神之錘III》里求平方根倒數(shù)的函數(shù)（快速平方根(倒數(shù))算法）

至于那個(gè)0x5f3759df，呃，我只能說，的確是一個(gè)超級(jí)的Magic Number。那個(gè)Magic Number是可以推導(dǎo)出來的，但我并不打算在這里討論，因?yàn)閷?shí)在太繁瑣了。簡單來說，其原理如下：因?yàn)镮EEE的浮點(diǎn)數(shù)中，尾數(shù)M省略了最前面的1，所以實(shí)際的尾數(shù)是1+M。如果你在大學(xué)上數(shù)學(xué)課沒有打瞌睡的話，那么當(dāng)你看到(1+M)^(-1/2)這樣的形式時(shí)，應(yīng)該會(huì)馬上聯(lián)想的到它的泰勒級(jí)數(shù)展開，而該展開式的第一項(xiàng)就是常數(shù)。下面給出簡單的推導(dǎo)過程：

對(duì)于實(shí)數(shù)R>0，假設(shè)其在IEEE的浮點(diǎn)表示中，指數(shù)為E，尾數(shù)為M，則： R^(-1/2)=(1+M)^(-1/2)* 2^(-E/2)將(1+M)^(-1/2)按泰勒級(jí)數(shù)展開，取第一項(xiàng)，得： 原式

=(1-M/2)* 2^(-E/2)= 2^(-E/2)(M/2)*2^(E/2)= C(R>>1)求出令到相對(duì)誤差最小的C值就可以了

上面的推導(dǎo)過程只是我個(gè)人的理解，并未得到證實(shí)。而Chris Lomont則在他的論文中詳細(xì)討論了最后那個(gè)方程的解法，并嘗試在實(shí)際的機(jī)器上尋找最佳的常數(shù)C。

而Chris Lomont則在他的論文中詳細(xì)討論了最后那個(gè)方程的解法，并嘗試在實(shí)際的機(jī)器上尋找最佳的常數(shù)C。

所以，所謂的Magic Number，并不是從N元宇宙的某個(gè)星系由于時(shí)空扭曲而掉到地球上的，而是幾百年前就有的數(shù)學(xué)理論。只要熟悉NR和泰勒級(jí)數(shù)，你我同樣有能力作出類似的優(yōu)化。

在GameDev.net上有人做過測試，該函數(shù)的相對(duì)誤差約為0.177585%，速度比C標(biāo)準(zhǔn)庫的sqrt提高超過20%。如果增加一次迭代過程，相對(duì)誤差可以降低到e-004的級(jí)數(shù)，但速度也會(huì)降到和sqrt差不多。據(jù)說在DOOM3中，Carmack通過查找表進(jìn)一步優(yōu)化了該算法，精度近乎完美，而且速度也比原版提高了一截。值得注意的是，在Chris Lomont的演算中，理論上最優(yōu)秀的常數(shù)（精度最高）是0x5f37642f，并且在實(shí)際測試中，如果只使用一次迭代的話，其效果也是最好的。但奇怪的是，經(jīng)過兩次NR后，在該常數(shù)下解的精度將降低得非常厲害（天知道是怎么回事?。?。經(jīng)過實(shí)際的測試，Chris Lomont認(rèn)為，最優(yōu)秀的常數(shù)是0x5f375a86。如果換成64位的double版本的話，算法還是一樣的，而最優(yōu)常數(shù)則為0x5fe6ec85e7de30da（又一個(gè)令人冒汗的Magic Number--b）。這個(gè)算法依賴于浮點(diǎn)數(shù)的內(nèi)部表示和字節(jié)順序，所以是不具移植性的。如果放到Mac上跑就會(huì)掛掉。如果想具備可移植性，還是乖乖用sqrt好了。但算法思想是通用的。大家可以嘗試推算一下相應(yīng)的平方根算法。