性能提高(3) (转载)经典重弹,不能不弹 QUAKE3(一款有名游戏)中使用的平方根算法

快速平方根倒数算法
本文介绍了一种在3D图形编程中用于快速计算平方根倒数的算法,该算法首次出现在QUAKE3中,利用牛顿迭代法并选取一个特殊常数进行优化,大幅提升了计算效率。

 

雷神之III》里求平方根倒数的函数(快速平方根(倒数)算法)

3D程中,常要求平方根或平方根的倒数,例如:求向量的度或将向量一化。C数学函数中的sqrt具有理想的精度,但3D程式来速度太慢。我希望能在保的精度的同提高速度。

Carmack
QUAKE3中使用了下面的算法,它第一次在公众合出候,几乎震住了所有的人。说该算法其并不是Carmack明的,它真正的作者是NvidiaGary Tarolli(未经证实)。

Java

  1.   
  2. //    
  3. // 算参数x的平方根的倒数    
  4. //    
  5. float InvSqrt (float x)    
  6. {    
  7. float xhalf = 0.5f*x;    
  8. int i = *(int*)&x;    
  9. i = 0x5f3759df - (i >> 1); // 算第一个近似根    
  10. x = *(float*)&i;    
  11. x = x*(1.5f - xhalf*x*x); // 迭代法    
  12. return x;    
  13. }   

 

//

// 算参数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 MethodNR),而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%,速度比Csqrt提高超20%。如果增加一次迭代程,相对误差可以降低到e-004数,但速度也会降到和sqrt差不多。据DOOM3中,Carmack过查找表步优化了算法,精度近乎完美,而且速度也比原版提高了一截。

得注意的是,在Chris Lomont的演算中,理上最秀的常数(精度最高)是0x5f37642f,并且在实际测试中,如果只使用一次迭代的,其效果也是最好的。但奇怪的是,经过两次NR后,在常数下解的精度将降低得非常害(天知道是怎回事!)。经过实际测试Chris Lomont认为,最秀的常数是0x5f375a86。如果64位的double版本的,算法是一的,而最常数则为0x5fe6ec85e7de30da(又一个令人冒汗的Magic Number - -b)。

个算法依于浮点数的内部表示和字节顺序,所以是不具移植性的。如果放到Mac上跑就会挂掉。如果想具可移植性,是乖乖用sqrt好了。但算法思想是通用的。大家可以尝试推算一下相的平方根算法。

下面CarmackQUAKE3中使用的平方根算法。CarmackQUAKE3的所有源代给开源了,所以大家可以放心使用,不用担心会收到律信。

Java

  1.   
  2. //    
  3. // CarmackQUAKE3中使用的算平方根的函数    
  4. //    
  5. float CarmSqrt(float x){    
  6. union{    
  7. int intPart;    
  8. float floatPart;    
  9. } convertor;    
  10. union{    
  11. int intPart;    
  12. float floatPart;    
  13. } convertor2;    
  14. convertor.floatPart = x;    
  15. convertor2.floatPart = x;    
  16. convertor.intPart = 0x1FBCF800 + (convertor.intPart >> 1);    
  17. convertor2.intPart = 0x5f3759df - (convertor2.intPart >> 1);    
  18. return 0.5f*(convertor.floatPart + (x * convertor2.floatPart));    
  19. }   

 

//

// CarmackQUAKE3中使用的算平方根的函数

//

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的作用。我没实际测试过该函数,所以劣无从评论,但估其精度应该会降低很多。

Java

  1.   
  2. float Faster_Sqrtf(float f)    
  3. {    
  4. float result;    
  5. _asm    
  6. {    
  7. mov eax, f    
  8. sub eax, 0x3f800000    
  9. sar eax, 1    
  10. add eax, 0x3f800000    
  11. mov result, eax    
  12. }    
  13. return result;    
  14. }   

 

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戏编程大技巧》,它使用一个多式来近似替代原来的度方程,但我搞不清楚作者使用的公式是怎出来的(如果你知道的话请我,谢谢)。

Java

  1.   
  2. //    
  3. // 个函数算从(0,0)(x,y)的距离,相对误3.5%    
  4. //    
  5. int FastDistance2D(int x, int y)    
  6. {    
  7. x = abs(x);    
  8. y = abs(y);    
  9. int mn = MIN(x,y);    
  10. return(x+y-(mn>>1)-(mn>>2)+(mn>>4));    
  11. }    
  12. //    
  13. // 函数(0,0,0)(x,y,z)的距离,相对误8%    
  14. //    
  15. float FastDistance3D(float fx, float fy, float fz)    
  16. {    
  17. int temp;    
  18. int x,y,z;    
  19. // 确保所有的值为    
  20. x = int(fabs(fx) * 1024);    
  21. y = int(fabs(fy) * 1024);    
  22. z = int(fabs(fz) * 1024);    
  23. // 排序    
  24. if (y < x) SWAP(x,y,temp)    
  25. if (z < y) SWAP(y,z,temp)    
  26. if (y < x) SWAP(x,y,temp)    
  27. int dist = (z + 11 * (y >> 5) + (x >> 2) );    
  28. return((float)(dist >> 10));    
  29. }   

 

//

// 个函数算从(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|))

求出向量v1v2
度, √(x^2+y^2) = (|v1|+|v2|)/2 * octagon(x,y)

到目前
止我都在讨论浮点数的方根算法,接下来到整数的方根算法。有人认为对整型数据求方根无任何意,因会得到99^(1/2)=9果。通常情况下确这样,但当我使用定点数的候(定点数仍然被用在很多系上面,例如任天堂的GBA的手持设备),整数的方根算法就得非常重要。整数平方的算法如下。我并不打算在这讨论它(事是我也没有仔考究,因在短期内都不会用到- -b),但你可以在文末James Ulery文中找到非常详细的推导过程。

Java

  1.   
  2. //    
  3. // 阅读的需要,我在下面的宏定中添加了行符    
  4. //    
  5. #define step(shift)    
  6. if((0x40000000l >> shift) + sqrtVal <= val)    
  7. {    
  8. val -= (0x40000000l >> shift) + sqrtVal;    
  9. sqrtVal = (sqrtVal >> 1) | (0x40000000l >> shift);    
  10. }    
  11. else    
  12. {    
  13. sqrtVal = sqrtVal >> 1;    
  14. }    
  15. //    
  16. // 32位整数的平方根    
  17. //    
  18. int32 xxgluSqrtFx(int32 val)    
  19. {    
  20. // Note: This fast square root function    
  21. // only works with an even Q_FACTOR    
  22. int32 sqrtVal = 0;    
  23. step(0);    
  24. step(2);    
  25. step(4);    
  26. step(6);    
  27. step(8);    
  28. step(10);    
  29. step(12);    
  30. step(14);    
  31. step(16);    
  32. step(18);    
  33. step(20);    
  34. step(22);    
  35. step(24);    
  36. step(26);    
  37. step(28);    
  38. step(30);    
  39. if(sqrtVal < val)    
  40. {    
  41. ++sqrtVal;    
  42. }    
  43. sqrtVal <<= (Q_FACTOR)/2;    
  44. return(sqrtVal);    
  45. }   

 

//

// 阅读的需要,我在下面的宏定中添加了行符

//

#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);

}




sqrt话题早在2003年便已在 GameDev.net上得到了广泛的讨论(可在非常火星了,当然不排除有其他尚在冥王星的人,嘿嘿)。而尝试探究该话题则完全是出于本人的趣和好奇心(话说就是无知)。其实现在随着FPU的提升和向量运算的硬件支持,大部分系上都提供了快速的sqrt实现。如果是理大批量的向量的,据最快的方法是使用SIMD(据而已,我根不懂),可同步计4个向量。


--------------------------------------------------------------------------------


克算法牛X的地方就是他了一个常数作起始。而个起始值让他只用一次迭代就了。

里看来的。QuakeIII自然就是奇高手卡克的杰作之一了。在有的CPU上,个函数比普通的(float)(1.0/sqrt(x)4倍!快的原因之一是用了一个神秘常数,0x5f3759df。普渡大学的Chris Lomont文里讨论个常数的意尝试格的方法推个常数(他提到有人认为这个函数是在NVidia工作Gary Tarolli写的)。Chris推出的常数是0x5f37642f),和Q_rsqrt里的稍有不同,而且实际也稍有不如。卡克到底怎推出个常数的就是了。这种高手不写在可惜。

Java

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

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值