魔力数字与快速平方根倒数算法

题图来自于世界名模KK的推文，图中的guy's code就是这篇文章要介绍的算法，它的作用是快速计算平方根倒数，在计算机图形学领域有着广泛的应用。其中"i=0.5f3759df-(i>>1)"是迄今为止我见过的最令我瞠目结舌的一行代码。

一.来源与背景

这段代码最早出现在《雷神之锤III竞技场》3D引擎的源码中，2002年左右被人发到了论坛上，并引发了广泛的讨论，下面是代码的原始版本（去掉了c预处理器指令，保留了源码注释）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

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;
}

二.IEEE754浮点数

要了解它的工作原理，需要先理解浮点数在计算机中的存储方式，大部分编程语言都遵循IEEE754浮点数规范，该规范规定了四种表示浮点数值的方式，其中比较常用的有单精度(32位)，双精度(64位)。c与java使用float指代单精度，double指代双精度，而js默认采用IEEE754双精度来存储数值类型。下表是一个32位浮点数实例：

下面简单介绍它们的存储格式
1.符号位(sign)
符号位使用1来代表负数，0代表正数。
2.指数位(exponent)
指数可以是正指数或是负指数，但其表示方式并不遵循传统的二进制补码规则，而是先将其视为无符号数，再让指数偏移固定值来代表实际值，以单精度浮点数来举例，IEEE754标准规定偏移值为127($2^{b-1}-1$，b是指数的位数，此例是8)，因此指数部分需要减去127从而得到实际值。比如上例中指数部分为11111111，转换成10进制就是255，再减去127得到128，这里128就是它的实际指数。再比如00000000表示的就是-127，按照这个规则，它能够表示-127~128范围内的数字。（假设采用二进制补码规则，11111111表示的是-1，能够表示的数值范围是-128~127）
3.尾数位(mantissa)
IEEE754中的尾数又称为有效数(significant)，它用来表示科学计数法的小数部分，同时它规定第一个有效数必定是1，因此1就不需要存储，例如单精度浮点数有23位用于存储有效数，加上最左边没有存储的1，一共就有24位。继续拿上表的例子来说，可以把它的尾数左边四位看成1010，因此它的有效数实际值就是$(2^{23}+2^{21})\times 2^{-23}=1+2^{21}\times 2^{-23}$。
最后，我们可以得到上表中的浮点数实际值为：$$(-1)^s\times (1+m\times 2^{-23})\times 2^{e-127}=(-1)^0\times(1+2^{21}\times 2^{-23})\times2^{255-127}=1.25\times 2^{128}$$

三.一行神奇的代码

下面来分析"i=0x5f3759df-(i>>1)"这行充满魔力的代码，它的作用是获取一个平方根倒数的近似值。其工作原理证明过程非常的复杂(参考论文Fast Inverse Square Root)，这里提供wikipedia上的证明版本（尽管这是我个人认为最简单的证明版本，但看起来依然惨不忍睹，不求甚解的读者可以跳过证明部分）：
我们要计算的值是$\frac{1}{\sqrt{x}}$，令$y=\frac{1}{\sqrt{x}}=x^{-\frac{1}{2}}$，等式两边同时取对数可以得到$logy=log(x^{-\frac{1}{2}})=-\frac{1}{2}logx$。
现在假设输入值x是一个n位的IEEE754浮点数，由于要求输入是正数，因此x的符号位一定是0，假设它有i个指数位，尾数位就是n-i-1。设它的指数位为$E_x$，表示的指数是$e_x$，尾数位为$M_x$，表示的有效数字为$1+m_x$，根据上面的对浮点数表示方式的定义，可以得出：$$e_x=E_x-B$$$$m_x=\frac{M_x}{L}$$
其中$B=(2^i-1)$，$L=2^{n-i-1}$，此时输入的x的值为$(1+m_x)\times 2^{e_x}$，同理y值也可以表示成$(1+m_y)\times 2^{e_y}$，将它们代入上面的对数方程可以得到：$$log((1+m_y)\times 2^{e_y})=-\frac{1}{2}log((1+m_x)\times 2^{e_x})$$$$\Rightarrow log(1+m_y)+e_y=-\frac{1}{2}log(1+m_x)-\frac{1}{2}e_x$$现在假设我们用$\sigma_a$来表示$log(1+m_y)$与$m_y$之间的误差，$\sigma_b$表示$log(1+m_x)$与$m_x$之间的误差，即$\sigma_a=log(1+m_y)-m_y$，$\sigma_b=log(1+m_x)-m_x$，将它们代入上式得到：$$m_y+\sigma_a+e_y=-\frac{1}{2}(m_x+\sigma_b)-\frac{1}{2}e_x$$$$\Rightarrow \frac{M_y}{L}+E_y-B+\sigma_a=-\frac{1}{2}(\frac{M_x}{L}+\sigma_b)-\frac{1}{2}(E_x-B)$$$$\Rightarrow M_y+E_yL=\frac{3}{2}BL-(\sigma_a+\frac{\sigma_b}{2})L-\frac{1}{2}(M_x+E_xL)$$
见证奇迹的时刻到了，回顾一下这张表：

由于$L=2^{n-i-1}$，因此$M_x+E_xL$，$M_y+E_yL$正好是x与y所代表的整数值，如果我们用R来代表原方程中$\frac{3}{2}BL-(\sigma_a+\frac{\sigma_b}{2})L$这一部分，会得到$I_y=R-\frac{I_x}{2}$，将除以2替换成计算机更喜欢的右移位后，就成了这行神奇的代码。现在只剩下最后一个问题了，R值是多少？

四.魔力数字

根据上面段落的推理过程，可以知道魔力数字R应该等于$\frac{3}{2}BL-(\sigma_a+\frac{\sigma_b}{2})L$，其中B与L都是常量，如果是32位浮点数，则B等于127，L等于$2^{23}$。而$\sigma_a=log(1+m_y)-m_y$，$\sigma_b=log(1+m_x)-m_x$。现在可以构建一个函数f(x)=log(1+x)-x，由于$m_x$，$m_y$均为尾数部分，因此自变量x的取值范围是[0,1]，到这里我的想法是利用积分计算出函数f(x)=log(1+x)-x在[0,1]区间的平均值，从而得出$\sigma$的近似值，即：
$$\int_0^1\left[log(1+x)-x\right]dx$$如果我没算错的话，结果应该是$2-\frac{1}{ln2}-\frac{1}{2}\approx0.05730495911$。(wikipedia对$\sigma$的计算结果是0.0450461875...，我没看懂它的计算方式，如果有数学高手看到这里希望能指导一下)。最后将该值带回原式：
$$1.5\times2^{23}\times(127 - 1.5\times0.05730495911)\approx1597308761=0x5f34ff59$$同样这里我的计算结果与源码(0x5f3759df)不同，我也不知道错在哪里，希望有大神指导。

五.牛顿法

至于代码中

1
2

	y  = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
//	y  = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed

六.总结

关于这段代码的出处至今仍有争议，相传是由天才程序员约翰·卡马克所开发（约翰·卡马克是编程领域畅销书DOOM启示录的主角），然而不管这段代码出自于谁，我都觉得是码农届的荣耀，毫不夸张的说我因为0x5f3759df这个数字在家花了一个月时间重温高数，还有里面的"evil floating point bit level hacking"，深刻的让我觉得这段代码的作者是用灵魂在撸代码。

参考链接：

• 平方根倒数速算法
• Fast Inverse Square Root