旭穹の陋室

与其他算术运算相比，除法在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的幂次时，我们可以用一个二进制移位操作来替换它，这个操作只需要一周期。二进制最大公约数算法就是这种技术的一个典型示例。

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

然后在运行时将计算：

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

Barrett Reduction

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

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

可以证明，总是存在这样的一对数字mmm和sss，编译器在遇到像x / y这样的常数除法的时候，会自动将其优化为 (x * m) >> s 这样的操作，其中 y 是常数，m 是根据 y 计算得出的魔数，s 是对应的二进制移位个数。 以下是一个例子，展示编译器如何将一个无符号长整数unsigned long long除以 (109+7)(10^9+7)(109+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”是因为它主要用于模运算，模运算可以通过下面这个公式进行替换，从而只需要一个除法、乘法和减法操作：

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

为什么行得通

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

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

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

让我们从mmm的边界开始：

现在，对于整个表达式：

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

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

在这种情况下，下边界为(2n−22n−1−1)(\frac{2^n-2}{2^n-1} - 1)(2n−12n−2​−1)，上边界为(2n−22n−1+2n−22s)(\frac{2^n-2}{2^n-1} + \frac{2^n-2}{2^s})(2n−12n−2​+2s2n−2​)。左边界非常接近一个整数（0），整个范围的长度1+x2s1 + \frac{x}{2^s}1+2sx​也是第二长的（xxx为第二大整数）。重点来了：如果s≥ns≥ns≥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"。

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

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

• 为了获取整数部分（29），我们可以对原数进行向下取整或在小数点前截断。

• 为了获取分数部分(56\frac{5}{6}65​)，我们只需要取小数点后的部分。

• 为了获取余数（5），我们可以将分数部分乘以除数。

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

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

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

因为我们这里所做的可以被看作三个浮点乘法的链式运算，总的相对误差为O(ϵ)O(\epsilon)O(ϵ)。由于ϵ=O(12s)\epsilon = O(\frac{1}{2^s})ϵ=O(2s1​)，s=2ns = 2ns=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;
}

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

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

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

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

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

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