快速平方根算法.docVIP

快速平方根算法(1) 默认分类 2010-09-03 06:42:49 阅读16 评论0 ??字号：大中小?订阅 快速平方根算法 在3D图形编程中，经常要求平方根或平方根的倒数，例如：求向量的长度或将向量归一化。C数学函数库中的sqrt具有理想的精度，但对于3D游戏程式来说速度太慢。我们希望能够在保证足够的精度的同时，进一步提高速度。 Carmack在QUAKE3中使用了下面的算法，它第一次在公众场合出现的时候，几乎震住了所有的人。据说该算法其实并不是Carmack发明的，它真正的作者是Nvidia的Gary Tarolli（未经证实）。 // // 计算参数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 Method，简称NR），而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)。现在就十分清晰了。语句i1其工作就是将指数除以2，实现2^(E/2)的部分。而前面用一个常数减去它，目的就是得到M^(1/2)同时反转所有指数的符号。 至于那个0x5f3759df，呃，我只能说，的确是一个超级的Magic Number。 那个Magic Number是可以推导出来的，但我并不打算在这里讨论，因为实在太繁琐了。简单来说，其原理如下：因为IEEE的浮点数中，尾数M省略了最前面的1，所以实际的尾数是1+M。如果你在大学上数学课没有打瞌睡的话，那么当你看到(1+M)^(-1/2)这样的形式时，应该会马上联想的到它的泰勒级数展开，而该展开式的第一项就是常数。下面给出简单的推导过程： 对于实数R0，假设其在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)正是(R1)， 而在IEEE表示中，指数的符号只需简单地加上一个偏移即可， 而式子的前半部分刚好是个常数，所以原式可以转化为： 原式 = C - (M/2)*2^(E/2) = C - (R1)，其中C为常数 所以只需要解方程： R^(-1/2) = (1+M)^(-1/2) * 2^(-E/2) = C - (R1) 求出令到相对误差最小的C值就可以了 上面的推导过程只是我个人的理解，并未得到证实。而Chris Lomont则在他的论文中详细讨论了最后那个方程