日前在書上看到一段使用多項式逼近計算平方根的代碼，至今都沒搞明白作者是怎樣推
算出那個公式的。但在嘗試解決問題的過程中，學到了不少東西，于是便有了這篇心
得，寫出來和大家共享。其中有錯漏的地方，還請大家多多指教。

的確，正如許多人所說的那樣，現(xiàn)在有有FPU，有3DNow，有SIMD，討論軟件算法好像不
合時宜。關(guān)于sqrt的話題其實早在2003年便已在?GameDev.net上得到了廣泛的討論（可
見我實在非常火星了，當然不排除還有其他尚在冥王星的人，嘿嘿）。而嘗試探究該話
題則完全是出于本人的興趣和好奇心（換句話說就是無知）。

我只是個beginner，所以這種大是大非的問題我也說不清楚（在GameDev.net上也有很多
類似的爭論）。但無論如何，Carmack在DOOM3中還是使用了軟件算法，而多知道一點數(shù)
學知識對3D編程來說也只有好處沒壞處。3D圖形編程其實就是數(shù)學，數(shù)學，還是數(shù)學。

文章原本是用HTML編排的，所以只截取了部分有比較有趣的東西放在這里。原文在我的
個人主頁上，同時也提供了2篇論文的下載：http:
//greatsorcerer.go2.icpcn.com/info/fastsqrt.html

=========================================================

在3D圖形編程中，經(jīng)常要求平方根或平方根的倒數(shù)，例如：求向量的長度或?qū)⑾蛄繗w一
化。C數(shù)學函數(shù)庫中的sqrt具有理想的精度，但對于3D游戲程式來說速度太慢。我們希望
能夠在保證足夠的精度的同時，進一步提高速度。

Carmack在QUAKE3中使用了下面的算法，它第一次在公眾場合出現(xiàn)的時候，幾乎震住了所
有的人。據(jù)說該算法其實并不是Carmack發(fā)明的，它真正的作者是Nvidia的Gary?Tarolli
（未經(jīng)證實）。
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
//
//?計算參數(shù)x的平方根的倒數(shù)
//
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;
}
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
該算法的本質(zhì)其實就是牛頓迭代法（Newton-Raphson?Method，簡稱NR），而NR的基礎(chǔ)則
是泰勒級數(shù)（Taylor?Series）。NR是一種求方程的近似根的方法。首先要估計一個與方
程的根比較靠近的數(shù)值，然后根據(jù)公式推算下一個更加近似的數(shù)值，不斷重復直到可以
獲得滿意的精度。其公式如下：
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
函數(shù)：y=f(x)
其一階導數(shù)為：y'=f'(x)
則方程：f(x)=0?的第n+1個近似根為
x[n+1]?=?x[n]?-?f(x[n])?/?f'(x[n])
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
NR最關(guān)鍵的地方在于估計第一個近似根。如果該近似根與真根足夠靠近的話，那么只需
要少數(shù)幾次迭代，就可以得到滿意的解。

現(xiàn)在回過頭來看看如何利用牛頓法來解決我們的問題。求平方根的倒數(shù)，實際就是求方
程1/(x^2)-a=0的解。將該方程按牛頓迭代法的公式展開為：
x[n+1]=1/2*x[n]*(3-a*x[n]*x[n])
將1/2放到括號里面，就得到了上面那個函數(shù)的倒數(shù)第二行。

接著，我們要設(shè)法估計第一個近似根。這也是上面的函數(shù)最神奇的地方。它通過某種方
法算出了一個與真根非常接近的近似根，因此它只需要使用一次迭代過程就獲得了較滿
意的解。它是怎樣做到的呢？所有的奧妙就在于這一行：
i?=?0x5f3759df?-?(i?>>?1);?//?計算第一個近似根

超級莫名其妙的語句，不是嗎？但仔細想一下的話，還是可以理解的。我們知道，IEEE
標準下，float類型的數(shù)據(jù)在32位系統(tǒng)上是這樣表示的（大體來說就是這樣，但省略了很
多細節(jié)，有興趣可以GOOGLE）：
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
bits：31?30?...?0
31：符號位
30-23：共8位，保存指數(shù)（E）
22-0：共23位，保存尾數(shù)（M）
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
所以，32位的浮點數(shù)用十進制實數(shù)表示就是：M*2^E。開根然后倒數(shù)就是：M^(-1/2)*2^
(-E/2)。現(xiàn)在就十分清晰了。語句i>?>1其工作就是將指數(shù)除以2，實現(xiàn)2^(E/2)的部分。
而前面用一個常數(shù)減去它，目的就是得到M^(1/2)同時反轉(zhuǎn)所有指數(shù)的符號。

至于那個0x5f3759df，呃，我只能說，的確是一個超級的Magic?Number。

那個Magic?Number是可以推導出來的，但我并不打算在這里討論，因為實在太繁瑣了。
簡單來說，其原理如下：因為IEEE的浮點數(shù)中，尾數(shù)M省略了最前面的1，所以實際的尾
數(shù)是1+M。如果你在大學上數(shù)學課沒有打瞌睡的話，那么當你看到(1+M)^(-1/2)這樣的形
式時，應(yīng)該會馬上聯(lián)想的到它的泰勒級數(shù)展開，而該展開式的第一項就是常數(shù)。下面給
出簡單的推導過程：
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
對于實數(shù)R>0，假設(shè)其在IEEE的浮點表示中，
指數(shù)為E，尾數(shù)為M，則：

R^(-1/2)
=?(1+M)^(-1/2)?*?2^(-E/2)

將(1+M)^(-1/2)按泰勒級數(shù)展開，取第一項，得：

原式
=?(1-M/2)?*?2^(-E/2)
=?2^(-E/2)?-?(M/2)?*?2^(-E/2)

如果不考慮指數(shù)的符號的話，
(M/2)*2^(E/2)正是(R>>1)，
而在IEEE表示中，指數(shù)的符號只需簡單地加上一個偏移即可，
而式子的前半部分剛好是個常數(shù)，所以原式可以轉(zhuǎn)化為：

原式?=?C?-?(M/2)*2^(E/2)?=?C?-?(R>>1)，其中C為常數(shù)

所以只需要解方程：
R^(-1/2)
=?(1+M)^(-1/2)?*?2^(-E/2)
=?C?-?(R>>1)
求出令到相對誤差最小的C值就可以了
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
上面的推導過程只是我個人的理解，并未得到證實。而Chris?Lomont則在他的論文中詳
細討論了最后那個方程的解法，并嘗試在實際的機器上尋找最佳的常數(shù)C。有興趣的朋友
可以在文末找到他的論文的鏈接。

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

在GameDev.net?上有人做過測試，該函數(shù)的相對誤差約為0.177585%，速度比C標準庫的
sqrt提高超過20%。如果增加一次迭代過程，相對誤差可以降低到e-?004?的級數(shù)，但速
度也會降到和sqrt差不多。據(jù)說在DOOM3中，Carmack通過查找表進一步優(yōu)化了該算法，
精度近乎完美，而且速度也比原版提高了一截（正在努力弄源碼，誰有發(fā)我一份）。

值得注意的是，在Chris?Lomont的演算中，理論上最優(yōu)秀的常數(shù)（精度最高）是
0x5f37642f，并且在實際測試中，如果只使用一次迭代的話，其效果也是最好的。但奇
怪的是，經(jīng)過兩次NR后，在該常數(shù)下解的精度將降低得非常厲害（天知道是怎么回
事！）。經(jīng)過實際的測試，Chris?Lomont認為，最優(yōu)秀的常數(shù)是0x5f375a86。如果換成
64位的double版本的話，算法還是一樣的，而最優(yōu)常數(shù)則為?0x5fe6ec85e7de30da（又一
個令人冒汗的Magic?Number?-?-b）。

這個算法依賴于浮點數(shù)的內(nèi)部表示和字節(jié)順序，所以是不具移植性的。如果放到Mac上跑
就會掛掉。如果想具備可移植性，還是乖乖用sqrt好了。但算法思想是通用的。大家可
以嘗試推算一下相應(yīng)的平方根算法。

下面給出Carmack在QUAKE3中使用的平方根算法。Carmack已經(jīng)將QUAKE3的所有源代碼捐
給開源了，所以大家可以放心使用，不用擔心會受到律師信。
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
//
//?Carmack在QUAKE3中使用的計算平方根的函數(shù)
//
float?CarmSqrt(float?x){
union{
int?intPart;
float?floatPart;
}?convertor;
union{
int?intPart;
float?floatPart;
}?convertor2;
convertor.floatPart?=?x;
convertor2.floatPart?=?x;
convertor.intPart?=?0x1FBCF800?+?(convertor.intPart?>>?1);
convertor2.intPart?=?0x5f3759df?-?(convertor2.intPart?>>?1);
return?0.5f*(convertor.floatPart?+?(x?*?convertor2.floatPart));
}