《微微一笑很傾城》中肖奈大神說的平方根倒數速演算法是什麼鬼？三十分鐘理解！
在電影《微微一笑很傾城》中，肖奈大神在玻璃上寫了一堆公式，提到平方根倒數速算演算法，這個到底是一個什麼演算法？筆者看電影的時候開啟手機學了一下，發現該演算法的作者真乃神人！今天有空，就把該演算法寫一寫。
在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的平方根倒數的首次近似值的整數表示值為：
，
其中
，
這個式子對應著原始碼中的這一行：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