《微微一笑很傾城》中肖奈大神說的平方根倒數速演算法是什麼鬼?三十分鐘理解!
在電影《微微一笑很傾城》中,肖奈大神在玻璃上寫了一堆公式,提到平方根倒數速算演算法,這個到底是一個什麼演算法?筆者看電影的時候開啟手機學了一下,發現該演算法的作者真乃神人!今天有空,就把該演算法寫一寫。
在3D圖形程式設計中,經常要求平方根的倒數,即1/Sqrt(x),如果用一般的程式碼(float)(1.0/sqrt(x)),,精度高,但是非常慢;我們需要一個快速,而又足夠高精度的演算法;著名遊戲《雷神之錘III》的程式碼在2002年左右被披露,人們發現了一段用於快速計算平方根倒數的程式碼,下面是整理後的程式碼(去掉了一些巨集定義)。
float InvSqrt (float x) { float xhalf = 0.5f*x; int i = *(int*)&x; // get bits for floating value i = 0x5f3759df - (i >> 1);// L1:gives initial guess y0 x = *(float*)&i; //l2:convert bits back to float x = x*(1.5f - xhalf*x*x); //l3:Newton step, repeating increases accuracy return x; }
該程式執行效率極高,經測試,基本是使用直接開根號求倒數程式的4倍速度!一時間驚為天人。那麼這段程式碼到底怎麼理解?為什麼中間出現了0x5f3759df這樣一個完全無厘頭的magic number?
--------------------------------------------------------------------------------------------------------------------------------------------------
該演算法的本質其實就是牛頓迭代法(Newton-Raphson Method,簡稱NR)。NR是一種求方程的近似根的方法。首先要估計一個與方程的根比較靠近的數值,然後根據公式推算下一個更加近似的數值,不斷重複直到可以獲得滿意的精度。其公式如下:
函式:y=f(x)
其一階導數為:y'=f'(x)
則方程:f(x)=0 的第n+1個近似根為
x[n+1] = x[n] - f(x[n]) / f'(x[n])
--------------------------------------------------------------------------------------------------------------------------------------------------
求一個數a的平方根的倒數,實際就是求方程f(x)=1/(x^2)-a=0的解;將該方程按牛頓迭代法的公式展開為:
x[n+1]=x[n]*(3/2-a/2*x[n]*x[n])
也就是上面程式碼藍色部分的第L3行;所以很明顯,這段程式碼就是進行一次迭代的牛頓迭代法。只進行一次迭代就結束程式,又要保證精度,也就是說,牛頓迭代法的初始值要非常精準,才能在一次迭代後完成計算。整個程式碼最精彩的就是L1行了,i = 0x5f3759df - (i >> 1);到底是什麼意思?為什麼初始值可以這樣算呢?瞭解該程式碼,需要先了解一些float point和fix point的表示方法[2]:
一個浮點數float是由32位二進位制位表示的有理數,分為三部分。其中符號佔1位,表示正負,記為Si;指數佔接下來的8位,表示經過偏移處理後的指數,即實際表示E(如圖中為124),需要偏移B(圖中為2的8次方減1,127。B為一個固定值),最後得指數值為E-B;有效數字(除最高位1以外,因為前面有指數,所以一個數肯定可以表示成1.xxxxx * 2^(kkk),只保留xxxxx)佔剩下的23位,記為m(0<m<1),圖中的。
所以浮點數的結構公式為:, 圖中
整數fix point的表示相對簡單,符號佔1位,數值佔剩下的31位。如果用上圖的浮點數字節序列來表示整數,那麼,即;平方根倒數函式僅能處理正數,所以符號位均為0。
對於同樣的32位二進位制數碼,若為浮點數表示時實際數值為,而若為整數表示時實際數值則為,其中,這裡n=32,b=8。式子中引入的新變數為:
------------------------------------等式1
,其中------------等式2
理解浮點數和整數的表示後,下面開始推導。
x的平方根倒數方程為:
兩邊取對數有:
因為浮點數可表示為:,所以也有,代入上式有:
由於,有,則在此可定義與x的關係為,表示誤差,所以將代入上式得:
將等式1,等式2代入到上述方程中,有:
移項整理得:
又因為浮點規格儲存的正浮點數x,若將其作為整數表示,則示值為:,所以x的平方根倒數的首次近似值的整數表示值為:
,
其中,,,n=32,b=8。
這個式子對應著原始碼中的這一行:i = 0x5f3759df - ( i >> 1 );,然後將整數表示值換回表示浮點數:x = * ( float * ) &i;。這樣就得到了浮點數的平方根倒數的近似值。
0x5f3759df 對應著R,即3/2(B-)L.當R為0x5f3759df時,有.
"現在不僅該演算法的原作者不明,人們也仍無法明確當初選擇這個“魔術數字”的方法。Chris Lomont在研究中曾做了個試驗:他編寫了一個函式,以在一個範圍內遍歷選取R值的方式將逼近誤差降到最小,以此方法他計算出了線性近似的最優R值0x5f37642f(與程式碼中使用的0x5f3759df相當接近),但以之代入演算法計算並進行一次牛頓迭代後,所得近似值與代入0x5f3759df的結果相比精度卻仍略微更低。……在Charles McEniry的論文中,他使用了一種類似Lomont但更復雜的方法來優化R值:他最開始使用窮舉搜尋,所得結果與Lomont相同;而後他嘗試用帶權二分法尋找最優值,所得結果恰是程式碼中所使用的魔術數字0x5f3759df"---維基百科
參考資料:
[1] CHRIS LOMONT, FAST INVERSE SQUARE ROOT
[2] http://blog.csdn.net/heloowird/article/details/21862251
[3] http://blog.chinaunix.net/uid-9255716-id-107951.html
[4] http://blog.csdn.net/xiaoguohaha/article/details/21652643