【回归本源】番外1-雷神之锤3的那段代码
这次来讲个小知识点，在计算机图形学中时常会用到的数值是 \frac{1}{\sqrt{x}} ，但是计算这个值并不是一个容易的事。我们知道计算机本质上只会做加法，减法也是用加法进行的。而乘法和除法的本质也是对应着加法与减法，也是因此相比与除法，乘法会更快一些。此前听说过《雷神之锤Ⅲ》(Quake Ⅲ Arena)中有这么一段代码让人直呼"WTF"就是用来快速计算这个值的，但直到昨天麦子才真正去看它，看完真的是不由得心生敬佩。
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
return y;
}感觉十分有必要写一下自己对这段代码的学习以及理解了。
浮点数的表达
在数学中如果我们考虑一个函数 f(x) = \frac{1}{\sqrt{x}} 时，这里的 x 自然是认为为实数（没考虑复数）。但是计算机时没有办法处理任意实数的，毕竟无限不循环这件事在物理层面是实现不了的。计算机中对实数有两种表达: double 和 float，两者主要差异是占用内存的大小区别。我们重点放在float，因为学会float理解double就容易了，而且在游戏计算中对计算速度要求很高，float也会更常用。
浮点数，全名应该是floating-point number，其实是一系列的规定。为了理解float，我们先思考计算机如何使用二进制表示数字。让我们一步一步来吧
正整数
举个例子，
(234)_{10} = 2 \times 10^2 + 3 \times 10^1 + 4 \times 10^0
我们对十进制的理解就是将数字从0数到9就进一位，因为没有一个数字或字母来直接表示10的含义，而二进制中只有0和1，也是因此 (1)_2 + (1)_{2} = (10)_{2} ，我们将每一个二进制的位称为比特 (bit) 。
但我们在写程序时也不是跟每一个bit打交道，由于历史以及某些硬件原因，我们一般与字节(byte) 打交道，也就是8个bits一起。那么每1个byte能代表的数字就是 0 ~ 255，即
(11111111)_{2} = 1 \times 2^7 + 1 \times 2^6 + 1 \times 2^5 + 1 \times 2^4 + 1 \times 2^3 + 1 \times 2^2 + 1 \times 2^1 + 1 \times 2^0 = (255)_{10}
整数
当然，事情自然不会这么简单，如果我们要表示负数咋办？如果你了解过如何用加法来进行减法，那就不难理解为啥我们会有一个咋一看很奇怪的规则：第一位bit的0、1来代表是否为负数，于是乎我们上面的那个例子就变成了
(11111111)_{2} = - (1 \times 2^6 + 1 \times 2^5 + 1 \times 2^4 + 1 \times 2^3 + 1 \times 2^2 + 1 \times 2^1 + 1 \times 2^0) = (-127)_{10}
1个byte能代表的数字也就变成了 -127 ~ 127
正负浮点数
好的，当你能理解上面的设定后，就可以来看看float的设定了，一般而言，如果不加以设定float一般时32 bits，也就是4 bytes的。其中
第1位bit是正负号；8个bits（1个byte）是幂 (Exponent，后面记为E)；后23个bits是小数位 (Mantissa, 后面记为M)
设定归设定，啥意思呢？还是先想十进制
(1.234)_{10} = 1\times10^0+2\times10^{-1} + 3\times10^{-2} + 4\times10^{-3}
好的，但是如果更大或更小的情况呢？比如 (1234000000)_{10} 或 (0.00001234)_{10} 呢，相信你可能就会用科学计数法了，即 (1.234 \times 10^9)_{10} 和 (1.234\times 10 ^{-5})_{10} 。这样就方便很多了，同样的，float的设定本质上算是科学计数法的拓展。
任意float可以通过下式转换为被表示为
\pm 2^{((E)_{10} - 127)}\times(1.0 + \frac{(M)_{10}}{2^{23}})
想看各种数如何变成二进制的可以去 IEEE-754 Floating Point Converter 试一下
看一个例子
(0.1)_{10} \approx (0,01111011,10011001100110011001101)_2
这里我用‘，’将符号、幂和小数位分开了
其中
(E)_{2} = (01111011)_{2} = (123)_{10} = (E)_{10}
(M)_{2} = (10011001100110011001101)_{2} = (5033165)_{10} = (M)_{10}
运用上面的公式我们得到了 0.10000000149011612 ，可以看出，float不能百分百还原我们要的小数，但是是一个不错的近似。
如果你不觉得这像是科学计数法，那我们就换一种写法
(10011001100110011001101)_{2} \\ = 1 \times 2^{-1} \\ + 0 \times 2^{-2}+ 0 \times 2^{-3} + 1 \times 2^{-4}+ 1 \times 2^{-5} \\ + 0 \times 2^{-6} + 0 \times 2^{-7}+ 1 \times 2^{-8} + 1 \times 2^{-9} \\ + 0 \times 2^{-10}+ 0 \times 2^{-11}+ 1 \times 2^{-12}+ 1 \times 2^{-13} \\ + 0 \times 2^{-14}+ 0 \times 2^{-15}+ 1 \times 2^{-16}+ 1 \times 2^{-17} \\ + 0 \times 2^{-18} + 0 \times 2^{-19} + 1 \times 2^{-20} + 1 \times 2^{-21} \\ + 0 \times 2^{-22}+ 1 \times 2^{-23} \\ = (0.6000000238418579)_{10}
再加一个1.0，即默认会有的 2^0 ，科学计数法是不以0做开头的（例如 0.23 \times 10^2 不是科学计数法，应该是 2.3 \times 10^1 ）。由此再看，我们这个就像是 +(1.M_1M_2M_3...M_{23})_2 \times 2^{-4} ，这里 M_i 指小数位 M 的第 i 位数，只不过每一位都是要乘2的对应幂，而不是10。
浮点数 -> 整数
刚我们看来如果你要将float从二进制换算为十进制，需要用将 E 和 M 按十进制读出来，然后再用公式
\pm 2^{((E)_{10} - 127)}\times(1.0 + \frac{(M)_{10}}{2^{23}})
但这只是规定而已，计算机中本质上还就是那些0和1，所以如果你非要将一个float的二进制读成一个二进制的整数，那也没人能拦着你，就像雷神三中的
i = * ( long * ) &y; 就是将float y读成了一个整数（long int）
y = * ( float * ) &i;然后又返回去了。能想着这么做真不是什么省油的灯。
如果你理解了前面整数的二进制表达部分，就不难理解这里 i 通过这么一系列变成了
F2I = (E)_{10} \times 2^{23} + (M)_{10}
这就是转换后会得到的整数。
对数函数近似
要看懂雷神三的代码还有两部分，一部分是牛顿法，另一部分就是一个对数函数的近似。牛顿法是很简单的，放到最后吧。有意思的是一个对数函数近似，
\log_2(1 + x) \approx x + k
虽然从图像上不难看出没啥问题，但是确实也挺没道理的。就算用泰勒展开，也不能直接写出关于x的表达（底数不是e）。
但不管怎样，从图像中，我们可以看出：当x在0到1之间，k都在0到0.1以内。
魔法开始
好了，有了全部这些东西，让我们来看看那段 What Technically Fabulous codes 的魔法吧
i = 0x5f3759df - ( i >> 1 ); 让我们先考虑
y = \frac{1}{\sqrt{x}} = x^{-\frac{1}{2}}
\Rightarrow \log_2(y) = -\frac{1}{2} \log_2(x)
将 y, x 认定为float，带入我们的float表达
\Rightarrow \log_2(2^{(E_y - 127)}\times(1.0 + \frac{M_y}{2^{23}})) = -\frac{1}{2} \log_2(2^{(E_x - 127)}\times(1.0 + \frac{M_x}{2^{23}})) \\ \Rightarrow (E_y - 127)+ \log_2(1.0 + \frac{M_y}{2^{23}}) = -\frac{1}{2} (E_x - 127) -\frac{1}{2} \log_2(1.0 + \frac{M_x}{2^{23}})
可以看到我们的对数函数内是符合之前那个近似的条件的，使用那个近似
\Rightarrow (E_y - 127)+ \frac{M_y}{2^{23}} + k = -\frac{1}{2} (E_x - 127) -\frac{1}{2} \times ( \frac{M_x}{2^{23}}+ k) \\ \Rightarrow E_y \times 2^{23} + M_y = \frac{3}{2}(127 - k)2^{23}-\frac{1}{2}(E_x \times 2^{23} + M_x)
这个形式是很有意思的，因为右边 E_x \times 2^{23} + M_x 是将需要进行操作的 x 的float强行看成整数时的表达，左边 E_y \times 2^{23} + M_y 是做完操作后的 y 看成整数时的表达。
所以代码中
i >> 1 就对应着
\frac{1}{2}(E_x \times 2^{23} + M_x)
位操作，向右移，差不多除以2，就不多细讲了，想明白整数的二进制表达就好。
而最后的magic number，0x5f3759df到底是什么鬼？
就应该是对应着我们的常数
\frac{3}{2}(127 - k)2^{23}
我们可以逆推一下，0x5f3759df是十六进制（想想二进制和十进制，你就应该能想明白），它换成十进制就是1597463007。那么k就应该是0.0450，与我们之前说k的范围是对应的。
i = * ( long * ) &y; // evil floating point bit level hacking
i = 0x5f3759df - ( i >> 1 ); // What the fuck?
y = * ( float * ) &i;再看一眼，也就是说这段代码就是将float y看成一个整数i，然后通过一个奇妙的常数去减去除以了2的i，这时就已经完成了 \frac{1}{\sqrt{y}} 的一个估计计算了，再让y等于将这个整数i看成float的情况。
妙啊，妙啊！
牛顿法
近似毕竟是近似，最后他们又用了一、两次牛顿法来得到一个更好的值，但这个就很好理解了。
x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}
以上就是牛顿法的形式了，这里
f(y) = \frac{1}{y^2} - x = 0
x 是输入，是一个定值，那么
f'(y) = - \frac{2}{y^3}
\Rightarrow y_{n+1} = y_n - \frac{\frac{1}{y_n^2} - x}{ - \frac{2}{y_n^3}} = y_n(\frac{3}{2} - \frac{x}{2} y_n^2)
好了，再看最后两句
y = y * ( threeHalfs - ( x2 * y * y ) ); // 1st iteration
y = y * ( threeHalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed牛顿法，朴实无华呐。
这段代码真的能看出一个程序员的功力啊，内心需要始终铭记计算机不是数字，而是数据。真的是让我大为震撼。
