平方根倒數速算法（卡馬克開方法）

平方根倒數速算法(Fast inverse square root)，常常和一個十六進制的常量 0x5f3759df聯繫起來。該算法被用來快速運算平方根倒數，速度是 float(1/sqrt(x)) 方法的4倍。該算法大概由上個世紀90年代的硅圖公司開發出來，後來出如今John Carmark的Quake III Arena的源碼中。

這是一個古老的算法，最先的討論見於2001年中國的CSDN論壇上。而且該段代碼可能已經不適用於當代的64bits機器，由於如今的64bits的機器上 long 型數據長度已經從4字節變成了8字節。可是我以爲，一代代新人對於一些老的經典的問題的討論是有實際應用價值之外的意義的。僅此而已，不作爭辯。

代碼總覽

下面的快速平方根倒數代碼來自Quake III Arena。

 1 float Q_rsqrt( float number )
 2 {
 3     long i;
 4     float x2, y;
 5     const float threehalfs = 1.5F;
 6 
 7     x2 = number * 0.5F;
 8     y  = number;
 9     i  = * ( long * ) &y;                       // evil floating point bit level hacking
10     i  = 0x5f3759df - ( i >> 1 );               // what the fuck? 
11     y  = * ( float * ) &i;
12     y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
13 //    y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed
14 
15     return y;
16 }

爲了求解平方根的倒數值，首先須要產生一個近似值。而後在利用數值計算的方法來減少近似值的偏差。近似值的選取很大程度上影響後面數值計算迭代次數。在這以前，這個近似值的選取是經過查表獲得的，利用所謂的lookup table。但時上面提供了比查表還要迅速的方法，而且速度是常規計算平方根倒數的四倍。雖然該方法在精度上有必定的損失，但時計算速度上彌補了這一點的不足。該算法是爲IEEE 754標準下的浮點數存儲規則設計的，Chris Lomot和Charles McEniry展現了其在其餘浮點存儲標準下的運用。

算法運行示例

這裏有一個例子，令x = 0.15625，咱們想計算的值。程序的第一步以下：

0011_1110_0010_0000_0000_0000_0000_0000  //x和i的各比特位表示
0001_1111_0001_0000_0000_0000_0000_0000  //右移一位後的結果: (i >> 1)
0101_1111_0011_0111_0101_1001_1101_1111  //magic number 0x5f3759df
0100_0000_0010_0111_0101_1001_1101_1111  //0x5f3759df - (i >> 1)的結果

 將上述的數若是是用IEEE 754標準存儲的浮點數的話，各個數值大小以下：

0_01111100_01000000000000000000000  //1.25 * 2^-3
0_00111110_00100000000000000000000  //1.125 * 2^-65
0_10111110_01101110101100111011111  //1.432430... * 2^+63
0_10000000_01001110101100111011111  //1.307430... * 2^+1

咱們神奇的發現最後的數字已經很接近準確值2.52982了！what the fuck, why??

算法原理

該算法經過下面幾個步驟計算的值：

1. 將32bits的浮點型變量reinterpret成整型變量，用於計算log₂(x)的近似值
2. 用上述的近似值計算log₂(1/√x)的近似值
3. 將上述的計算結果強轉回32bits浮點型變量，獲得迭代初值
4. 通過一次簡單的牛頓迭代，獲得足夠高精度的結果

• 將浮點數型數據強制轉換成整型，近似求解對數值

前面一篇博文專門介紹了IEEE 754標準下的浮點數表示方法，下面簡單的回顧一下。爲了可以存儲非零的浮點數，第一步就是將浮點數x寫成標準化的二進制形式：

這裏e_x是整數，m_x ∈ [0, 1)顯而易見。咱們把1.b₁b₂b₃...稱爲尾數1 + m_x的二進制表示，其中小數點左邊的1稱爲前導數位。很顯然，前導數位始終爲1，因此在浮點數的存儲過程當中將其忽略。

咱們如今來計算三個數：

• S_x，符號位。當x>0，等於0，x<0，等於1；
  
• E_x = e_x + B，偏移指數。爲了既能存儲正指數，又能存儲負指數。IEEE 754標準在真實的指數上加一個偏移量來做爲存儲指數，對於單精度浮點數來講B=127；
• M_x = m_x × L，L = 2²³
  

 首先一個被開方數確定是個非負數，因此符號位必定爲0，那麼，對其取以2爲底的對數可得。由於m_x ∈ [0, 1)，因此有。其中σ是用來調節近似的精確度的，當σ ≈ 0.0430357時能夠得到最優的近似結果。因此咱們有

若是咱們把浮點數強制轉換成整數的話，偏移指數變成了這個整數的高8位，尾數部分變成了整數的低23位。

這個整數的數值爲：

 從上式咱們能夠看出，整數值I_x實際上是對log2(x)平移和伸縮後的結果，而且很容易把log2(x)從新表示出來：

• 迭代初值的計算

好的，到如今爲止咱們知道了將浮點數強轉成整數發生了什麼。卡馬克開方法最重要的一步就是找出了一個迭代初值，這個初值神接近真實結果，使得迭代收斂速度奇快。那麼咱們接下來將要看到，所謂的Magic Number的來歷，以及這個迭代初值是怎麼求得的。

令y = 1/√x ，而後就有

根據上面剛剛獲得的結論有：

整理後獲得：

咱們發現了什麼？將平方根倒數和被開方數強轉成整數居然有這種近似關係！因此所謂的Magic Number就是3/2L(B-σ)的十六進制寫法，因此就有了以下的代碼：

i  = 0x5f3759df - ( i >> 1 );

第二項是經過右移操做計算½ I_x。

• 牛頓迭代法

經過上述的操做後，該算法將求得的整數從新強轉回浮點數，這樣就獲得了迭代的初始值。下面再來看看迭代是如何進行的。

令

y就是下面方程的解：

牛頓迭代法給出了下面的方法來迭代求解方程的根：

其中是前一次迭代的結果，對於第一次迭代來講就是迭代初值。因此上式寫成代碼的形式以下：

y = y*(threehalfs - x2*y*y);

重複上面的步驟，將每次計算的結果做爲輸入，就會不斷逼近真實結果。

因爲咱們經過前面的方法選取的初始值已經特別接近真實值了，沒有必要進行更屢次的迭代，所以源碼中甚至將第二次迭代的過程都註釋掉了。

 

Reference:

[1] https://en.wikipedia.org/wiki/Fast_inverse_square_root