首页 > 代码库 > 速求平方根倒数
速求平方根倒数
在游戏3D建模方面很多时候要用到求平方根的倒数，而本文章打算介绍的算法会比正常算法快上4倍左右。这对于产品性能将是一个大幅度的提高。
那我们要从哪里开始呢？首先不得不提一提 idsoftware。这是一个创建之初只有13个人的小公司，但它推出的毁灭战士（DOOM）系列游戏可以说改变了游戏世界，极大地推动了游戏产业的发展，因为在当时贫瘠的电脑性能的支撑下，一个开发者能够在游戏中加入一段流畅的动画都会让人惊叹不已，而我们所说的idsoftware，在那个年代就已经做出了画面远超同代其余作品的游戏，像之前提到的DOOM以及Quake都曾是煊赫一时的产品，这个仅有13人的小公司，在1995年一年的纯利润收入就达到了1500多万美元，堪称业界神话，那么这和我们要讲的快速算法有什么关系呢？
我们要知道idsoftware这家公司不止制作游戏，他们还在做引擎，像大名鼎鼎的半条命（Half-Life）以及反恐精英（CS）就是用这家公司的Quake引擎制作完成的。依据GPL协议，idsoftware公司在2005年08月22日正式对外公开了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 -x。x为我们为程序输入的数，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轴的交点即为一个零点的近似值，我们称之为x1。x1一定比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,，e和m部分的值变为原来的一半。
现在我们设输入的浮点型数的各位数的值分别为 0 E M。Magic 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，老哥怎能善罢甘休，使用了暴力穷举法，找到了最好的constant即0x5f375a86。文末老哥还表示对这段代码的创始人很感兴趣，想问问他，原始数字到底是算出来的还是一个个测试出来的。
参考文献：
FAST INVERSE SQUARE ROOT CHRIS LOMONT
http://www.matrix67.com/data/InvSqrt.pdf
---恢复内容结束---
速求平方根倒数