侧耳聆听

译者：zhfkt

zhfkt@hotmail.com

原创译文欢迎转载，转载时请注明译者信息和文章原始出处，如有错误尽请指正。

原文：FAST INVERSE SQUARE ROOT

WIKI上更详尽的解释http://en.wikipedia.org/wiki/Fast_inverse_square_root

 

 

如何快速求解平方根的倒数 

CHRIS LOMONT

 

摘要：在许多应用程序里，计算逆平方根都是必不可缺的：比如在3D游戏中计算单位向量的时候。通常，由于计算速率的大幅度提高，一些精度的损失是可以接受的。这篇文章详细说明了在一些在线的代码库里所发现的快速求解平方根的倒数源代码，并且提供了思路，让读者可以从类似的方法推出其它函数快速计算代码。

 

1.引言

当我在www.gamedev.net的“数学编程”论坛的时候，遇到了一个快速求解平方根倒数的有趣算法。这个代码本质上用C语言描述就是：

 

 

float InvSqrt(float x)

{

float xhalf = 0.5f*x;

 

int i = *(int*)&x; // 取得浮点数的每一个bit

i = 0x5f3759df - (i>>1); // 猜测一个初始的y0

x = *(float*)&i; // 把每一个bit重新转换成浮点数

x = x*(1.5f-xhalf*x*x); // 运用牛顿迭代法反复计算提高精度

return x;

}

 

 

     最有意思的地方就是常数0x5f3759df：这个常数从哪里来？为什么这些代码会起作用？在Visual C++.NET中的某些快速测试中表明，上述的代码大约比原生的(float)(1.0/sqrt(x))快4倍，并且遍历所有浮点数后最大的误差只有0.00175228，因此这种方法看起来非常有效。而紧接着的三个问题就是：

 

（1）       这个算法是如何工作的？

（2）       这个算法能够被改进吗？

（3）       是哪一位“位运算牛人”(bit master)设计出了如此不可思议的算法？

 

 

     在快速搜索了0x5f3759df之后Google给予了一些提示，关于这段代码最相关的参考是2002年1月9号在comp.graphics.algorithms上的一个帖子，它是由D. Eberly所得出的一个错误但已经非常接近的解释，不过他的阐释很有启发性。而更进一步的挖掘也没有发现关于这段代码的其它正确解释。这段代码最初是出现在Quake 3的源代码里，由传奇的游戏编程大师John Carmack所编写，但是在网上有人把它归于曾经在Nvidia的Gary Tarolli（我没有找到出处）。有谁能够确认作者的身份吗？这段代码也在Crystal Space的源代码 the Titan Engine和the Fast Code Library中出现过，尽管看起来每一段代码都似乎来源于Quake 3.

 

     尝试这个算法的动机在Eberly的文档中有更加清楚的解释，他假定这种变换的目的是对于求逆平方根的线性插值。注意有多种方法来加快这一代码，但本文不打算进一步优化。当然也还有更快的方法，或许通过查表的技巧可以达到，但其中大多数都没有这种方法来的精确。

 

     本文假定PC是32位或是相似的架构。特别是浮点型的表示方法在IEEE 754-1985中有叙述。这段C源码已经被当成是字节序中立的（据称是在Macintosh上进行的测试）。我并没有去验证这一点，因为这个方法在32位情况下似乎出人意料的和字节序无关，所以很容易运用本文所提供的想法把这段代码扩展到别的情况和位数长度下（比如double型）。无论如何，让我们现在着手开始解决问题1.2.和3吧。

 

2.背景知识

 

     浮点数是作为32位比特(bits)被存储在PC中的：

 

 

     S是符号位（1代表的是负数），E是8位阶数(exponent)，M是23位尾数(mantissa)。阶数位可以减去127而得到正指数或是负指数，并且尾数M不储存大于1的值。因此考虑到M作为一个完全在小数点在右边的二进制数，它的值域为I = [0, 1)。而x所代表的值为：

     这些比特(bits)能够被看成真实数值(real number)的浮点表示，而单单考虑比特自身也能够被看作整数(integer)。所以M可以通过前后关系被当做一个小数或是整数，并且M被当作真实数值的时候可以通过M作为整数时除以2^23求得。

 

3.算法

     算法的主要思想是通过牛顿迭代法得到近似，那个神奇的常量被用作初始迭代时的一个精确猜测。考虑到浮点数X>0，我们希望计算的是。定义，我们想要的值就是f(x)的正根。牛顿求根的方法给出一个合适的近似根，并通过迭代找到更加近似的结果：

 

     由于f(y)已知，这个式子可以化简成，所对应的代码就是x = x*(1.5f-xhalf*x*x)。x是一个初始的猜测，之后将被称作y0。

 

     代码i = 0x5f3759df - (i>>1)计算了这个初始的y0，大致是x的-1/2次方，接着通过筛选比特（picking bits）将误差减少至最低。i>>1是C语言中的右移运算符，会将整数的位数向右移一位。这种方法通过去掉最低有效位，高效的将原数折半除2。Eberly的解释是这是一种线性逼近，然而他的观点却是错误的。我们将看到这是一种分段的线性(piecewise linear)，该函数的近似并非在所有情况下都一样。不过，我猜他的解释是原作者创建该算法的灵感。

 

 

4.神奇的0x5f3759df

       因此现在我们只剩下寻找一个初始的猜测y0的过程了。假设我们有一个浮点数X>0，所对应的32位就是

     指数e=E−127，接着令，而我们想要的数，在之中要把e和M当作真实所能表示的值（real values）而非整数(integers)。对于一般情况下我们取得了神奇的常数R，并且分析y0 = R−(i>>1)。这里的减法是对于整数，i是以比特(bits)表示的X，然而我们将y0视为真实所能表示的值。R作为浮点数的整形表示方法可以看作

 

     即R1是在阶数域而R2是在尾数域。当我们将i向右移一位的时候，我们可能将阶数域比特移入尾数域比特，因此我们要分析两种情况。

 

     在本文剩下的部分对于实数a来说，[a]代表不大于a的最大整数。 

 

4.1 当阶数域E为偶数

 

       在这种情况下，当我们使用i>>1时，没有比特(bits)从阶数域移入尾数域，并且[E/2]= E/2。而对于真实的指数e=E-127是一个奇数，因此令e=2d+1，d是一个整数。这些比特所代表的初始猜测给予了一个新的阶数域。

       我们要求这个式子的最后结果为正，如果是负的符号位将会为1，这个方法将不能成功的返回一个正数。如果最终结果为0，那么尾数域将不会从阶数域借位。我们在下面将会看到这个条件非常必要的。因为这对于任何的偶数来说，我们都需要R1>=128。

 

       由于在i>>1的情况下没有比特(bit)从阶数域移入尾数域，因此新的尾数部分作为整数来说就是R2-[M/2]。假设R2>=[M/2]，那么初始的猜测y0就是

 

 

       可以看到M/2替换了[M/2] ，因为作为实数部分尾数最终结果的误差最多是，相比于这种方法中的其他误差来说是微不足道的。

 

       如果R2>1)作为整数相减时，在尾数域中的比特（bit）将不得不从阶数域中借位（这就是为什么我们需要新的阶数域为正）。尾数域中的比特就会变成(1+R2)-[M/2]，因此如果R2

 

 

      在这里我们把y0的阶数域写的与R2>=[M/2]情况时一样

 

   如果我们定义

 

 

      我们就能把这两个与y0相关的等式用一个统一的式子来表示。

      注意到除了在R=M/2的时候，都是连续且可导的。把e替换成e=2d+1，近似后的值就是 

 

 

      相对误差是我们用来衡量这种方法好坏的标准

 

 

      化简成

 

 

      注意到现在这个式子不再与d的取值有任何关系，并且M可以取为I=[0,1)之间的任何值。这个公式仅当E为偶数时才有效，因此也和E有隐藏的依赖关系。

      假设我们需要找到一个仅使最大相对误差<1/8的R2，那么

   扩展出来后

 

      最后一步用到了这个事实，将左右两边同时取以2为底的log后加上190.5就得到190.7>R1>189.2，因此R1=190=0xbe对于偶数E来说是一个保证任何结果的相对误差小于0.125的独一无二的整数。在以浮点数表示中的第24位至第30位，这将使(0xbe>>1)=0x5f成为神奇常数R的一部分（当然还需要第23位为0）。

 

      为了说明这一点，图1显示了在R2=0.43 d=2的情况下作为实线画出的y0和作为虚线函数 所近似的结果。很显然，y0是一个非线性逼近，然而，用非线性的方法却更加合适。图2说明了当与时相对误差的值。显然常数R2在整个定义域中最小的极大误差大约为R2=0.45。

 

4.2 当阶数域E为奇数时

      在已经有了分析前一种情况的经验下，我们来分析一种更加困难的情况。不同之处就在于阶数域E中的最低位比特(bit)会由于i>>1而移入尾数域中的最高位，这样移入后尾数域真正的值就会比原先多加0.5，即[M/2]+1/2。这个式子是对尾数本身所代表的整数进行取整运算，而对尾数域所代表的小数进行加法运算。e=E-127=2d是偶数。和上面的结论相似新的阶数域为

   像之前的情况一样这里必须是正数。

 

      由于同样地尾数最终结果的误差最多是，可以把[M/2]替换成M/2。当R2>=(M+1)/2时，新的尾数域所代表的小数为R2-(M/2+1/2)。这要求尾数域不能从阶数域借位。

      我们为y0选择了与E是偶数时相同的阶数域。如果R2<(M+1)/2。减法运算R−(i<<1)需要从正阶数域借一个比特，而使y0/2。而尾数域所代表的小数也会变成1+R2-(M/2+1/2),且依然在定义域I内。因此如果R2<(M+1)/2，我们得到的初始猜测y0为：

 

      与 相似，定义

      我们可以把这两个y0立即写成

      注意到除了在2R2=M+1之外都是连续且可导的，因此对于y0来说也是如此。令e=2d，我们把y0近似表示成

      需要注意的是这不是和我们在E为偶数时所得到的相同的函数，这种情况少了一个的因子

      我们想得到的最小相对误差像上面一样化简成

 

      又一次与d的取值无关（不过像上面一样对E有着隐藏的依赖关系）

      和之前一样如果对于奇数E我们想让近似结果的最大相对误差小于1/8，我们可以对M=0或是M=1的情况进行分析，但是因为我们要找的常数R1既要适合奇数情况下的E，又要适合偶数情况下的E，我们可以取由上一种情况分析得到对两种情况都可行的R1=190。注意到这个R1的值满足本篇开头时所要求的R1>=128。因此对于本文余下的部分假设R1=190=0xbe，并且重新定义两个误差函数。

      这两个误差函数仅仅和M与R2有关，且定义域皆为I。

 

      到目前为止我们已经完全理解了算法的原理，并且计算机测试证实了模型中近似步骤的准确性，所以，我们已经达到了目标1，即“这个算法是如何工作的？”。

 

4.2.1 例子

 

   当X为16时，E=4+127=131，M=0，阶数域当中的一个bit被移入了尾数域当中。因此在i = 0x5f3759df - (i>>1)这一行之后，初始的猜测 y0= 0x3E7759DF，新的阶数域就变成190 − [131/2] = 125 = 0x7d，对应于e=125-127=-2。新的尾数域需要从阶数域借位一个bit，使e=-3，并且得到 0xb759df - 0x4000000= 0x7759DF，对应于M= 0.932430148。因此，这已经与=0.25的真实结果相当的接近了。

 

 未完待续