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 (Chris Lomont) 试图弄清楚原始算法是如何工作的：

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 的 sqrt 的数学“高级”版本（达到几乎相同的常数）证明不如 Gary 最初开发的版本。他无法解释为什么 id 的 iirc 如此出色。

Answer 2

当然，这些天来，它比仅使用 FPU 的 sqrt 慢得多（尤其是在 360/PS3 上），因为在 float 和 int 寄存器之间交换会导致 load-hit-store，而浮点单元可以在硬件中做倒数平方根。

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

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 会去而是提前停止，从而提高计算速度。

Answer 4

我很想知道作为浮点数的常量是什么，所以我简单地编写了这段代码并在谷歌上搜索了弹出的整数。

long i = 0x5F3759DF;
float* fp = (float*)&i;
printf("(2^127)^(1/2) = %f\n", *fp);
//Output
//(2^127)^(1/2) = 13211836172961054720.000000

看起来这个常数是“2^127 平方根的整数近似值，其浮点表示的十六进制形式更广为人知，0x5f3759df”https://mrob.com/pub/math/numbers-18.html

在同一个网站上，它解释了整个事情。 https://mrob.com/pub/math/numbers-16.html#le009_16

Answer 5

根据to this nice article不久前写的...

代码的魔力，即使你 跟不上，突出为 i = 0x5f3759df - (i>>1);线。简化， Newton-Raphson 是一个近似值 从猜测开始 通过迭代对其进行细化。服用 32 位 x86 特性的优势 处理器，i，一个整数，是 最初设置为 你要取的浮点数 的平方反比，使用 整数转换。然后我设置为 0x5f3759df，减去自身移一 向右一点。正确的转变 删除 i 的最低有效位， 基本上减半。

真的很好读。这只是其中的一小部分。