【轉】代碼優化之-優化開根（卡馬克算法）

          最原始的版本不是求開方，而是求開方倒數，也即。爲啥這樣，緣由有二。首先，開方倒數在實際應用中比開方更常見，例如在遊戲中常常會執行向量的歸一化操做，而該操做就須要用到開方倒數。另外一個緣由就是開方倒數的牛頓迭代沒有除法操做，於是會比先前的牛頓迭代（ 從X_i-1=1開始迭代）開方要快。

           

        由這個公式咱們就很清楚地明白代碼y=y*(threehalfs-(x2*y*y))的含義，這其實就是執行了單次牛頓迭代。爲啥只執行了單次迭代就完事了呢？由於單次迭代的精度已經達到至關高的程度。 

         爲何單次迭代就能夠達到精度要求呢？根據以前的分析咱們能夠知道，最根本的緣由就是選擇的初值很是接近精確解。而估計初始解的關鍵就是下面這句代碼：

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

    正是因爲這句代碼，特別是其中的「magic number」使算法的初始解很是接近精確解。具體的原理是地址強轉：首先將float類型的數直接進行地址轉換轉成int型（代碼中long在32位機器上等價於int），而後對int型的值進行一個神奇的操做，最後再進行地址轉換轉成float類型就是很精確的初始解。

          float型浮點數和對應的int型整數之間的關係給出一個公式：

         

     其中，Ix表示 float型浮點數地址強轉後的int型整數，L=2^23，x是原始的浮點數（還沒有表示成float類型），B=127，是一個無窮小量。化簡一下上述公式咱們獲得：

 

             

        有了這個公式咱們就能夠推導初始解的由來了。要求 ，咱們能夠將其等價轉化成 ，而後代入上面的公式咱們就獲得：

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

          算法的最終目的是要對浮點數開平方，該算法性能很是高，並且精度也很高，三次迭代精度就和系統函數同樣，可是速度只有系統函數sqrtf的十分之一不到，至關了得。

 

#include "stdio.h"

#include "conio.h"

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 煩人的浮點位級處理 */

    i  = 0x5f3759df - ( i >> 1 ); /* what the fuck? 0x5f3759df or 0x5f375a86 什麼該死的？ 卡馬克算法 - dong - 北風寒*/

    y  = * ( float * ) &i;  /* 取長整型數i的地址，將其存儲單元轉換成浮點型，而後再把轉換後的數取出來*/

    y  = y * ( threehalfs - ( x2 * y * y ) );   /* 1st iteration 第一次迭代*/   

    y  = y * ( threehalfs - ( x2 * y * y ) );   /*  2nd iteration, this can be removed 第二次迭代，可以移除*/ 

    return y;  

}  

int main()

{

    float n,z=1.0;

    printf("請輸入一個須要求其平方根的數：");

    scanf("%f",&n);

    z=Q_rsqrt(n);

    printf("平方根爲%f\n",1.0/z);

    getch(); 

  return 0;
}

 

 舉例：X=2^e(1+f)=5.125=2^2(1+0.28125)

  Ix=EL+F=L(e+B+f)=2^23（2+127+0.28125）=2^23*10000001.01001=0(符號)10000001(階碼)

01001000000000000000000(尾數)（8388608*129.28125=1084489728）

 Ix表示浮點數的整數表示，E=e+B表示IEEE階碼值，L=表示階碼的起始位置，F=Lf表示尾數的整數表示

=12582912*(127-0.0450466)-1/2*1084489728=1597463007-542244864=1055218143=01111101 11001010101100111011111

y=*(float*)&i=2^(-2)*1.791805148124694824≈  0.447951287

y₁ =y(1.5-2.5625y^2)≈ 0.441593890

 

 

參考來源：夏風習習的博客->卡馬克算法，http://blog.163.com/lxd007_2005/blog/static/405618252015112410210140/

 

 

 

 

附：整數開根處理函數

 

/**************************************************************
 * 函數名：Q_rsqrt
 * 描  述：開根號函數（整型）
 * 輸  入：uint32_t radicand，被開方數
 * 輸  出：uint32_t result_sqrt，開方結果
 * 說  明：返回開方後的整數結果，運行16次。
**************************************************************/
uint32_t Q_rsqrt(uint32_t radicand)
{
    uint32_t i;    
    uint32_t result_sqrt;            // 開方結果
    uint32_t radicand_top,radicand_rmng;    
    if(radicand == 0)
    {
        return 0;
    }
    result_sqrt = 0;
    radicand_top = (radicand >> 30);
    radicand <<= 2;
    if(radicand_top > 1)
    {
        result_sqrt++;
        radicand_top -= result_sqrt;
    }
    for(i=15; i>0; i--)
    {
        result_sqrt <<= 1;
        radicand_top <<= 2;
        radicand_top += (radicand >> 30);
        radicand_rmng = result_sqrt;
        radicand_rmng = (radicand_rmng << 1) + 1;
        radicand <<= 2;
        if(radicand_top >= radicand_rmng)
        {
            radicand_top -= radicand_rmng;
            result_sqrt++;
        }
    }    
    return result_sqrt;
}