一、寫在前面

1. 提示

請朋友們先行閱讀Fast inverse square root algorithm(下稱FISR)演算法，本文內容適用于已經了解此演算法基本原理的朋友們，還不了解的讀者請移步以下參考資料（推薦程度按照以下順序）:

(1) 中文且強推: 魔力數字與快速平方根倒數演算法
(2) Wikipedia官方資料(需要翻墻) : Wikipedia Fast inverse square root
(3) 中文且做補充資料: 數學之美：平方根倒數速演算法中的神奇數字 0x5f3759df
(4) 英文paper原文: InvSqrt.pdf(CHRIS LOMONT)

2. 背景與前情

大學非常厲害的舍友面試阿里的后端時，面試官讓他不用函式庫，實作一個快速開平方根的演算法，當舍友把這個問題拋給我的時候，我滿腦子想的都是：二分法！
顯然這是大家都能想得到的演算法，其時間復雜度為 O ( l o g N ) O(logN) O(logN)，更準確點，其實需要做的操作是 l o g 2 ( X e r r o r ) log_2(\cfrac{X}{error}) log2?(errorX?)次，X為被開平方的數，error是預設的可允許誤差，

But！ 據說面試官聽到舍友的二分法后，略帶失望的搖了搖頭(有演繹的成分)，說：再想想，經過5分鐘的煎熬，最后面試官妥協了：“那用二分法做吧，其實這個題有線性時間解法”，

這個線性時間的解法一下子就讓我充滿了興趣，開平方根，怎么可能有線性時間的演算法？？？第一個反應是Impossible，第二個反應是竭盡腦汁思考除了二分法以外的其他方法，在經過長達10分鐘的思考后，我終于…
明白了一句話：終日而思矣，不如須臾之所學也！好吧，上網搜得了！

最終的最終，經過一番學習以后，我現在腦子里對此題有了三種解法，

• 二分法(正常學過演算法的都應該秒想到的)
• 牛頓迭代法(數值計算學過，應該屬于正常人能想得到的非常好用而且有效的方法)
• Fast inverse sqaure root后再取倒數，也就是本文的核心

接下來，我試圖從自己理解的角度，闡述一下對于FISR的理解，

二、正文

先把FISR代碼貼在這兒，注意注釋為what the fuck的那一行：

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

	return y;
}

1. 需求分析

問題: 求 y = 1 x y=\cfrac{1}{\sqrt{x}} y=x ?1?
對此問題，我們變變形，方便以下計算 y = 1 x = x ? 1 / 2 = x 0 × x ? 1 / 2 (1) y=\cfrac{1}{\sqrt{x}}=x^{-1/2}=x^0 \times x^{-1/2} \tag{1} y=x ?1?=x?1/2=x0×x?1/2(1)

已有工具: 因為不讓用函式庫，因此我們有的就是 + , ? , × , ÷ +, -, \times, \div +,?,×,÷四種基礎工具，
聯想: 如果能把平方這種級別的運算“降維”變成 + , ? +, - +,?就好了，怎么做呢？
? \longrightarrow ?: 數學上最常用的手段就是取對數 l o g log log，比如 l o g ( a b × a c ) = b l o g a + c l o g a = ( b + c ) l o g a log(a^b \times a^c)=bloga + cloga = (b+c) loga log(ab×ac)=bloga+cloga=(b+c)loga，
那么，如果我們對(1)式取對數，就降維成了
l o g y = l o g x 0 ? 1 2 l o g x (2) logy=logx^0-\cfrac{1}{2}logx \tag{2} logy=logx0?21?logx(2) 那么，開平方根倒數的計算就變成了 ? 1 2 - \cfrac{1}{2} ?21?的計算，那不就大功告成了！！！
事實上，觀察(2)，它跟以下的方程很像 I ( y ) = R ? 1 2 I ( x ) (3) I(y) = R - \cfrac{1}{2}I(x) \tag{3} I(y)=R?21?I(x)(3) 或者是那一行代碼非常之像
i = 0 x 5 f 3759 d f ? ( i > > 1 ) (4) i = 0x5f3759df - ( i >> 1 ) \tag{4} i=0x5f3759df?(i>>1)(4) 注意(4)中等是左邊的 i i i其實就是 I ( y ) I(y) I(y), 神奇數字 0 x 5 f 3759 d f 0x5f3759df 0x5f3759df其實就是 R R R，最后 ? ( i > > 1 ) -(i>>1) ?(i>>1)其實就是 ? 1 2 i - \cfrac{1}{2} i ?21?i，

是不是感覺有點思路了？有了上面非常相似的(2)(3)式，我們現在關心的有兩點：1.是求 I I I和 a r g I argI argI，因為我們需要把 x x x變為 I ( x ) I(x) I(x)以及把 I ( y ) I(y) I(y)還原為 y y y，2.是求R所對應的神奇數字，

2. 必備工具之IEEE-754浮點數表示方法

相應的Wikipedia的介紹和例子放在下面(英文不好的同學們去別的地方搜一下吧，這部分的基礎原理不在本文介紹之中)，

注意到，我們被開平方根的數字是正數，所以符號位S=0在本文中永遠成立，接下來不再討論，此外，下文中簡稱IEEE-754的32 bits表述一個浮點數的方法為SEM法(S符號位：一位二進制數，E指數位：八位二進制數，M尾數位：二十三位二進制數)，

下面的表格對比了：對于同一個浮點數正數x的不同表示方法，方便接下來討論，注意，沒有理解此表格的朋友請理解好了再往下走，切勿囫圇吞棗，

為了以下方便表示，我們引入兩個常量
B = 127 , L = 2 23 B=127, L=2^{23} B=127,L=223 至于原因，是因為IEEE 32位浮點數表示中，E有8位，所以 B = 2 8 ? 1 = 127 B=2^8-1=127 B=28?1=127， 而M有23位，所以 L = 2 23 L=2^{23} L=223，

3. 同一儲存單元32bits的兩種不同意義

以上的表格中，我們詳細的對比了對于同一個浮點數x，1.人類或是2.計算機是如何理解、儲存及還原它的，
然而，不僅x的理解方式有兩種，同樣有兩種理解方式的還有一個32 bits的數字，當我們在計算機的存盤單元中發現了一個32位的二進制數，在沒有被告知它的型別是1.整形或是2.浮點型時，我們是并不能判斷次32 bits的數字到底代表什么意義的，

下面的表格就對比了同一個32 bits的數字的兩種不同釋義，注意，我們可以把這個32 bits數字劃分為三部分：S是這個32bits的從左起第0位，E是第1-第8位，M是第9-第31位，E, M的意義完全跟上面相同，同樣的還有B,L，我們都延續第2部分的概念，
S 0 E 1 E 2 . . . E 7 E 8 M 9 M 10 . . . M 30 M 31 S_0 E_1 E_2 ...E_7 E_8 M_{9} M_{10} ...M_{30}M_{31} S0?E1?E2?...E7?E8?M9?M10?...M30?M31?

注意到， 1. 這里的 F l o a t = a r g I Float = argI Float=argI，即二者的函式(function，也可理解為功能)是一樣的，2. 注意整數的 I n t ( E , M ) Int(E, M) Int(E,M)和浮點數 F l o a t ( E , M ) Float(E, M) Float(E,M)的形式都很相似，都是 ( a + M L ) × b (a+\cfrac{M}{L}) \times b (a+LM?)×b的形式，

有了上面注意點里面的1、2兩點內容，我們就會思考，這是巧合么？我們又能利用這兩點巧合做到什么么？ 那么接下來的重頭戲即將呈現！！！不過在此之前，我們停下來回顧一下上面的準備作業都準備了哪些：

1. 一個任意的浮點數x都可以被Encoding成為兩種形式，一種是我們人腦作為encoder記住其e, m；另一種是計算機作為encoder把此浮點數按SEM方法存盤成一個32bits的二進制數存入存盤單元， 注意，這兩種Encoding的形式都是一一對應(雙射)的！！！
2. 一個任意的位于存盤單元中的32bits 二進制數字t，我們都可以將其Decoding為兩個不同的數字，一種方法是告訴我們x是整數，那我們就要按照其為整數進decoding；另一種是告訴我們t是浮點數，那我們就按照其為浮點數進行decoding，同樣，這兩種Decoding的方法也都是一一對應(雙射)的！！！

4. 公式推導

掌握以上知識，我們試圖著手將輸入的自變數x從乘方級運算(開平方根就是乘 ? 1 2 -\cfrac{1}{2} ?21?次方)通過映射變為 I ( x ) I(x) I(x)后變為加減乘除級別運算，此處的x在計算機是以SEM方法存盤，因此通過回顧第2部分的表格 x = a r g I ( E , M ) = ( 1 + M L ) × 2 E ? B (5) x = argI(E, M) = (1 + \cfrac{M}{L}) \times 2^{E - B} \tag{5} x=argI(E,M)=(1+LM?)×2E?B(5)
以2為底取對數 l o g 2 log_2 log2?
l o g 2 ( x ) = l o g 2 [ ( 1 + M L ) × 2 E ? B ] log_2(x) = log_2[(1 + \cfrac{M}{L}) \times 2^{E - B}] log2?(x)=log2?[(1+LM?)×2E?B] = ( E ? B ) + l o g 2 ( 1 + M L ) =(E-B) + log_2(1+\cfrac{M}{L}) =(E?B)+log2?(1+LM?) = [ E + l o g 2 ( 1 + M L ) ] ? B (6) =[E + log_2(1+\cfrac{M}{L})] - B \tag{6} =[E+log2?(1+LM?)]?B(6)
OK，稍稍稍微又卡住了，又推不下去了，那么我們需要一點近似，如果把 M L \cfrac{M}{L} LM?當做一個整體 m m m, m就是剛才第二部分的m，其意義是二進制科學計數法的小數部分, m ∈ [ 0 , 1 ) m\isin [0, 1) m∈[0,1)，那么 F 1 ( m ) = l o g 2 ( 1 + m ) F_1(m)=log_2(1+m) F1?(m)=log2?(1+m)與 F 2 ( m ) = m F_2(m)=m F2?(m)=m兩條直線在[0, 1]上及其相似：

1. F 1 ( 0 ) = F 2 ( 0 ) = 0 F_1(0) = F_2(0) = 0 F1?(0)=F2?(0)=0且 F 1 ( 1 ) = F 2 ( 1 ) = 1 F_1(1) = F_2(1) = 1 F1?(1)=F2?(1)=1，意味著兩條直線在[0, 1]兩個端點重合
2. 兩條直線的最大差 M a x ( F 1 ( m ) ? F 2 ( m ) ) = 0.086071 , m ∈ [ 0 , 1 ) Max(F_1(m) - F_2(m)) = 0.086071 , m\isin [0, 1) Max(F1?(m)?F2?(m))=0.086071,m∈[0,1)
3. 兩條直線的平均差 ∫ 0 1 F 1 ( m ) ? F 2 ( m ) d x / ( 1 ? 0 ) = 0.057304959 , m ∈ [ 0 , 1 ) \int_0^1 F_1(m) - F_2(m) dx / (1 - 0)=0.057304959, m\isin [0, 1) ∫01?F1?(m)?F2?(m)dx/(1?0)=0.057304959,m∈[0,1)

注意， 上述第2條、第3條可以看我另一篇博客：數學公式之求 log2(1+x)-x的積分
因此，由于 F 1 ( m ) = l o g 2 ( 1 + m ) F_1(m)=log_2(1+m) F1?(m)=log2?(1+m)與 F 2 ( m ) = m F_2(m)=m F2?(m)=m兩條直線在[0, 1]上的相似性，我們可以用 m m m來做 l o g 2 ( 1 + m ) log_2(1+m) log2?(1+m)的近似，只不過需要加一個誤差變數 α \alpha α, 即 l o g 2 ( 1 + m ) = m + α , m ∈ [ 0 , 1 ) (7) log_2(1+m) = m + \alpha \tag{7}, m\isin[0, 1) log2?(1+m)=m+α,m∈[0,1)(7)
既然有了(7), 將其帶入(6)，我們又可以繼續前進了！！！
l o g ( x ) = [ E + l o g 2 ( 1 + M L ) ] ? B log(x)=[E + log_2(1+\cfrac{M}{L})] - B log(x)=[E+log2?(1+LM?)]?B = [ E + M L + α ] ? B =[E+\cfrac{M}{L} + \alpha] - B =[E+LM?+α]?B = ( E + M L ) × L L + α ? B (8) =\cfrac {(E+\cfrac{M}{L} ) \times L}{L} + \alpha- B \tag{8} =L(E+LM?)×L?+α?B(8)
我的天！我們看到了什么啊！一個似曾見過的數字 ( E + M L ) × L (E+\cfrac{M}{L} ) \times L (E+LM?)×L， 這個東西好像在哪里見過，趕緊去回顧一下上面第3部分的表格，發現它的意義是32bits 二進制數decoding為整數，

4. 本文核心

有了(8)，我們便得到了本文的核心，考慮到其重要性，我單開了一個文章標題來介紹它
l o g ( 小 數 x ) = ( 小 數 先 按 S E M 方 法 存 儲 為 32 b i t s ， 再 按 照 整 數 方 法 還 原 ) × 1 L + α ? B log(小數x) = (小數先按SEM方法存盤為32bits，再按照整數方法還原) \times \cfrac{1}{L} + \alpha - B log(小數x)=(小數先按SEM方法存儲為32bits，再按照整數方法還原)×L1?+α?B
即
l o g ( x ) = F l o a t ( I E ( x ) , I M ( x ) ) × 1 L + α ? B (9) log(x) = Float(I_E(x), I_M(x)) \times \cfrac{1}{L} + \alpha - B \tag{9} log(x)=Float(IE?(x),IM?(x))×L1?+α?B(9)
那么(9)式就是我們本文的核心，如果能看懂(9)，就說明你明白了本文在想干什么了，注意： L , B L, B L,B都是常量，而 α \alpha α是一個變數； I E , I M I_E, I_M IE?,IM?在第2部分表格有介紹， F l o a t ( E , M ) Float(E, M) Float(E,M)在第3部分表格有介紹，

5. 開平方根倒數

關鍵性的技術突破了，那接下來就很容易了，對 x x x進行inverse square root l o g 2 ( x ? 1 / 2 ) = ? 1 2 × l o g 2 [ ( 1 + M L ) × 2 E ? B ] log_2(x^{-1/2}) = - \cfrac{1}{2} \times log_2[(1 + \cfrac{M}{L}) \times 2^{E - B}] log2?(x?1/2)=?21?×log2?[(1+LM?)×2E?B] = ? 1 2 ( α ? B ) ? 1 2 × F l o a t ( I E ( x ) , I M ( x ) ) L (10) = - \cfrac{1}{2}(\alpha - B) - \cfrac{1}{2} \times \cfrac{Float(I_E(x), I_M(x))}{L} \tag{10} =?21?(α?B)?21?×LFloat(IE?(x),IM?(x))?(10)

什么？ 對小數開平方根的倒數 取對數后，只需要把其化為IEEE-754標準encoding存入存盤單元，再按照整數標準decoding把此32bits 二進制數還原為整數，那么對這個整數只需要進行加減乘除法的運算就相當于對原數進行乘方開方的運算，

換句話說，我們徹底完成了需求分析中的，乘方到加減乘除的降維！！！
現在重新來看我們的問題(1)式
y = x ? 1 / 2 y=x^{-1/2} y=x?1/2 兩邊取對數 l o g y = ? 1 2 l o g x logy=-\cfrac{1}{2}logx logy=?21?logx 把(8)帶入上式
( E y + M y L ) × L L + α y ? B = ? 1 2 × ( ( E x + M x L ) × L L + α x ? B ) \cfrac {(E_y+\cfrac{M_y}{L} ) \times L}{L} + \alpha_y- B = -\cfrac{1}{2} \times (\cfrac {(E_x+\cfrac{M_x}{L} ) \times L}{L} + \alpha_x- B) L(Ey?+LMy??)×L?+αy??B=?21?×(L(Ex?+LMx??)×L?+αx??B) 整理移項得 M y + E y L = [ 3 2 B L ? ( α a + α b 2 ) L ] ? 1 2 ( M x + E x L ) (11) M_y+E_yL=[\cfrac{3}{2}BL - (\alpha_a + \cfrac{\alpha_b}{2})L] - \cfrac{1}{2}(M_x + E_xL) \tag{11} My?+Ey?L=[23?BL?(αa?+2αb??)L]?21?(Mx?+Ex?L)(11)
如果把 M + E L M+EL M+EL當做一個整體，它是一個function I I I；且把(11)式[]方括號中心的全部內容看作R的話，那么(11)式可以變成什么？
I ( y ) = R ? 1 2 I ( x ) I(y) = R - \cfrac{1}{2}I(x) I(y)=R?21?I(x)
回頭一看，此式正是(3)式！！！對于此 I I I，我們可再熟悉不過了，它就是第2部分表中的encoder I I I，再幫大家回憶一下： E x = I E ( x ) = e + B E_x = I_E(x) = e + B Ex?=IE?(x)=e+B M x = I M ( x ) = m L M_x = I_M(x) = \cfrac{m}{L} Mx?=IM?(x)=Lm?這樣，有了x，我們就可以通過 I I I求出為三十二位bits二進制數的 I ( x ) I(x) I(x)，在對 I ( x ) I(x) I(x)進行乘以 ? 1 2 -\cfrac{1}{2} ?21?(或者是- I ( x ) I(x) I(x)>>1)的基礎操作后，我們給其加一個常數R就能得到 I ( y ) I(y) I(y)，之后，我們需要利用同樣是第2部分表中的decoder
y = a r g I ( E y , M y ) = ( 1 + M y L ) × 2 E y ? B y = argI(E_y, M_y) = (1 + \cfrac{M_y}{L}) \times 2^{E_y - B} y=argI(Ey?,My?)=(1+LMy??)×2Ey??B 就可以將三十二位bits的二進制數還原為一個浮點數，而這個浮點數，恰恰就是我們要求的 y y y，也就是 x ? 1 / 2 x^{-1/2} x?1/2，

由此，我們算是講完了Fast inverse square root演算法中最核心的思想，這是一個從開平方到加減乘除質的飛躍！！！

剩下的一共還有兩項作業，

1. 是求
   R = [ 3 2 B L ? ( α a + α b 2 ) L ] (12) R = [\cfrac{3}{2}BL - (\alpha_a + \cfrac{\alpha_b}{2})L] \tag{12} R=[23?BL?(αa?+2αb??)L](12)
   注意到上式 B , L B, L B,L都是確定的常數，而 α a , α b \alpha_a, \alpha_b αa?,αb?是變數，因此我們要做的是試圖用一個確定的常量 來近似所有的 α \alpha α(也包括 α a , α b \alpha_a, \alpha_b αa?,αb?)，使得(12)大致成立！！！
2. 由于1.的近似R，因此求出的 y y y只是 x ? 1 / 2 x^{-1/2} x?1/2的近似解，我們需要用牛頓迭代加深這個初始解 y y y的精確度，（根據源代碼，1次牛頓迭代就夠了，換句話說，用了FISR演算法以后只需要一次牛頓迭代，這個數值解就已經近似于決議解了，當然你可以做n次都可以，可以但沒必要好吧）

6. 作業掃尾

針對遺留下來的2個小問題，第1個問題在我的另一篇文章兩條曲線相似度的探究(MSE推廣) 有所提及，(博主為寫這一篇文章容易嘛，不僅花了一個禮拜絞盡腦汁理解演算法，還寫了兩篇輔助作業的博客，那能看到這里的人都是人才了，樓主深表欣慰，也希望你們點個贊評個論呀~)
由于 α = E r r o r ( x ) = l o g 2 ( 1 + x ) ? x , x ∈ [ 0 , 1 ) \alpha =Error(x) =log_2(1+x) - x , x\isin[0, 1) α=Error(x)=log2?(1+x)?x,x∈[0,1)
因此這里有兩種計算 α \alpha α的方法，
1.是取Error最大值得一半， α = M a x ( E r r o r ( x ) ) 2 = 0.0430357 \alpha = \cfrac{Max(Error(x))}{2}=0.0430357 α=2Max(Error(x))?=0.0430357
2.是取Error函式在[0,1]上面積分的平均， α = ∫ 0 1 E r r o r ( x ) d x 1 ? 0 = 0.0573049 \alpha =\cfrac{\int_0^1Error(x)dx}{1-0}=0.0573049 α=1?0∫01?Error(x)dx?=0.0573049

然而，這些都只是猜測與理論，經過實驗，最終源代碼的作者最終采用了 α = 0.0450465 (13) \alpha=0.0450465 \tag{13} α=0.0450465(13)這個數值（我也不知道怎么來的，應該是按照性能和準確率做實驗做出來的），
那么，將(13)結果帶入(12):
R = [ 3 2 B L ? ( α + α 2 ) L ] R = [\cfrac{3}{2}BL - (\alpha + \cfrac{\alpha}{2})L] R=[23?BL?(α+2α?)L] = 1.5 × L × ( B ? 1.5 α ) =1.5 \times L \times(B - 1.5 \alpha) =1.5×L×(B?1.5α) = 1.5 × 2 23 × ( 127 ? 1.5 × 0.0450465 ) =1.5 \times 2^{23} \times(127 - 1.5 \times 0.0450465) =1.5×223×(127?1.5×0.0450465) ≈ 1 , 597 , 463 , 007 \approx 1,597,463,007 ≈1,597,463,007 = 0 x 5 f 3759 d f (14) =0x 5f3759df \tag{14} =0x5f3759df(14)
至此，我們就解決了遺留作業1. R = 0 x 5 f 3759 d f R=0x 5f3759df R=0x5f3759df，這也就是代碼

	i  = 0x5f3759df - ( i >> 1 );               // what the fuck? 

中的那個magic number，
懂了么？明白這個magic number怎么來的了么？！

至于第2個遺留作業就不是本文的討論內容了，那是牛頓老人家的研究成果，詳情自行google或百度即可，

三、總結

額，太累了，還沒想好，下次再說吧，就五一快樂！(8點終于可以回家了)