序

這是一個神奇的演算法！

一、介紹

起源於一篇《改變計算技術的偉大演算法》文章，知道這個演算法，然後google一下，維基講的還不錯，本文權當自己理清下思路。先貼原始碼，為《雷神之錘III競技場》原始碼中的應用例項，剝離了C語言前處理器的指令，並附上了原有的註釋。

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（對浮點數的邪惡位級hack）
	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（第二次迭代，可以刪除）
 
	return y;
}
 

演算法的思路：

1、計算出該浮點數的平方根倒數的首次近似值，見原始碼中的8-11行；

2、利用牛頓法迭代以加強精度，得到要求精度內的值（迭代次數根據精度要求調整，原始碼中一次迭代就滿足精度要求）。

演算法的巧妙之處在於程式碼中的四行藍色程式碼，大致過程是：將表示浮點數的位元組序列用來表示整數（i  = * ( long * ) &y;），然後通過一個巧妙的整數運算得到一個新的整數（i  = 0x5f3759df - ( i >> 1 );），最後將表示整數的位元組序列換回表示浮點數。所以弄清演算法關鍵障礙是：在計算機中是如何表示浮點數和整數的、整數運算又怎能算出浮點數的平方根倒數的近似值、0x5f3759df怎麼來的。 

二、浮點數和整數的表示方法

浮點數和整數儲存位數一定是相同的。這裡浮點數和整數都佔4位元組，32位。

一個浮點數是由32位二進位制位表示的有理數，分為三部分。其中符號佔1位，表示正負，記為Si；指數佔接下來的8位，表示經過偏移處理後的指數，即實際表示E（如圖中為124），需要偏移B（圖中為2的8次方減1,127。B為一個固定值），最後得指數值為E-B；有效數字（除最高位以外）佔剩下的23位，記為m（0<m<1）,圖中的。

所以浮點數的結構公式為：，  圖中

整數的表示相對簡單，符號佔1位，數值佔剩下的31位。如果用上圖的浮點數字節序列來表示整數，那麼，即.平方根倒數函式僅能處理正數，所以符號位均為0。

小結：對於同樣的32位二進位制數碼，若為浮點數表示時實際數值為，而若為整數表示時實際數值則為，其中，這裡n=32，b=8。式子中引入的新變數為：

   ------------------------------------等式1

，其中------------等式2

三、浮點數的平方根倒數近似值

理解浮點數和整數的表示後，下面開始推導。

平方根倒數方程為：

兩邊取對數有：

因為浮點數可表示為：，所以也有，代入上式有：

再度引入新數描述與近似值R間的誤差：由於，有，則在此可定義與x的關係為，其中介於0到1/3，所以將代入上式得：

將第二部分小結中等式1，等式2代入到上述方程中，有：

移項整理得：

又因為浮點規格儲存的正浮點數x，若將其作為長整型表示，則示值為：，所以x的平方根倒數的首次近似值的整數表示值為：

，其中，，，

n=32，b=8。

這個式子對應著原始碼中的這一行：i  = 0x5f3759df - ( i >> 1 );，然後將整數表示值換回表示浮點數：y  = * ( float * ) &i;。這樣就得到了浮點數的平方根倒數的近似值。

四、神祕的0x5f3759df 

由第三部分可知：0x5f3759df 對應著R，即3/2(B-)L.當R為0x5f3759df時，有.

"現在不僅該演算法的原作者不明，人們也仍無法明確當初選擇這個“魔術數字”的方法。Chris Lomont在研究中曾做了個試驗：他編寫了一個函式，以在一個範圍內遍歷選取R值的方式將逼近誤差降到最小，以此方法他計算出了線性近似的最優R值0x5f37642f（與程式碼中使用的0x5f3759df相當接近），但以之代入演算法計算並進行一次牛頓迭代後，所得近似值與代入0x5f3759df的結果相比精度卻仍略微更低。……在Charles McEniry的論文中，他使用了一種類似Lomont但更復雜的方法來優化R值：他最開始使用窮舉搜尋，所得結果與Lomont相同；而後他嘗試用帶權二分法尋找最優值，所得結果恰是程式碼中所使用的魔術數字0x5f3759df"---維基百科

五、結束

至於最後一步牛頓法，自行google。平方根倒數速演算法的神奇之處在於：1、充分利用了浮點數和整數在計算機中的表示，然後以兩次轉換表示和一次整數運算替換複雜的浮點數計算，最後通過牛頓法加強精度；2、R的取值。

參考資料