算法趣谈(2)——快速平方根逆
float inv_sqrt(float number) {
uint32_t i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = *(uint32_t *) &y;
i = 0x5f3759df - (i >> 1);
y = *(float *) &i;
y = y * (threehalfs - (x2 * y * y));
// y = y * (threehalfs - (x2 * y * y));
return y;
}
这是一段奇怪却又极富盛名的代码，通过一些使用令人费解的操作之后和奇怪的 0x5f3759df，返回了 number 的平方根的倒数，即 \(\frac1{\sqrt{number}}\). 这段代码就是 Quake III Arena ，作者为 John Carmack.
本文将介绍这段代码的原理。
IEEE 754
要想了解这个算法，我们得首先看看浮点数在计算机下的表示方法，即 IEEE 754 标准。
假设数字 \(x = 0.15625 = {0.00101}_2 = 1.01 \times 2^{-3}\)，那么在 IEEE 754 中，其表示为：
- 符号位：\(0\)，表示正数.
- 指数位：\(E = 01111100\)，为 \(-3 + 127 = 124\) 得到，其目的是能够表示正负指数。
- 尾数位：\(M = 01000000000000000000000\)，其真实值为 \(1 + \frac{M}{2^{23}} = 1 + 0.25 = 1.25 = 1.01_2\). 这里的 1 是因为 IEEE 754 中默认尾数位为 \(1.M\)，所以可以节省一位。
我们可以得到
\[x = (1 + \frac{M}{2^{23}}) \times 2^{E - 127}\]对数的魔法
接下来我们要做的就是看看我们能否尽可能的利用等价变形化简这个值并进行近似。 看到指数时一个自然的想法就是取对数，即
\[\begin{align*} x &= (1 + \frac{M}{2^{23}}) \times 2^{E - 127}\\ \log_2 x &= \log_2 (1 + \frac{M}{2^{23}}) + E - 127\\ &\approx \frac{M}{2^{23}} + \mu + E - 127\\ &= \frac{1}{2^{23}}(M + 2^{23} \times E) + \mu - 127 \end{align*}\]注意到 \(M + 2^{23} \times E\) 恰好 是 \(x\) 二进制下的表示。这一洞察为我们提供了一个近似的方法。
主算法
有了以上的铺垫，我们就可以开始看看这个算法了。
Bit Hack
y = number;
i = *(uint32_t *) &y;
这是利用指针和地址将 \(y\) 作为无符号整型来得到该浮点数的二进制表示，不多赘述，后面的代码也是同理，只不过是将整型转为浮点型。
Magic Number
i = 0x5f3759df - (i >> 1); // 这里的 i 即为 M + 2^23 x E
为了计算 \(z = \frac{1}{\sqrt{y}}\)，我们还是取对数，即 \(\log_2 z = -\frac12 \log_2 y\)，而根据上面的推导，我们有
\[\frac{1}{2^{23}}(M_z + 2^{23} \times E_z) + \mu - 127 = -\frac12(\frac{1}{2^{23}}(M_y + 2^{23} \times E_y) + \mu - 127)\]不难得到
\[\begin{align*} M_z + 2^{23} \times E_z &= \frac32 2^{23}(127-\mu) - \frac12(M_y + 2^{23} \times E_y)\\ &= Magic - (i >> 1); \end{align*}\]注意 \(\mu\) 是上文近似对数函数得到的，对于 \(\mu\) 值的选择也就决定了 magic number 的值。
枚举还是计算？ 在以前听到这个算法的时候有个最大比较大的争议是 0x5f3759df 这个神奇的数字究竟是纯粹数学计算出来的还是枚举出来的。通过维基百科资料我的感觉是这是在一个范围内枚举出来的，因为这个值本身即使考虑牛顿迭代法也不是最优的。所以一个可能的结果是在通过计算得到一个范围，例如 0x5f37642f 这个在某些分布下表现较好的值，然后再基于这个值上下枚举出一个实际使用最合适的。
牛顿迭代法
y = y * (threehalfs - (x2 * y * y));
最后一步较为简单，是利用牛顿迭代法来进行近似，回顾一下牛顿迭代法的原理：
\[y_{n+1} = y_n - \frac{f(y_n)}{f'(y_n)}\]而对应于 \(f(y) = \frac{1}{y^2} - x\)，我们有
\[\begin{align*} f'(y) &= -\frac2{y^3}\\ y_{n+1} &= y_n - \frac{\frac1{y_n^2} - x}{-\frac2{y_n^3}}\\ &= y_n + \frac{y_n - xy_n^3}{2}\\ &= y_n(\frac32 - \frac{x}{2}y_n^2) \end{align*}\]也就对应了上面的代码。
意义与应用
快速平方根逆算法的应用或者说是来源主要是图形学中的单位向量计算。如 \(v = (v_1, v_2, v_3)\), \(\Vert v \Vert = \sqrt{v_1^2 + v_2^2 + v_3^2}\), 则单位向量为 \(\frac{v}{\Vert v\Vert}\), 也就要用到平方根逆算法。
但随着硬件的发展，这个算法已经没那么重要了，例如在 SSE 中有一个指令 rsqrtss 就是专门用来计算平方根逆的，速度和准度都比这个算法要好。当然这些并不影响这个算法的魅力。