快速平方根倒數算法深度了解
快速平方根倒數算法是什麼？
簡單來說這個算法避開了開方和除法運算快速實作了
y = 1 x y= \frac{1}{\sqrt x} y=x
1
快速平方根倒數算法首次出現是在著名的雷神之錘3的圖形計算優化代碼中
話不多說先上c代碼
一.代碼
小白的我以為這麼寫的：
float y = 1/sqrt(x);
開發遊戲的大佬是這麼寫的：
float Q_rsqrt(float number)
{
long i;
float x2,y;
const float threehalfs = 1.5F;
x2= number * 0.5F;
y= number;
i= * (long * ) &y;
i= 0x5f3759df-(i>>1);
y= * (float * ) &i;
//上述代碼給牛頓疊代法提供了一個較好的初始值
y=y * (threehalfs - (x2*y*y)); //牛頓疊代一次
//y=y * (threehalfs - (x2*y*y)) 可持續疊代
}
經大佬測試代碼比(float)(1/sqrt(x))快4倍。
初步觀看代碼：
①定義了32位整型變量i
②定義了兩個32位浮點數
③将1.5存入一個常量
④将形參number除2存入x2
⑤将形參number存入y
之後的四行代碼乍一看？？？莫慌下面開始正文
二.了解算法前的知識鋪墊
1.IEEE 754 （二進制浮點數算術标準）
32位分為三個部分：符号位 指數部分
①符号位：0為正數 1為負數
②8位指數部分：8位可以表示0到255所有的數字，但我們還需要得到負指數。
是以我們對0-255減去127就可以得到-127到128之間所有的整數。
舉個例子如果指數部分我想要4
那這八位必須表示4+127=131 即10000011
這樣通過移碼我們可以将正負指數都變成0-255之間的整數
③23位尾數：利用23位我們可以表示0到2^23-1之間的整數
對于十進制科學計數法，尾數是1到10
對于二進制科學計數法，尾數是1-2
11000 = 1.1 ∗ 2 4 11000=1.1*{2}^{4} 11000=1.1∗24
0.0101 = 1.01 ∗ 2 − 3 0.0101=1.01*{2}^{-3} 0.0101=1.01∗2−3
在科學計數法中，第一位按照定義不可以是0，對于二進制這個數隻能是1。是以我們不存儲這一位。
說了一堆理論的我們舉一個例子看一下浮點數是怎麼存儲在計算機中的：
我們拿 8.25 舉例子
Ⅰ.寫成二進制：1000.01
Ⅱ.變成科學計數法： 1.00001 ∗ 2 3 = ( 1 + 0.00001 ) ∗ 2 3 1.00001 *{2^3}=(1 + 0.00001) * {2^3} 1.00001∗23=(1+0.00001)∗23
Ⅲ.解釋
指數部分我們用 8bit 表示 3+127
3對應的就是2^3中的3
(1 + 0.00001)：加号前面的1我們丢掉不存
對于0.00001中小數點後的數我們用23位尾數去存
羅嗦了這麼久現在我們可以輕松了解以下内容：
對于一個浮點數x(十進制或者二進制)有：
x = ( − 1 ) s ∗ ( 1 + m ) ∗ 2 e x=(-1)^s*(1+m)*{2}^{e} x=(−1)s∗(1+m)∗2e
對于下式 ( 1 + 0.00001 ) ∗ 2 3 (1 + 0.00001) * {2^3} (1+0.00001)∗23
m : 00001 m:00001 m:00001
e : 3 e:3 e:3
我們将存取之後的8位指數部分用E表示
将存取之後的23位尾數部分用M表示
E : 3 + 127 E:3+127 E:3+127
M : 00001 ⋯ （ 補 18 個 0 填 滿 23 位 ） M:00001\cdots（補18個0填滿23位） M:00001⋯（補18個0填滿23位）
我們可以得出E e和M m的關系式
E = e + 127 E=e+127 E=e+127
M = 2 23 ∗ m M={2}^{23}*m M=223∗m
三.數學推導
y = 1 x y= \frac{1}{\sqrt x} \quad y=x
1
取對數：
① log 2 y = − 0.5 ∗ log 2 ①\log_2 y=-0.5* \log_2 ①log2y=−0.5∗log2
将y和x分别用二進制科學計數法表示：
② y = ( 1 + m y ) ∗ 2 e y ②y=(1+m_y)*{2}^{e_y} ②y=(1+my)∗2ey
③ x = ( 1 + m x ) ∗ 2 e x ③x=(1+m_x)*{2}^{e_x} ③x=(1+mx)∗2ex
聯合上述三個式子化簡：
l o g 2 ( 1 + m y ) + e y = − 0.5 ∗ ( l o g 2 ( 1 + m x ) + e x ) log_2 (1+m_y)+e_y=-0.5*(log_2(1+m_x)+e_x) log2(1+my)+ey=−0.5∗(log2(1+mx)+ex)
通過一次線性近似：
m y + b + e y = − 0.5 ∗ ( m x + b + e x ) m_y+b+e_y=-0.5*(m_x+b+e_x) my+b+ey=−0.5∗(mx+b+ex)
m y = M y 2 23 m_y= \frac{M_y}{{2}^{23}} \quad my=223My
m x = M x 2 23 m_x= \frac{M_x}{{2}^{23}} \quad mx=223Mx
E y = e y + 127 E_y=e_y+127 Ey=ey+127
E x = e x + 127 E_x=e_x+127 Ex=ex+127
上述五個式化簡：
M y + 2 23 ∗ E y = 1.5 ∗ 2 23 ∗ ( 127 − b ) − 0.5 ∗ ( M x + 2 23 ∗ E x ) M_y+{2}^{23}*E_y=1.5*{2}^{23}*(127-b)-0.5*(M_x+{2}^{23}*E_x) My+223∗Ey=1.5∗223∗(127−b)−0.5∗(Mx+223∗Ex)
令
I y = M y + 2 23 ∗ E y I_y=M_y+{2}^{23}*E_y Iy=My+223∗Ey
I x = M x + 2 23 ∗ E x I_x=M_x+{2}^{23}*E_x Ix=Mx+223∗Ex
對于 1.5 ∗ 2 23 ∗ ( 127 − b ) 1.5*{2}^{23}*(127-b) 1.5∗223∗(127−b)我們令其為K，可以知道通過調節b的大小就可以改變K的值
前輩們給出了當
b = 0.0450465 b=0.0450465 b=0.0450465時有K=
K = 0 x 5 f 3759 d f K=0x5f3759df K=0x5f3759df
這剛好對應代碼中：
通過一系列不怎麼複雜的數學推導我們究竟得到了什麼?
我們知道了一個輸出結果的整形和輸入整形之間的關系，即：
I y = K − 0.5 ∗ I x I_y=K-0.5*I_x Iy=K−0.5∗Ix
四.代碼解釋
①
long i;
float x2,y;
const float threehalfs = 1.5F;
x2= number * 0.5F;
y= number;
i= * (long * ) &y;
我們将數字number存入y中準備對其進行位操作，但是衆所周知浮點數沒有位操作的語句。
但是長整型數可以進行位操作！
如果執行下面這一行代碼：
i=long(y);
y=4.456時i=4
但是我們想要的是将y逐位轉化成長整型，換句話說按照上面的代碼我們改變了二進制逐位的值。
是以我們采取下面代碼：
這句話幹了什麼呢？
簡單來說它轉換了y對應的記憶體位址而不是y這個數字本身。
或者說這句話告訴cpu把y指向的記憶體資料當成長整形來讀取
&y 取到浮點數y的位址
(long * ) &y 将該位址轉化成長整型
*(long * ) &y 再取其值
6666指針确實是個神奇的東西，我們繼續向下
通過這句話浮點數y的二進制表示存入到了i中。
②
這句的了解在第三大部分可以輕松找到。
簡單來說我們通過放縮系數和平移操作代替了除法和開方的過程，這大大提高了代碼運作效率。
③
和①類似我們再次通過這種操作将二進制位轉換成浮點數。
通過上述三步我們得到了一個近似值。是以我們還需要進行牛頓疊代法繼續使結果精确。
④
下面進行牛頓疊代法的數學推導：
y = 1 x y= \frac{1}{\sqrt x} y=x
1
f ( y ) = 1 y 2 − x f(y)=\frac{1}{y^2}-x f(y)=y21−x
y n + 1 = y n − 1 y n 2 − x − 2 y n 3 = y n ∗ ( 3 2 − 1 2 ∗ x ∗ y n 2 ) y_{n+1}=y_n-\frac{\frac{1}{y_n^2}-x}{\frac{-2}{y_n^3}}=y_n*(\frac{3}{2}-\frac{1}{2}*x*y_n^2) yn+1=yn−yn3−2yn21−x=yn∗(23−21∗x∗yn2)
五.總結
短短的幾行代碼背後是強大的數學支撐，萬萬想不到之前學過的高數會在算法中這樣展現。
第一次試寫涉及到數學公式推導的文章，如有疑問請大家多多留言！
![查找算法之二分查找查找算法之二分查找[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)