125

John Carmack 在 Quake III 源代码中有一个特殊的函数，它计算一个浮点数的反平方根，比常规快 4 倍(float)(1.0/sqrt(x))，包括一个奇怪的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;
}
4

5 回答 5

80

供参考。卡马克没有写。Terje Mathisen 和 Gary Tarolli 都对它进行了部分（并且非常谦虚）的归功于它，以及归功于其他一些来源。

神话常数是如何得出的，这是一个谜。

引用加里·塔罗利的话：

它实际上是以整数进行浮点计算 - 花了很长时间才弄清楚它是如何工作的以及为什么工作，我已经不记得细节了。

由专家数学家 (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 如此出色。

于 2009-08-28T21:52:23.207 回答
56

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

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

于 2009-08-28T22:01:40.863 回答
39

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 会去但更早停止，因此计算速度的提高。

于 2017-02-13T09:51:42.473 回答
25

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

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

于 2018-01-19T23:05:17.970 回答
23

根据不久前写的这篇不错的文章...

代码的魔力，即使你无法理解，也可以从 i = 0x5f3759df - (i>>1); 中脱颖而出。线。简而言之，Newton-Raphson 是一种近似，它从猜测开始，然后通过迭代对其进行细化。利用 32 位 x86 处理器的特性，i 是一个整数，最初使用整数转换设置为要取其平方反比的浮点数的值。然后将 i 设置为 0x5f3759df，减去自身向右移动一位。右移会降低 i 的最低有效位，基本上将其减半。

这真是一本好书。这只是其中的一小部分。

于 2009-08-28T21:57:52.667 回答