歡迎來到Linux教程網
Linux教程網
Linux教程網
Linux教程網
您现在的位置: Linux教程網 >> UnixLinux >  >> Linux編程 >> Linux編程

快速平方根算法

在3D圖形編程中,經常要求平方根或平方根的倒數,例如:求向量的長度或將向量歸一化。C數學函數庫中的sqrt具有理想的精度,但對於3D游戲程式來說速度太慢。我們希望能夠在保證足夠的精度的同時,進一步提高速度。

Carmack在QUAKE3中使用了下面的算法,它第一次在公眾場合出現的時候,幾乎震住了所有的人。據說該算法其實並不是Carmack發明的,它真正的作者是Nvidia的Gary Tarolli(未經證實)。

//

// 計算參數x的平方根的倒數

//

float InvSqrt (float x)

{

      float xhalf = 0.5f*x;

      int i = *(int*)&x;

      i = 0x5f3759df - (i >> 1);    // 計算第一個近似根

      x = *(float*)&i;

      x = x*(1.5f - xhalf*x*x);      // 牛頓迭代法

      return x;

}

該算法的本質其實就是牛頓迭代法(Newton-Raphson Method,簡稱NR),而NR的基礎則是泰勒級數(Taylor Series)。NR是一種求方程的近似根的方法。首先要估計一個與方程的根比較靠近的數值,然後根據公式推算下一個更加近似的數值,不斷重復直到可以獲得滿意的精度。其公式如下:

函數:y=f(x)

其一階導數為:y""=f""(x)

則方程:f(x)=0 的第n+1個近似根為

x[n+1] = x[n] - f(x[n]) / f""(x[n])

NR最關鍵的地方在於估計第一個近似根。如果該近似根與真根足夠靠近的話,那麼只需要少數幾次迭代,就可以得到滿意的解。


現在回過頭來看看如何利用牛頓法來解決我們的問題。求平方根的倒數,實際就是求方程1/(x^2)-a=0的解。將該方程按牛頓迭代法的公式展開為:

x[n+1]=1/2*x[n]*(3-a*x[n]*x[n])

將1/2放到括號裡面,就得到了上面那個函數的倒數第二行。


接著,我們要設法估計第一個近似根。這也是上面的函數最神奇的地方。它通過某種方法算出了一個與真根非常接近的近似根,因此它只需要使用一次迭代過程就獲得了較滿意的解。它是怎樣做到的呢?所有的奧妙就在於這一行:

i = 0x5f3759df - (i >> 1);    // 計算第一個近似根

超級莫名其妙的語句,不是嗎?但仔細想一下的話,還是可以理解的。我們知道,IEEE標准下,float類型的數據在32位系統上是這樣表示的(大體來說就是這樣,但省略了很多細節,有興趣可以GOOGLE):

bits:31 30 ... 0

31:符號位

30-23:共8位,保存指數(E)

22-0:共23位,保存尾數(M)

所以,32位的浮點數用十進制實數表示就是:M*2^E。開根然後倒數就是:M^(-1/2)*2^(-E/2)。現在就十分清晰了。語句i>>1其工作就是將指數除以2,實現2^(E/2)的部分。而前面用一個常數減去它,目的就是得到M^(1/2)同時反轉所有指數的符號。


至於那個0x5f3759df,呃,我只能說,的確是一個超級的Magic Number。


那個Magic Number是可以推導出來的,但我並不打算在這裡討論,因為實在太繁瑣了。簡單來說,其原理如下:因為IEEE的浮點數中,尾數M省略了最前面的1,所以實際的尾數是1+M。如果你在大學上數學課沒有打瞌睡的話,那麼當你看到(1+M)^(-1/2)這樣的形式時,應該會馬上聯想的到它的泰勒級數展開,而該展開式的第一項就是常數。下面給出簡單的推導過程:

對於實數R>0,假設其在IEEE的浮點表示中,

指數為E,尾數為M,則:

R^(-1/2)

= (1+M)^(-1/2) * 2^(-E/2)

將(1+M)^(-1/2)按泰勒級數展開,取第一項,得:

原式

= (1-M/2) * 2^(-E/2)

= 2^(-E/2) - (M/2) * 2^(-E/2)

如果不考慮指數的符號的話,

(M/2)*2^(E/2)正是(R>>1),

而在IEEE表示中,指數的符號只需簡單地加上一個偏移即可,

而式子的前半部分剛好是個常數,所以原式可以轉化為:

原式 = C - (M/2)*2^(E/2) = C - (R>>1),其中C為常數

所以只需要解方程:

R^(-1/2)

= (1+M)^(-1/2) * 2^(-E/2)

= C - (R>>1)

求出令到相對誤差最小的C值就可以了

上面的推導過程只是我個人的理解,並未得到證實。而Chris Lomont則在他的論文中詳細討論了最後那個方程的解法,並嘗試在實際的機器上尋找最佳的常數C。有興趣的朋友可以在文末找到他的論文的鏈接。


所以,所謂的Magic Number,並不是從N元宇宙的某個星系由於時空扭曲而掉到地球上的,而是幾百年前就有的數學理論。只要熟悉NR和泰勒級數,你我同樣有能力作出類似的優化。


在GameDev.net上有人做過測試,該函數的相對誤差約為0.177585%,速度比C標准庫的sqrt提高超過20%。如果增加一次迭代過程,相對誤差可以降低到e-004的級數,但速度也會降到和sqrt差不多。據說在DOOM3中,Carmack通過查找表進一步優化了該算法,精度近乎完美,而且速度也比原版提高了一截(正在努力弄源碼,誰有發我一份)。


值得注意的是,在Chris Lomont的演算中,理論上最優秀的常數(精度最高)是0x5f37642f,並且在實際測試中,如果只使用一次迭代的話,其效果也是最好的。但奇怪的是,經過兩次NR後,在該常數下解的精度將降低得非常厲害(天知道是怎麼回事!)。經過實際的測試,Chris Lomont認為,最優秀的常數是0x5f375a86。如果換成64位的double版本的話,算法還是一樣的,而最優常數則為0x5fe6ec85e7de30da(又一個令人冒汗的Magic Number - -b)。


這個算法依賴於浮點數的內部表示和字節順序,所以是不具移植性的。如果放到Mac上跑就會掛掉。如果想具備可移植性,還是乖乖用sqrt好了。但算法思想是通用的。大家可以嘗試推算一下相應的平方根算法。


下面給出Carmack在QUAKE3中使用的平方根算法。Carmack已經將QUAKE3的所有源代碼捐給開源了,所以大家可以放心使用,不用擔心會收到律師信。

//

// Carmack在QUAKE3中使用的計算平方根的函數

//

float CarmSqrt(float x){

      union{

            int intPart;

            float floatPart;

      } convertor;

      union{

            int intPart;

            float floatPart;

      } convertor2;

      convertor.floatPart = x;

      convertor2.floatPart = x;

      convertor.intPart = 0x1FBCF800 + (convertor.intPart >> 1);

      convertor2.intPart = 0x5f3759df - (convertor2.intPart >> 1);

      return 0.5f*(convertor.floatPart + (x * convertor2.floatPart));

}

另一個基於同樣算法的更高速度的sqrt實現如下。其只是簡單地將指數除以2,並沒有考慮尾數的方根。要看懂該代碼的話必須知道,在IEEE浮點數的格式中,E是由實際的指數加127得到的。例如,如果實數是0.1234*2^10,在浮點表示中,E(第23-30位)的值其實為10+127=137。所以下面的代碼中,要處理127偏移,這就是常數0x3f800000的作用。我沒實際測試過該函數,所以對其優劣無從評論,但估計其精度應該會降低很多。

float Faster_Sqrtf(float f)

{

      float result;

      _asm

      {

            mov eax, f

            sub eax, 0x3f800000

            sar eax, 1

            add eax, 0x3f800000

            mov result, eax

      }

      return result;

}

除了基於NR的方法外,其他常見的快速算法還有多項式逼近。下面的函數取自《3D游戲編程大師技巧》,它使用一個多項式來近似替代原來的長度方程,但我搞不清楚作者使用的公式是怎麼推導出來的(如果你知道的話請告訴我,謝謝)。

//

//    這個函數計算從(0,0)到(x,y)的距離,相對誤差為3.5%

//

int FastDistance2D(int x, int y)

{

      x = abs(x);

      y = abs(y);

      int mn = MIN(x,y);

      return(x+y-(mn>>1)-(mn>>2)+(mn>>4));

}

//

// 該函數計算(0,0,0)到(x,y,z)的距離,相對誤差為8%

//

float FastDistance3D(float fx, float fy, float fz)

{

      int temp;

      int x,y,z;

      // 確保所有的值為正

      x = int(fabs(fx) * 1024);

      y = int(fabs(fy) * 1024);

      z = int(fabs(fz) * 1024);

      // 排序

      if (y < x) SWAP(x,y,temp)

      if (z < y) SWAP(y,z,temp)

      if (y < x) SWAP(x,y,temp)

      int dist = (z + 11 * (y >> 5) + (x >> 2) );

      return((float)(dist >> 10));

}

還有一種方法稱為Distance Estimates(距離評估?),如下圖所示:

紅線所描繪的正八邊形上的點為:

octagon(x,y) = min((1/√2) * (|x|+|y|), max(|x|,|y|))

求出向量v1和v2的長度,則:

√(x^2+y^2) = (|v1|+|v2|)/2 * octagon(x,y)

到目前為止我們都在討論浮點數的方根算法,接下來輪到整數的方根算法。也許有人認為對整型數據求方根無任何意義,因為會得到類似99^(1/2)=9的結果。通常情況下確實是這樣,但當我們使用定點數的時候(定點數仍然被應用在很多系統上面,例如任天堂的GBA之類的手持設備),整數的方根算法就顯得非常重要。對整數開平方的算法如下。我並不打算在這討論它(事實是我也沒有仔細考究,因為在短期內都不會用到- -b),但你可以在文末James Ulery的論文中找到非常詳細的推導過程。

//

// 為了閱讀的需要,我在下面的宏定義中添加了換行符

//

#define step(shift)

if((0x40000000l >> shift) + sqrtVal <= val)

{

      val -= (0x40000000l >> shift) + sqrtVal;

      sqrtVal = (sqrtVal >> 1) | (0x40000000l >> shift);

}

else

{

      sqrtVal = sqrtVal >> 1;

}

//

// 計算32位整數的平方根

//

int32 xxgluSqrtFx(int32 val)

{

      // Note: This fast square root function

      // only works with an even Q_FACTOR

      int32 sqrtVal = 0;

      step(0);

      step(2);

      step(4);

      step(6);

      step(8);

      step(10);

      step(12);

      step(14);

      step(16);

      step(18);

      step(20);

      step(22);

      step(24);

      step(26);

      step(28);

      step(30);

      if(sqrtVal < val)

      {

            ++sqrtVal;

      }

      sqrtVal <<= (Q_FACTOR)/2;

      return(sqrtVal);

}

sqrt2003 GameDev.net,,FPU,sqrt,SIMD,,4 

Copyright © Linux教程網 All Rights Reserved