揭祕·變態的平方根倒數算法

時間 2019-11-06

標籤揭祕變態平方根倒數算法简体版

原文原文鏈接

本文首發於個人博客 轉載請註明出處算法

神的時代已離去
神的故事卻化爲傳說
流落凡間
供凡人傳頌、膜拜app

這是什麼

在上世紀 90 年代，出現過一款難以想象的遊戲——雷神之錘（Quake series）。除了優秀的情節設定和精美的畫面，最讓人稱道的莫過於它的運行效率——要知道在那個計算機配置低下的時代，一段小動畫都是一個奇蹟，但 Quake 卻能流暢地運行於各類配置的電腦上。性能

直至 2005 年，當 Quake Engine 開源時，Quake 系列的祕密才被揭開。在代碼庫中，人們發現了許多堪稱神來之筆的算法。它們以極其變態的高效率，壓榨着計算機的性能，進而才支撐起了 90 年代 3D 遊戲的傳奇。其中的某些算法，甚至比系統原生的實現還要快！優化

咱們今天的主角——快速平方根倒數算法（Fast Inverse Square Root）就是其中一個。動畫

這是一個用於計算 $1 \over {\sqrt x}$ 的算法。這是一步重要的運算，由於在 3D 繪圖中，計算機須要大量求解一個矢量的方向矢量，平方根倒數爲其中的一步（$(v_x,v_y, v_z) \over {\sqrt{v_x^2+v_y^2+v_z^2}})$），並且爲最麻煩的一步。若是能在此做優化，渲染效率無異會極大地提升。this

然而，高效並非它的惟一特色，吸引人們的更是其中的一個神祕的常數——$(5f3759df)_{16}$。spa

且看代碼：翻譯

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

翻譯成數學語言就是：code

設輸入數爲 $number$，令 $ x_2=number/2, y=number $。blog
上式中，$y$ 是一個 32 位浮點數（能夠理解爲數學中的實數），如今咱們將其「看做」一個整數，並賦給 $i$。
最神奇的一步出現了——用一個十六進制數 $5f3759df$ 減去 $i$ 自身右移一位的結果，並將結果賦給 $i$。
將整數 $i$ 從新「看做」是一個浮點數，並賦值給 $y$。
令 $y=y(1.5-x_2 \times y^2)$。
重複一遍操做 5（事實上這一步能夠忽略）。

在此須要向無計算機基礎的朋友解釋一下：

右移（Right Shifting）指的是：將一個數表示成二進制以後總體向右移動一位，並抹去溢出的末位。如 $(4)_{10} >> 1=(100)_2 >> 1=(10)_2=(2)_{10}$。其效果等價爲整除 2，但因爲有 CPU 指令直接支持，速度比整除快若干個數量級。
不管是實數仍是整數，在計算機中只不過是一串二進制序列。所以一串二進制序列，既能夠被看做一個實數，也能夠被看做一個整數。

這個算法獲得的只是一個近似值。對於 $[0.01, max\_float]$ 內的全部浮點數，最大偏差爲 $0.175%$（見 Accuracy）。但它卻比內置算法 1.0/sqrt(x) 快 4 倍，可謂瑕不掩瑜。

那麼，它到底是怎麼實現的呢？

牛頓法

咱們從後面開始分析。

從代碼中能夠看出：算法最後有兩步相同的操做，像是在對一個數進行某種迭代。而其中的第二步被註釋掉了，彷佛是由於和性能損耗相比對結果的二次迭代意義不大，也說明 一次迭代的結果在偏差容許範圍內——這讓人想到了牛頓法。

牛頓法是什麼？

牛頓法是一種經常使用的求方程數值解的方法。其敘述以下：

若在某個區間 $I$ 中，$f(x)$ 連續可導，且有惟一零點 $x_0$，則任取 $x_1 \in I$，定義數列 $x_{n+1}=x_n-{f(x_n) \over {f'(x_n)}}$，則有 $lim_{n \to \infty}{x_n}=x_0$。

用牛頓法進行迭代，能夠完成對解任意精度的數值逼近。下面咱們嘗試寫出 $x=1 \over {\sqrt {a}}$的迭代式。

令 $f(x)={ {1 \over x^2}-a}$，則有

$$ x_{n+1}=x_n-{f(x_n) \over {f'(x_n)}} $$
$$ =x_n-{ { {1 \over x_n^2} -a} \over {-2 \over {x_n^3} } } $$
$$ ={3 \over 2}x_n - { {ax_n^3} \over 2} $$
$$ ={x_n}(1.5-{a \over 2}x_n^2)$$

若是咱們將$x_{n+1}$、$x_n$替換成$y$，將$a \over 2$替換成$x2$，能夠發現和算法的最後一步是吻合的。由此可知：算法的最後確實採用了牛頓法。

也許你注意到：能解出 $1 \over {\sqrt a}$ 的方程不止這一條，迭代式的形式有不少。事實上，做者有意選擇了這條方程——由於只有從這條方程得出的迭代式是 不用除法的。除法的性能十分糟糕，應該儘可能避免。

可是，牛頓法不是須要迭代屢次的嗎？怎麼在這裏只進行一次就足夠精確了呢？

牛頓法的收斂過程依賴於初值$x_1$的選取。若想一步到位，除非初值自己已經足夠精確了。

初值是什麼？

就是那神奇的第 3 步獲得的結果。

神奇的 0x5f3759df

這是整個算法的奧妙所在。咱們再來回顧一下：

上式中，$y$ 是一個 32 位浮點數（能夠理解爲數學中的實數），如今咱們將其「看做」一個整數，並賦給 $i$。
最神奇的一步出現了——用一個十六進制數 $5f3759df$ 減去 $i$ 自身右移一位的結果，並將結果賦給 $i$。
將整數 $i$ 從新「看做」是一個浮點數，並賦值給 $y$。

有了上一節的分析，幾乎能夠確定：這是爲了獲得 $1 \over {\sqrt {number} }$ 的一個粗略值，即應該有 $y \approx {1 \over {\sqrt {number} } }$。

爲了進一步的論證，咱們首先要了解一個知識點：

IEEE 754 （浮點數儲存標準）

IEEE 754: Standard for Binary Floating-Point Arithmetic

IEEE floating point - Wikipedia

這是計算機內實數的儲存標準。

在本算法中，咱們使用的是 32 位浮點數（即用32個二進制位表示一個實數）。儲存方式以下：

對任意一個實數 $x$，總能夠將其分解爲以下形式：
$$ x=2^E \times (1+M)\ (E \in Z,M \in [0, 1))$$
則 32 個二進制位安排以下：
$$ S,EEEEEEEE,MMMMMMMMMMMMMMMMMMMMMM $$

首位爲符號位，取 0 爲正數，取 1 爲負數

接下來 8 位爲 帶符號 指數。根據帶符號數的儲存方式，該數減去 127 才爲真實的指數 $E$

接下來 23 位爲底數。是 $M$ 左移 23 位再取整的結果

舉個例子：將 $0.15625=(0.00101)_2$ 換算成浮點數：

$$ 0.15625=+2^{-3}\times(1+0.25) $$

則有：

符號位爲 $0$
$E=-3$，換算成帶符號正數爲 $-3+127=124=(01111100)_2$
$M=0.25=(0.01)_2$，左移 23 位後爲 $(01000000000000000000000)_2$

從而浮點表示爲：$0,01111100,01000000000000000000000$

回到原題

如今，讓咱們思考一下：

如何只用加、減、乘和位運算神出鬼沒地快速逼近 $1 \over {\sqrt x}$？

若想回答這個問題，得看咱們對 $x$ 瞭解多少。

從上一小節的浮點數標準能夠看出：計算機看到的 x 和咱們接觸的 x 結構不一樣，早在程序的編譯期，$x$ 就被拆成了指數和底數兩部分，並被打包存好——這一步是不耗時間的，但卻給咱們提供了海量的信息。

等等！指數和底數！彷佛在暗示着什麼！

咱們，是否是能夠經過對數進行運算？

設想，若是咱們能利用以上信息，輕易地轉換 $x$ 和 $\log_2 x$，由 $\log_2 {1 \over {\sqrt x}}={-{1 \over 2}}{\log_2 x}$ 就能夠求得 $1 \over {\sqrt x}$ 的值了。

那麼，$\log_2 x$ 該如何表示呢？且看如下變換。

根據 IEEE 754，對於一個 $x=2^E \times (1+M)>0$，若是咱們將其 32 位浮點表示看做一個 32 位整數 $I$，則有

$$I=2^{23}\times (E+127)+2^{23}\times M$$
$$=2^{23}\times (E+127+M)$$

這個 $I$ 爲已知量。經過已知量 $I$，咱們能夠獲得已知量 $E+M$。

而另外一方面，

$$ \log_2 x=\log_2 {2^E \times (1+M)} $$
$$ = E + \log_2(1+M) $$

記

$$ F(M) = M $$
$$ G(M) = \log_2(1+M) $$

注意到：當 $M\in[0,1)$ 時，$G(M)\in[0,1)$。同時觀察圖像能夠發現 $G(M)$ 在此區間接近線性。

所以，咱們能夠經過上下移動 $F(M)$ 對 $G(M)$ 進行線性擬合。

注意一下這個例子中的偏差衡量標準：咱們在儘可能減少偏差的同時，也要保證偏差分佈儘可能均勻，從而最大偏差要儘可能小。

綜合以上考慮，問題最終化簡爲：

找到一個$\sigma$，對$\Delta(M)=\log_2(1+M)-(M+\sigma)$，使$|\Delta(M)|$在 $[0,1)$ 上的最大值儘可能小。

易知：$|\Delta(0)|=|\Delta(1)|=\sigma$，另外一個可能的極值點$M_0$知足$\Delta'(M_0)=0$，解得$M_0={1\over \ln 2}$

當 $max\{|\Delta(M)|\}$ 最小時，必有：$\Delta(M_0)=\sigma$。解得 $\sigma=0.0430356660279671$。

從而咱們有：$ \log_2 (1+M) \approx M+\sigma $。進而：

$$ I=2^{23}\times(E+M+127)$$
$$ = 2^{23}\times(127+E+M+\sigma-\sigma)$$
$$ \approx 2^{23}\times(127-\sigma+(E+\log_2(1+M)))$$
$$ = 2^{23}\times(127-\sigma+\log_2 x)$$

則有：

$$ {I \over 2^{23}} + \sigma - 127 \approx \log_2 x $$

記 $r={1 \over {\sqrt x} }$ 對應的 $I$ 值爲 $I_r$。則：

$$ I_r \approx 2^{23}\times(127-\sigma+\log_2 m)$$
$$ = 2^{23}\times(127-\sigma-{1 \over 2}{\log_2 x})$$
$$ \approx 2^{23}\times(127-\sigma-{1 \over 2}({I \over 2^{23}} + \sigma - 127))$$
$$ = (127-\sigma)\times 2^{22} \times 3-{I \over 2}$$
$$ = (5f37bcb6)_{16}-{I \over 2}$$
$$ = (5f37bcb6)_{16}-(I >> 1)$$

這即是步驟 3 的那個式子，神祕的常數終於出現了。

這裏算得的常數和原算法有一些不一樣，主要是 $\sigma$ 取值差別形成的。若是咱們取 $\sigma= 0.0450466$ 就會獲得原算法的常數。做者選取了後者，應該是綜合考慮了牛頓法迭代產生的偏差。我沒有深刻研究——事實上，在容許範圍內微調 $\sigma$ 不會對精度和速度產生顯著影響。

Lomont 的另外一種詮釋

儘管 Quake Engine 的源碼直到 2005 年才被公開，這個算法卻在 2000 年左右即爲人們所知。2003 年，數學家 Chris Lomont 曾寫過一篇論文討論過這一常數的由來。他採用了另外一種不一樣的方法。

Lomont 也發現 i = 0x5f3759df - ( i >> 1 ) 實爲 $1 \over {\sqrt x}$ 的粗略逼近，但他沒有深究這個式子背後的數學含義，而是將此常數泛化爲一個變量 $R$，經過最小化其與 $1 \over {\sqrt x}$ 的偏差，反解出常數 $R$，其中運用的一些技巧值得品味。

有趣的是，他在論文中提供了另外一個他認爲最佳的常數 5f375a86，並經過實測證明了他的想法。但此常數只是以毫釐之差險勝原算法（一次迭代 $0.175228$ vs $0.175124$，二次迭代 $4.66 \times 10^{-4}$ vs $4.65437 \times 10^{-4}$）。Lomont 在論文中表達了他的驚訝，這個實驗同時也使那位不知名的做者（傳聞說是 Quake 的做者 Carmack，然而他本身不認可）更加撲朔迷離。