快速浮點開方運算，浮點開方運算

快速浮點開方運算，浮點開方運算

代碼下載：開根號的幾種演算法實現

在之前的部落格中我們介紹了資料類型的地址轉換，利用它我們可以將一個float型的值直接看成一個int類型。這種地址轉換到底有什麼意義，或者說有什麼用途呢？今天，給大家展示一個執行個體—快速浮點開方運算，讓大家更加明白地址轉換的含義和它們之間的對應關係。

1 二分法

          浮點開方也就是給定一個浮點數x，求。這個簡單的問題有很多解，我們從最簡單最容易想到的二分開始講起。利用二分進行開平方的思想很簡單，就是假定中值為最終解。假定下限為0，上限為x，然後求中值；然後比較中值的平方和x的大小，並根據大小修改下限或者上限；重新計算中值，開始新的迴圈，直到前後兩次中值的距離小於給定的精度為止。需要注意的一點是，如果x小於1，我們需要將上限置為1，原因你懂的。代碼如下：

float SqrtByBisection(float n){float low,up,mid,last; low=0,up=(n<1?1:n); mid=(low+up)/2; do{if(mid*mid>n)up=mid; else low=mid;last=mid;mid=(up+low)/2; }while(fabsf(mid-last) > eps);return mid; }

這種方法非常直觀，也是面試過程中經常會問到的問題，不過這裡有一點需要特別注意：在精度判別時不能利用上下限而要利用前後兩次mid值，否則可能會陷入死迴圈！這是因為由於精度問題，在迴圈過程中可能會產生mid值和up或low中的一個相同。這種情況下，後面的計算都不會再改變mid值，因而在達不到精度內時就陷入死迴圈。但是改為判斷前後兩次mid值就不會有任何問題（為啥自己想）。大家可以找一些例子試一下，這可以算是二分法中的一個trick。二分雖然簡單，但是卻有一個非常大的問題：收斂太慢！也即需要迴圈很多次才能達到精度要求。這也比較容易理解，因為往往需要迭代3到4次才能獲得一位準確結果。為了能提升收斂速度，我們需要採用其它的方法。

2 牛頓迭代法

          原理也比較簡單，就是將中值替換為切線方程的零根作為最終解。原理可以利用解釋（from matrix67）：

圖一 牛頓迭代法求開方

假設現在要求的值（圖中a=2），我們將其等價轉化為求函數與x軸大於0的交點。為了獲得該交點的值，我們先假設一個初始值，在圖一中為。過的直線與交於一點，過該點做切線交x軸於，則是比好的一個結果。重複上述步驟，過直線與交於一點，過該點做切線交x軸於，則是比更好的一個結果……
很明顯可以看出該方法斜著逼近目標值，收斂速度應該快於二分法。但是如何由獲得呢，我們需要獲得一個遞推公式。看圖中的陰影三角形，豎邊的長度為，如果我們能求得橫邊的長度l，則很容易得到。因為三角形的斜邊其實是過的切線，所以我們可以很容易知道該切線的斜率為，然後利用正切的定義就可以獲得l的長度為，由此我們得到遞推公式為：

後面我們需要做的就是利用上面的公式去迭代，直到達到精度要求，代碼如下：

float SqrtByNewton(float x){float val=x;//初始值float last;do{last = val;val =(val + x/val) / 2;}while(fabsf(val-last) > eps);return val;}

對上述代碼進行測試，結果確實比二分法快，對前300萬的所有整數進行開方的時間分別為1600毫秒和1000毫秒，快的原因主要是迭代次數比二分法更少。雖然牛頓迭代更快，但是還有進一步最佳化的餘地：首先牛頓迭代的代碼中有兩次除法，而二分法中只有一次，通常除法要比乘法慢個幾倍，因而會導致單次迭代速度的下降，如果能消除除法，速度還能提高不少，後面會介紹沒有除法的演算法；其次我們選擇原始值作為初始估值，這其實不是一個好的估計，這就導致需要迭代多次才能達到精度要求。當然二分法也存在這個問題，但是上下限不容易估計，只能採用最保守的方式。而牛頓迭代則可以任意選擇初始值，所以就存在選擇的問題。

          我們分析一下為什麼牛頓迭代法可以任意選擇初值（當然必須要大於0）。由公式

我們可以得出幾個結論：

所以，牛頓迭代存在一個初值選擇的問題，選擇得好會極大降低迭代的次數，選擇得差效率也可能會低於二分法。我們先給出一個採用新初值的代碼：

float SqrtByNewton(float x){int temp = (((*(int *)&x)&0xff7fffff)>>1)+(64<<23);float val=*(float*)&temp;float last;do{last = val;val =(val + x/val) / 2;}while(fabsf(val-last) > eps);return val;}

對上述代碼重複之前的測試，已耗用時間由1000毫秒降為240毫秒，效能提升了接近4倍多！為啥改用上面複雜的兩句代碼就能使速度提升這麼多呢？這就需要用到我們之前部落格介紹的IEEE浮點數表示。我們知道，IEEE浮點標準用的形式來表示一個數，將該數存入float類型之後變為：

現在需要對這個浮點數進行開方，我們看看各部分都會大致發生什麼變化。指數E肯定會除以2,127保持不變，m需要進行開方。由於指數部分是浮點數的大頭，所以對指數的修改最容易使初始值接近精確值。幸運的是，對指數的開平方我們只需要除以2即可，也即右移一位。但是由於E+127可能是奇數，右移一位會修改指數，我們將先將指數的最低位清零，這就是& 0xff7fffff的目的。然後將該轉換後的整數右移一位，也即將指數除以2，同時尾數也除以2（其實只是尾數的小數部分除以2）。由於右移也會將127除以2，所以我們還需要補償一個64，這就是最後還需要加一個(64<<23)的原因。

          這裡大家可能會有疑問，最後為什麼加(64<<23)而不是(63<<23)，還有能不能不將指數最後一位清零？答案是都可以，但是速度都沒有我上面寫的快。這說明我上面的估計更接近精確值。下面簡單分析一下原因。首先假設e為偶數，不妨設e=2n，開方之後e則應該變為n,127保持不變，我們看看上述代碼會變為啥。e+127是奇數，會清零，這等價於e+126，右移一位變為n+63，加上補償的64，指數為n+127，正是所需！再假設e為奇數，不妨設e=2n+1，開方之後e應該變為n+1（不精確），127保持不變，我們看看上述代碼會變為啥。e+127是偶數等於2n+128，右移一位變為n+64，加上補償的64，指數為n+1+127，也是所需！這確實說明上述的估計比其他方法更精確一些，因而速度也更快一些。

          雖然最佳化之後的牛頓迭代演算法比二分快了很多，但是速度都還是低於庫函數sqrtf，同樣的測試sqrtf只需要100毫秒，效能是最佳化之後牛頓迭代演算法的3倍！庫函數到底是如何?的！這說明我們估計的初始值還不是那麼精確。不要著急，我們下面介紹一種比庫函數還要快的演算法，其效能又是庫函數的10倍！

3卡馬克演算法

          這個演算法是99年被人從一個遊戲源碼中扒出來的，作者號稱是遊戲界的大神卡馬克，但是追根溯源，貌似這個演算法存在的還要更久遠，原始作者已不可考，暫且稱為卡馬克演算法。啥都不說，先上代碼一睹為快：

float SqrtByCarmack( float number ){int i;float x2, y;const float threehalfs = 1.5F;x2 = number * 0.5F;y  = number;i  = * ( int * ) &y;     i  = 0x5f375a86 - ( i >> 1 ); y  = * ( float * ) &i;y  = y * ( threehalfs - ( x2 * y * y ) ); y  = y * ( threehalfs - ( x2 * y * y ) );  y  = y * ( threehalfs - ( x2 * y * y ) ); return number*y;}

掃一眼上面的代碼會有兩個直觀的感覺：這代碼居然沒有迴圈！這代碼居然沒有除法！第一眼見到該代碼的人都會被震撼！下面對該演算法進行解釋，要解釋需要看最原始的版本：

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 hackingi  = 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 removedreturn y;}

 圖2 牛頓迭代法求開方倒數

          最原始的版本不是求開方，而是求開方倒數，也即。為啥這樣，原因有二。首先，開方倒數在實際應用中比開方更常見，例如在遊戲中經常會執行向量的歸一化操作，而該操作就需要用到開方倒數。另一個原因就是開方倒數的牛頓迭代沒有除法操作，因而會比先前的牛頓迭代開方要快。但是上面的代碼貌似很難看出牛頓迭代的樣子，這是因為函數變了，由變為，因而求解公式也需改變，但是遞推公式的推導不變，二所示。按照之前的推導方式我們有：

由這個公式我們就很清楚地明白代碼y  =y*(threehalfs-(x2*y*y));  的含義，這其實就是執行了單次牛頓迭代。為啥只執行了單次迭代就完事了呢？因為單次迭代的精度已經達到相當高的程度，代碼也特別註明無需第二次迭代（達到遊戲要求的精度）。圖三給出了對從0.01到10000之間的數進行開方倒數的誤差（from維基百科），可以看出誤差很小，而且隨著數的增大而減小。

圖三 卡馬克演算法的誤差

          為什麼單次迭代就可以達到精度要求呢？根據之前的分析我們可以知道，最根本的原因就是選擇的初值非常接近精確解。而估計初始解的關鍵就是下面這句代碼：

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

正是由於這句代碼，特別是其中的“magic number”使演算法的初始解非常接近精確解。具體的原理又用到前面部落格介紹的地址強轉：首先將float類型的數直接進行地址轉換轉成int型（代碼中long在32位機器上等價於int），然後對int型的值進行一個神奇的操作，最後再進行地址轉換轉成float類型就是很精確的初始解。

          在前面的部落格中，我們曾經針對float型浮點數和對應的int型整數之間的關係給出一個公式：

其中，表示float型浮點數地址強轉後的int型整數，，x是原始的浮點數（尚未表示成float類型），B=127，是一個無窮小量。化簡一下上述公式我們得到：

有了這個公式我們就可以推導初始解的由來了。要求，我們可以將其等價轉化成，然後代入上面的公式我們就得到：

這個公式就是神奇操作的數學表示，公式中只有是未知量，其它都已知。的值沒有好的求解方法，數學家通過暴力搜尋加實驗的方法求得最優值為0.0450466，此時第一項就對應0x5f3759df。但是後來經過更仔細的實驗，大家發現用0x5f375a86可以獲得更好的精度，所以後來就改用此數。

          演算法的最終目的是要對浮點數開平方，而原始的卡馬克演算法求的是開方倒數，所以我們最初的代碼返回的結果是原始值乘以開方倒數。該演算法效能非常高，而且精度也很高，三次迭代精度就和系統函數一樣，但是速度只有系統函數sqrtf的十分之一不到，相當了得。

4 改進的牛頓迭代          卡馬克演算法也啟發我們能不能對原始的牛頓迭代開方演算法進行類似的修改。之前我們已經提供了一種方法去估計初始值，而且獲得了不錯的效能提升，我們希望通過按照卡馬克演算法的思路修改初始值的估計來獲得更大的效能提升。要求，我們可以將其等價轉化成，然後再代入上面的公式我們就得到：

新公式和卡馬克公式非常相似，只是係數發生了變化。這也很好理解，本來開方和開方倒數的對數表示只差一個負號。這個公式也只有是未知量，其它都已知。我沒有精力去暴力搜尋它的最優值，只能估計幾個：可以選擇等於0，也可以選擇和卡馬克演算法中一樣，這樣得到的magic number分別為0x1fc00000和0x1fbd1e2d。分別用這兩個值去計算，效果差不多，和我們之前的初始值估計相比效能又提升了大約25%，但是還比庫函數sqrtf慢一倍。

          這樣的結果讓人感到沮喪，按理說應該會比卡馬克演算法慢，但是依舊差20多倍就不能理解了，難道初始解選擇的還不好？一怒之下，我將迴圈去掉，也改成只迭代三次，得到的結果和系統函數得到的結果一樣，只是速度上慢了一倍而已！這個結果很令人吃驚，這說明do迴圈的開銷其實很大。這個結論可以通過在卡馬克演算法中也添加do迴圈得到：如果在卡馬克演算法中添加for迴圈之後，運行速度立刻降了10倍！為什麼do迴圈會如此慢呢？一個原因可能是只需fabsf的原因，其他原因還不詳。

          通過速度慢一倍我們能得到什麼結論呢？原始的牛頓迭代用了兩次除法（加法忽略），而卡馬克演算法則用了三次乘法（在返回結果時還有一次），都是迭代三次，它們的效能相差一倍，我們可以推出除法的已耗用時間大約是乘法的三倍多。這個結論啟示我們，最佳化掉除法是提速的一個重要途徑。

          在猜測的值時，我們實驗了兩個很隨意的值，但是結果卻很好，是否會存在更好的呢？答案是肯定的，但是它只有在一次迭代的時候才會有影響（和原始的卡馬克演算法相比），如果迭代三次，則的值將影響不大，在某個區間裡面的值都會得到同樣的結果，已耗用時間也一樣，因為結果已經足夠精確。

          如果將我最開始估計初始值的do迴圈也去掉，則三次迭代精度達不到要求，說明我自己臆想出來的初始值還是太差，初始值估計確實是一門學問。總結一下牛頓迭代和卡馬克演算法，我們能得到什麼經驗教訓呢？首先，為了獲得最好的效能，代碼中盡量不要有迴圈，這也就是迴圈展開存在的意義；其次，兩個演算法最後的對比完全是除法和乘法的對比，盡量通過數學變換消除代碼中的除法，這也會帶來不少的效能提升；深入瞭解浮點數在電腦中的儲存結構很重要，在不少問題上會給我們帶來很多極致效能的解法。

5 SSE彙編指令

          最佳化是永無止境的。在進階語言領域，卡馬克演算法是目前最快的開方演算法，但是在機器指令層面，邪惡的 Intel 提供了這樣一條指令RSQRTSS，從硬體上支援卡馬克演算法。RSQRTSS是一條SSE指令，SSE（Streaming SIMD Extensions）指令也即單指令多資料流式擴充指令能夠有效增強CPU浮點運算的能力。通常編譯器也會提供SSE指令的高層實現，從而允許使用者在C++代碼中不用編寫彙編代碼就可直接使用SSE指令的功能。下面給出用該指令的代碼（from尋找更快的平方根倒數演算法）：

float SqrtByRSQRTSS(float a){float b=a;__m128 in = _mm_load_ss(&b);__m128 out = _mm_rsqrt_ss(in);_mm_store_ss(&b, out);return a*b;}

由於sse指令用到的寄存器是128位，我們無法直接將float類型轉化為__m128，而是調用sse專門的load和store指令。至於效能，在TIMINGSQUARE ROOT中有一個結果（四），但是我沒有複現。而且當不做任何設定的時候，上述代碼的運行速度比卡馬克演算法慢3倍，這說明編譯器預設不啟用SSE指令最佳化。當指定參數/arch:SSE之後效能匹配卡馬克演算法，但是也沒有達到圖四的三倍效能差距。但是Intel肯定不可能實現一個效能不如進階語言的指令，對SSE指令的更詳細使用規則我不是很清楚，應該還是我沒有掌握SSE編程的基本配置規則，因而圖四的結果還是可信的。此外，還有一個稍微嚴重的問題，RSQRTSS的結果不精確，誤差是 ±1.5*2^-12，對精度要求比較嚴格的應用不能使用該指令。

圖四SSE指令與卡馬克演算法的效能對比

          在圖四中第二行還顯示了SQRTSS指令的效能，和改進的牛頓迭代法類似，效能也是慢於RSQRTSS指令，而且差距還不小。

6 總結          本部落格對浮點開方常用的方法進行了一一介紹，希望能對大家產生些積極的影響。最後，將我的所有實驗結果貼出來，供大家參考。我的CPU型號為雙核32位酷睿2 T5750，記憶體2G，測試的內容是對1到300w內的整數進行開方運算，已耗用時間如下：

後來又在公司的伺服器上重測了一遍，直接傷心了。伺服器CPU為24核64位至強E5-2630，記憶體128G，編譯器為花錢購買的icc編譯器。測試的內容是對1到1000w內的整數進行開方運算，已耗用時間如下：

看到這個結果，我只能說系統函數無敵了（沒有算錯）！上面寫的全部作廢，以編譯器實測為準……

目前電腦的速度最快可達到__次浮點運算

天河一號的峰值運算速度為每秒4700萬億次。天河一號運算1小時，相當於全國13億人同時計算340年以上；“天河一號”運算1天，相當於1台雙核的高檔案頭電腦運算620年以上。
 
對於一個數開方的簡便，快速運算

2的5*3/5次方=8

請採納!!!!!!!!