Edison.X. Blog

這是個讓 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，其相對誤差如下表所列

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. 測試結果與結論

(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