John Carmack的不同寻常的快速反平方根(Quake III)

Question

John Carmack在Quake III源代码中有一个特殊的函数，它可以计算浮点数的平方根的倒数，比常规的(float)(1.0/sqrt(x))快4倍，其中包括一个奇怪的0x5f3759df常量。请参阅下面的代码。有人能逐行解释一下这里到底是怎么回事吗?为什么它比常规实现快这么多？

float Q_rsqrt( float number )
{
  long i;
  float x2, y;
  const float threehalfs = 1.5F;

  x2 = number * 0.5F;
  y  = number;
  i  = * ( long * ) &y;
  i  = 0x5f3759df - ( i >> 1 );
  y  = * ( float * ) &i;
  y  = y * ( threehalfs - ( x2 * y * y ) );

  #ifndef Q3_VM
  #ifdef __linux__
    assert( !isnan(y) );
  #endif
  #endif
  return y;
}

Answer 1

仅供参考。不是卡马克写的。Terje Mathisen和Gary Tarolli都将其部分(非常适度)归功于此，以及归功于其他一些来源。

这个虚构的常数是如何推导出来的，这是一个谜。

引用加里·塔罗利的话：

实际上是在做整数形式的浮点计算-它花了很长时间才弄清楚它是如何工作的，为什么它是这样的，我已经记不住细节了。

一个稍微好一点的常量，developed by an expert mathematician (克里斯·洛蒙特饰)试图弄清楚原始算法是如何工作的：

float InvSqrt(float x)
{
    float xhalf = 0.5f * x;
    int i = *(int*)&x;              // get bits for floating value
    i = 0x5f375a86 - (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;
}

尽管如此，他最初尝试的id's sqrt‘s sqrt的数学上的’高级‘版本(几乎是相同的常数)被证明比Gary最初开发的版本要差得多，尽管它在数学上要’纯粹‘得多。他无法解释为什么id的iirc如此优秀。

Answer 2

当然，如今，它比仅仅使用FPU的sqrt慢得多(尤其是在360/PS3上)，因为在浮点和整型寄存器之间交换会导致加载-命中-存储，而浮点单元可以在硬件中进行倒数平方根。

它只是展示了随着底层硬件的性质变化，优化必须如何发展。

Answer 3

Greg Hewgill和IllidanS4给出了一个极好的数学解释。对于那些不想深入细节的人，我将试着在这里总结一下。

任何数学函数，除了某些例外，都可以用多项式求和来表示：

y = f(x)

可以精确地将转换为：

y = a0 + a1*x + a2*(x^2) + a3*(x^3) + a4*(x^4) + ...

其中a0，a1，a2，...是常量。问题是，对于许多函数，比如平方根，对于精确值，这个和有无穷多个成员，它不会在某个x^n结束。但是，如果我们停在某个x^n，我们仍然会得到一个达到一定精度的结果。

所以，如果我们有：

y = 1/sqrt(x)

在这种特殊情况下，他们决定丢弃秒以上的所有多项式成员，可能是因为计算速度：

y = a0 + a1*x + [...discarded...]

现在任务归结为计算a0和a1，以便y与精确值的差异最小。他们计算出最合适的值是：

a0 = 0x5f375a86
a1 = -0.5

所以，当你把这个放入等式中时，你会得到：

y = 0x5f375a86 - 0.5*x

这与您在代码中看到的代码行相同：

i = 0x5f375a86 - (i >> 1);

编辑:实际上在这里y = 0x5f375a86 - 0.5*x和i = 0x5f375a86 - (i >> 1);是不一样的，因为将浮点数作为整数进行移位不仅会被二除，还会将指数除以二，从而导致一些其他的伪像，但它仍然需要计算一些系数a0，a1，a2……。

在这一点上，他们发现这个结果的精度不足以达到目的。因此，他们只做了牛顿迭代的一步来提高结果的准确性：

x = x * (1.5f - xhalf * x * x)

他们可以在循环中进行更多的迭代，每次迭代都会改善结果，直到达到所需的精度。这正是它在CPU/FPU!中的工作方式，但似乎只有一次迭代就足够了，这也是速度的福音。CPU/FPU根据需要进行尽可能多的迭代，以达到存储结果的浮点数的精度，并且它具有更通用的算法，适用于所有情况。

因此，简而言之，他们所做的是：

使用(几乎)与CPU/FPU相同的算法，针对1/sqrt(x)的特殊情况利用初始条件的改进，并且不会一直计算到精度CPU/FPU会提前停止，从而提高计算速度。