首页 > 代码库 > 速求平方根倒数

速求平方根倒数

 

 在游戏3D建模方面很多时候要用到求平方根的倒数,而本文章打算介绍的算法会比正常算法快上4倍左右。这对于产品性能将是一个大幅度的提高。

  那我们要从哪里开始呢?首先不得不提一提 idsoftware。这是一个创建之初只有13个人的小公司,但它推出的毁灭战士(DOOM)系列游戏可以说改变了游戏世界,极大地推动了游戏产业的发展,因为在当时贫瘠的电脑性能的支撑下,一个开发者能够在游戏中加入一段流畅的动画都会让人惊叹不已,而我们所说的idsoftware,在那个年代就已经做出了画面远超同代其余作品的游戏,像之前提到的DOOM以及Quake都曾是煊赫一时的产品,这个仅有13人的小公司,在1995年一年的纯利润收入就达到了1500多万美元,堪称业界神话,那么这和我们要讲的快速算法有什么关系呢?

  我们要知道idsoftware这家公司不止制作游戏,他们还在做引擎,像大名鼎鼎的半条命(Half-Life)以及反恐精英(CS)就是用这家公司的Quake引擎制作完成的。依据GPL协议,idsoftware公司在20050822日正式对外公开了Quake3引擎的源代码。引得业界群众普大喜奔,大家早就想看看这些能够压榨电脑性能到极限的代码是什么样子。

  很快人们发现了下面这段代码,也就是我们今天的主角,速求平方根倒数算法:

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#ifndef Q3_VM
#ifdef __linux__
assert( !isnan(y) ); // bk010122 - FPE?
#endif
#endif
return y;
}  

  了解了什么是牛顿迭代法之后,我们再回过头来看我们的问题,由我们的问题可以抽象出一个函数f(x) = 1/y2 -xx为我们为程序输入的数,y为我们所要求的结果。由牛顿迭代法可知: y1 = y0(1.5 - 0.5xy02) ,这也就是x = x*(1.5f-xhalf*x*x); 这一步的数学表达式。也就是说我们所看不懂的“magic code到,这个算法的主要实现思想就是牛顿迭代法,但是这种方法可以在两次迭代以内就能得到精度极高的结果,并且第一次迭代所用方法是基于位数直接操作的,对于当时很多没有加载乘法功能的cpu来说,极大地提高了运算的速度,这也是这种算法能够比正常算法快上四倍的原因所在。  而这个算法最神奇也最令人费解的地方莫过于i = 0x5f3759df - ( i >> 1 );这一句了。它甚至被标以what the fuck?这种注释,足以见得其本身之神奇。有一天,数学家Chris Lomont意外的发现了这段代码,他打算亲自下场,算一算这个语句是怎么来的,他使用了各种高端方法来进行线性拟合,也求出了一个常数0x5f37642f。然而实际使用时发现,这个常数虽然在线性拟合上的匹配度很高,但是在一次牛顿迭代后,精度反不及原常数0x5f3759df。怎么办呢,Lomont老哥一计不成又生一计,我暴力破解啊,于是Lomont老哥在几经尝试之后,找到了一个最优常数0x5f375a86。然后他据此写了一篇论文,本文后半部分就将结合Lomont老哥的论文,解释一下这个神奇的常数是怎么来的。  首先给大家放上Lonmont老哥的精简代码版本:float InvSqrt(float x) {   float xhalf = 0.5f*x;  

   int i = *(int*)&x; // get bits for floating value   

   i = 0x5f3759df - (i>>1); // gives initial guess y0   

  x = *(float*)&i; // convert bits back to float   

  x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy return x;

}

   Lomont老哥这个版本基本上就是主要步骤的集合,也就是研究的目标。那我们从哪里开始呢?从牛顿迭代法开始。什么是牛顿迭代法,比如说你想求一个函数f(x)的零点在哪,我们可以先猜一个x0 ,然后求得f’(x0),即求得一条在x0点与函数f(x)相切的直线,直线与x轴的交点即为一个零点的近似值,我们称之为x1x1一定比x0更接近零点(如果x0很显然,x1要比x0更接近零点,在数次迭代后,我们就可以得到误差足够小的近似解。易得公式x1 = x0 - f(x0/f’(x0)其实是在为这一步迭代给予一个 initial guess。那么我们现在就可以集中精力去算算,这个0x5f3759df到底是从哪里蹦出来的了。

 

  作为基础知识储备,先给大家说一下浮点型的存储(32位)。浮点型的存储分为3个部分,表示数字符号的s部分(31号位),表示数字指数大小的e部分(30-23号位)以及表示数字位数大小的m部分(22-0号位)。即一个浮点型可以表示为 x = (-1)s(1+m)2e-127。我们可以看到,在进行magic code部分之前,我们将这个数强制转换为了int类型,如此一来 i>>1 (位运算符大家还没忘吧︿( ̄︶ ̄)︿)就意味着对整个数除以2,em部分的值变为原来的一半。

 

  现在我们设输入的浮点型数的各位数的值分别为 0 E MMagic constant部分的数的各部分值位 0 R1 R2

 

  在我们对输入的浮点型数除以2时,会发生两种情况,第一种情况,E是偶数,第二种情况E是奇数,区别在于会不会有一位属于E的位数进入到M中。为什么我们不讨论M是不是奇数呢,因为哪怕在位移中丢了一个1,那也将只是2-24级别的误差,对于整体的精度影响太小,故不以考虑。准备工作做好了,现在我们就开始对E进行分情况讨论。

 

E为奇数:

 

  此时我们设实际意义上的指数e = 2d + 1,则有E= 128 + 2d,那么

 

R1 - [e/2] = R1 - 64 - d (1)

 

[ a]表示不超过a的最大整数)因为负数对我们没有意义,所以我们要保证(1)式结果为正,又因为E属于[0,254],则我们要求R1 >= 128。然而对于R2我们还有两种情况要予以讨论,第一种,R2大于[M/2],此时没什么变化,正常做减法,y = (1 + R2 - M/2)2R1 - 64 - d - 127。第二种,R2小于[M/2],此时R2要向前一位借1。所以结果为

 

y = ( 1 + ( 1 + R2 - [M/2] ) )2R1 - 64 - d - 127 - 1

 

  因为函数处于零点时会有 1 / y2 = x 。所以由上述方程可得

 

Y  = (1 / (1 + M)0.5) 2-e/2,

 

若想使我们的magic code 的结果y与实际结果Y最接近则有

 

相对误差 H = ( Y - y)/Y  = 1 - (2 + A )2R1 - 192

 

A为两种情况下位数的值)

 

若使误差最小,即H最小,则有R1 必属于[189.2,190.7],也就是说,作为一个整数,R1只有在取190 的时候,误差最小。190则可以用0xbe来表示。而Lomont老哥在论文里原话是这样的:In bit positions 24 through 30,this gives leading (0xbe >>1) = 0x5f part of the magic constant R( as well as requiring bit 23 to be 0)。也就是说我们有了“magic constant”的前一部分。

 

  E为偶数的时候算法是相似的,不再赘述,现在我们最优化尾数部分。由于R1已经确定了,所以原函数成为了M R的方程,即f(M,R),在Lomont老哥的高超算法技巧以及数学工具的帮助下,解得 R2 = 0x37642f。(详情请见文末原论文链接)。与前文连起来即可得0x5f37642f,属于Lomont老哥自己的 magic constant。很明显这跟源代码所提供的magic costant 有所区别,Lomont 老哥就对此进行了测试,测试之后发现,虽然自己的costant更符合原理,但在第一次迭代之后效果反不如原始constant,老哥怎能善罢甘休,使用了暴力穷举法,找到了最好的constant0x5f375a86。文末老哥还表示对这段代码的创始人很感兴趣,想问问他,原始数字到底是算出来的还是一个个测试出来的。

 

 

 

 

 

参考文献:

 

                                                                                                                     FAST INVERSE SQUARE ROOT      CHRIS LOMONT

 

                                                                                                                          http://www.matrix67.com/data/InvSqrt.pdf

 

---恢复内容结束---

速求平方根倒数