《微微一笑很傾城》中肖奈大神說的平方根倒數速算法是什麼鬼？三十分鐘理解！

時間 2019-12-12

標籤微微一笑很傾城平方根倒數速算什麼三十分 30分理解简体版

原文原文鏈接

在電影《微微一笑很傾城》中，肖奈大神在玻璃上寫了一堆公式，提到平方根倒數速算算法，這個究竟是一個什麼算法？筆者看電影的時候打開手機學了一下，發現該算法的做者真乃神人！今天有空，就把該算法寫一寫。html

在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？編程

--------------------------------------------------------------------------------------------------------------------------------------------------app

該算法的本質其實就是牛頓迭代法（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的解；將該方程按牛頓迭代法的公式展開爲：ui

x[n+1]=x[n]*(3/2-a/2*x[n]*x[n])
spa

也就是上面代碼藍色部分的第L3行；因此很明顯，這段代碼就是進行一次迭代的牛頓迭代法。只進行一次迭代就結束程序，又要保證精度，也就是說，牛頓迭代法的初始值要很是精準，才能在一次迭代後完成計算。整個代碼最精彩的就是L1行了，i = 0x5f3759df - (i >> 1);究竟是什麼意思？爲何初始值能夠這樣算呢？瞭解該代碼，須要先了解一些float point和fix point的表示方法[2]：.net

一個浮點數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）,圖中的 $\scriptstyle m=1\times 2^{-2}=0.250$ 。

因此浮點數的結構公式爲： $\scriptstyle x=(-1)^{\mathrm{Si}}\cdot(1+m)\cdot 2^{(E-B)}$ ，圖中 $\scriptstyle x=(1+0.250)\cdot 2^{-3}=0.15625$

整數fix point的表示相對簡單，符號佔1位，數值佔剩下的31位。若是用上圖的浮點數字節序列來表示整數，那麼 $\scriptstyle I=E\times 2^{23}+M$ ，即 $I=124\times 2^{23} + 2^{21}$ ；平方根倒數函數僅能處理正數，因此符號位均爲0。

對於一樣的32位二進制數碼，若爲浮點數表示時實際數值爲 $\scriptstyle x=(1+m_x)2^{e_x}$ ，而若爲整數表示時實際數值則爲 $\scriptstyle I_x=E_xL+M_x$ ，其中 $\scriptstyle L=2^{n-1-b}$ ，這裏n=32，b=8。式子中引入的新變量爲：

$m_x=\frac{M_x}{L}$ ------------------------------------等式1

$e_x=E_x-B$ ，其中 $B=2^{b-1}-1$ ------------等式2

理解浮點數和整數的表示後，下面開始推導。

x的平方根倒數方程爲：

$y=\frac{1}{\sqrt{x}}$

兩邊取對數有：

$\log_2{(y)}=-\frac{1}{2}\log_2{(x)}$

由於浮點數可表示爲： $\scriptstyle x=(1+m_x)2^{e_x}$ ，因此也有，代入上式有：

$\log_2(1+m_y)+e_y=-\frac{1}{2}\log_2{(1+m_x)}-\frac{1}{2}e_x$

因爲 $\scriptstyle 0 \le x < 1$ ，有 $\scriptstyle \log_2{(1+x)}\approx {x}$ ，則在此可定義 $\sigma$ 與x的關係爲 $\scriptstyle \log_2{(1+x)}\cong x+\sigma$ ， $\sigma$ 表示偏差，因此將 $\scriptstyle \log_2{(1+x)}= x+\sigma$ 代入上式得：

$m_y+\sigma+e_y=-\frac{1}{2}m_x-\frac{1}{2}\sigma-\frac{1}{2}e_x$

將等式1，等式2代入到上述方程中，有：

$M_y+(E_y-B)L=-\frac{3}{2}\sigma{L}-\frac{1}{2}M_x-\frac{1}{2}(E_x-B)L$

移項整理得：

$E_yL+M_y=\frac{3}{2}(B-\sigma)L-\frac{1}{2}(E_xL+M_x)$

又由於浮點規格存儲的正浮點數x，若將其做爲整數表示，則示值爲： $\scriptstyle I_x=E_xL+M_x$ ，因此x的平方根倒數的首次近似值的整數表示值爲：

$I_y=E_yL+M_y=R-\frac{1}{2}(E_xL+M_x)=R-\frac{1}{2}I_x$ ，

其中 $R=\frac{3}{2}(B-\sigma)L$ ， $B=2^{b-1}-1$ ， $\scriptstyle L=2^{n-1-b}$ ，n=32，b=8。

這個式子對應着源代碼中的這一行：i = 0x5f3759df - ( i >> 1 );，而後將整數表示值換回表示浮點數：x = * ( float * ) &i;。這樣就獲得了浮點數的平方根倒數的近似值。

0x5f3759df 對應着R，即3/2(B- $\sigma$ )L.當R爲0x5f3759df時，有 $\scriptstyle \sigma=0.0450461875791687011756$ .

"如今不只該算法的原做者不明，人們也仍沒法明確當初選擇這個「魔術數字」的方法。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