魔術數字

概述

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 = 餘數