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;
}
顺便提一下，Carmack并没有编写它。Terje Mathisen和Gary Tarolli都在部分（非常谦虚的）认领了它，并给出了其他一些来源的贡献。
这个神秘的常数是如何得出的还有点不明确。
引用Gary Tarolli的话：
实际上是在整数中执行浮点计算——花费了很长时间才弄清楚为什么以及如何工作，我已经记不起细节了。
一个稍微更好的常数，由一位专业数学家（Chris Lomont）开发，他试图弄清楚原始算法的工作方式，它的链接是：http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf
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;
}
当然，在当今这些日子里，它的速度比仅使用FPU的sqrt要慢得多（特别是在360/PS3上），因为在浮点寄存器和整数寄存器之间切换会导致load-hit-store，而浮点单位可以通过硬件执行倒数平方根。
这只是展示了随着底层硬件性质的变化，优化必须不断地发展进化。
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)
我很好奇将常量转换为浮点数后得到的结果是多少，于是我写了这段代码，并在谷歌上搜索弹出的整数。
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.htmlhttps://mrob.com/pub/math/numbers-16.html#le009_16sqrt(2^127)，而是一个接近的近似值（最高有效位数为两位）。sqrt(2^127) = 1.3043x10^19 而 0x5F3759DF = 1.3211x10^19。 - Loves Probabilityi = * ( long * ) &y;
i = 0x5f3759df - ( i >> 1 );
y = * ( float * ) &i;
第一行将浮点数y的表示形式作为整数i处理。第二行将i向左移动一位并从一个神秘的常量中减去它。第三行将得到的数字转换回标准的float32格式。那么这是为什么呢？
假设g是一个将浮点数映射到其浮点表示形式（读作整数）的函数。上面的第一行设置i = g(y)。
以下好的近似g存在(*)：g(y) ≈ Clog_2 y + D，其中C和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是一个在相当广泛的范围内产生良好结果的常数之一。
有关此过程的更详细描述，请参见维基百科：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 + σ，其中 sigma 是一个小常数。另外，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。正的零点是 y = 1/sqrt(num)，这是我们想要计算的。y_n+1 = y_n - f(y_n)/f'(y_n)。y_n+1 = y_n - (-y_n+y_n^3*num)/2 = y_n * (3/2 - num/2 * y_n * y_n)。上面代码行所做的就是这个。
i = * ( long * ) &y;中，为什么要将 y 的地址作为长整型指针来获取，然后再将其解除引用？ - Nubcakey是一个float，这里将其强制转换为整数。这是不安全的，因为它违反了C语言的严格别名规则。在C99中使用union，或在C89 / C++中使用memcpy可以遵循语言规则做同样的事情，并且至少在现代优化编译器中编译结果相同。 - Peter Cordes