【魔术数字 Magic Number】神奇的平方根倒数快速估算法

引言

我们把一些出现在代码中的不明所以却又发挥重要作用的数字称作魔术数字(Magic Number)。

这通常是因为程序员不加注释造成的后果，数字的意义只有本人清楚（更准确地说，可能只有本人在写代码的时候还记得清楚）。

我们不妨学习一下费马先生，解释一句：我已经想到一个绝妙的数字用来解决这个问题，但这里的地方太小，写不下完整的注释。

而在算法界流传着一个著名的魔术数字案例：用来快速估算平方根倒数的0x5f3759df，第一个提出这种算法的神人真的可以为自己辩驳一句：真的不是几句注释能解释得明白的呀。

问题引入

首先，先问问是什么？所谓平方根倒数，就是给定$y$，要计算下式：
$$x=\frac{1}{\sqrt{y}}$$
再问问为什么需要这种算法呢？在图形学，尤其是游戏计算场景中，经常需要将某个向量的进行归一化（标准化），即Matlab中的normalize操作，将某向量的各个值除以向量的模，使得该向量变成单位方向向量。如向量$(x,y,z)$执行归一化，即变为：
$$(\frac{x}{\sqrt{x^2+y^2+z^2}},\frac{y}{\sqrt{x^2+y^2+z^2}},\frac{z}{\sqrt{x^2+y^2+z^2}})$$
这样一来，每次标准化都涉及到大量的平方根倒数计算问题。众所周知，计算机在计算加法和乘法时速度很快，而无论是计算平方根还是计算倒数（即除法），速度都是要大打折扣的。这也就是为什么我们需要一个算法来快速计算平方根倒数。

神奇的代码

接下来，让我们看看这个算法的代码是怎么做得呢？
据说这个算法很早就被提出了，作者已经不可考证（据说不是卡马克大神，被本人否认了）。但是它被大众所熟知的原因，是曾经风靡一时的《雷神之锤3竞技场》的源代码中大量使用了这一段代码。这段代码长这个样子：

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

短短几行代码，很好地诠释了什么叫大道至简。代码一经公布，便让大家都摸不着头脑了——0x5f3759df? 这是什么玩意儿？为什么能速算平方根倒数呢？第一次看到或许内心真的要爆一句注释里的那句粗口了。

不着急，要解释清楚的话还真的需要花点时间，先来学习/复习一下前置知识。

前置知识

移码

类似于原码、反码、补码，移码也是一种二进制数字的表示方式，通常用于浮点数的阶码的表示（这个概念后续再说），也称作增码/偏置码。

从形式上说，其实就是补码的符号位取了个反。考虑这样一种做法：若有8个二进制位，原码表示范围为$[0,255]$，若要表示负数，我们可以偏移128，这样表示范围变成了$[-128,127]$，这就是移码的表示方式。

128即对应8位二进制位中的最高位，所以最高位可以理解成符号位，只不过这个符号位的01与通常意义上的相反，0表示负数，1表示非负数，这样的好处是更容易比较大小（由于是直接偏移得到的，所以不管是负数、零、还是正数，直接比较就行了）。

举个例子，3转换为8位移码为 $10000011$ ，0表示为 $10000000$ ，-3表示为 $01111101$ 。
简单来说，就是把原数字加上128转成原码就行了，反过来就是减128。

浮点数表示法（IEEE 754）

引言

在计算机中，小数的表示方法是一个较为麻烦的事情，更准确的说，是浮点数的表示。
浮点数这一概念相对的是定点数，这里的点指的是小数点的位置，例如整数就是定点数，小数点的位置固定在个位后，再例如如果规定一种表示金额的数，xx.xx元，精确到分，那么其实这个小数也是定点数，因为我们知道小数点总是在倒数第二位。

而对于一般的小数，小数点的位置是浮动的，应该如何表示呢？一种直观的想法是，例如有32个二进制位，我们可以规定第1位为符号位，第2-16位为小数点前的部分，后16位为小数点后的部分，但是这样32个二进制位表示的浮点数范围只能是$[-2^{15}+2^{-16},2^{15}-2^{-16}]$，能表示的范围过小。这是因为这种表示法其实也是强行规定了小数点的位置，算是定点数表示。

考虑一下科学计数法，我们在表示小数时也可以用科学计数法，我们通常更关心小数的有效数字位数。若将小数转化为二进制，如$3.25$，二进制表示为$11.01$，类比科学计数法，可以写作：
$$11.01_{(2)}=1.101\times 2^1$$
再深入分析一下更普适的情况，任意一个小数转成二进制后可以表示为以下形式：
$$N\times 2^e$$
由于二进制中只有0和1，所以如果该浮点数不是0，写成这种形式的话其中的$N$一定是$1.X$这种形式，可以表示为：
$$1.M\times 2^e$$
至于$0$的情况，我们可以特殊处理，这样可以节约一个bit位，默认第一位为1。
这便是IEEE 754的一个基本思想，了解了这些，接下来我们来正式说一说IEEE 754的标准。

IEEE二进制浮点数算术标准（IEEE 754）

以单精度浮点数为例，有32个二进制位，规定第1位为符号位（0正1负），之后8位表示指数阶码，最后23位表示尾数，即上式中的$M$。注：双精度浮点数也是同理的，只不过是有64个二进制位，依次划分为1+11+52，另外还有延伸单精确度与延伸双精确度，这里不再展开。

指数阶码，这里IEEE 754规定，阶码偏移量为127，即原来的8个bit位表示的$[0,255]$偏移到$[-127,128]$（实际是$[-126,127]$，因为-127和128被用于表示特殊数字了，后面会详细展开）。我们用小$e$表示实际的指数，大$E$表示阶码，即$e=E-127$。这里为了方便叙述，后文提到的指数默认指代实际的指数，阶码默认指代指数阶码。

这里需要注意几点问题，首先一些特殊数字的处理，比如$0$的表示。为了解决这个问题，IEEE 754规定有以下几种类型的数字：

数字类型	表示方法
零	尾数M为0且阶码为0（即e=-127），根据符号位的不同有$\pm 0$的区别
规则化数	即常规的数字的表示，阶码范围在$[1,254]$
非规格化数	阶码为0但是尾数M不为0，表示一些绝对值过小的小数（即M前的位数默认不再为1而为0）
无穷	阶码为255且尾数M为0，根据符号位的不同有$\pm \infty$的区别
非数NaN	阶码为255且尾数M不为0

至于为什么偏移量是127不是128，个人理解是因为虽然这两个数字都可以保证正负表示范围的“对称”，但是偏移127的表示范围$[-126,127]$，偏移128的表示范围为$[-127,126]$，偏移127的表示范围更大。

除此以外还有一些溢出、计算等规则，这里也不再展开，我们最关心的是规格化数的表示，我们可以写出通式，后面会用到：
$$(-1{)}^S\times 1.M\times 2^{E-127}$$

牛顿迭代法

牛顿迭代法（Newton’s method）是一种求解非线性方程（或方程组）近似根的迭代算法。它基于局部线性逼近的方法，通过迭代过程逐步改进对根的估计。牛顿迭代法的收敛速度通常很快，特别是在初始估计值接近实际根的情况下。

假设我们要求解非线性方程$f(x)=0$的根。牛顿迭代法的基本思想是从一个初始的根的估计值$x_0$开始，然后通过以下迭代公式逐步改进对根的估计：
$$x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)}$$
其中，$x_n$ 是第$n$次迭代的根的估计值，$f^{\prime }(x_n)$ 是 $f(x)$ 在 $x_n$ 处的导数值。

代码剖析

强制以整型方式操作浮点数

    i = * ( long * ) &y;                       // evil floating point bit level hacking
    i = 0x5f3759df - (i >> 1);                 // what the fuck?
    y = * ( float * ) &i;

先来看看这段代码的第1行和第3行。平方根倒数快速估算法，一个重要的部分是利用IEEE 754的规则，强制对浮点数进行位运算。而C语言中并不允许我们这么做，所以代码中用了很巧妙的一招，将浮点型的指针转成了整型的指针，让编译器将那些表示浮点数的32个bit位强制看成了整型，第1行操作实际上没有修改任何bit位的值，只是改变了编译器的“认知”，让我们可以直接以整型位运算的方式去操作浮点数的bit位。之后第3行再转回来，获取计算后的浮点数。

魔术数字 0x5f3759df

有了前序的铺垫，我们来看看这个算法的核心部分。

首先，我们的任务是计算$y=\frac{1}{\sqrt{x}}$，我们把$x$表示为IEEE 754的规格化形式（默认$x$为正数）：
$$x=1.M_x\times 2^{E_x-127}=(1+\frac{M_x}{2^{23}})\times 2^{E_x-127}$$
我们考虑直接从这个式子入手，去计算平方根的倒数，即：
$$y=\frac{1}{\sqrt{x}}=x^{-\frac{1}{2}}=[(1+\frac{M_x}{2^{23}})\times 2^{E_x-127}{]}^{-\frac{1}{2}}$$
数学上的直觉告诉我们，应该取个对数试一试，且取完对数之后可以发现，有一个$\log 2$，所以我们可以取以2为底数的对数，那么有：
$${\log }_{2}y=-\frac{1}{2}[{\log }_2(1+\frac{M_x}{2^{23}})+E_x-127]$$
接下来估算的部分来了，我们可以发现$\frac{M_x}{2^{23}}$其实表示的$1.M_x$的小数部分，所以$\frac{M_x}{2^{23}}$的值域为$[0,1)$，在这个区间上，$f(x)={\log }_2(1+x)$和$g(x)=x$是近似相等的，如下图所示：

如果再加上某个固定的常数μ，就更加接近了，如上图的中0.035，所以我们可以用$h(x)=x+\mu$来近似估算代替${\log }_2(1+x)$，于是有：
$${\log }_{2}y=-\frac{1}{2}(\frac{M_x}{2^{23}}+E_x+\mu -127)$$
再然后，又是一步巧妙的操作，我们把式子中的$\frac{1}{2^{23}}$提出来，有：
$${\log }_{2}y=-\frac{1}{2}\times \frac{1}{2^{23}}\times [(M_x+E_x\times 2^{23})+(\mu -127)\times 2^{23}]$$
离谱的地方来了，我们可以惊讶地发现，这个$M_x+E_x\times 2^{23}$，其实就是IEEE 754规定的浮点的表示形式（$E_x$在前，$M_x$在后），于是乎我们可以直接把内存中的32个bit位直接看成一个整体进行操作，这也就是为什么我们前面提到了要强制以整型方式操作浮点数。个人认为IEEE 754在设计的时候就考虑到了这个取log的问题，所以把E放在M前面了。

我们把最后的结果$y$也用这种方式类似操作，左式改写后，可以得到下式：
$$\frac{1}{2^{23}}\times [(M_y+E_y\times 2^{23})+(\mu -127)\times 2^{23}]=-\frac{1}{2}\times \frac{1}{2^{23}}\times [(M_x+E_x\times 2^{23})+(\mu -127)\times 2^{23}]$$
上式左边的$M_y+E_y\times 2^{23}$即为我们要求的答案，化简即得：
$$M_y+E_y\times 2^{23}=\frac{3}{2}(127-\mu )\times 2^{23}-\frac{1}{2}(M_x+E_x\times 2^{23})$$
上式的右式中，减号前为一个可以预处理的常值，我们预先设定一个合理的μ，然后计算$\frac{3}{2}(127-\mu )\times 2^{23}$，四舍五入再转为十六进制，即为代码中的魔术数字。所以这个魔术数字的取值与μ的选值有关系，据说理论上算得的最优值实际表现并不如这个代码里的魔数0x5f3759df，最后是通过枚举相对最优而选定的魔数的值。
这样一来就明了了，上面推得的公式即对应代码中这一行位运算：

    i = 0x5f3759df - (i >> 1);                 // what the fuck?

两次牛顿迭代

由于是利用μ进行的估算，已经非常接近于实际的准确值，我们可以用两次牛顿迭代来进一步逼近准确解。

本场景中，$y=\frac{1}{\sqrt{x}}$，其中$x$已知，我们要计算的实际上是$f(y)=y-\frac{1}{\sqrt{x}}=0$的根，而我们当然不希望去计算$\frac{1}{\sqrt{x}}$（废话，要算的就是这个啊哈哈），改写一下变为：
$$g(y)=x-\frac{1}{y^2}=0$$
对$g(y)$求个导：
$$g^{\prime }(y)=\frac{2}{y^3}$$
代入牛顿迭代法的公式中：
$$y_{n+1}=y_n-\frac{g(y_n)}{g^{\prime }(y_n)}=y_n(\frac{3}{2}-\frac{x}{2}\cdot y_n^2)$$
对应代码中的这一段（做了两次牛顿迭代）：

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

小结

现在，我们已经可以直接用1/sqrt(y)来计算平方根倒数，因为现在的计算机已经集成了更快更准的算法了。但是这段代码依然是有研究意义的。
神奇的代码，背后隐藏的是精妙的数学分析过程和对底层知识的透彻理解，这也提醒我们应当脚踏实地，打牢基础，同时培养数学推理和应用的能力。

关键词

平方根倒数速算 0x5f3759df 魔术数字 IEEE754 Magic Number 魔法数字 魔数 卡马克

参考文献

FAST INVERSE SQUARE ROOT