數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

JLee發表於2016-09-03

在伯樂小組的這個討論帖《你見過或寫過的最複雜的 C 語言程式是?》中,就提到了「平方根倒數速演算法」中的神奇數字 0x5f3759df。Christian 在本文中探討了該演算法諸多有趣的地方,解釋了它的原理,對指數在 -1 到 1 的情況進行了擴充,還有一些數學相關的新思路。

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

0x5f3759df

這篇文章講述一個神奇的常數 0x5f3759df 和一個非常簡潔的演算法,平方根倒數速演算法,也是這個常數的來源。

來看看平方根倒數速演算法:

這段程式碼可以快速計算

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

的高精度近似值。

現如今這個函式很出名,而它的流行是從出現在 2005 年的 Quake III Arena 的原始碼中開始的,起初歸功於 John Carmack,但後來發現早在 Quake 之前,80 年代中期,Ardent Computer 公司就已使用,比 SGI 和 3dfx 還早,已能追溯到的最初作者是 Greg Walsh。上面具體的程式碼是 Quake 原始碼的改編版(Quake 原始碼是所有討論的源頭)。

這篇文章探討了這個演算法諸多有趣的地方,解釋了它的原理,對指數在 -1 到 1 的情況進行了擴充,還有一些數學相關的新思路。

(本文確實包含許多數學知識,你可以把這些方程看做註釋,無需閱讀即可掌握本文要點,但如果你想要完整的閱讀本文,驗證我的觀點的正確性,則需要關注這些方程。)

為什麼?

為何需要計算倒平方根,甚至於值得實現一個怪異的演算法來加速計算?因為它一直是 3D 程式設計計算中的一部分。在 3D 圖形中,你使用平面法線,長度為 1 的三座標向量,來表示光線和反射。你會使用很多平面法線,計算它們時需要對向量進行標準化。如何標準化一個向量呢?每個座標分別除以向量的長度,因此,每個座標需乘上

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

計算 數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df 相對來說開銷很小,計算平方根並做除法開銷很大。進入FastInvSqrt

如何做到的?

這個函式究竟做了什麼計算出了結果?它有 4 個主要步驟。首先它把浮點型別輸入的二進位制表示重定義為一個整形數字的二進位制表示。

運用得到的結果進行整形算術計算,得到我們所求值的近似值:

儘管結果不是近似值本身,但如果你把這個整形數字的二進位制表示重定義為一個浮點數,那麼就會得到近似值。所以程式碼做了與步驟一相反的變換回到了浮點數:

最終進行一次牛頓法迭代來提高精度。

這便得到 x 倒平方根的非常好的近似值。最後一步牛頓法相對比較直觀,我便不多花時間講解。關鍵步驟是第二步:對原始浮點數轉化來的整形數進行算術運算,得到一個有意義的結果,下面重點關注這一步。

怎麼回事?

這部分解釋第二步背後的數學原理(下面推導的第一部分,至常數值的計算,最早由 McEniry 發現的)。

在進入這一有趣的部分之前,我先快速解釋下標準浮點數編碼,我只講解我要用的部分,完整的背景知識請參見維基百科。浮點數由三部分組成:符號,指數和尾數。這是單精度(32 位)浮點數二進位制表示:

最高位是符號,接下來的 8 位是指數,底部 23 位是尾數。由於將要計算的平方根只對正數定義,所以從現在起假設符號位是 0 。

當把一個浮點數看作一堆 0 、1 時,指數和底數便只是普通的正整數,沒什麼特別的。把它們記為 EM(會經常使用)。另一方面,當把這些 0、1 解釋為浮點值時,尾數則是介於 0 到 1 之間的值,全 0 意味著 0,全 1 是非常接近但略微小於 1 的數。並非將指數當做 8 位無符號整形使用,而是減去一個偏量 B,使之成為介於 -127 到 128 之間的有符號整形。把這些值的浮點解釋記為 em。總之,我按照 McEnry 把整形解釋相關的值用大寫字母表示,把浮點解釋相關的值用小寫字母表示。

這兩種解釋之間的轉換很直接:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

對 32 位浮點數來說, L 是 223B 是 127。給定 e m,浮點數的值可計算得到:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

相應地整形解釋的數是

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

現在解釋這個演算法的所有基礎幾乎都有了,因此可以開始了,遺漏的部分會在解釋過程中引入。給定輸入 x,想要計算的是它的倒平方根,或者

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

我們先對等式兩邊取以 2 為底的對數,至於原因很快就會知道:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

既然我們進行計算的值都是浮點數,則可以把 x 和 y 用各自的浮點組成替換:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

對數比較麻煩,幸運的是可以很輕鬆的擺脫它們,但在這之前,需暫停一下來解釋其原理。

在方程的兩邊都含有這樣的項:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

其中 v 介於 0 到 1 之間,也僅僅當 v 在 0 到 1 之間時,這個函式和一條直線很接近:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

或者方程形式:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

其中 σ 是一個可選常數,儘管不能完美匹配但可以調節 σ 使得二者非常接近。通過它我們可以把上述包含對數的方程近似為一個線性方程,計算起來更輕鬆:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

到達關鍵的部分了!這時,用整形解釋下的指數和尾數替代浮點解釋下的表示進行計算就很方便了:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

通過少許步驟進行移項合併,就會得到很熟悉的形式(細節很乏味,可以跳過):

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

走完最後一步之後,有趣的事情發生了:在這些雜亂的項中,方程兩邊都有整形解釋的值:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

換句話說,y 的整形解釋等於某個常數減去 x 的整形解釋的一半,或者用 C 程式碼表示:

K 看起來很熟悉對吧?

剩下的工作就是找到這個常數。我們已經知道 B L 的值,σ 的值還不知道。記住,σ 是對對數函式近似的調節因子,所以它的取值有些自由。我選取最初應用過的一個值,0.0450465,得到:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

想知道這個值的 16 進製表示?0x5f3759df。(當然理應是它,因為選的特定的 σ 得到的),這個常數不是一個位模式,這點你可能通過它被寫成 16 進位制形式而推知,它只是一個普通計算的結果化為整形的形式。

但是正如 Knuth 所說,目前為止我們只是證明這個方法應該有效,還沒檢測。為了對這個近似的精確程度有直觀的認識,用它的圖線和精確的倒平方根的圖線進行對比:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

這裡輸入的值介於 0 到 100 之間。很精確對嗎?理應如此。這不僅僅是神奇,正如上面介紹的其推導,它的計算也很奇怪,但都是清晰、有意義的操作——浮點和整形之間的位轉換

等等,還有更多!

這個操作的推導告訴你的不只是一個常數的值,你會注意到這個推導不依賴於任何項的具體值,它們只是用於變換的常數,這意味著即使我們改變它們,推導也成立。

首先,這個計算不關心 LB 具體是什麼,它們通過浮點表示給出,這意味著對 64 位和 128 位浮點數可以使用同樣的方法,所有要做的只是重新計算常數,這是唯一依賴它們的部分。

第二,它不關心 σ 的取值,使得對數函式和 x+σ 之間差異最小的 σ 可能不會得到最精確的近似,事實確實如此。這是浮點數的舍入和牛頓法聯合造成的,σ 取值本身是一個有趣的課題,McEniry 和 Lomont 有過論述。

最後,不僅僅只能計算倒平方根,這裡指數碰巧是 -1/2,但是對於 -1 到 1 之間的指數,推導依然成立。我們把指數記為 p(因為 e 已經使用),替代掉 -1/2 作同樣的推導,得到:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

來試試其他一些指數,首先 p=0.5,普通平方根計算:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df
數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

或者程式碼形式,

同樣有效嗎?當然:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

這也許是近似計算平方根的一個很著名的方法,但 Google 和維基百科的一個粗略調查顯示它並不是。

甚至對奇數次開方也適用,如立方根

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

相應地:

由於這個奇數因子,我們無法使用移位來代替乘法。同樣近似很接近:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

這時,你可能會注意到,當指數改變時,我們只需做非常簡單的工作:只是通過改變一個線性因子來調整常數值,改變和輸入的整形解釋相乘的因子。這些操作開銷都不大,因此在執行時進行是可行的,而不需要重新計算。如果我們事先對其他兩個因子的乘積進行計算:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

則不用提前知道指數,也可計算:

對項進行整理,還可以減少一步乘法計算:

這便是對 -1 到 1 之間任意指數進行快速指數計算的神奇部分,現在我們還差一部分,就是對牛頓近似進行一般化,使得快速指數計算函式對所有指數都適用,且和之前的倒平方根函式一樣準確。這部分我沒深入,就交給其他部落格文章了(很可能是其他人的而不是我)。

上面表示式包含了一個新的神奇的常數,  0x3f7a3bea。即使它在某種程度上比最初的常數更神奇,最初的常數可由它演變而來,但它依賴於 σ 的主觀選取,所以它也不具通用性。把這個常數記為 Cσ ,馬上我們會更進一步探究它。

但首先,我們可以對 p=0 時的公式進行合理性檢驗。對於 p 取零時,結果應該始終為 1,因為 x0 等於 1,與 x 無關。確實,第二項消失了,因為它乘以了 0 ,我們得到簡潔的:

確實是常數,當解釋為浮點值時等於 0.977477,近似為 1,所以合理性檢驗成立。這也告訴我們一些資訊: 當 Cσ  轉化為浮點解釋時,是一個有意義的值,為 1 或非常接近 1。

這很有趣,更進一步,Cσ 的整形解釋是:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

這很像但不是一個浮點數的形式,唯一的問題是這裡是減去而不是加上第二項,這個問題可以很輕易解決:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

現在它看起來像一個浮點數的整形解釋,為了具體來看,我們先分辨出指數和尾數,然後計算浮點數值,cσ, 這是指數:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

這是尾數:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

所以常數的浮點數值是:

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

確實如果你先用 σ 作除法,0.0450465 除以 2 得 0.02252325,用 1 減去它的 0.97747675 或者剛才的“近似為 1”。這讓我們從另一種角度來看待 Cσ,作為一個浮點數的整形解釋。用程式碼計算它:

注意,對於固定的 σ , 這些皆為常數,編譯器應該能夠優化整個計算過程。結果為 0x3f7a3beb,不是之前的 0x3f7a3bea ,但只有一位元組的不同(最次要的一個位元組),這對浮點數計算很正常。要得到最初的常數,本文的標題,只需把結果乘以 1.5。

我們已經足夠深入了,至少我很滿足了,沒有什麼可以繼續探究的了。對於我來說,通過這個練習,主要學到的是整形和浮點型之間的位轉換操作不是沒有意義的,它很怪異但卻很經濟的數字操作,在計算中很有用。我認為它還有更多有待發現的用處。

打賞支援我翻譯更多好文章,謝謝!

打賞譯者

打賞支援我翻譯更多好文章,謝謝!

任選一種支付方式

數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df 數學之美:平方根倒數速演算法中的神奇數字 0x5f3759df

相關文章