约翰 · 卡马克的非同寻常的平方根倒数速算法(《地震3》)

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;
}
72284 次浏览
小开

仅供参考。不是卡马克写的。泰耶•马蒂森和加里•塔罗利(Gary Tarolli)都为此获得了部分(也是非常有限的)功劳，同时还获得了其他一些来源的功劳。

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

引用 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 是如此优秀的 IRC。

小开

根据 这篇好文章前阵子写的..。

代码的魔力，即使你 跟不上它，因为 i = 显得很突出 0x5f3759df-(i > > 1) ; 行, 牛顿-拉斐逊是一个近似值 以猜测开始 通过迭代来提炼它 32位 x86性质的优势 处理器，i，一个整数，是 初始设置为 要取的浮点数 反平方，使用 然后将 i 设置为 0x5f3759df 减去本身的位移1 往右一点，右转 掉了最不重要的一点, 实际上是减半。

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

小开

当然，现在它比仅仅使用 FPU 的 sqrt (特别是在360/PS3上)要慢得多，因为浮点寄存器和 int 寄存器之间的交换会导致负载命中存储，而浮点寄存器可以在硬件中执行倒数平方根。

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

小开

Greg Hewgill 和 伊利丹 S4给出了一个极好的数学解释。 我在这里总结一下，对于那些不想过多谈论细节的人来说。

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

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);是不一样的，因为把浮点数作为整数不仅要除以2，而且还要除以指数，还会产生一些其他的伪影，但是它还是要计算一些系数 a0，a1，a2..。

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

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

他们可以在一个循环中进行更多的迭代，每一次迭代都会改进结果，直到达到所需的准确性为止。但似乎只有一次迭代就足够了，这也是速度的福音。CPU/FPU 根据需要进行多次迭代，以达到存储结果的浮点数的准确性，它具有更通用的算法，适用于所有情况。

简而言之，他们所做的是:

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

小开

我很想知道浮点数的常量是什么所以我写了这段代码然后在谷歌上搜索出来的整数。

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”

在同一个网站上，它解释了整个事情

小开

代码由两个主要部分组成。第一部分计算1/sqrt (y)的近似值，第二部分采用牛顿方法的一次迭代得到一个更好的近似值。

计算1/sqrt (y)的近似值

i  = * ( long * ) &y;
i  = 0x5f3759df - ( i >> 1 );
y  = * ( float * ) &i;

第1行获取 y 的浮点表示，并将其视为一个整数 i。第2行将 i 移过一位，并从一个神秘的常量中减去它。第3行获取结果数字并将其转换回标准 float32。为什么这个有用？

设 g 是一个函数，它将一个浮点数映射到它的浮点表示，读取为一个整数。上面的第一行设置的是 i = g(y)。

下面是 g 存在的最佳近似(*) : 一些常量 C 和 D 的 g(y) ≈ Clog_2 y + D。关于为什么存在如此好的近似的一个直觉是 y 的浮点表示在指数中大致是线性的。

第2行的目的是将 g (y)映射到 g (1/sqrt (y)) ，然后第3行可以使用 g ^-1将该数字映射到1/sqrt (y)。使用上面的近似值，我们得到了 g(1/sqrt(y)) ≈ Clog_2 (1/sqrt(y)) + D = -C/2 log_2 y + D。我们可以使用这些公式来计算从 g(y)到 g(1/sqrt(y))的映射，即 g(1/sqrt(y)) ≈ 3D/2 - 1/2 * g(y)。在第2行，我们有 0x5f3759df ≈ 3D/2和 i >> 1 ≈ 1/2*g(y)。

常数0x5f3759df 略小于给出 g(1/sqrt(y))最佳近似值的常数。这是因为这一步不是单独完成的。由于牛顿的方法容易错过的方向，使用稍微小一点的常数往往会产生更好的结果。在此设置中使用的确切最佳常数取决于 y 的输入分布，但0x5f3759df 就是这样一个常数，它可以在相当宽的范围内提供良好的结果。

这个过程的更详细的描述可以在 Wikipedia: https://en.wikipedia.org/wiki/Fast_inverse_square_root#Algorithm上找到

(*)更明确地说，让 y = 2^e*(1+f)。取两边的对数，我们得到 log_2 y = e + log_2(1+f)，对于一个小常数 σ，它可以近似为 log_2 y ≈ e + f + σ。另外，表示为整数的 y 的 Float32编码是 g(y) ≈ 2^23 * (e+127) + f * 2^23。结合这两个方程，我们得到了 g(y) ≈ 2^23 * log_2 y + 2^23 * (127 - σ)。

用牛顿的方法

y  = y * ( threehalfs - ( x2 * y * y ) );

考虑函数 f(y) = 1/y^2 - num，f 的正零是 y = 1/sqrt(num)，这正是我们感兴趣的计算。

牛顿方法是一种迭代算法，用于对函数 f 的零点取近似 y _ n，并用以下方程计算更好的近似 y _ n + 1: y_n+1 = y_n - f(y_n)/f'(y_n)。

计算函数 f 的形状，得到如下方程: y_n+1 = y_n - (-y_n+y_n^3*num)/2 = y_n * (3/2 - num/2 * y_n * y_n)。这正是上面这行代码正在做的事情。

你可以在这里了解更多关于牛顿方法的细节: https://en.wikipedia.org/wiki/Newton%27s_method


提交答案