魔数 — 平方根倒数速算
1 | float Q_rsqrt( float number ) |
“魔数”通常指代码中凭空出现的数字。在游戏中，渲染物理效果时，会涉及平方根倒数的运算。1999年发布的《雷神之锤3》中有这样一段C语言代码，只用了4行，就可以快速的计算平方根倒数，读过的人都会说 what the fuck !
$$f(x) = \frac 1 {\sqrt{x}}$$
现在的计算机会对这种计算做特殊优化，但是在20年前使用正常的计算方式会涉及到浮点数的除法运算，在二进制计算机中，加法和乘法运算非常快，但是减法和除法运算相对较慢，这个算法巧妙的避免了除法运算，并且在数据的存储方面，浮点数比较特殊，它存储方式不利于位运算，这个算法使用C语言的未定义规则对浮点数进行了位运算，进一步加快了算法的速度。
代码分别涉及到了:
- IEEE二进制浮点数算术标准(IEEE 754)
- C语言UB(Undefined Behavior)
- 平方根倒数估算
- 牛顿迭代法
IEEE 754
**IEEE754**是在20世纪80年代使用最广泛的浮点数标准。
符号位 sign：通常表示浮点数的正负。
阶码 exponent：其中阶码通过移码来表示，移码就是实际的编码值加上一个固定的值，在IEEE754中该固定的值为$2^{n-1}-1$。
尾数码 fraction：该部分用原码表示，也称为有效位，尾数码的长度决定了该浮点数的精度。
举个例子
- 用IEEE754来表示 9.625
确定符号位
- 9.625是一个正数，因此符号位为0。
转换为二进制
- 整数部分9的二进制形式是1001。
- 小数部分0.625转换为二进制可以通过不断乘以2并取整数部分的方式得到，直到小数部分为0或者达到指定精度为止。由此得到小数部分的二进制形式是0.101。
- 合并整数部分和小数部分，得到完整的二进制表示形式是1001.101。
规格化
- 上面得到二进制数继续用规范的浮点数表达应为 1.001101×2^3
确定阶码
- 这里的指数为 3，加上单精度浮点数的偏移量127，得到真正的阶码值：3+127=130，即二进制的 1000 0010。
填充尾数
- 有效数字省略掉小数点左侧的1之后为001101，然后在右侧用零补齐。
组合32位浮点数
- 符号位（1位）：0
- 阶码位（8位）：1000 0010
- 尾数位（23位）：001 1010 0000 0000 0000 0000
1 | 0 | 10000010 | 00110100000000000000000 |
现在对这个二进制数进行特殊表示，设阶码的真值为E，
$$
E = 1000,0010_{(2)} = 130_{(10)}
$$
设尾数为M，则有
$$
M =001 1010 0000 0000 0000 0000_{(2)} = 0.001101 * {2^{23}}_{(2)}
$$
此时抛开IEEE754浮点数规则，设这个二进制数为L，L的值就等于，阶码乘上 2的23次方 再加上尾码的值
$$
L = E × 2^{23} +M
$$
设 F 为 L 的十进制，则可以用下列公式表示F
$$
F = (1+{\frac {M} {2^{23}}}) ×2^{E-127}
$$
现在对F取对数
$$
\log_2 F = \log_2(1+{\frac {M} {2^{23}}}) + E -127
$$
还记得M是什么吗，它是浮点数的尾数，取值范围在[0,1)
结合$y=log_2(1+x)$和$y = x$图像可以得出
$$
\lim\limits_{x \to 0+} log_2 (1+x) = x
$$
$$
\lim\limits_{x \to 1-} log_2 (1+x) = x
$$
由此可知在[0,1]，取系数$\mu$得
$$
log_2 (1+x) ≈ x + \mu
$$
现在对$log_2F$用上述公式进行简化
$$
\begin{align*} \log_2 F &= \log_2(1+{\frac {M} {2^{23}}}) + E -127 \ &≈ (\frac M {2^{23}}+\mu) + E -127 \ &= ({\frac M {2^{23}}} +{\frac {E×2^{23}} {2^{23}}}) + \mu - 127 \ &= \frac 1 2^{23}(M + E×2^{23}) + \mu -127 \end{align*}
$$
还记得之前的二进制数L吗？
$$
L = E × 2^{23} +M
$$
结合两个公式可以得出
$$
log_2F = \frac 1 2^{23} L + \mu -127
$$
由此得知浮点数F的对数，可以通和它的二进制数L进行变化得出，此时对于$y = \frac 1 {\sqrt{x}}$两边同时取对数，用上述公式简化
$$
\begin{align*} log_2 y &= log_2 (\frac 1 {\sqrt x}) \ log_2y &= -\frac 1 2 log_2x \ \frac 1 2^{23} L_y + \mu -127 &= -\frac 1 2(\frac 1 2^{23} L_x + \mu -127) \ L_y &= 3 × 2^{22}(127 -\mu) - \frac 1 2 L_x \end{align*}
$$
此时当这个值$\mu=0.450465$，$3 × 2^{22}(127 -\mu)$ 就等于1597463007，它的16进制就是我们说的魔数0x5f3759df
$$
\begin{align*} L_y &= 3 × 2^{22}(127 -\mu) - \frac 1 2 L_x \ &= 0x5f3759df - \frac 1 2 L_x\end{align*}
$$
此时终于能看出点源码的样子了0x5f3759df - ( i >> 1 )还需要计算一个除法，看起来很简单啊，对一个数除2就对它二进制数字算数右移1位，但是别忘了此时的i是一个浮点数，无法进行位运算
那就把它转为整形吧，普通的浮点数转整形会舍去小数点后的值，如果y=3.14执行long i = (long)y 后i的值为3 。看看作者如何实现的吧
1 | i = * ( long * ) &y; //将存放在y中的浮点数（float）用长整形（long）取出 |
众所周知，每个变量会有一个十六进制的地址。 &y 可以获取变量y的地址，通过 * 可以获取该地址中存放的值。( long * ) &y这段代码对y变量地址进行了强制类型转换，让编译器认为这个浮点地址是长整形的，最后通过* 完整的取出了该地址中的二进制值，并赋值给变量i 。
到此为止变量i中就存放了一串浮点数的二进制数字，再通过一样的方法把它转回浮点数，就完成了对浮点数的位运算
1 | y = * ( float * ) &i; |
最后通过牛顿迭代可以求出更精确的值
1 | y = y * ( threehalfs - ( x2 * y * y ) ); |
总结
20年后的今天再欣赏这段代码，我仍然忍不住要说一句what the fuck，即便现在计算机可能自带更快速的计算方式，但是段代码涉及的知识和它对数据的处理方法值得我们研究和学习。
