平方根倒數速算法

<image>

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

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

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

## 演算法的切入點

<image>

浮點數的平方根倒數常用於計算正規化向量Blinn。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競技場》使用了平方根倒數速算法，而後來採用現場可編程邏輯門陣列的頂點著色器也應用了此演算法Middendorf。

## 代碼概覽

下列代碼是《雷神之錘III競技場》原始碼中平方根倒數速算法之實例。示例去除了C預處理器的指令，但附上了原有的注釋（括號內為注釋的翻譯）：

```c
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;
}
```

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

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

### 未定義行為的問題

從 C99 和 C++98 標準開始，由於嚴格別名規則，形如 `i = * (long*) &y;` 的寫法是未定義行為。想要正確做到類型雙關，一種同時符合標準 C 和 C++ 的方式是使用 `memcpy`，但這不是唯一一種方法。

下面是更符合標準 C 和標準 C++，可移植性更高的代碼：

```c
#include <string.h>
#include <stdint.h>

float Q_rsqrt(float number)
{
	float x = number;

	// 第一步：将 x 当作整数，进行移位和减法操作，变成近似的结果
	uint32_t i;
	memcpy(&i, &x, sizeof i); // 编译器往往是会优化的，不必担心时间和空间开销
	i = 0x5f3759dfu - (i >> 1);
	memcpy(&x, &i, sizeof x);

	// 第二步：对 x 进行牛顿迭代，使之更接近结果
	const float xhalf = number * 0.5f;
	x = x * (1.5f - (xhalf * x * x));
	//x = x * (1.5f - (xhalf * x * x)); 第二次迭代可以去掉了

	return x;
}
```

但下文仍然會用第一個版本的代碼進行講解。

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

<image>

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

0｜0｜0｜0｜0｜0｜1｜0｜=｜2

0｜0｜0｜0｜0｜0｜0｜1｜=｜1

0｜0｜0｜0｜0｜0｜0｜0｜=｜0

1｜1｜1｜1｜1｜1｜1｜1｜=｜−1

1｜1｜1｜1｜1｜1｜1｜0｜=｜−2

1｜0｜0｜0｜0｜0｜0｜1｜=｜−127

1｜0｜0｜0｜0｜0｜0｜0｜=｜−128

如上所述，一個有符號正整數在二進位補碼系統中的表示中首位為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所除Hennessey。

### 「魔術數字」

S(ign,符號)｜E(xponent,指數)｜M(antissa,尾數)

1 位｜b位｜(n-1-b)位

n位McEniry

對猜測平方根倒數速算法的最初構想來說，計算首次近似值所使用的常數0x5f3759df也是重要的線索。為確定程式設計師最初選此常數以近似求取平方根倒數的方法，Charles McEniry首先檢驗了在代碼中選擇任意常數R所求取出的首次近似值的精度。回想上一節關於整數和浮點數表示的比較：對於同樣的32位二進位數碼，若為浮點數表示時實際數值為$\textstyle x=(1+m_x)2^{e_x}$，而若為整數表示時實際數值則為$\textstyle I_x=E_xL+M_x$此處「浮點數」所指為標準化浮點數，也即有效數字部分必須滿足$\textstyle 0 \le m_x <1)$，可參見，其中$\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。在此McEniry再度引入新數$\sigma$描述$\textstyle \log_2{(1+x)}$與近似值R間的誤差Lomont 2003確定R的方式則有所不同，他先將R分解為$\textstyle R_1$與$\textstyle R_2$，其中$\textstyle R_1$與$\textstyle R_2$分別代表R的有效數字域和指數域Lomont。由於$\textstyle 0 \le x < 1$，有$\textstyle \log_2{(1+x)}\approx {x}$，則在此可定義$\sigma$與x的關係為$\textstyle \log_2{(1+x)}\cong x+\sigma$，這一定義就能提供二進位對數的首次精度值（此處$0\le\sigma\le\tfrac{1}{3}$；當R為0x5f3759df時，有$\textstyle \sigma=0.0450461875791687011756$McEniry）。由此將$\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_yL+M_y=\frac{3}{2}(B-\sigma)L-\frac{1}{2}(E_xL+M_x)$$

如上所述，對於以浮點規格存儲的正浮點數x，若將其作為長整型表示則示值為$\textstyle I_x=E_xL+M_x$，由此即可根據x的整數表示導出y（在此$\textstyle y=\frac{1}{\sqrt{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$.

最後導出的等式$\textstyle I_y=R-\frac{1}{2}I_x$即與上節代碼中`i = 0x5f3759df - (i>>1);`一行相契合，由此可見，在平方根倒數速算法中，對浮點數進行一次移位操作與整數減法，就可以可靠地輸出一個浮點數的對應近似值McEniry。到此為止，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$Hennessey，這時用R（在此即為「魔術數字」0x5f3759df）減之即可使指數的最低有效數位轉入有效數字域，之後重新轉換為浮點數時，就能得到相當精確的平方根倒數近似值。在這裡對常數R的選取亦有所講究，若能選取一個好的R值，便可減少對指數進行除法與對有效數字域進行移位時可能產生的錯誤。基於這一標準，0xbe即是最合適的R值，而0xbe右移一位即可得到0x5f，這恰是魔術數字R的第一個字節Lomont。

### 精確度

<image>

如上所述，平方根倒數速算法所得的近似值驚人的精確，右圖亦展示了以上述代碼計算（以平方根倒數速算法計算後再進行一次牛頓法迭代）所得近似值的誤差：當輸入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}$）作為二次迭代的輸入）正是為了進一步提高結果的精度McEniry，但由於雷神之錘III引擎的圖形計算中並不需要太高的精度，所以代碼中只進行了一次迭代，二次迭代的代碼則被注釋Eberly。

## 歷史與考究

<image>

《雷神之錘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；因此Lomont將目標改為尋找在進行1－2次牛頓迭代後能得到最大精度的R值，在暴力搜尋後得出最優R值為0x5f375a86，以此值代入演算法並進行牛頓迭代，所得的結果都比代入原始值（0x5f3759df）更精確Lomont，於是他的論文最後提到「如果可能我想詢問原作者，此速算法是以數學推導還是以反覆試錯的方式求得？」Lomont。在論文中，Lomont亦指出，64位的IEEE754浮點數（即雙精度類型）所對應的魔術數字是0x5fe6ec85e7de30da，但後來的研究表明，代入0x5fe6eb50c7aa19f9的結果精確度更高（McEniry得出的結果則是0x5FE6EB50C7B537AA，精度介於兩者之間，Matthew Robertson給出的精度更高的值是0x5FE6EB50C7B537A9）。在Charles McEniry的論文中，他使用了一種類似Lomont但更複雜的方法來優化R值：他最開始使用窮舉搜索，所得結果與Lomont相同McEniry；而後他嘗試用調日法尋找最優值，所得結果恰是代碼中所使用的魔術數字0x5f3759df，因此，McEniry認為，這一常數最初或許便是以「在可容忍誤差範圍內使用調日法」的方式求得McEniry。

## 外部連結

- Origin of Quake3's Fast InvSqrt , Beyond3D.com
- Quake III Arena source code , id Software