平方根倒數速算法

**反平方根快速演算法**（英語：Fast Inverse Square Root，亦常以「Fast InvSqrt ()」或其使用的十六進位常數 **0x5f3759df** 代稱）是用於快速計算 $\textstyle x^{-1/2}$ （即 $\textstyle x$ 的平方根的倒數，在此 $\textstyle x$ 需取符合 IEEE 754 標準格式的 32 位元浮點數）的一種演算法。這一演算法的優勢在於減少了求平方根倒數時浮點運算操作帶來的巨大的運算耗費，而在電腦圖學領域，若要求取照明和投影的波動角度與反射效果，就常需計算平方根倒數。

此演算法首先接收一個 32 位元帶符浮點數，然後將之作為一個 32 位元整數看待，以將其向右進行一次邏輯移位的方式將之取半，並用在浮點數規格代表√2127 近似值的十六進位「魔術數字」0x5f3759df 減之，如此即可得對輸入的浮點數的平方根倒數的首次近似值；而後重新將其作為浮點數，以牛頓法反覆迭代，以求出更精確的近似值，直至求出符合精確度要求的近似值。在計算浮點數的平方根倒數的同一精度的近似值時，此演算法比直接使用浮點數除法要快四倍。

此演算法最早被認為是由約翰・卡馬克於 90 年代前期在 SGI Indigo 的開發中使用，後來則於 1999 年在《雷神之鎚 III 競技場》的原始碼中應用，但直到 2002－2003 年間才在 Usenet 一類的公共論壇上出現。後來的調查顯示，該演算法在這之前就於電腦圖學的硬體與軟體領域有所應用，如 SGI 和 3dfx 就曾在產品中應用此演算法，但至今為止仍未能確切知曉演算法中所使用的特殊常數的起源。

## 演算法的切入點

浮點數的平方根倒數常用於計算正規化向量。3D 圖形程式需要使用正規化向量來實現光照和投影效果，因此每秒都需做上百萬次平方根倒數運算，而在處理坐標轉換與光源的專用硬體裝置出現前，這些計算都由軟體完成，計算速度亦相當之慢；在 1990 年代這段代碼開發出來之時，多數浮點數操作的速度更是遠遠滯後於整數操作，因而針對正規化向量演算法的最佳化就顯得尤為重要。下面陳述計算正規化向量的原理：

要將一個向量標準化，就必須計算其歐幾里得範數，以求得向量長度，為此便需對向量的各分量的平方和求平方根；而當求取到其長度，並以之除該向量的每個分量後，所得的新向量就是與原向量同向的單位向量，若以公式表示：

　
: $\|{\boldsymbol {v}}\|={\sqrt {v_{1}^{2}+v_{2}^{2}+v_{3}^{2}}}$ 可求得向量 v 的歐幾里得範數，此演算法正類如對歐幾里得空間的兩點求取其歐幾里得距離，

　
: 而 ${\boldsymbol {\hat {v}}}={\boldsymbol {v}}/\|{\boldsymbol {v}}\|$ 求得的就是標準化的向量，若以 ${\boldsymbol {x}}$ 代表 $v_{1}^{2}+v_{2}^{2}+v_{3}^{2}$ ，則有 ${\boldsymbol {\hat {v}}}={\boldsymbol {v}}/{\sqrt {x}}$ ，

可見標準化向量時，對向量分量計算平方根倒數實為必需，所以，對平方根倒數計算演算法的最佳化對計算正規化向量也大有裨益。

為了加速圖像處理單元計算，《雷神之鎚 III 競技場》使用了反平方根快速演算法，而後來採用現場可程式化邏輯閘陣列的頂點著色器也應用了此演算法。

## 代碼概覽

下列代碼是《雷神之鎚 III 競技場》原始碼中反平方根快速演算法之實例。範例去除了 C 預處理器的指令，但附上了原有的注釋：

為計算平方根倒數的值，軟體首先要先確定一個近似值，而後則使用某些數值方法不斷計算修改近似值，直至達到可接受的精度。在 1990 年代初（也即該演算法發明的大概時間），軟體開發時通用的平方根演算法多是從尋找表中取得近似值，而這段代碼取近似值耗時比之更短，達到精確度要求的速度也比通常使用的浮點除法計演算法快四倍，雖然此演算法會損失一些精度，但效能上的巨大優勢已足以補償損失的精度。由代碼中對原資料的變數類型聲明為 float 可看出，這一演算法是針對 IEEE 754 標準格式的 32 位元浮點數設計的，不過據 Chris Lomont 和後來的 Charles McEniry 的研究看，這一演算法亦可套用於其他類型的浮點數上。

反平方根快速演算法在速度上的優勢源自將浮點數轉化為長整型以作整數看待，並用特定常數 0x5f3759df 與之相減。然而對於代碼閱讀者來說，他們卻難以立即領悟出使用這一常數的目的，因此和其它在代碼中出現的難以理解的常數一樣，這一常數亦被稱為「魔術數字」。如此將浮點數當作整數先位移後減法，所得的浮點數結果即是對輸入數字的平方根倒數的粗略估計值，而後再進行一次牛頓迭代法，以使之更精確後，代碼即執行完畢。由於演算法所生成的用於輸入牛頓法的首次近似值已經相當精確，此演算法所得近似值的精度已可接受，而若使用與《雷神之鎚 III 競技場》同為 1999 年發布的 Pentium III 中的 SSE 指令 rsqrtss 計算，則計算平方根倒數的收斂速度更慢，精度也更低。

### 將浮點數轉化為整數

要理解這段代碼，首先需了解浮點數的儲存格式。一個浮點數以 32 個位元表示一個有理數，而這 32 位元由其意義分為三段：首先首位為符號位，如若是 0 則為正數，反之為負數；接下來的 8 位元表示經過偏移處理（這是為了使之能表示 - 127－128）後的指數；最後 23 位表示的則是有效數字中除最高位以外的其餘數字。將上述結構表示成公式即為 $\textstyle x=(-1)^{\mathrm {Si} }\cdot (1+m)\cdot 2^{(E-B)}$ ，其中 $\textstyle m$ 表示有效數字的尾數（此處 $\textstyle 0\leq m<1$ ，偏移量 $\textstyle B=127$ ，而指數的值 $\textstyle E-B$ 決定了有效數字（在 Lomont 和 McEniry 的論文中稱為「尾數」（_mantissa_））代表的是小數還是整數。以上圖為例，將描述帶入有 $\textstyle m=1\times 2^{-2}=0.250$ ），且 $\textstyle E-B=124-127=-3$ ，則可得其表示的浮點數為 $\textstyle x=(1+0.250)\cdot 2^{-3}=0.15625$ 。

如上所述，一個有符號正整數在二補數系統中的表示中首位為 0，而後面的各位則用於表示其數值。將浮點數取別名儲存為整數時，該整數的數值即為 $\textstyle I=E\times 2^{23}+M$ ，其中 E 表示指數，M 表示有效數字；若以上圖為例，圖中樣例若作為浮點數看待有 $\textstyle E=124$ ， $M=1\cdot 2^{21}$ ，則易知其轉化而得的整數型號數值為 $I=124\times 2^{23}+2^{21}$ 。由於平方根倒數函式僅能處理正數，因此浮點數的符號位（即如上的 Si）必為 0，而這就保證了轉換所得的有符號整數也必為正數。以上轉換就為後面的計算帶來了可行性，之後的第一步操作（邏輯右移一位）即是使該數的長整形式被 2 所除。

### 「魔術數字」

對猜測反平方根快速演算法的最初構想來說，計算首次近似值所使用的常數 **0x5f3759df** 也是重要的線索。為確定程式設計師最初選此常數以近似求取平方根倒數的方法，Charles McEniry 首先檢定了在代碼中選擇任意常數 R 所求取出的首次近似值的精度。回想上一節關於整數和浮點數表示的比較：對於同樣的 32 位元二進位數位，若為浮點數表示時實際數值為 $\textstyle x=(1+m_{x})2^{e_{x}}$ ，而若為整數表示時實際數值則為 $\textstyle I_{x}=E_{x}L+M_{x}$ ，其中 $\textstyle L=2^{n-1-b}$ 。以下等式引入了一些由指數和有效數字匯出的新元素：

　
: $m_{x}={\frac {M_{x}}{L}}$

　
: $e_{x}=E_{x}-B$ ，其中 $B=2^{b-1}-1$

再繼續看 McEniry 2007 里的進一步說明：

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

對等式的兩邊取二進位對數（ $\textstyle \log _{2}$ ，即函式 $\textstyle f(n)=2^{n}$ 的反函式），有

　
: $\log _{2}{(y)}=-{\frac {1}{2}}\log _{2}{(x)}$

　
: $\log _{2}(1+m_{y})+e_{y}=-{\frac {1}{2}}\log _{2}{(1+m_{x})}-{\frac {1}{2}}e_{x}$

以如上方法，就能將浮點數 x 和 y 的相關指數消去，從而將乘方運算化為加法運算。而由於 $\textstyle \log _{2}{(x)}$ 與 $\textstyle \log _{2}{(x^{-1/2})}$ 線性相依，因此在 $\textstyle x$ 與 $\textstyle y_{0}$ （即輸入值與首次近似值）間就可以線性組合的方式建立方程式。在此 McEniry 再度引入新數 $\sigma $ 描述 $\textstyle \log _{2}{(1+x)}$ 與近似值 R 間的誤差：由於 $\textstyle 0\leq x<1$ ，有 $\textstyle \log _{2}{(1+x)}\approx {x}$ ，則在此可定義 $\sigma $ 與 x 的關係為 $\textstyle \log _{2}{(1+x)}\cong x+\sigma $ ，這一定義就能提供二進位對數的首次精度值（此處 $0\leq \sigma \leq {\tfrac {1}{3}}$ ；當 R 為 0x5f3759df 時，有 $\textstyle \sigma =0.0450461875791687011756$ ）。由此將 $\textstyle \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}$

參照首段等式代入 $M_{x}$ ， $E_{x}$ ， $B$ 與 $L$ ，有：

　
: $M_{y}+(E_{y}-B)L=-{\frac {3}{2}}\sigma {L}-{\frac {1}{2}}M_{x}-{\frac {1}{2}}(E_{x}-B)L$

移項整理得：

　
: $E_{y}L+M_{y}={\frac {3}{2}}(B-\sigma )L-{\frac {1}{2}}(E_{x}L+M_{x})$

如上所述，對於以浮點規格儲存的正浮點數 x，若將其作為長整型表示則示值為 $\textstyle I_{x}=E_{x}L+M_{x}$ ，由此即可根據 x 的整數表示匯出 y（在此 $\textstyle y={\frac {1}{\sqrt {x}}}$ ，亦即 x 的平方根倒數的首次近似值）的整數表示值，也即：

　
: $I_{y}=E_{y}L+M_{y}=R-{\frac {1}{2}}(E_{x}L+M_{x})=R-{\frac {1}{2}}I_{x}$ ，其中 $R={\frac {3}{2}}(B-\sigma )L$ .

最後匯出的等式 $\textstyle I_{y}=R-{\frac {1}{2}}I_{x}$ 即與上節代碼中 i = 0x5f3759df - (i>>1);一行相契合，由此可見，在反平方根快速演算法中，對浮點數進行一次移位元運算與整數減法，就可以可靠地輸出一個浮點數的對應近似值。到此為止，McEniry 只證明了，在常數 R 的輔助下，可近似求取浮點數的平方根倒數，但仍未能確定代碼中的 R 值的選取方法。

關於作一次移位與減法操作以使浮點數的指數被 - 2 除的原理，Chris Lomont 的論文中亦有個相對簡單的解釋：以 $\textstyle 10000=10^{4}$ 為例，將其指數除 - 2 可得 $\textstyle 10000^{-1/2}=10^{-2}=1/100$ ；而由於浮點表示的指數有進行過偏移處理，所以指數的真實值 e 應為 $\textstyle e=E-127$ ，因此可知除法操作的實際結果為 $\textstyle -e/2+127$ ，這時用 R（在此即為「魔術數字」0x5f3759df）減之即可使指數的最低有效數位轉入有效數字域，之後重新轉換為浮點數時，就能得到相當精確的平方根倒數近似值。在這裡對常數 R 的選取亦有所講究，若能選取一個好的 R 值，便可減少對指數進行除法與對有效數字域進行移位時可能產生的錯誤。基於這一標準，0xbe 即是最合適的 R 值，而 0xbe 右移一位即可得到 0x5f，這恰是魔術數字 R 的第一個位元組。

### 精確度

如上所述，反平方根快速演算法所得的近似值驚人的精確，右圖亦展示了以上述代碼計算（以反平方根快速演算法計算後再進行一次牛頓法迭代）所得近似值的誤差：當輸入 0.01 時，以 C 語言標準庫函式計算可得 10.0，而 InvSqrt () 得值為 9.9825822，其間誤差為 0.017479，相對誤差則為 0.175%，且當輸入更大的數值時，絕對誤差不斷下降，相對誤差也一直控制在一定的範圍之內。

### 牛頓法提高精度

在進行了如上的整數操作之後，範例程式再度將被轉為長整型的浮點數迴轉為浮點數（對應 x = \*(float\*)&i;），並對其進行一次浮點運算操作（對應 x = x\*(1.5f - xhalf\*x\*x);），這裡的浮點運算操作就是對其進行一次牛頓法迭代，若以此例說明：

　
: $y={\frac {1}{\sqrt {x}}}$ 所求的是 y 的平方根倒數，以之構造以 y 為自變數的函式，有 $f(y)={\frac {1}{y^{2}}}-x=0$ ，

　
: 將其代入牛頓法的通用公式 $y_{n+1}=y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}$ （其中 $\,y_{n}$ 為首次近似值），

　
: 有 $y_{n+1}={\frac {y_{n}(3-xy_{n}^{2})}{2}}$ ，其中 $f(y)={\frac {1}{y^{2}}}-x$ ， $f'(y)={\frac {-2}{y^{3}}}$ 。

　
: 整理有 $\,y_{n+1}={\frac {y_{n}(3-xy_{n}^{2})}{2}}=y_{n}(1.5-{\frac {xy_{n}^{2}}{2}})$ ，對應的代碼即為 y = y\*(threehalfs - xhalf\*y\*y);。

在以上一節的整數操作產生首次近似值後，程式會將首次近似值作為參數送入函式最後兩句進行精化處理，代碼中的兩次迭代（以一次迭代的輸出（對應公式中的 $y_{n+1}$ ）作為二次迭代的輸入）正是為了進一步提高結果的精度，但由於雷神之鎚 III 引擎的圖形計算中並不需要太高的精度，所以代碼中只進行了一次迭代，二次迭代的代碼則被注釋。

## 歷史與考究

《雷神之鎚 III》的代碼直到 QuakeCon 2005 才正式放出，但早在 2002 年（或 2003 年）時，反平方根快速演算法的代碼就已經出現在 Usenet 與其他論壇上了。最初人們猜測是卡馬克寫下了這段代碼，但他在回覆詢問他的郵件時否定了這個觀點，並猜測可能是先前曾幫 id Software 最佳化雷神之鎚的資深組譯程式設計師 Terje Mathisen 寫下了這段代碼；而在 Mathisen 的郵件裡，他表示，在 1990 年代初，他只曾作過類似的實作，確切來說這段代碼亦非他所作。現在所知的最早實作是由 Gary Tarolli 在 SGI Indigo 中實作的，但他亦坦承他僅對常數 R 的取值做了一定的改進，實際上他也不是作者。在向以發明 MATLAB 而聞名的 Cleve Moler 查證後，Rys Sommefeldt 則認為原始的演算法是 Ardent Computer 公司的 Greg Walsh 所發明，但他也沒有任何決定性的證據能證明這一點。

不僅該演算法的原作者不明，人們也仍無法確定當初選擇這個「魔術數字」的方法。Chris Lomont 曾做了個研究：他推算出了一個函式以討論此速演算法的誤差，並找出了使誤差最小的最佳 R 值 0x5f37642f（與代碼中使用的 0x5f3759df 相當接近），但以之代入演算法計算並進行一次牛頓迭代後，所得近似值之精度仍略低於代入 0x5f3759df 的結果；因此 Lomont 將目標改為尋找在進行 1－2 次牛頓迭代後能得到最大精度的 R 值，在暴力搜尋後得出最佳 R 值為 0x5f375a86，以此值代入演算法並進行牛頓迭代，所得的結果都比代入原始值（0x5f3759df）更精確，於是他的論文最後提到「如果可能我想詢問原作者，此速演算法是以數學推導還是以反覆試錯的方式求得？」。在論文中，Lomont 亦指出，64 位元的 IEEE754 浮點數（即雙精度類型）所對應的魔術數字是 0x5fe6ec85e7de30da，但後來的研究表明，代入 0x5fe6eb50c7aa19f9 的結果精確度更高（McEniry 得出的結果則是 0x5FE6EB50C7B537AA，精度介於兩者之間，英文 wiki 給出的精度更高的值是 0x5FE6EB50C7B537A9）。在 Charles McEniry 的論文中，他使用了一種類似 Lomont 但更複雜的方法來最佳化 R 值：他最開始使用窮舉搜尋，所得結果與 Lomont 相同；而後他嘗試用帶權二分法尋找最佳值，所得結果恰是代碼中所使用的魔術數字 0x5f3759df，因此，McEniry 認為，這一常數最初或許便是以「在可容忍誤差範圍內使用二分法」的方式求得。

## 注釋

## 參考

### 註腳

### 參考文獻

## 延伸閱讀

* Kushner, David. The wizardry of Id. IEEE Spectrum. Aug 2002, **39** (8): 42–47. doi:10.1109/MSPEC.2002.1021943.
* 顧森（Matrix67）. 旧闻一则：神秘的 0x5f3759df 不可思议的 Quake III 源码. 2007-11-24 [2015-08-01]. （原始內容存檔於 2015-08-04）.

## 外部連結

* Origin of Quake3's Fast InvSqrt () （頁面存檔備份，存於網際網路檔案館）, Beyond3D.com
* Quake III Arena source code （頁面存檔備份，存於網際網路檔案館）, id Software
* Margolin, Tomer. Magical Square Root Implementation In Quake III. CodeMaestro. The Coding Experience. 2005-08-27. （原始內容存檔於 2012-04-14）.