猜猜我要讲什么
奇怪的数字
我们有一个很简单的函数，将输入的参数除以三并返回：
int div(int num) {
return num / 3;
}
我们将其编译为汇编：
div(int):
movsx rax, edi
imul rax, rax, 1431655766
shr rax, 32
sar edi, 31
sub eax, edi
ret
-
movslq %edi, %rax：函数接收的第一个参数需要存储在 edi 中。这条指令将其扩展为 64 bit 然后移动到 rax 中； -
imulq $1431655766, %rax, %rax：rax = rax * 1431655766。这步中两个 32 bit 的数相乘，存储在 64 bit 寄存器中，不会产生溢出； -
shrq $32, %rax：rax 逻辑右移了 32 位； -
sarl $31, %edi：edi 算数右移了 32 位，也就是判断了正负。如果为正，则为全零，否则为全一； -
subl %edi, %eax：eax = eax - edi； -
ret：函数返回。
上面的过程相当于下面的函数：
int div(int v){
int magic = 1431655766;
int result = ((long long)v * (long long)magic) >> 32;
if (v < 0) result += 1;
return result;
}
？？？
这里根本没有涉及到除法，怎么就完成了除三的操作呢？这个奇怪的 1431655766 又是什么呢？
证明
对于整数除法来说，由于返回值也是整数，故需要做截断。具体的为：
\[n\div d=\left\{\begin{array}{ll}\lfloor n/d\rfloor&nd\geq0\\\lceil n/d\rceil&nd<0\end{array}\right.\]那么这个奇怪的 1431655766 呢？
\[1431655766=\frac{2^{32}+2}{3}\]那么 result = ((long long)v * (long long)magic) >> 32 这一步便可以写成
当 \(n\geq0\) 时，由于
\[0\leq\frac{2n}{3\cdot2^{32}}<\frac{1}{3}\]故有
\[q=\left\lfloor\frac{n}{3}\right\rfloor\]接下来证明 \(n<0\) 的情况。
\[\begin{array}{lll}q&=&\left\lfloor\frac{2^{32}+2}{3}\frac{n}{2^{32}}\right\rfloor+1\\ &=&\left\lfloor\frac{2^{32}n+2n+3\cdot2^{32}}{3\cdot2^{32}}\right\rfloor\\ &=&\left\lfloor\frac{2^{32}n+2n+1}{3\cdot2^{32}}\right\rfloor\\ &=&\left\lfloor\frac{n}{3}+\frac{2n+1}{3\cdot2^{32}}\right\rfloor\\ \end{array}\]由于
\[-\frac{1}{3}+\frac{1}{3\cdot2^{32}}\leq\frac{2n+1}{3\cdot2^{32}}\leq\frac{1}{3\cdot2^{32}}\]故有
\[q=\left\lceil\frac{n}{3}\right\rceil\]证毕。
根据以上过程，我们也很容易得到，对于除以 \(d\)，我们只需要把那个神奇的数字改为 \(\frac{2^{32}+2}{d}\)。
比较
我们来比较下面两段汇编：
; div.asm
section .text
global div31
global div32
div31:
mov eax, edi
movsx rdx, edi
shr rdx, 32
mov ecx, 3
idiv ecx
ret
div32:
movsx rax, edi
imul rax, 1431655766
shr rax, 32
sar edi, 31
sub eax, edi
ret
其中 div31 为正常的汇编除法写法，div32 则是利用了我们这个数字。
我们测试一下这两种写法的差别。
div31: 1675 div32: 1954
div31: 1995 div32: 1565
div31: 1754 div32: 1581
div31: 1913 div32: 1818
div31: 1776 div32: 1900
div31: 1634 div32: 1849
div31: 1954 div32: 1958
div31: 1652 div32: 1765
div31: 1703 div32: 2160
div31: 2540 div32: 1573
呃呃……好像没啥差别……速度不相上下……
为什么呢？我也不知道……反正做 gcc 的那帮人肯定比我聪明。
雷神之锤
谈到高效运算中的 magic number，肯定绕不过 Quake III 的浮点开方：
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;
}
？？？？？？
这又是什么？
我们先挑简单的看。
i = * ( long * ) &y; 把 y 看作整数赋值给 i；
然后 i = 0x5f3759df - ( i >> 1 ); 一通操作；
y = * ( float * ) &i; 把 i 再赋值给 y。
最后一句 y = y * ( threehalfs - ( x2 * y * y ) ); 实际上是牛顿法。
对于函数
\[f(y)=\frac{1}{y^2}-x=0\]根据牛顿法的求解方法，有
\[y_{n+1}=y_{n}-\frac{\frac{1}{y_n^2}-x}{-\frac{2}{y_n^3}}=y_n\left(\frac{3}{2}-\frac{x}{2}y_n^2\right)\]注释掉的那句呢？则是第二次牛顿法迭代，可以让结果更加准确。
现在，就剩了中间这句 i = 0x5f3759df - ( i >> 1 ); 实在费解。根据上面牛顿法的思路，这个地方应该是估计出了 \(\frac{1}{\sqrt{y}}\) 的大概值。这是怎么做到的？
对于一个 float，它可以被表示为
\[F=\pm2^{((E)_{10}-127)}\times\left(1+\frac{(M)_{10}}{2^{23}}\right)\]而当我们把它转换成整数时，显然，它变成了
\[I=(E)_{10}\times2^{23}+(M)_{10}\]我们想要解
\[y=x^{-\frac{1}{2}}\]可以等价为
\[\log_2(y)=-\frac{1}{2}\log_2(x)\]带入上面浮点数的式子，可以得到
\[(E_y-127)+\log_2\left(1+\frac{M_y}{2^{23}}\right)=-\frac{1}{2}(E_x-127)-\frac{1}{2}\log_2\left(1+\frac{M_x}{2^{23}}\right)\]然后，利用对数函数的一个神秘的近似
\[(E_y-127)+\left(k+\frac{M_y}{2^{23}}\right)=-\frac{1}{2}(E_x-127)-\frac{1}{2}\left(k+\frac{M_x}{2^{23}}\right)\]其中，\(k\) 是一个需要试的数字。
移相得到
\[E_y\times2^{23}+M_y=\frac{3}{2}(127-k)2^{23}-\frac{1}{2}(E_x\times2^{23}+M_x)\]这里可以一眼看出来了
\[0x5f3759df=\frac{3}{2}(127-k)2^{23}\]作者在这里取了 \(k\approx0.045047\)
事实上，根据后人计算，取 0x5f375a86 时效果会更好。
这就是魔法！
参考资料
- Lomont, Chris. “Fast inverse square root.” Tech-315 nical Report 32 (2003).
- https://en.wikipedia.org/wiki/Fast_inverse_square_root
- https://cjting.me/2021/03/16/the-missing-div-instruction-part1/
- https://www.zhihu.com/question/26287650
Comments