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

本文轉載自查看原文 2016-08-14 08:20 3972 卡馬克開方法/ C++/ Algorithm/ IEEE 754

平方根倒數速算法(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??

算法原理

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

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

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

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

${\begin{aligned}x&=\pm 1.b_{1}b_{2}b_{3}\ldots \times 2^{e_{x}}\\&=\pm 2^{e_{x}}(1+m_{x})\end{aligned}}$

這里 $e x 是整數， m x \in [0, 1) 顯而易見。我們把 1.b 1 b 2 b 3 ... 稱為尾數 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 \times L ， L = 2 23$