math 之 sqrt 0x5f3759df

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

    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;                        // evil floating point bit level hacking
    i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//    y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

#ifndef Q3_VM
#ifdef __linux__
    assert( !isnan(y) ); // bk010122 - FPE?
#endif
#endif
    return y;
}

以上代码来自一款老游戏（quake3-1.32b-source）的源码，实现sqrt函数，比C库里面的sqrt函数快了4倍，是最快的sqrt版本；

0x5f3759df 神奇地完成了核心运算，想知道是如何得来的。

2022/8/25更新

在网上看到一则具有出色数学解释：

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

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);不同，因为将float转换为整数不仅会除以2，还会将指数除以2，并导致其他一些伪像，但仍然归结为一些系数a0，a1， a2... 。

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

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

他们可以在一个循环中完成更多迭代，每次迭代都会提高结果，直到满足所需的精度为止。这正是它在CPU / FPU中的工作方式！但是似乎只有一次迭代就足够了，这对于速度也有好处。 CPU / FPU会根据需要进行尽可能多的迭代，以达到存储结果的浮点数的准确性，并且它具有适用于所有情况的更通用的算法。

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

使用(几乎)与CPU / FPU相同的算法，针对1 / sqrt(x)的特殊情况利用初始条件的改善，而不是一路计算出精确的CPU / FPU将会到达但停止得更早，因此提高计算速度。

相关讨论

将指针强制转换为long类型是log_2(float)的近似值。将其投射回去大约2 ^长。这意味着您可以使比率近似线性。