求平方根倒数的算法

下面这个求$1/\sqrt{x}$ 的函数号称比直接调用 sqrt 库函数快 4 倍，来自游戏 Quake III 的源代码。

float InvSqrt (float x){
    float xhalf = 0.5f*x;
    int i = *(int*)&x;
    i = 0x5f3759df - (i>>1);
    float y = *(float*)&i;
    y = y*(1.5f - xhalf*y*y);
    return y;
}

我们这里分析一下它的原理（指程序的正确性，而不是解释为何快）。

分析程序之前，我们必须解释一下 float 数据在计算机里的表示方式。一般而言，一个 float 数据$x$ 共 32 个 bit ，和 int 数据一样。其中前 23 位为有效数字$M_x$ ，后面接着一个 8 位数据$E_x$ 表示指数，最后一位表示符号，由于这里被开方的数总是大于 0 ，所以我们暂不考虑最后一个符号位。此时

如果我们把计算机内的浮点数$x$ 看做一个整数$I_x$ ，那么

现在开始逐步分析函数。这个函数的主体有四个语句，分别的功能是：

1. int i = *(int*)&x; 这条语句把$x$ 转成$i=I_x$ 。
2. i = 0x5f3759df - (i>>1); 这条语句从$I_x$ 计算$I_{1/\sqrt{x}}$ 。
3. y = *(float*)&i; 这条语句将$I_{1/\sqrt{x}}$ 转换为$1/\sqrt{x}$ 。
4. y = y\*(1.5f - xhalf\*y\*y); 这时候的 y 是近似解；此步就是经典的牛顿迭代法。迭代次数越多越准确。

关键是第二步 i = 0x5f3759df - (i>>1); 这条语句从$I_x$ 计算$I_{1/\sqrt{x}}$ ，原理:

令$y=1/\sqrt{x}$ ，用$x=(1+m_x)2^{e_x}$ 和$y=(1+m_y)2^{e_y}$ 带入之后两边取对数，再利用近似表示$\log_2(1+z)\sim z+\delta$ ，算一算就得到

若取$\delta=0.0450465679168701171875$ ，$\frac{3}{2}(127-\delta)2^{23}$ 就是程序里所用的常量 0x5f3759df。至于为何选择这个$\delta$ ，则应该是曲线拟合实验的结果。

Q. E. D.