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

 快速平方根（平方根倒數）演算法

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

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

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

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

在3D圖形程式設計中，經常要求平方根或平方根的倒數，例如：求向量的長度或將向量歸一化。C數學函式庫中的sqrt具有理想的精度，但對於3D遊戲程式來說速度太慢。我們希望能夠在保證足夠的精度的同時，進一步提高速度。

Carmack在QUAKE3中使用了下面的演算法，它第一次在公眾場合出現的時候，幾乎震住了所有的人。據說該演算法其實並不是Carmack發明的，它真正的作者是Nvidia的Gary Tarolli（未經證實）。

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

//

// 計算引數x的平方根的倒數

//

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;

}

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

該演算法的本質其實就是牛頓迭代法（Newton-Raphson Method，簡稱NR），而NR的基礎則是泰勒級數（Taylor Series）。NR是一種求方程的近似根的方法。首先要估計一個與方程的根比較靠近的數值，然後根據公式推算下一個更加近似的數值，不斷重複直到可以獲得滿意的精度。其公式如下：

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

函式：y=f(x)

其一階導數為：y'=f'(x)

則方程：f(x)=0 的第n+1個近似根為

x[n+1] = x[n] - f(x[n]) / f'(x[n])

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

 NR最關鍵的地方在於估計第一個近似根。如果該近似根與真根足夠靠近的話，那麼只需要少數幾次迭代，就可以得到滿意的解。

現在回過頭來看看如何利用牛頓法來解決我們的問題。求平方根的倒數，實際就是求方程1/(x^2)-a=0的解。將該方程按牛頓迭代法的公式展開為：

x[n+1]=1/2*x[n]*(3-a*x[n]*x[n])

 將1/2放到括號裡面，就得到了上面那個函式的倒數第二行。

接著，我們要設法估計第一個近似根。這也是上面的函式最神奇的地方。它通過某種方法算出了一個與真根非常接近的近似根，因此它只需要使用一次迭代過程就獲得了較滿意的解。它是怎樣做到的呢？所有的奧妙就在於這一行：

i = 0x5f3759df - (i >> 1); // 計算第一個近似根

超級莫名其妙的語句，不是嗎？但仔細想一下的話，還是可以理解的。我們知道，IEEE標準下，float型別的資料在32位系統上是這樣表示的（大體來說就是這樣，但省略了很多細節，有興趣可以GOOGLE）：

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

bits：31 30 ... 0

31：符號位

30-23：共8位，儲存指數（E）

22-0：共23位，儲存尾數（M）

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

所以，32位的浮點數用十進位制實數表示就是：M*2^E。開根然後倒數就是：M^(-1/2)*2^(-E/2)。現在就十分清晰了。語句i> >1其工作就是將指數除以2，實現2^(E/2)的部分。而前面用一個常數減去它，目的就是得到M^(1/2)同時反轉所有指數的符號。

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

那個Magic Number是可以推匯出來的，但我並不打算在這裡討論，因為實在太繁瑣了。簡單來說，其原理如下：因為IEEE的浮點數中，尾數M省略了最前面的1，所以實際的尾數是1+M。如果你在大學上數學課沒有打瞌睡的話，那麼當你看到(1+M)^(-1/2)這樣的形式時，應該會馬上聯想的到它的泰勒級數展開，而該展開式的第一項就是常數。下面給出簡單的推導過程：

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

對於實數R>0，假設其在IEEE的浮點表示中，

指數為E，尾數為M，則：

R^(-1/2)

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

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

原式

= (1-M/2) * 2^(-E/2)

= 2^(-E/2) - (M/2) * 2^(-E/2)

如果不考慮指數的符號的話，

(M/2)*2^(E/2)正是(R>>1)，

而在IEEE表示中，指數的符號只需簡單地加上一個偏移即可，

而式子的前半部分剛好是個常數，所以原式可以轉化為：

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

所以只需要解方程：

R^(-1/2)

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

= C - (R>>1)

求出令到相對誤差最小的C值就可以了

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

上面的推導過程只是我個人的理解，並未得到證實。而Chris Lomont則在他的論文中詳細討論了最後那個方程的解法，並嘗試在實際的機器上尋找最佳的常數C。有興趣的朋友可以在文末找到他的論文的連結。

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

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

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

這個演算法依賴於浮點數的內部表示和位元組順序，所以是不具移植性的。如果放到Mac上跑就會掛掉。如果想具備可移植性，還是乖乖用sqrt好了。但演算法思想是通用的。大家可以嘗試推算一下相應的平方根演算法。

下面給出Carmack在QUAKE3中使用的平方根演算法。Carmack已經將QUAKE3的所有原始碼捐給開源了，所以大家可以放心使用，不用擔心會受到律師信。

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

//

// Carmack在QUAKE3中使用的計算平方根的函式

//

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));

}