文章出處

概述

        平方根倒數速算法，是用于快速計算1/Sqrt(x)的值的一種算法，在這里x需取符合IEEE 754標準格式的32位正浮點數。讓我們先來看這段代碼：

 1 float Q_rsqrt( float number )
 2 {
 3     long i;
 4     float x2, y;
 5     const float threehalfs = 1.5F;
 6 
 7     x2 = number * 0.5F;
 8     y  = number;
 9     i  = * ( long * ) &y;                        // evil floating point bit level hacking
10     i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
11     y  = * ( float * ) &i;
12     y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
13 //    y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed
14 
15 #ifndef Q3_VM
16 #ifdef __linux__
17     assert( !isnan(y) ); // bk010122 - FPE?
18 #endif
19 #endif
20     return y;
21 }

        2010年，John Carmack公開了Quake III Arena的源代碼，這段代碼便出自其中，這種算法在游戲《雷神之錘III競技場》（Quake III Arena）被廣泛應用，作為《雷神之錘》3D游戲引擎的作者，該算法理所應當的被認為該是Carmark所創，但是后來Carmack在對Beyond3D上關于該段代碼的作者的討論回信表示這段代碼并不是出自他之手，也許是其開發小組的另一位成員所寫，而至今該算法的最終起源仍然是個謎。

魔數0x5f3759df

        此算法最為牛叉之處在于注釋了what the fuck?的那一行代碼，此處將浮點數作為整數右移一位，通過一個神秘的十六進制數0x5f3759df進行減法運算，便完成了牛頓迭代需要的第一次猜測值，而這次猜測值已經很接近滿足精度的結果了，因此只通過了一次牛頓迭代就完成了所需精度的結果，其實就相當于直接進行計算就得到了所需結果。

32位浮點數在計算機中的存儲方式

        如果要理解該段代碼除了魔數一行的代碼，就必須理解浮點數在計算機中的存儲方式。任何數據在計算機內存中都是以二進制形式存在的，浮點數也不例外， 32位的浮點數在計算機中的存儲方式如下圖，分為三個部分：

 

  1. 標志位（Sign），第31位為標志位，1代表負數，0代表正數；
  2. 指數部分（Exponent），指數部分共占8位，而 float類型規定其偏移量為127，由于指數可正可負，對于8位二進制數其表示范圍為[-128, 127]；
  3. 尾數部分（Mantissa），表示有效數字位，由于尾數部分的整數部分恒為1，該位將不被存儲，因此原本23位的尾數部分可以表示的精度就變成了24位；

因此在二進制科學計數法中，我們可以將任意的32位浮點數表示為：

而平方根倒數速算法中的第9行的目的就是將浮點數轉化為32位整型數的表現形式。

牛頓迭代法

        牛頓迭代法又稱為牛頓-拉菲森法（Newton-Raphson method），這種算法的原理就是一個求近似解的方法，該算法正是采用該方法來迭代平方根倒數的。

        要求F(x)=0的解，首先選取一個接近函數F(x)零點的點(x0, F(x0))，過該點做切線求得其與X軸交點的x的值記為x1，該值通常會比x0更接近方程的解，然后不斷使用新的點進行迭代至滿足的精度為止。

        在(x0, F(x0))點的斜率為，過該點的切線方程為，因此求得該切線與x軸的交點橫坐標值，因此簡化后的迭代公式為：

 以求該平方根倒數求值為例，通過該方法轉化為求方程的解，申明，根據牛頓法，第一次迭代結果為：

而這也正好是算法中的第12行代碼對該算法的應用。

Lomont的逆向數學推理

        對于該算法最讓人難以理解的就是魔數0x5f3759df，對于這個問題，Purdue大學的數學家Chris Lomont在論壇上得知后覺得非常有趣，于是決定研究一下它的工作原理。

        Lomont在假設該算法成立的前提下（事實上已被廣泛證實該算法的正確性），希望通過數學推理推導出該值。下面簡要介紹一下Lomont的推導過程：

        由于是求平方根倒數，因此假設所有浮點數的標志位均為0，浮點數二進制科學計數公式簡化為：

我們假設所求的魔數為R，那么第一次猜測值就可以表示為，其中i表示浮點數x的整數形式，R為整型，如下圖：

在對i進行按位右移的過程中，因為i的二進制指數位可能會被右移到尾數部分，所以必須按兩種情況分類討論：

1. 當i的指數部分為偶數，那么指數部分的最低位為0，指數位將不會被右移到尾數部分中，而真正的指數e = E - 127就是奇數，令e = 2d + 1，那么在經過運算后，y0的指數部分的值為：

由于我們所求的平方根倒數大于0，所以新的指數部分也必須大于0，且E的取值在[0, 254]之間，所以R1 ≥ 128；因為E是偶數，所以指數位并沒有移動到尾數中，所以運算之后尾數y0的尾數部分為：，

如果，那么初始猜測值為：

如果，二進制減法正如十進制中的減法運算，需向指數部分借位，此時，

定義

那么

因為e = 2d + 1，那么我們要求解的平方根倒數的近似值為：

所以y和y0的相對誤差則為：

那么根據要求相對誤差需小于0.125，則可以推導出R2的取值范圍(189.2, 190.7)，所以R2=190=0xbe；

2. 當i的指數部分為奇數時，同理推導出：

如果我們希望以上兩種情況的相對誤差均小于0.125， R1確定且與E的奇偶無關，那么R1=190=0xbe，此時重新定義以上兩個誤差函數：

最終通過當M=0，1，以及R2=M/2的幾個臨界點獲取相對誤差最大值函數，然后在所有函數中獲取使得這些誤差函數值都最小的R2的值，這個值便是我們要找的結果，分類討論一共得到了9個關于R2(取值范圍是[0, 1))的函數：

 

最終這個求值過程交給Mathematica來完成，也不知為何我最后求得的結果和Lomont略有差異：

r0 = 0.432744889959443195468521587014

我求出的結果：

r0 = 0.432744834277619894180588744347

其他討論

        Lomont求解魔數0x5f3759df的推導過程的前提是假設算法中

i  = 0x5f3759df - ( i >> 1 );

這行代碼成立，通過反推來求解魔數的，但是這行代碼最初是如何得出的，依據是什么，Lomont并沒有給出更多的解釋或猜測，好在其他的論壇也有相關的討論。

        在一則slashdot.org關于Origin of Quake3's Fast InvSqrt()的討論中，kent.dickey的回復令人印象深刻，而我個人也非常贊同他的思路，非常簡單的推導，也許他的這種思路才是該算法的最初雛形。

        他在回復中說：對于浮點數平方根倒數的求解，指數部分的初始猜測值就是對指數部分減求反（暫時忽略尾數部分，因為對于2以內的求解，一步就可以得到結果），比如求16的平方根倒數：

但是由于在內存中浮點數特殊的存儲方式，因此在內存中指數的整型值為：

NewExpIntValueInMemory=(ActualExp+127) << 23

我們想要的指數結果是：

 NewActualExp=-ActualExp/2

由于在內存中：Exp=ActualExp+127，那么ActualExp=Exp-127，所以：

NewExp-127=-(Exp-127)/2

那么：

NewExp=127-(Exp-127)/2=(127×3)/2-Exp/2

那 么指數在內存中指數的整型值為：

所以在完全忽略尾數部分的前提下，初始猜測近似值為：

i=0x5F400000-i≫1

隨后kent.dickey試圖使用類似方法找到對于尾數部分的類似特定模式，遺憾的是沒能找到，也許作者當初在有了這個思路（如果該思路真的就是原算法的雛形）之后，也采用了類似Lomont的暴力求解來獲得更精確的初始猜測。

參考文獻


文章列表


不含病毒。www.avast.com
arrow
arrow
    全站熱搜
    創作者介紹
    創作者 大師兄 的頭像
    大師兄

    IT工程師數位筆記本

    大師兄 發表在 痞客邦 留言(0) 人氣()

    E-mail轉寄