整数除法在除数为常量时的优化原理

5/9/2024

Compiler

搬运自我的知乎回答：编译器是如何实现 32 位整型的常量整数除法优化的？

题主说的应该是将除以编译期常量的除法变换为一次乘法和一堆加减法和按位右移的优化吧：

unsigned test(unsigned x) {
  return x / 10;
}

被编译器优化为：

test:
  mov   eax, edi
  mov   edx, 3435973837
  imul  rax, rdx
  shr   rax, 35
  ret

这是一个很经典的优化，在 1994 年的论文 Division by Invariant Integers using Multiplication 中提出了明确的构造方法及其正确性证明。

现代 CPU 内部的硬件除法器通常延迟很高。以 x86 架构为例，uops.info 上给出了 Alder Lake 架构下 32 位整数除法指令 idiv 的延迟为 10-15 个周期左右，相比之下 32 位整数乘法指令 imul 的延迟为 3-4 个周期。因此，将除法操作替换为若干其他算术运算操作是有效的优化。

无符号整数除法

记被除数为，除数为，且是一个常数。我们试图计算。

这个优化的基本原理其实非常简单。我们知道，在计算整数除法时，如果除数是 2 的整数次幂，那么这个除法操作可以被一个非常高效的按位右移操作取代。因此我们可以考虑把除数换为 2 的某个整数次幂，并将原来的除法操作替换为一次乘法和一个按位右移操作：

进一步，这里虽然不是一个整数，但是由于最后我们会把乘法的结果按位右移位得到原除法操作的结果，因此我们可以尝试用一个与之足够接近的整数近似替代使得不精确的乘法只会影响乘法结果的低位：

因此，只要能找到这样一组使得上式中的约等号成为等号即可完成这一优化。

假设被除数和除数都是位无符号整数，其数值范围为。前述论文中的定理 4.2 证明了，对于整数和，如果有：

那么一定有：

证明也不难，只需要证明

即可。感兴趣的读者可以在论文中找到具体证明细节。

有了这个结论就很好构造了。根据条件 (1)，我们需要在区间内找到一个的倍数。这个区间的长度为，显然为了使区间内一定存在这样的倍数，我们可以令。在确保区间内一定存在的倍数后，我们就可以很容易构造出来。在构造出来和后，原无符号整数除法便可以被变换为一次无符号整数乘法操作和一次按位逻辑右移操作。

从优化原理也不难看出，对于同一个优化前的除法操作，可能存在多种组合满足优化条件，因此不同的编译器可能会将同一个除法操作转换为不同的优化代码。本文接下来只介绍其中的几个具体优化实现方法。

上面的理论优化方案在实际应用过程中会遇到问题。我们有：

且当为 2 的整数次幂时取等号。这意味着计算可能会发生溢出。当是一个很大的无符号位整数时，的计算结果可能会超出位无符号整数的范围。

因此，编译器一般会首先尝试令并依次确认 (1) 中的区间内是否存在的倍数。如果存在，那么相应的可以被塞入一个位无符号整数中参与乘法运算。此时的结果不会超出位无符号整数的范围，能够直接通过 CPU 乘法指令进行计算。

实例：对于以下函数 test1：

unsigned test1(unsigned n) {
  return n / 10;
}

目前最新的 gcc 和 clang 编译器都把它优化为了：

unsigned test1(unsigned n) {
  // mov   ecx, $n
  // mov   eax, 3435973837
  // imul  rax, rcx
  // shr   rax, 35
  // ret
  return (uint64_t)n * 3435973837 >> 35;
}

这里，根据上面的优化方法，当时 (1) 中的区间为，这个区间内刚好存在一个 10 的倍数 34359738370 。因此，，分别对应于优化后的代码中的两个魔数。

如果对于所有都没有在 (1) 中的区间内找到的倍数，但是是一个偶数，那么我们还可以进一步尝试将原除法约分。记且是一个奇数，那么有：

此时和就成为新的被除数和除数，他们都是位无符号整数。并且由于除数变小了，我们也就更有可能可以在 (1) 中的区间内找到其倍数了。我们可以再按照上述的流程，依次令并分别尝试为在区间内寻找倍数。如果成功找到了倍数，就可以计算相应的并进行优化变换。

实例：对于以下函数 test2：

unsigned test2(unsigned n) {
  return n / 42;
}

目前最新的 clang 编译器把它优化为了：

unsigned test2(unsigned n) {
  // shr   $n
  // imul  rax, rdi, 818089009
  // shr   rax, 34
  // ret
  return (uint64_t)(n >> 1) * 818089009 >> 34;
}

当时 (1) 中的区间内均不存在 42 的倍数。但是 42 是一个偶数，因此 n / 42 可以被变换为 (n >> 1) / 21 ，因此优化后的代码有一个将 x 逻辑右移 1 位的操作。新的除数为，约分后的被除数和除数都是位无符号整数。当时，区间内刚好有一个的倍数 17179869189 。因此有，，分别对应于优化后代码中的另外两个魔数。与前一种情况相比，这种情况会多出一个按位逻辑右移操作。

在尝试了上面所有办法后，如果 (1) 中的区间内仍然不存在的倍数，我们就只能令。此时有：

为了解决溢出的问题，我们令，显然，并且有：

这样就避免了乘法计算过程中的溢出。

记，那么有，因此可以由一个位无符号整数表示。在按照上面的优化方案计算时，虽然乘法操作不会溢出，但的结果仍然可能溢出位无符号整数的上界。考虑到，我们可以进一步优化上面的计算方法消除这样的溢出情况：

其中这一部分是计算两个整数的平均值，利用经典的无溢出整数平均值计算算法可计算为：

实例：对于以下函数 test3：

unsigned test3(unsigned n) {
  return n / 31;
}

目前最新的 clang 编译器把它优化为了：

unsigned test3(unsigned n) {
  // mov   eax, $n
  // imul  rax, rax, 138547333
  // shr   rax, 32
  // sub   $n, eax
  // shr   $n
  // add   eax, $n
  // shr   eax, 4
  // ret
  unsigned k = (uint64_t)n * 138547333 >> 32;
  return (((n - k) >> 1) + k) >> 4;
}

对于所有，区间内都不存在 31 的倍数，且 31 并不是偶数，因此只能令。此时有，，即计算 k 时用到的魔数。k值计算得到后再直接套用 (3) 并带入就可以得到原除法运算结果。

有符号整数除法

针对有符号整数除以常数的优化方案可以基于无符号整数的优化方案进行针对性调整。假设被除数为，除数为，且是一个常数。我们希望计算，其中的含义是取的整数部分（即朝 0 的方向把取整到最近的整数）。和都是位有符号整数，他们的数值范围为且。

从前述论文中的定理 5.1 可以引申出定理 (2) 在有符号整数上的情况。如果有：

那么一定有：

其中若成立，的值为 1，否则其值为 0 。也就是说，对于有符号整数的除法，可以首先按照无符号整数除法的优化方法进行优化，随后判断除法结果是否为负数，如果结果为负则将除法结果加一即可。

实例：对于以下函数：

int test(int n) {
  return n / -10;
}

目前最新的 clang 编译器把它优化为了：

int test(int n) {
  // movsxd rax, $n
  // imul   rax, rax, -1717986919
  // mov    rcx, rax
  // shr    rcx, 63
  // sar    rax, 34
  // add    eax, ecx
  // ret
  int64_t k = (int64_t)n * -1717986919;
  return (k >> 34) + (k < 0);
}

优化后代码的最后一步就是判断 k 是否为负并在其为负时将结果加 1 。对于优化后代码中的两个魔数，则是完全按照无符号除法的优化方法计算得到的。我们尝试分别令并依次判断区间内是否存在的倍数。当时，区间内恰好有一个 -10 的倍数 17179869190 。因此令，，分别对应于优化后的代码中的两个魔数。

与无符号整数除法类似，如果当时区间内均不存在的倍数，那么便只能令，。此时有。为处理乘法溢出问题，当时，效仿无符号整数除法的处理情况，可以将优化后的操作变换如下：

其中表示算术右移操作。当时也可以类似地得出相应的计算方法。

实例：对于以下函数：

int test(int n) {
  return n / 15;
}

目前最新的 clang 编译器把它优化为了：

int test(int n) {
  // movsxd rax, $n
  // imul   rcx, rax, -2004318071
  // shr    rcx, 32
  // add    eax, ecx
  // mov    ecx, eax
  // shr    ecx, 31
  // sar    eax, 3
  // add    eax, ecx
  // ret
  int t1 = (uint64_t)((int64_t)n * -2004318071) >> 32;
  int t2 = n + t1;
  return (t2 >> 3) + (n < 0);
}

同样，优化后代码的最后一步加法就是判断除法结果是否为负并在其为负时将结果加 1 。由于当时区间内均不存在 15 的倍数，因此只能令，因此需要将 t2 右移位。此时，，对应于优化后代码中的乘法操作的魔数操作数。