2013年10月26日 星期六

[轉貼]一個Sqrt函數引發的血案

轉自http://kb.cnblogs.com/page/189867/
講到當年Quake-III中一段開平方根的函數,竟然比系統內建sqrt()還要快!
我覺得它的設計原理可以給重視運算速度的programmer有很大的啟發 :)


    好吧,我承認我標題黨了,不過既然你來了,就認真看下去吧,保證你有收穫。
  我們平時經常會有一些數據運算的操作,需要調用sqrt,exp,abs等函數,那麼時候你有沒有想過:這個些函數係統是如何實現的?就拿最常用的sqrt函數來說吧,系統怎麼來實現這個經常調用的函數呢?
  雖然有可能你平時沒有想過這個問題,不過正所謂是“臨​​陣磨槍,不快也光”,你“眉頭一皺,計上心來”,這個不是太簡單了嘛,用二分的方法,在一個區間中,每次拿中間數的平方來試驗,如果大了,就再試左區間的中間數;如果小了,就再拿右區間的中間數來試。比如求sqrt(16)的結果,你先試(0+16)/2=8,8*8=64,64比16大,然後就向左移,試(0+8)/2=4, 4*4=16剛好,你得到了正確的結果sqrt(16)=4。然後你三下五除二就把程序寫出來了:

float SqrtByBisection(float n) //用二分法
{
    if(n<0) //小於0的按照你需要的處理 
    return n; 
    float mid,last; 
    float low,up; 
    low=0, up=n; 
    mid=(low+up)/2; 
    do
    {
        if(mid*mid>n)
            up=mid; 
        else 
            low=mid;
        last=mid;
        mid=(up+low)/2; 
    }while(abs(mid-last) > eps);//精度控制
    return mid; 
}

   然後看看和系統函數性能和精度的差別(其中時間單位不是秒也不是毫秒,而是CPU Tick,不管單位是什麼,統一了就有可比性)
      從圖中可以看出,二分法和系統的方法結果上完全相同,但是性能上整整差了幾百倍。為什麼會有這麼大的區別呢?難道系統有什麼更好的辦法?難道。。。。哦,對了,回憶下我們曾經的高數課,曾經老師教過我們“牛頓迭代法快速尋找平方根”,或者這種方法可以幫助我們,具體步驟如下:
求出根號a的近似值:首先隨便猜一個近似值x,然後不斷令x等於x和a/x的平均數,反覆運算個六七次後x的值就已經相當精確了。
例如,我想求根號2等於多少。假如我猜測的結果為4,雖然錯的離譜,但你可以看到使用牛頓反覆運算法後這個值很快就趨近於根號2了:
(4 + 2/4) / 2 = 2.25
(2.25 + 2/2.25) / 2 = 1.56944..
(1.56944.. + 2/1.56944..) / 2 = 1.42189..
(1.42189.. + 2/1.42189..) / 2 = 1.41423..
....
這種演算法的原理很簡單,我們僅僅是不斷用(x,f(x))的切線來逼近方程x^2-a=0的根。根號a實際上就是x^2-a=0的一個正實根,這個函數的導數是2x。也就是說,函數上任一點(x,f(x))處的切線斜率是2x。那麼,x-f(x)/(2x)就是一個比x更接近的近似值。代入 f(x)=x^2-a得到x-(x^2-a)/(2x),也就是(x+a/x)/2。
  相關的代碼如下:

float SqrtByNewton(float x)
{
    float val = x; //最终
    float last; //保存上一個計算的值
    do
    {
        last = val;
        val = (val + x/val) / 2;
    }while(abs(val-last) > eps);
    return val;
}

  然後我們再來看下性能測試:
  哇塞,性能提高了很多,可是和系統函數相比,還是有這麼大差距,這是為什麼呀?想啊想啊,想了很久仍然百思不得其解。突然有一天,我在網上看到一個神奇的方法,於是就有了今天的這篇文章,廢話不多說,看代碼先:

float InvSqrt(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
    x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
    x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
    return 1/x;
}

  然後我们最後一次来看下性能測試:
  這次真的是質變了,結果竟然比系統的還要好。。。哥真的是震驚了!!!哥吐血了!!!一個函數引發了血案!!!血案,血案。。。
 到現在你是不是還不明白那個“鬼函數”,到底為什麼速度那麼快嗎?不急,先看看下面的故事吧:
Quake-III Arena (雷神之錘3)是90年代的經典遊戲之一。該系列的遊戲不但畫面和內容不錯,而且即使電腦配置低,也能極其流暢地運行。這要歸功於它3D引擎的開發者約翰·卡馬克(John Carmack)。事實上早在90年代初DOS時代,只要能在PC上搞個小動畫都能讓人驚歎一番的時候,John Carmack就推出了石破天驚的Castle Wolfstein,然後再接再勵,doom, doomII, Quake...每次都把3-D技術推到極致。他的3D引擎代碼資極度高效,幾乎是在壓榨PC機的每條運算指令。當初MS的Direct3D也得聽取他的意見,修改了不少API。
最近,QUAKE的開發商ID SOFTWARE 遵守GPL協議,公開了QUAKE-III的原代碼,讓世人有幸目睹Carmack傳奇的3D引擎的原碼。這是QUAKE-III原代碼的下載位址:
http://www.fileshack.com/file.x?fid=7547
(下麵是官方的下載網址,搜索 “quake3-1.32b-source.zip” 可以找到一大堆中文網頁的。ftp://ftp.idsoftware.com/idstuff/source/quake3-1.32b-source.zip)
我們知道,越底層的函數,調用越頻繁。3D引擎歸根到底還是數學運算。那麼找到最底層的數學運算函數(在game/code/q_math.c), 必然是精心編寫的。裡面有很多有趣的函數,很多都令人驚奇,估計我們幾年時間都學不完。在game/code/q_math.c裡發現了這樣一段代碼。它的作用是將一個數開平方並取倒,經測試這段代碼比(float)(1.0/sqrt(x))快4倍:
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?
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
    //y  = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration,
     this can be removed

    #ifndef Q3_VM
    #ifdef __linux__
         assert( !isnan(y) ); // bk010122 - FPE?
    #endif
    #endif
    return y;
}  
函數返回1/sqrt(x),這個函數在影像處理中比sqrt(x)更有用。 
注意到這個函數只用了一次疊代!(其實就是根本沒用疊代,直接運算)。編譯,實驗,這個函數不僅工作的很好,而且比標準的sqrt()函數快4倍!要知道,編譯器自帶的函數,可是經過嚴格仔細的彙編優化的啊!
這個簡潔的函數,最核心,也是最讓人費解的,就是標注了“what the fuck?”的一句:
 i = 0x5f3759df - ( i >> 1 );
再加上:
y  = y * ( threehalfs - ( x2 * y * y ) ); 
兩句話就完成了開方運算!而且注意到,核心那句是定點移位元運算,速度極快!特別在很多沒有乘法指令的RISC結構CPU上,這樣做是極其高效的。
演算法的原理其實不複雜,就是牛頓反覆運算法,用x-f(x)/f'(x)來不斷的逼近f(x)=a的根。
沒錯,一般的求平方根都是這麼迴圈反覆運算算的,但是卡馬克(quake3作者)真正牛B的地方是他選擇了一個神秘的常數 0x5f3759df 來計算那個猜測值,就是我們加注釋的那一行,那一行算出的值非常接近1/sqrt(n),這樣我們只需要2次牛頓反覆運算就可以達到我們所需要的精度。好吧如果這個還不算NB,接著看:
普渡大學的數學家Chris Lomont看了以後覺得有趣,決定要研究一下卡馬克弄出來的這個猜測值有什麼奧秘。Lomont也是個牛人,在精心研究之後從理論上也推導出一個最佳猜測值,和卡馬克的數字非常接近,0x5f37642f。卡馬克真牛,他是外星人嗎?
傳奇並沒有在這裡結束。Lomont計算出結果以後非常滿意,於是拿自己計算出的起始值和卡馬克的神秘數字做比賽,看看誰的數字能夠更快更精確的求得平方根。結果是卡馬克贏了... 誰也不知道卡馬克是怎麼找到這個數字的。
最後Lomont怒了,採用暴力方法一個數位一個數位試過來,終於找到一個比卡馬克數字要好上那麼一丁點的數字,雖然實際上這兩個數字所產生的結果非常近似,這個暴力得出的數字是0x5f375a86。
Lomont為此寫下一篇論文,"Fast Inverse Square Root"。 論文下載地址:
http://www.math.purdue.edu/~clomont/Math/Papers/2003/InvSqrt.pdf
http://www.matrix67.com/data/InvSqrt.pdf
參考:<IEEE Standard 754 for Binary Floating-Point Arithmetic><FAST INVERSE SQUARE ROOT>
最後,給出最精簡的1/sqrt()函數:
float InvSqrt(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;
}  
大家可以嘗試在PC機、51、AVR、430、ARM、上面編譯並實驗,驚訝一下它的工作效率。
前兩天有一則新聞,大意是說 Ryszard Sommefeldt 很久以前看到這麼樣的一段 code (可能出自 Quake III 的 source code):
float InvSqrt (float x) 
{
    float xhalf = 0.5f*x;
    int i = *(int*)&x;
    i = 0x5f3759df - (i>>1);
    x = *(float*)&i;
    x = x*(1.5f - xhalf*x*x);
    return x;
}
他一看之下驚為天人,想要拜見這位前輩高人,但是一路追尋下去卻一直找不到人;同時間也有其他人在找,雖然也沒找到出處,但是 Chris Lomont 寫了一篇論文 (in PDF) 解析這段 code 的演算法 (用的是 Newton’s Method,牛頓法;比較重要的是後半段講到怎麼找出神奇的 0x5f3759df 的)。 
PS. 這個 function 之所以重要,是因為求“開根號倒數”這個動作在 3D 運算 (向量運算的部份) 裡面常常會用到,如果你用最原始的 sqrt() 然後再倒數的話,速度比上面的這個版本大概慢了四倍吧… XD
PS2. 在他們追尋的過程中,有人提到一份叫做 MIT HACKMEM 的檔,這是 1970 年代的 MIT 強者們做的一些筆記 (hack memo),大部份是 algorithm,有些 code 是 PDP-10 asm 寫的,另外有少數是 C code (有人整理了一份列表)
  好了,故事就到這裡结束了,希望大家能有收獲:)

p.s.感謝有人提供另一位文章,做了各種sqrt的速度及正確率比較:
http://www.codeproject.com/Articles/69941/Best-Square-Root-Method-Algorithm-Function-Precisi

2 則留言: