快速平方根算法.docxVIP

快速平方根算法(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