速求平方根倒数

Posted 2020-09-05 望山海

 

 在游戏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/y² -x。x为我们为程序输入的数，y为我们所要求的结果。由牛顿迭代法可知: y₁ = y₀(1.5 - 0.5xy₀²) ，这也就是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)的零点在哪，我们可以先猜一个x₀ ，然后求得f’(x₀)，即求得一条在x₀点与函数f(x)相切的直线，直线与x轴的交点即为一个零点的近似值，我们称之为x₁。x₁一定比x₀更接近零点（如果x₀很显然，x₁要比x₀更接近零点，在数次迭代后，我们就可以得到误差足够小的近似解。易得公式x₁ = x₀ - f(x₀）/f’(x₀)其实是在为这一步迭代给予一个 initial guess。那么我们现在就可以集中精力去算算，这个0x5f3759df到底是从哪里蹦出来的了。

 

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

 

  现在我们设输入的浮点型数的各位数的值分别为 0 E M。Magic constant部分的数的各部分值位 0 R₁ R₂。

 

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

 

E为奇数：

 

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

 

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

 

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

 

y = ( 1 + ( 1 + R₂ - [M/2] ) )2^{R1 - 64 - d - 127 - 1}。

 

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

 

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

 

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

 

相对误差 ： H = ( Y - y)/Y  = 1 - (2 + A )2^{R1 - 192}。

 

（A为两种情况下位数的值）

 

若使误差最小，即H最小，则有R₁ 必属于[189.2,190.7]，也就是说，作为一个整数，R₁只有在取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为偶数的时候算法是相似的，不再赘述，现在我们最优化尾数部分。由于R₁已经确定了，所以原函数成为了M 与R的方程，即f(M,R)，在Lomont老哥的高超算法技巧以及数学工具的帮助下，解得 R₂ = 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

 

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