【轉】神一樣的算法，一個數開平方並取倒，0x5f3759df

gameplayer

http://www.douban.com/note/93460299/

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的根。

簡單來說比如求平方根,f(x)=x^2=a ,f'(x)= 2*x,f(x)/f'(x)=x/2,把f(x)代入

x-f(x)/f'(x)後有(x+a/x)/2，現在我們選a=5,選一個猜測值比如2，
那麼我們可以這麼算
5/2 = 2.5; (2.5+2)/2 = 2.25; 5/2.25 = xxx; (2.25+xxx)/2 = xxxx ...
這樣反覆迭代下去，結果必定收斂於sqrt(5)，沒錯，一般的求平方根都是這麼算的
但是卡馬克(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
點擊此處下載 ourdev_714253IQ3KPV.pdf(文件大小:148K) (原文件名: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、上面編譯並實驗，驚訝一下它的工作效率。


==============增加================
看matrix67大俠的解釋
http://www.matrix67.com/blog/archives/362

這樣的代碼速度肯定飛快，我就不用多說了；但算法的原理是什麼呢？其實說穿了也不是很神，程序首先猜測了一個接近1/sqrt(number)的值，然後兩次使用牛頓迭代法進行迭代。根號a的倒數實際上就是方程1/x^2 - a = 0的一個正實根，它的導數是-2/x^3。運用牛頓迭代公式x' = x - f(x)/f'(x)，式子化簡為x' = x * (1.5 - 0.5a * x^2)。迭代幾次後，x的值將趨於1/sqrt(a)。
    但這段代碼真正牛B的是那個神秘的0x5f3759df，因為0x5f3759df - (i >> 1)出人意料地接近根號y的倒數。人們都不知道這個神秘的常數是怎麼來的，只能把它當作神來膜拜。這個富有傳奇色彩的常數到底咋回事，很少有人說得清楚。這篇論文（就是上面的附件）比較科學地解釋了這個常數。