Ch3nyang's blog

home

home

person

about

hive

project

collections_bookmark

mindclip

rss_feed

rss

汇编中的除法

calendar_month 2021-08
archive 汇编
tag assembly

猜猜我要讲什么

奇怪的数字

我们有一个很简单的函数,将输入的参数除以三并返回:

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 这一步便可以写成

\[\begin{array}{lll}q&=&(1431655766n)>>32\\ &=&\left\lfloor\frac{2^{32}+2}{3}\frac{n}{2^{32}}\right\rfloor\\ &=&\left\lfloor\frac{n}{3}+\frac{2n}{3\cdot2^{32}}\right\rfloor\\ \end{array}\]

当 \(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 时效果会更好。

这就是魔法!

参考资料

Comments

Share This Post