Skip to content

Latest commit

 

History

History
280 lines (201 loc) · 6.38 KB

2021-09-01-汇编中的除法.md

File metadata and controls

280 lines (201 loc) · 6.38 KB
layout title date categories tags comments mathjax copyrights recommend
post
汇编中的除法
2021-09-01 00:00:00 +0800
汇编
assembly
true
true
原创
true

猜猜我要讲什么

奇怪的数字

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

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&lt;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 时效果会更好。

这就是魔法！

参考资料