Quake 3 中的 InvSqrt 算法
亦闻楚歌起
·
2023-07-13 12:50:53
·
个人记录
[toc]
0. 前言
在一些图形化界面，我们很经常需要取单位向量。这就需要乘以模长的逆元。具体地，我们需要使用函数 f(x)=x^{-0.5} 。因为这个函数将要被调用很多次，除了给它加上 inline 之外，我们还希望它能有一定程度的优化。我们允许的代价是一定范围内的精度误差（微小的误差在图形化时可以忽略）。
广为流传的快速 InvSqrt 代码据说是来自 John Carmack 的 Quake 3，传世版本大概如下：
float invSqrt(float x) {
float xhalf = x * .5f;
int i = *(int*)&x; // get bits for float type.
i = 0x5f3759df - (i >> 1); // **Words having been harmonised**
x = *(float*)&i; // return to float type.
x = x * (1.5f - xhalf * x * x); // Newtown Step.
// x = x * (1.5f - xhalf * x * x); // Repetition of Newtown Step enhances accuracy.
return x;
}这段代码据说比 1/std::sqrt(x) 快 4 倍，而精度误差也不算大。这个方法最迷人的地方就是第 4 行的操作。对应的注释是一句被和谐化的语言。对于这行代码的操作，有一篇英文论文 InvSqrt.dvi (lomont.org) 专门研究了这行代码。
然而，当我以不写暑假作业为代价读完了这篇论文之后，回顾网上的中文解释，似乎都不尽人意。
所以就有了这篇文章。
1. 前置知识
为了不影响阅读体验，在阅读主菜之前您应该先了解以下知识：
1.1 牛顿迭代法求函数零点
对于连续曲线，如何求曲线与 x 轴的交点？
首先我们有个原始猜测 x_0 （initial guess）。然后作曲线在 x_0 处的切线。切线与 x 轴的交点作为更精准的猜测 x_1 。
然后再作 x_1 处的切线，与 x 轴交点为 x_2 。以此类推
比如说，下面这张图，蓝色是给定曲线，绿色是切线。
原始猜测 x_0=0.8
得到交点 x_1=0.6603 ，然后作下一条切线。
得到交点 x_2=0.6213 。
可以发现，x_0,x_1,x_2,\dots 这些交点越来越接近零点 0.618 。无论原始猜测是什么，牛顿迭代的次数越多，得到的结果就越接近零点。
一般地，对于曲线 y=f(x) ，给定初始猜测 x_0 ，使用牛顿迭代法得到第 n+1 次猜测时所做切线为：
y-f(x_n)=f'(x_n)(x-x_n)
代入 y=0 得到其与 x 轴的交点。
x_{n+1}=x_n-\dfrac{f(x_n)}{f'(x_n)}\quad(n\in N)
其中，N 为自然数集。上述式子即为牛顿迭代法的递推式。
1.2 浮点数的二进制表示
计算机是二进制的，所以计算机能表示的数也只能是关于 x=2 的多项式，对次数还有一定的要求。
比如，9=1\times 2^3+0\times 2^2+0\times 2^1+1\times 2^0 。
可以证明，所有的自然数都可以表示为形如 \sum_{i=0}^{+\infty}a_i2^i,a_i\in\{0,1\} 的多项式。当然，计算机的内存不可能去到正无穷，所以 unsigned int 类型仅能表示所有形如 \sum_{i=0}^{31}a_i2^i 的多项式（对于 signed 类型还需要考虑符号位，在这里不加以讨论）计算机中保存的就是多项式的系数 a_i 。
对于小数部分（即小数点后的部分），二进制同样只能表示形如 \sum_{i=1}^{+\infty}b_i2^{-i},b_i\in\{0,1\} 的多项式。
比如，0.6875=1\times2^{-1}+0\times2^{-2}+1\times2^{-3}+1\times2^{-4} 。
同样，由于计算机内存不可能去到正无穷，所以对于 i 的最大值也会有限制。一些小数不可能精确表示，比如 0.1 就不能表示为有限个 2 的若干次幂相加。如果设置 i 的最大值为 24 ，那么 0.1 的二进制表示将会是：
0.000110011001100110011001
转换为十进制就是：
0.09999996423721313
它非常接近 0.1 ，所以在实际处理中，计算机就会直接把他当做 0.1 输出。这是一种近似处理，在有效位数 7 位的情况下，它仍然可以被视为 0.1 。而当有效位数是 8 位时，他就是 0.09999996 了。所以计算机表示的小数会存在有效位数。
那么我们保存多项式的系数 b_i ？不行，对于一个数，可以同时拥有整数部分和小数部分。这样一来我们或许还要保存小数点的位置。如果钦定整数部分和小数部分的范围，会造成浪费：比如，我要表示 (0,1) 之间的小数，那么整数部分就非常鸡肋；或者，我要表示整数部分很大的一位小数，小数部分又会非常鸡肋。
为了协调需求，计算机中一般采用 IEEE 754-1985 表示浮点数。这是一种科学记数法。
当然，在一些场合还是会使用定点数，比如 ttf 字体文件中字形顶点的位置。
十进制中科学计数法大概是这样的：对于任意正实数 x ，可以表示为：
x=m\times 10^e,\quad m\in[1,10)
类似的，在二进制中：
x = m\times2^e,\quad m\in[1,2)
其中，m,e 都是实数，e 不表示自然对数的底数，而表示 2 的指数。
这样，我们只要保存一个符号位以及 m,e 就可以了。只要协调好 m,e 之间的内存分配问题。对于 32 位有符号浮点数 float，它是这样分配的：
首位为符号位，保存符号 s 。
接下来 8 位为指数位，保存与指数相关的数 E 。
最后 23 位为系数位，保存与系数相关的数 M 。
为了表示 e 为负数的情况，在保存时会给其加上 127 。所以 E=e+127 。
注意到 m 的整数部分一定是 1 ，可以不保存。所以 M=m-1 。
这里 E 一定是整数，但是 M 是实数，未必可以用有限位二进制来保存。所以这里保存的是其 23 位近似值。
总的来说，一个 float 类型变量 x 代表的数是：
x = (-1)^s\cdot(1+M)2^{E-127}
举个例子，13.5=(-1)^0\cdot(1+0.6875)2^{130-127} ，所以其 32 位二进制表示是这样的：
0 10000010 10110000000000000000000
2. 算法基础
给定 a ，要求 x=a^{-0.5} 的近似解。两边平方取倒数，即 \frac 1{x^2}=a 。也就是说，我们要求 f(x)=\dfrac1{x^2}-a 的零点。
根据牛顿迭代法的递推式，给定初始猜测 x_0 ，可知 x_1=x_0(1.5-0.5ax_0^2) 。
代码的大意是：通过第 4 行计算出一个比较准确的初始猜测 x_0 ，然后第 6 行用一次牛顿迭代修正。因为 x_0 已经比较精确了，所以最后的结果也会较为精确。
现在我们就要弄清楚，第 4 行究竟干了什么？
3. 神奇猜测
我们已经知道浮点数在计算机中是使用科学记数法表示的。但是，对于浮点数 13.5 ，如何得到其位数？其实很简单，只需要将其地址强制转换为 int 类型的地址，然后计算机就会按照 int 的规则去读取。
比如下面代码：
float x = 13.5f;
int i = *(int*)&x;然后输出看一下 i 的值是多少？应该输出的是 1096286208。因为其二进制表示就是 0 10000010 10110000000000000000000，恰好是 13.5 按照 IEEE 754-1985 法则储存的二进制序列。
所以，将其地址类型进行偷梁换柱就可以得到储存的二进制序列。
其中，有 8 位储存的是指数相关的数 E 。现在我们希望把指数减小一半并取相反数，那么一定会有的操作就是
指数右移一位
取反
其中，取反操作并不能很好地给出较为准确的初始猜测，于是我们选择用某个数 R 减去指数，这样也能起到类似取反的效果。
当然，我们不可能只对 E 进行操作。哪怕我们把二进制序列拿了出来存储在 int 里面，我们也需要对整个 int 进行操作。不妨把 R 切分成 3 部分：符号位，接下来 8 位记为 R_1 ，最后 23 位记为 R_2 。这样切分就和 IEEE 754-1985 法则对应上了，相减的时候就是 R_1 减 E ，R_2 减 M 。还有借位的情况，这些稍后再讨论。
为了使算法正常工作，我们必须时刻保持符号位为 0 （实数的开平方根操作不涉及负数），所以 R 的符号位也是 0 。
接下来就要分析，R_1 和 R_2 分别取什么值能使得误差最小。
3.1 指数部分
指数部分 E 先右移一位，这在整个序列中会出现一个问题：如果 E 是奇数，那么就会有一个 1 被右移到系数部分。简单起见我们先讨论偶数情况。
3.1.1 指数为奇数
当指数 e=2d+1 为奇数时，对应的 E 就会是偶数，这时右移对系数部分没有任何影响。
首先看一下右移再相减对指数部分的效果：
\begin{aligned}
E'&=R_1-\frac E2\\
&=R_1-\frac{e+127}{2}\\
&=R_1-\frac{2d+128}{2}\\
&=R_1-64-d
\end{aligned}
无论把上面的字母看做二进制序列还是一个具体的数，这个等式都是成立的。因为在对应法则的约束下，二进制序列和具体的数一一对应，他们的减法有同构关系。
请记住这个结果 必须是正数 。当它是负数时，二进制序列相减就会向符号位借位，然后你将会得到一个负数作为算术平方根，这很明显是不对的。当它是 0 时，就无法应对来自系数部分可能得借位。所以这个结果 必须是正数 。
偶数 E\in[0, 254] ，根据这个范围可以得出 d 的范围，从而 R_1 必须大于 127 ，即 R_1\ge 128 。
当 E 为偶数时，指数部分的右移对系数部分没有影响，所以系数部分正常运算，只需要考虑是否借位。
仍然先考虑简单的，若无需借位，则直接相减。
下面等式中，使用 \frac M2 来计算右移。因为右移舍弃的末位即使是 1 ，导致的误差在 IEEE 754-1985 法则中也只是 2^{-24} ，对于算法本身的误差来说可以忽略不计。
M'=R_2-\frac M2
那么我们将要得到的初始猜测 x_0 就会是：
\begin{aligned}
x_0&=(1+M')\times 2^{E'-127}\\
&=(1+R_2-\frac M2)\times 2^{R_1-64-d-127}\\
&=(2+2R_2-M)\times 2^{R_1-192-d}
\end{aligned}
如果需要借位，那么也会对指数部分产生一点影响（借走了 1 ）。这是系数就是
M'=1+R_2-\frac M2
初始猜测将会是：
\begin{aligned}
x_0&=(1+M')\times2^{E'-1-127}\\
&=(2+R_2-\frac M2)\times 2^{R_1-192-d}
\end{aligned}
对比一下两个式子，只有中间一点点不同。为了将它们写在一起，记
\beta_{M}^{R_2}=\begin{cases}2R_2-M&R_2\ge \frac M2\\R_2-M&R_2<\frac M2\end{cases}
则
x_0=(2+\beta_M^{R_2})\times 2^{R_1-192-d}
注意到 x_0 被表示成了一个关于 R_2 的分段函数（把 M 看作常量），而且函数图像连续，在每一段上可导，连接处不可导。说直接一点，就是因为 \beta_M^{R_2} 的图像是一条折线。
现在让我们来做误差分析：设真实值 a^{-0.5}=x ，因为 a 可以表示成 (1+M)\times 2^e ，那么对应的真实值 x 就应该是 \frac{1}{\sqrt{1+M}}\cdot2^{-0.5e} ，代入 e=2d+1 可以得到：
x = \frac{1}{\sqrt{1+M}}\cdot 2^{-d-\frac 12}=\frac{1}{\sqrt{2+2M}}\cdot2^{-d}
相对误差为 |\frac{x-x_0}{x}| ，代入展开即：
\begin{aligned}
&\left|\frac{\frac{1}{\sqrt{2+2M}}\cdot2^{-d} - (2+\beta_M^{R_2})\cdot2^{R_1-192-d}}{\frac{1}{\sqrt{2+2M}}\cdot 2^{-d}}\right|\\\\
=&|1-\sqrt{2+2M}(2+\beta_M^{R_2})2^{R_1-192}|
\end{aligned}
记为 |\varepsilon_0(M,R_1,R_2)| 。我们希望其最大值最小。若其最大值为 A_0 ，则必有 A_0 \ge |\varepsilon_0(0,R_1,R_2)| 。如果想要找到 R 使得 A_0 < \frac 18 ，则必有：
\begin{aligned}
\frac 18 &> A_0\\
&\geq|\varepsilon(0,R_1,R_2)|\\
&=\left|1-\sqrt2(2+2R_2)\cdot2^{R_1-192}\right|\\
&=|1-(1+R_2)\cdot2^{R_1-190.5}|
\end{aligned}
拆绝对值，即
-\frac 18<1-(1+R_2)\cdot2^{R_1-190.5}<\frac18
注意到 1\le1+R_2<2 ，可以得到
\frac 7{16}<2^{R_1-190.5}<\frac 98
取对数，估算出 189.2<R_1<190.7 。又因为 R_1 必须是整数，所以 R_1 = 190 就是唯一一个有机会把最大误差控制在 \frac 18 以内的数。
3.1.2 指数为偶数
这个时候，指数部分右移就会导致有一个 $1$ 右移到系数部分的最高位，对应地，就是给系数的实际值增加了 $2^{-1}=0.5$。这是后话。
让我们同样先从指数部分开始分析。
$$
\begin{aligned}
E'&=R_1-\left\lfloor\frac E2\right\rfloor\\
&=R_1-\left\lfloor\frac {e+127}2\right\rfloor\\
&=R_1-d-63
\end{aligned}
$$
然后计算系数部分，当不需要借位时：
$$
\begin{aligned}
x_0&=\left[1+R_2-\left(\frac M2+\frac 12\right)\right]\cdot 2^{R_1-63-d-127}\\
&=(2+4R_2-2M)\cdot 2^{R_1-192-d}
\end{aligned}
$$
当需要借位时：
$$
\begin{aligned}
x_0&=\left[1+1+R_2-\left(\frac M2 + \frac 12\right)\right]\cdot2^{R_1-63-d-1-127}\\
&=(3+2R_2-M)\cdot 2^{R_1-192-d}
\end{aligned}
$$
同样，记
$$
\gamma_M^{R_2}=\begin{cases}4R_2-2M&R_2\ge\frac M2+\frac12\\1+2R_2-M&R_2<\frac M2+\frac 12\end{cases}
$$
有
$$
x_0=(2+\gamma_m^{R_2})\cdot 2^{R_1-192-d}
$$
同样地，误差分析函数：
$$
\varepsilon_1(M,R_1,R_2)=1-\sqrt{1+M}(2+\gamma_M^{R_2})\cdot2^{R_1-192}
$$
相比于第一种情况少了一个 $\sqrt 2$。
同样设其最大值为 $A_1$，必有 $A_1\ge\varepsilon_1(1,R_1,R_2)$。虽说 $M\in[0,1)$ 不包含 $1$，但是函数连续且表达式在 $1$ 处有意义，为了不分类讨论就直接取 $M=1$ 了。
若最大误差小于 $\frac 18$，必有：
$$
\begin{aligned}
\frac 18 &> A_1\\
&\ge|\varepsilon_1(1,R_1,R_2)|\\
&=|1-\sqrt2(2+2R_2)\cdot 2^{R_1-192}|
\end{aligned}
$$
非常眼熟，因为和上面的一模一样！
所以得出结论：$R_1=190$ 是唯一一个有机会把最大误差控制在 $\frac 18$ 以内的指数部分的取值。
## 3.2 系数部分
现在有必要提醒您：接下来的内容涉及比较多暴力成分（基本没有技巧可言）。
既然 $R_1=190$，那么我们就可以重新定义误差函数了：
$$
\varepsilon_0(M,R_2)=1-\frac{\sqrt{2+2M}}4\cdot(2+\beta_M^{R_2})
$$
$$
\varepsilon_1(M,R_2)=1-\frac{\sqrt{1+M}}4\cdot(2+\gamma_M^{R_2})
$$
我们希望找到一个 $R_2$，使得 $\max\{\varepsilon_0(M,R_2),\varepsilon_1(M,R_2)\}$ 最小。
让我们先固定 $R_2$，然后枚举 $M$ 取什么值可能使得误差最大。
由于误差函数都是连续的、分段可导的，所以最值只会在段端点或极值点取到。
先看 $\epsilon_0$，最大值可能是下面的 5 种情况：
1. 端点 $M=0$。记此时最大值为 $f_1(R_2)=\varepsilon_0(0,R_2)=1-\dfrac{1+R_2}{\sqrt2}
端点 M=1 。记此时最大值为 f_2(R_2)=\varepsilon_0(1,R_2)=\begin{cases}\frac{1-2R_2}4&R<\frac 12\\\frac{1-2R_2}2&R\ge\frac 12\end{cases}
端点 M=2R_2 。记此时最大值为 f_3(R_2)=\varepsilon_0(2R_2,R_2)=1-\sqrt{\dfrac{1+2R_2}2}
M\le 2R_2$ 时的极值点。此时 $M=\frac 23R_2$。记最大值为 $f_4(R_2)=\varepsilon_0(\frac 23R_2,R_2)=1-\dfrac{(1+\frac 23R_2)^{\frac 32}}{\sqrt 2}
M > 2R_2$ 时的极值点。此时 $M=\frac 23(1+R_2)$。注意到 $M>2R_2$ 所以 $0\le R_2 < \frac 12$。记最大值为 $f_5(R_2) = \begin{cases}\varepsilon_0(\frac 23(1+R_2),R_2)=1-\dfrac{(5+2R_2)^{\frac 32}}{6\sqrt 6}&R_2 < \frac 12\\0&R_2\ge\frac 12\end{cases}
对于 \varepsilon_1 ，同样分为 5 种情况：
端点 3 种：
\begin{aligned}
&f_6(R_2)=\varepsilon_1(0, R_2)=\begin{cases}\frac{1-2R_2}4&R_2<\frac 12\\\frac 12-R_2&R_2 \ge \frac 12\end{cases}\\
&f_7(R_2)=\varepsilon_1(1, R_2)=1-\dfrac{1+R_2}{\sqrt 2}\\
&f_8(R_2)=\begin{cases}0&R_2<\frac 12\\\varepsilon_1(2R_2-1, R_2)=1-\sqrt{2R_2}&R_2\ge \frac 12\end{cases}
\end{aligned}
极值点 2 种：
于是乎合起来写：
f_9(R_2)=\begin{cases}
\varepsilon_1(\frac{2R_2+1}3,R_2)=1-\sqrt2(\frac{2+R_2}3)^{\frac 32}&R_2 < \frac 12\\
\varepsilon_1(\frac{2R_2-1}3, R_2)=1-[\frac 23(1+R_2)]^{\frac 32}&R_2\ge \frac 12
\end{cases}
现在，最大值所有可能得情况都分析完了。如果 R_2 能使得 \max\{\varepsilon_0(M, R_2),\varepsilon_1(M, R_2)\} 最小，那么 R_2 也必定使得 \max_{i=1}^9f_i(R_2) 最小。
这个时候让我们把问题交给 matplotlib。
可以发现，在 0.4 附近的时候最大值最小。
让我们放大一下。
可以发现，是在 0.43275 附近 f_3(R_2) 和 f_4(R_2) 的交点。
端点 M=2R_2 。记此时最大值为 f_3(R_2)=\varepsilon_0(2R_2,R_2)=1-\sqrt{\dfrac{1+2R_2}2}
M\le 2R_2$ 时的极值点。此时 $M=\frac 23R_2$。记最大值为 $f_4(R_2)=\varepsilon_0(\frac 23R_2,R_2)=1-\dfrac{(1+\frac 23R_2)^{\frac 32}}{\sqrt 2}
显然它们应该都是单调递减的。而黄色线却在上升，就是绝对值起了作用。去绝对值后应该有 f_3(R_2) + f_4(R_2) = 0 。展开，即
(1+\frac23R_2)^{\frac 32}+(1+2R_2)^{\frac 12}=2^{\frac 32}
所以 R_2 就是方程 (1+\frac 23x)^{\frac 32}+(1+2x)^{\frac 12}-2^{\frac 32}=0 在 [0,1] 内的根。
很好啊，三个根式，两次平方走不脱了。
左边函数是单调的，所以方程至多一个解。求解过程我们希望尽量精确，所以使用二分法（二分法可以使得近似解能用 IEEE 754-1985 法则表示）。为了减少误差我们希望方程中都是整数。
先直接平方
(1+\frac 23x)^3+(1+2x)+2(1+\frac23x)\sqrt{(1+\frac 23x)(1+2x)}=8
然后乘以 27 去分母，移项
(3+2x)^3+(27+54x)-216=-2(9+6x)\sqrt{12x^2+24x+9}
接着平方、展开、化简，最终得到
64x^6+576x^5+1296x^4-3456x^3-15552x^2-46656x+23328=0
（科技真好用）
然后对左边的式子轻轻在脑子里求导。导函数中，常数项 -46656 ，一次项系数 -31104 ，次数再高一点的对应系数都不过万，然后 x\in[0,1] 就导致它更小，所以整个导函数在 [0,1] 内都是负的，原函数单调递减。
所以就只有一个零点，取左端点为 0 ，右端点为 1 ，要求精度为 2^{-23} ，进行二分求解。
然而，如果在小数范围内求解仍然会面临精度问题，所以将等式两边同时乘以 (2^{23})^6 ，然后换元令 t=2^{23}x 。只要将 t 保留到整数即可，左端点为 0 ，右端点为 2^{23} 对下面方程进行二分。
64t^6M^0+576t^5M^1+1296t^4M^2-3456t^3M^3-15552t^2M^4-46656tM^5+23328M^6=0
其中，M=2^{23} 是一个常数。
很明显计算的过程需要用到高精度，所以我选择 python。
vec = [23328, -46656, -15552, -3456, 1296, 576, 64]
M = 2**23
def func(x):
res = 0
for i in range(len(vec)):
res += (vec[i] * M ** (6-i)) * x ** i
return res
l, r = 0, M
while l + 1 < r:
mid = (l+r) // 2
if func(mid) > 0:
l = mid
else:
r = mid
print(hex(l), func(l))
print(hex(r), func(r))输出结果是：
0x37642f 628423524624020249592480082211712497728
0x376430 -19274642075805850948756907632004956160000x37642f 的结果更接近 0 。事实上，这两个数都可以，他们对整个方法的误差不会太大。
在 IEEE 754-1985 规定的表示方法中，保存系数 R_2 的二进制序列，如果把这个序列看做一个整数，它其实就是 R_2 对应的二进制小数乘上 2^{23} 。因为他们都是二进制下的多项式，系数对应的指数都恰好相差了 23 ，所以这里就取 0x37642f 作为 R_2 的二进制序列。
然后前面提到 R_1 = 190 ，它被保存在符号位之后的 8 位。
所以使得初始猜测 i = R - (i >> 1); 误差最小的常数 R 在 IEEE 754-1985 法则下的二进制表示就是
0 10111110 01101110110010000101111四位一分就是：
0101 1111 0011 0111 0110 0100 0010 1111转成十六进制即为 0x5f37642f。
你可能会发现：这个值怎么和开头代码给的值不太一样？实际上，这个值得到的初始猜测确实更加接近真实值，但是在都进行一次牛顿迭代后，0x5f3759df 反而能够做到更优。
可以说，0x5f37642f 在初始猜测上更贴近，但是牛顿迭代对它不怎么友好。所以 Chris Lomont 在 0x5f37642f 附近进行一通暴力搜索后，找到了一个最“牛迭友好”的常数 0x5f375a86。
一般来说，这个方法都会进行两次牛顿迭代。因为双迭和单迭时间上相差极小，而精确度却提升了将近三个数量级。这是一笔亏小赚大的买卖。
（完）
4. 后记
翻译自开始时提到的论文。
这篇文章对代码中出现的神奇常数加以解释，并且尝试理解它以达到可能的优化。新常数 0x5f375a86 比原来的常数表现好一点点。然而这两者都是近似估计，所以两个常数都是可以的。如果可以的话，我希望找到原作者，询问一下原来的常数是算出来的还是枚举出来的。
这篇文章的功能是梳理这个方法的原理，为这个方法迁移到其他平台或其他数据类型提供理论基础。当然，任何理论都必须在实际中加以检验，因为硬件未必符合所有理论（有些平台可能不使用 IEEE 754-1985 法则）。