快速平方根演算法

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

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

//
// 計算引數x的平方根的倒數
//float invsqrt (float 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) convertor;
union 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));
}

float faster_sqrtf(float f)
return result;
}

除了基於nr的方法外，其他常見的快速演算法還有多項式逼近。下面的函式取自《3d遊戲程式設計大師技巧》，它使用乙個多項式來近似替代原來的長度方程，但我搞不清楚作者使用的公式是怎麼推導出來的（如果你知道的話請告訴我，謝謝）。

//
//  這個函式計算從(0,0)到(x,y)的距離，相對誤差為3.5%
//int fastdistance2d(int x, int y)
//// 該函式計算(0,0,0)到(x,y,z)的距離，相對誤差為8%
//float fastdistance3d(float fx, float fy, float fz)

還有一種方法稱為distance estimates（距離評估？），如下圖所示：

紅線所描繪的正八邊形上的點為：

octagon(x,y) = min((1/√2) * (|x|+|y|), max(|x|,|y|))

求出向量v1和v2的長度，則：

√(x^2+y^2) = (|v1|+|v2|)/2 * octagon(x,y)

到目前為止我們都在討論浮點數的方根演算法，接下來輪到整數的方根演算法。也許有人認為對整型資料求方根無任何意義，因為會得到類似99^(1/2)=9的結果。通常情況下確實是這樣，但當我們使用定點數的時候（定點數仍然被應用在很多系統上面，例如任天堂的gba之類的手持裝置），整數的方根演算法就顯得非常重要。對整數開平方的演算法如下。我並不打算在這討論它（事實是我也沒有仔細考究，因為在短期內都不會用到- -b），但你可以在文末james ulery的**中找到非常詳細的推導過程。

//
// 為了閱讀的需要，我在下面的巨集定義中新增了換行符
//#define step(shift)
if((0x40000000l >> shift) + sqrtval <= val)
else
//// 計算32位整數的平方根
//int32 xxglusqrtfx(int32 val)
sqrtval <<= (q_factor)/2;
return(sqrtval);
}

關於sqrt的話題早在2023年便已在 gamedev.net上得到了廣泛的討論（可見我實在非常火星了，當然不排除還有其他尚在冥王星的人，嘿嘿）。而嘗試**該話題則完全是出於本人的興趣和好奇心（換句話說就是無知）。其實現在隨著fpu的提公升和對向量運算的硬體支援，大部分系統上都提供了快速的sqrt實現。如果是處理大批量的向量的話，據說最快的方法是使用simd（據說而已，我壓根不懂），可同步計算4個向量。

快速平方根演算法

快速計算平方根

計算平方根的演算法

平方根問題

快速平方根演算法

快速計算平方根

計算平方根的演算法

平方根問題

相關推薦