這是個讓 CS 領域驚豔的計算方式,同時裡面使用的 magic number -  0x5f3759df 拿去 google 可得到一狗票的東西,本文並不進行 magic number 之推導,整個推導最詳盡的,應是 這篇pdf 。一開始的原始碼長得像這樣 (轉自 wiki)

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 [sic]
        i  = 0x5f3759df - ( i >> 1 );               // what the fuck? [sic]
        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. 原作

原作者是誰似乎已不詳。日前確定這段碼是在發展 3D 遊戲時所衍生。因 3D 遊戲 ( 沒記錯的話是 Quake-III Arena,雷神之錘3 )裡使用 1 / sqrt(x) 此動作頻繁,但速度很慢,在考慮沒必要那麼精確及速度上之問題,於是推了這個 magic number。至於作者是誰?至少確定不是上述連結 pdf 之原作 - Chris Lomont。

2. 原由

關於相關原由,可至 wiki 看,這裡只做簡介。

Quake III (雷神之錘3) 在 2005 前原始碼並未 release,但在 2002~2003 時就在論壇中出現。 John Carmack 看完之後為之震驚,便想尋訪拜會這位高人。一開始其他人建議去找 Terje Mathisen,Terje Mathisen 在 1990 年後期寫了許多相關類似的 bit code。此時 Jim Blinn 亦提出了其它文獻類似之方法,此法見於 1997 之 IEEE Computer Graphics and Applications,但最後卻沒結論得到誰是推導出這段碼之原作。

而後,Chris Lomont 發展後,求出當時他認為是最小誤差之 magic number - 0x5f37642f,趨近於 0x5f3759df,但這新之值卻較比一次迭代後之準確性為低,使用的是牛頓法之二次迭代。隨後 Lomont 找到的最佳數字 (應是用硬暴法) 為0x5f375a86(最大之相對誤差為 0.00176),這比先前找到的二個數字誤差都還要小!

隨後,在 64bits IEEE 754 上,Chris Lomont 也給了另一個 magic number:0x5fe6ec85e7de30da,已非常接近由 Charles McEniry 所找出之最佳值:0x5fe6eb50c7aa19f9。一開始  Charles McEniry 找到之數值與 Lomont 相同,之後再才 weighted bisection 方式找出了更好的數值。

3. 精度問題

對於 IEEE754 - 32bits,有給出 3 個 magic number,其中效果最好的為 0x5f375a86,並非最有名的 0x5f3759df,其相對誤差如下表所列

Value Predicted Tested 1 Iteration(%) 2 Iterations(%)
0x5f3759df 3.43758 3.43756 0.175228 4.66E-04
0x5f37642f 3.42128 3.42128 0.177585 4.78E-04
0x5f375a86 3.43655 3.43652 0.175124 4.65E-04

4. 程式碼

依原版之原始碼,程式碼之要分二部份撰寫,一部份用 pointer 硬轉型,一部份用 union 方式做運算,而 union 依合併方式也分二種,另再以 poniter 多寫一份二次迭代,因進行二次迭代後精度會大符提昇,再拿來測時比較需多耗之時間,以方便 coder 決定是否要再進行二次迭代。

 

(4.1) floating - pointer

float InvSqrt32_pointer(float x)
{
    float xhalf = 0.5f*x;
    int i = *(int*)&x; // get bits for floating value
    i = 0x5f375a86 - (i>>1); // gives initial guess y0
    x = *(float*)&i; // convert bits back to float
    x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
    return x;
}

 

 (4.2) floating - pointer - iterator 2

float InvSqrt32_pointer_itera2(float x)
{
    long i;
    float x2, y=x;
    const float threehalfs = 1.5f;

    x2 = x * 0.5f;
    i  = * ( long * ) &y;  // evil floating point bit level hacking
    i  = 0x5f37a86 - ( i >> 1 );
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
    y  = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration
    return y;
}

 

(4.3) floating - union1

float InvSqrt32_union1 (float x)
{    
    union{        
       int intPart;
       float floatPart;
    } convertor;    
    convertor.floatPart = x;    
    convertor.intPart = 0x5f378a86 - (convertor.intPart >> 1);    
    return 0.5f*convertor.floatPart*(3.0f - x*convertor.floatPart*convertor.floatPart);
}

 

 (4.4) floating - union2

float InvSqrt32_union2(float x)
{
    union{        
       int intPart;
       float floatPart;    
    } convertor;    
    convertor.floatPart = x;    
    convertor.intPart = 0x5f375a86 - (convertor.intPart >> 1);    
    return convertor.floatPart* \
(1.5f - 0.5f*x*convertor.floatPart*convertor.floatPart); }

 

(4.5) double - pointer

float InvSqrt32_union2(float x)
{
    union{        
       int intPart;
       float floatPart;    
    } convertor;    
    convertor.floatPart = x;    
    convertor.intPart = 0x5f375a86 - (convertor.intPart >> 1);    
    return convertor.floatPart* \
(1.5f - 0.5f*x*convertor.floatPart*convertor.floatPart); }

 (4.6) double - pointer - iterator2

double InvSqrt64_pointer_itera2(double x)
{
    long long int i;
    double x2, y=x;
    const double threehalfs = 1.5;

    x2 = x * 0.5;
    i  = * ( long long int* ) &y;  // evil floating point bit level hacking
    i  = 0x5fe6ec85e7de30da - ( i >> 1 );
    y  = * ( double * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
    y  = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration
    return y;
}

 

 (4.7) double - union1

double InvSqrt64_union1 (double x)
{    
    union{        
       long long int intPart;
       double floatPart;
    } convertor;    
    convertor.floatPart = x;    
    convertor.intPart = 0x5fe6ec85e7de30da - (convertor.intPart >> 1);    
    return 0.5*convertor.floatPart* \
(3.0 - x*convertor.floatPart*convertor.floatPart); }

 

(4.8) double - union2

double InvSqrt64_union2(double x)
{
    union{        
       long long int intPart;
       double floatPart;    
    } convertor;    
    convertor.floatPart = x;    
    convertor.intPart = 0x5fe6ec85e7de30da - (convertor.intPart >> 1);    
    return convertor.floatPart* \
(1.5 - 0.5*x*convertor.floatPart*convertor.floatPart); }

 

 5. 測試結果與結論

AMD Athlom™ II   VC2008   DevC-4.9.9.2  
X2 245 2.91G sum(n=1..10^9) elaspe sum(n=1..10^9) elaspe
1.0f/sqrt(float) 8192.000000 23051 8192.000000 86875
InvSqrt32_pointer 8192.000000 14922 8192.000000 16328
InvSqrt32_pointer_itera2 8192.000000 31234 8192.000000 27172
InvSqrt32_union1 8192.000000 14156 8192.000000 10719
InvSqrt32_union2 8192.000000 14610 8192.000000 10484
1.0/sqrt(double) 63244.092833 15891 63244.092833 81156
InvSqrt64_pointer 63184.364968 22359 63184.364968 31937
InvSqrt64_pointer_itera2 63243.975838 32797 63243.945838 45297
InvSqrt64_union1 63184.364968 25656 63184.364968 28391
InvSqrt64_union2 63184.364968 23438 63184.364968 29734

 

(5.1) 一般而言,的確是以 union 較 pointer 為佳,但 VC 於 double 之表示為 pointer 較 union 為佳,推測應為對於union 之 long long int 處理速度慢。

(5.2) 以 VC 言,在浮點數 double 之運算已有相當之速度,此公式唯於 floating 上表示較佳,約省下 33% 之運算時間;但若要再求進一步之精度使用二次之迭代,其運算時間卻又比 1/sqrt 方式為差。

(5.3) 以 dev-c  而言,都是以 union 為佳。雖 dev-c 內建之數學函式速度慢是事實,但值得注意的是,dev-c 之 floating-union 速度竟為最快!此為令吾人所驚。

 

6. Lemma

magic number 推導過程實為繁複,數學不夠強的也一定看不懂,不強的話有二條路可走:

(1) 暴力破解法:一個數字一個數字慢慢帶進去試,float 要試 2^32 次 ,double 要試 2^64 (若是 80bits 之 double 要試更多) 次;對 double 而言1、2天應該跑不掉。

(2) GA algorithm :考慮用 GA algorithm 建立模型,目前非常確定最小差誤是 float 之 0x5f375f86,以此下去試。若收斂性、準確性均佳,可考慮以此模型尋其餘公式之 magic number,提醒為,GA algorithm 找到並非保證最佳解,而是「近似」最佳解 (近似程度視於參數設定、模型建構方式,有空實做看看結果)。

 

7. 參考資料

wiki - methods of computing square roots

wiki - fast inverse squart root

edisonx 發表在 痞客邦 PIXNET 留言(0) 人氣()