0x5f3759df這個快速開方中的常數的數學依據是什麼？

01-14

看了關於卡馬克和雷神之錘3的故事，網上的幾篇文章並沒有打消我的疑問
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; }
代碼中標識WTF的地方神奇地完成了核心運算，這個是靠猜的嗎？如果不是猜的依據在哪裡？如果是猜的為何會這麼准？

實名反對一樓 @Yong He 的答案。

單精度的浮點數結構是這樣：

那麼現在，每個正浮點數 y 可以用尾數和指數的形式寫成 $2^e(1+m)$ ，其中 m 是尾數部分，取值範圍是 $[0, 1)$ ；e 是指數部分，一個整數。每個浮點數所對應的「整數形式」則可以用 $EL+M$ 表示，其中 L 是指數部分需要的位移次數（用 2 的冪表示），E 和 M 是指數部分和小數部分的整數版本。兩者之間的關係是

$left{ egin{array}{ll} E = e + B \ M = mL end{array} ight.$

對單精度浮點而言， $L=2^{23},B=127$ 。

考慮對數 $log_2 y=e+log_2(1+m)$ ，由於 $minleft[0, 1 ight)$ ， $log_2(1+m)inleft[0,1 ight)$ ，取近似 $log_2(1+m)approx m+delta$ ，可以算出整體「偏差」最小的 $delta=0.0430357$ ，此時兩者基本相當。因此我們可以說

$log_2 yapprox e+m+delta$ ……………… (1)

那麼，對於 y 的整數形式 Y 而言，展開並帶代入 (1) 有：

$Y = EL+M\ = L(e+B+m) \ approx L(log_2 y + B - delta)$

即

$log_2 yapprox {Y over L}-(B-delta)$ ………………(2)

那麼對於 $y$ 來說， $log_2{y$ ，代入 (2) 得：

${Y$

解得

$Y$ 。

這個就是代碼

i = 0x5F3759DF - ( i &>&> 1 );

的秘密所在。0x5F3759DF 的具體取值是根據 δ 變的，至於卡馬克為啥用 0x5F3759DF 而不是其他相近的值，估計是做實驗測的吧……

牛頓迭代法

0x5f3759df？這是個什麼東西？學過數值分析就知道，演算法裡面求平方根一般採用的是無限逼近的方法，比如牛頓迭代法，抱歉當年我數值分析學的太爛，也講不清楚。簡單來說比如求5的平方根，選一個猜測值比如2，那麼我們可以這麼算


5/2 = 2.5; 2.5+2/2 = 2.25; 5/2.25 = xxx; 2.25+xxx/2 = xxxx ...

這樣反覆迭代下去，結果必定收斂於sqrt(5)，沒錯，一般的求平方根都是這麼算的

。而卡馬克的不同之處在於，他選擇了一個神秘的猜測值0x5f3759df作為起始，使得

整個逼近過程收斂速度暴漲，對於Quake III所要求的精度10的負三次方，只需要一

次迭代就能夠得到結果。 
好吧，如果這還不算牛b，接著看。 
普渡大學的數學家Chris Lomont看了以後覺得有趣，決定要研究一下卡馬克弄出來的

這個猜測值有什麼奧秘。Lomont也是個牛人，在精心研究之後從理論上也推導出一個

最佳猜測值，和卡馬克的數字非常接近, 0x5f37642f。卡馬克真牛，他是外星人嗎？ 
傳奇並沒有在這裡結束。Lomont計算出結果以後非常滿意，於是拿自己計算出的起始

值和卡馬克的神秘數字做比賽，看看誰的數字能夠更快更精確的求得平方根。結果是

卡馬克贏了... 誰也不知道卡馬克是怎麼找到這個數字的。

最後Lomont怒了，採用暴力方法一個數字一個數字試過來，終於找到一個比卡馬克數字要好上那麼一丁點的數字，雖然實際上這兩個數字所產生的結果非常近似，這個暴力得出的數字是0x5f375a86。

假設原浮點數為0EM,

E為8位指數部分，M為23位尾數部分。

位操作的主要目的是為了讓新的浮點數的e = -(e0)/2

我們知道E = e0+127

而我們想讓位操作後的E"=-e0/2+127

右移一位可以得到：

E_s = E/2 = e0/2 + 127/2

如果我們用190去減E_s, 就可以得到

E_1 = 190-E_s=190-e0/2-63=-e0/2+127

而190正是這個魔法數的第31-23位。

因此這個魔法數幫助我們得到了正確的新浮點數的指數部分。小數部分則比較複雜，因為在原指數部分為奇數的情況下右移一位會無端端給尾數增加1/2，且相減後尾數可能會借用指數部分導致指數變小，因此通過調整尾數來修正。如需詳細了解可參考這篇paper: http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf

需要指出現在的Intel CPU已經提供快速rsqrt指令，所以無需自己去實現這個演算法。

Geometry tools 上面有論文，包括詳細的推導過程和數值分析，以及magic number的選取方法

首先，豆瓣和果殼已經有過一些解答了。網址：

http://m.guokr.com/post/90718/?page=2

還有篇論文解釋過了：http://pan.baidu.com/share/link?shareid=1180931030uk=2989683554

下面，說說個人的見解，其實數學之中有許多巧合，比如∑1/n2=？，這在當時難倒了很多數學家，只好強行運算，算到小數後面第七位的時候是1.6449341，發現和π2/6很像,於是就猜測是不是就等於π2/6，後來歐拉證明了確實如此，但是據說歐拉給出的證明方法是錯的。

我就想問那個發現和π2/6很像的人，你沒事記著π2/6的值幹嘛？！太巧了！

有的時候不得不相信巧合。

記得Purdue的Chris Lomont後來暴力算了一個比這個稍微好一點的數發了篇paper 總之就是最開始找出這個數太無解了。。。。

我記得matrix67博客上討論過，好像是二元擬合出來的

這個值好像最早是在cray超級計算機中出現的