整数除法

与其他算术运算相比，除法在x86等普通计算机上的执行效率非常差。在硬件层面上实现整数和浮点数除法是非常困难的。这涉及到算术逻辑单元（ALU）的设计，除法运算电路会占据大量的空间，同时计算过程也需要很多步骤，因此导致完成div除法运算通常需要10-20个计算周期。对于处理小数据类型的除法运算，其延迟稍微少一些。

x86上的除法和取模

没有人愿意为了独立的模数运算复制所有的这些乱七八糟的东西，所以div指令既用于除法运算也用于取模运算。要执行32位的整数除法运算，你需要特定地将被除数（dividend）放在eax寄存器中，并调用div指令，除数（divisor）作为唯一操作数。这样操作后，商数（quotient）将被存储在eax寄存器中，余数（remainder）将被存储在edx寄存器中。

需要注意的是，被除数实际上需要被存储在eax和edx两个寄存器中：这种机制使得64位除以32位或者甚至128位除以64位的除法运算成为可能，这类似于128位乘法的工作方式。当执行常规的32位除以32位有符号的除法运算时，我们需要将eax进行符号扩展到64位，并存储其高位部分在edx中。

div(int, int):
    mov  eax, edi
    cdq
    idiv esi
    ret

对于无符号除法，您可以将edx设置为零，这样它就不会干扰运算：

div(unsigned, unsigned):
    mov  eax, edi
    xor  edx, edx
    div  esi
    ret

在这两种情况下，除了eax寄存器中的商，我们还可以获得edx中的余数：

mod(unsigned, unsigned):
    mov  eax, edi
    xor  edx, edx 
    div  esi
    mov  eax, edx
    ret

你还可以将128位整数（存储在rdx:rax中）除以64位整数：

div(u128, u64):
    ; a = rdi + rsi, b = rdx
    mov  rcx, rdx
    mov  rax, rdi
    mov  rdx, rsi
    div  edx 
    ret

被除数的高部分必须小于除数，否则会发生溢出。由于这种限制，编译器很难自动生成上述代码。编译器一般会进行额外的检查来避免这个问题，但这种检测可能实际上是不必要的。

常数除法

整数除法的运算速度非常慢，即使在全硬件实现的情况下也是如此。但是，如果除数是一个常数，我们在某些情况下可以避免使用除法运算。一个众所周知的例子是当除以2的幂次时，我们可以用一个二进制移位操作来替换它，这个操作只需要一周期。二进制最大公约数算法就是这种技术的一个典型示例。

更一般的情况下，有一些巧妙的技巧可以将除法运算替换为乘法运算，但需要产生一些预先计算出的数值。所有这些技巧都基于下述想法：考虑这样一个问题，我们要将一个浮点数 $x$ 除以另一个已知的浮点数 $y$ 。我们可以做的是提前计算一个常数，这样我们就可以用乘法来替代除法了：

$d \approx y^{-1}$

然后在运行时将计算：

$x/y = x \cdot y^{-1} \approx x \cdot d$

对于 $\frac{1}{y}$ 这种运算，最大的误差是ϵ。而对于 $x⋅d$ 这种运算，因为 $d$ 本身可能已经含有误差 $d(1+\epsilon)$ ，所以整个运算的误差将是这两种误差的叠加，即 $\frac{x \cdot (d(1+\epsilon)) \cdot (1+\epsilon) - x \cdot d}{x \cdot d} = 2ϵ+ϵ^2$ 。这里当 $ϵ$ 比较小的时候， $ϵ^2$ 将更小，所以可以忽略不计，于是总的误差控制在 $O(ϵ)$ 以内。这在大多数浮点运算的场景下是可以接受的。

Barrett Reduction

如何将这个技巧推广到整数呢？似乎计算int d = 1/y并不可行，因为这会得到0。我们尽力能做的就是：

$d = \frac{m}{2^s}$

然后找一个魔数 $m$ 和一个二进制移位个数 $s$ ，使得对所有范围内的x有x / y == (x * m) >> s 。

$⌊x / y⌋ = ⌊x \cdot y^{-1}⌋ = ⌊x \cdot d⌋ = ⌊x \cdot \frac{m}{2^s}⌋$

可以证明，总是存在这样的一对数字 $m$ 和 $s$ ，编译器在遇到像x / y这样的常数除法的时候，会自动将其优化为 (x * m) >> s 这样的操作，其中 y 是常数，m 是根据 y 计算得出的魔数，s 是对应的二进制移位个数。以下是一个例子，展示编译器如何将一个无符号长整数unsigned long long除以 $(10^9+7)$ 的操作，转换成乘法和位移操作的底层汇编代码：

;  input (rdi): x
; output (rax): x mod (m=1e9+7)
mov    rax, rdi
movabs rdx, -8543223828751151131  ; load magic constant into a register
mul    rdx                        ; perform multiplication
mov    rax, rdx
shr    rax, 29                    ; binary shift of the result

... ; more operation to do mod

Barrett reduction 是一种用于优化模运算（取余数）的技术。这种技术称为“reduction”是因为它主要用于模运算，模运算可以通过下面这个公式进行替换，从而只需要一个除法、乘法和减法操作：

$r = x - ⌊x / y⌋ \cdot y$

这种方法需要一些预先计算，包括实际执行一次除法。因此，这种方法只有在执行不只一次除法，而且所有的除数都是同一个固定的常数时，才会有所益处。

为什么行得通

我们目前不是很清楚为什么总是存在这样的 $m$ 和 $s$ ，更不用说如何找到它们了。但是如果给定 $s$ ，直觉告诉我们 $m$ 应当接近 $2^s/y$ 以便能够约掉 $2^s$ 。所有有两个自然而然的选择： $⌊2^s/y⌋$ 和 $⌈2^s/y⌉$ 。

第一个选择是不可行的，因为你将其带入

$⌊\frac{x\cdot⌊2^s/y⌋}{2^s}⌋$

可以看到对于任何可以整除的 $\frac{x}{y}$ ，若 $y$ 不为偶数，则上述结果一定会小于真实结果。这就只剩下另一种情况， $m = ⌈2^s/y⌉$ 。现在，让我们尝试导出计算结果的下限和上限：

$⌊x/y⌋ = ⌊\frac{x \cdot m}{2^s}⌋ = ⌊\frac{x \cdot ⌈2^s/y⌉}{2^s}⌋$

让我们从 $m$ 的边界开始：

$2^s/y ≤ ⌈2^s/y⌉ < 2^s/y + 1$

现在，对于整个表达式：

$x/y - 1 < ⌊\frac{x \cdot ⌈2^s/y⌉}{2^s}⌋ < x/y + x/s^s$

我们可以结果落在一个长度为 $1 + \frac{x}{2^s}$ 的范围内，如果对所有的 $x/y$ ，这个范围内总是有一个整数，那么该算法则能保证输出正确的答案。事实证明，我们总是可以将 $s$ 设置的足够大，使得上述要求得到满足。

这里最糟糕的情况是什么？怎么挑选 $x$ 和 $y$ 的值，使得 $(x/y - 1,x/y + x/s^s)$ 范围内有两个整数？我们可以看到， $x$ 和 $y$ 的比率为整数其实没关系，因为左边为开区间，如果假设 $x/s^s < 1$ ，则只有 $x/y$ 自己在这个范围内。最糟糕的情况其实是 $x/y$ 接近1但是不超过1。对于n位整数，这种情况可能是由第二大整数除以第一大整数导致的，即：

$x = 2^n - 2 \\ y = 2^n - 1$

在这种情况下，下边界为 $(\frac{2^n-2}{2^n-1} - 1)$ ，上边界为 $(\frac{2^n-2}{2^n-1} + \frac{2^n-2}{2^s})$ 。左边界非常接近一个整数（0），整个范围的长度 $1 + \frac{x}{2^s}$ 也是第二长的（ $x$ 为第二大整数）。重点来了：如果 $s≥n$ ，则范围内的整数只有1，算法总会返回该值（不是很理解，左边界理论上是小于0，则范围内应该还有0啊，除非作者认为左边界就是0）。

Lemire Reduction

Barrett reduction有点复杂，而且由于它是间接计算的，所以进行模运算时会生成一长串的指令序列。2019年 Daniel Lemire提出的一种新方法——“Lemire reduction”。这种方法消除了Barrett reduction的一些缺点。首先，Lemire算法比Barrett reduction简单。其次，在一些情况下，Lemire reduction实际上比Barrett reduction速度更快。尽管这种新的方法还没有一个广为人知的名称，作者仍决定暂时将其命名为"Lemire reduction"。

下面是该算法的核心思想。考虑某个整数分数的浮点表示：

$\frac{179}{6} = 11101.1101010101... = 29\frac{5}{6} \approx 29.83$

我们应当如何“解剖”它，以获得我们需要的“零部件”呢？

为了获取整数部分（29），我们可以对原数进行向下取整或在小数点前截断。
为了获取分数部分( $\frac{5}{6}$ )，我们只需要取小数点后的部分。
为了获取余数（5），我们可以将分数部分乘以除数。

对于32位整数，我们可以设定 $s$ 值为64，然后查看我们在乘法和移位方案中所做的计算：

$⌊x/y⌋ = ⌊\frac{x \cdot m}{2^s}⌋ = ⌊\frac{x \cdot ⌈2^s/y⌉}{2^s}⌋$

我们在这里实际做的是首先将整数 $x$ 乘以一个浮点数常量 $(x \cdot m)$ ，然后将得到的结果进行截断 $(⌊\frac{\cdot}{2^s}⌋)$ 。

如果我们取低比特位而不是高比特位呢？这将对应于分数部分——如果我们进一步将其乘上 $y$ 并且截断结果，那我们将得到余数：

$r = ⌊\frac{(x \cdot ⌈2^s/y⌉ \mod 2^s) \cdot y}{2^s}⌋$

因为我们这里所做的可以被看作三个浮点乘法的链式运算，总的相对误差为 $O(\epsilon)$ 。由于 $\epsilon = O(\frac{1}{2^s})$ ， $s = 2n$ ，所以误差总是小于1，因此结果将是准确的。

uint32_t y;

uint64_t m = uint64_t(-1) / y + 1; // ceil(2^64 / y)

uint32_t mod(uint32_t x) {
    uint64_t lowbits = m * x;
    return ((__uint128_t) lowbits * y) >> 64; 
}

uint32_t div(uint32_t x) {
    return ((__uint128_t) m * x) >> 64;
}

我们还可以只使用一次乘法来判断 $x$ 是否能被 $y$ 整除，这利用了当且仅当分数部分（ $m \cdot x$ 的低64位）小于 $m$ 时，余数为0（否则当分数部分乘上 $y$ 并右移64位后，结果将不为0）。

bool is_divisible(uint32_t x) {
    return m * x < m;
}

这个方法唯一的的缺点是在进行乘法运算时，它需要的整数类型是原来的四倍大小，而其他的reduction方法只需要双倍大小就可以。

此外，还有一种通过精心地分割和处理中间结果的来执行 64位整型mod运算的方法；这种算法的实现作为一个练习留给读者。

拓展阅读

可以查看 libdivide 和 GMP ，以了解优化整数除法的更通用实现。

同样值得一读的是 Hacker’s Delight ，它有一整章专门讨论整数除法。