魔术数字

概述

cpu做除法div比起其他指令耗费大量的性能，如下表

    Operands  µops   µops each port   Latency    Reciprocal throughput
DIV   r64      36    p0 p1 p5 p6      35-88      21-83
MUL   r64      2     p1 p6            3          1
SHR   r,i      1     p06              1          0.5

(本表引用自 https://agner.org/optimize/instruction_tables.pdf 的Intel Coffee Lake章节)

用乘法取代除法

乘法和其它简单操作的效率远高于一个div，而魔术数字(magic number)的功能就是用乘法取代除法，假设除数固定，需要事前算出一个魔术数

m = \lfloor \frac{2^n + d - 1}{d} \rfloor

其中：

$d$ ：除数
$n$ ：次方数， $2^n$ 需要 $\ge$ 被除数

然后事先计算出常数m，比如 n = 64 要除以10的话 m = 1844674407370955162
接着在程式用被除数乘以m，再右移n位数

q = (a \cdot m) \gg n

其中：

$q$ ：商
$a$ ：被除数
$m$ ：上面算出的常数
$n$ ：上面使用的次方数
$\gg n$ ：数学上除以 $2^n$ ，但>>右移是上面的shr非常快的指令

实际上程式做的运算就只有mul和shr得到商
但是需要右移64位的表示商刚好在mul后的rdx
所以mul直接取rdx就是除法得到的商

取得余数

前言

以下是经由rax取得余数的方式，若用原本方法将商rdx乘以除数再用被除数减掉更快的话请用原本方法

r_{\text{(余数)}} = a_{\text{(被除数)}} - \text{rdx}_{\text{(商)}} \cdot d_{\text{(除数)}}

2560 ÷ 10

我们观察一个示例

$a$ ：2560
$d$ ：10
$n$ ：64

事先计算出常数 m = 1844674407370955162
然后直接将被除数 $a$ 乘以 m

\text{hex}(a \cdot m) = 0\text{x} \overbrace{100}^{\text{rdx}} \overbrace{0000000000000400}^{\text{rax}}

以十六进制来看，原本要右移8 bytes等于后面16位数都是之前盖掉的
前面三位数0x100是商256，那后面0x400是1024，也就是商数乘以4
我们将乘数4以 $\delta$ 代称， $\delta$ 的来源是

\delta = (- 2^n) \bmod d

也就是 $(- 2^{64}) \bmod 10 = 4$ 这个同样可以在事前计算出来

2562 ÷ 10

将 2562 乘以 m 得到

\text{hex}(a \cdot m) = 0\text{x} \overbrace{100}^{\text{rdx}} \overbrace{3333333333333734}^{\text{rax}}

我们关注低64位rax的部分，这是 $m \cdot 2 + q \cdot \delta$ 的结果，当没有余数的时候只有商数 q 乘以 $\delta$ (4) 的0x400，现在加上了两份m，余数 r = 2，将算式移向得到

r = \frac{\text{rax} - \delta \cdot q}{m}

这边又到了div问题，我们将除数带入最上面的常数m的算式会得到m刚好等于原本的除数10
也就是这次 m = 除数 d，将 rax 做完上面分子的运算之后乘以除数 d 再右移 64 位就是余数 r

r = ((\text{rax} - \delta \cdot q) \cdot d) \gg 64

2_{\text{new rdx}} = (\text{rax} - 4 \cdot \text{rdx}) \cdot 10

这里的 $\delta$ 如果是2的次方数请用最快的shl左移 shl rdx, 2

65535 ÷ 2

经过事前运算得到 m = 9223372036854775808 (0x8000000000000000)
```
(2**64+d-1)//d  # python代碼
```
除数 65535 乘以 m 得到 $\text{hex}(65535 \cdot m) = 0\text{x} \overbrace{\text{7fff}}^{\text{rdx}} \overbrace{8000000000000000}^{\text{rax}}$
得到商 rdx = 0x7fff = 32767
事前计算 $\delta$ $δ$ 得到 0
```
-(2**64) % 2
# 0
```
$(\text{rax} \cdot d) \gg 64 = (\text{rax} \ll 1) \gg 64 = \text{rax} \gg 63 = 1$ 得到余数 1

示例代码

global _start

section .text
_start:
        ; 2562 ÷ 10
        mov     rax, 2562  ; 被除数
        mov     rbx, 1844674407370955162  ; magic number
        mul     rbx
        ; rdx = 商
        shl     rdx, 2     ; 事前计算得到δ是4
        sub     rax, rdx
        mov     rbx, 10
        mul     rbx
        ; rdx = 余数