为什么 GCC 在实现整数除法时使用乘以一个奇怪的数字？
我一直在阅读有关 div 和 mul 汇编操作的文章，我决定通过用 C 编写一个简单的程序来了解它们的实际效果：
文件分割.c
#include <stdlib.h>
#include <stdio.h>
int main()
{
size_t i = 9;
size_t j = i / 5;
printf("%zu\n",j);
return 0;
}
然后生成汇编语言代码：
gcc -S division.c -O0 -masm=intel
但是查看生成的 division.s 文件，它不包含任何 div 操作！相反，它使用位移和幻数执行某种黑魔法。下面是计算 i/5 的代码片段：
mov rax, QWORD PTR [rbp-16] ; Move i (=9) to RAX
movabs rdx, -3689348814741910323 ; Move some magic number to RDX (?)
mul rdx ; Multiply 9 by magic number
mov rax, rdx ; Take only the upper 64 bits of the result
shr rax, 2 ; Shift these bits 2 places to the right (?)
mov QWORD PTR [rbp-8], rax ; Magically, RAX contains 9/5=1 now,
; so we can assign it to j
这里发生了什么？为什么 GCC 根本不使用 div？它是如何产生这个幻数的，为什么一切都有效？
-3689348814741910323 转换为 CCCCCCCCCCCCCCCD 作为 uint64_t 或大约 (2^64)*4/5。
-O0 处发出 div 指令。 （抄送@克利福德）
整数除法是您可以在现代处理器上执行的最慢的算术运算之一，延迟高达几十个周期，吞吐量很差。 （对于 x86，请参阅 Agner Fog's instruction tables and microarch guide）。
如果您提前知道除数，则可以通过将其替换为具有等效效果的一组其他操作（乘法、加法和移位）来避免除法。即使需要几个操作，它通常仍然比整数除法本身快得多。
以这种方式实现 C / 运算符而不是使用涉及 div 的多指令序列只是 GCC 进行常量除法的默认方式。它不需要跨操作进行优化，即使是调试也不会改变任何东西。 （不过，对小代码量使用 -Os 确实让 GCC 使用 div。）使用乘法逆而不是除法就像使用 lea 而不是 mul 和 add
因此，如果除数在编译时未知，您只会在输出中看到 div 或 idiv。
有关编译器如何生成这些序列的信息，以及让您自己生成它们的代码（几乎肯定没有必要，除非您使用的是 Braindead 编译器），请参阅 libdivide。
除以 5 与乘以 1/5 相同，这又与乘以 4/5 并右移 2 位相同。相关的值是十六进制的 CCCCCCCCCCCCCCCD，如果放在十六进制点之后，它是 4/5 的二进制表示（即五分之四的二进制是 0.110011001100 重复出现 - 原因见下文）。我想你可以从这里拿走它！您可能想查看 fixed point arithmetic（尽管请注意它在末尾四舍五入为整数）。
至于为什么，乘法比除法快，而当除数固定时，这是一条更快的路线。
有关其工作原理的详细说明，请参阅 Reciprocal Multiplication, a tutorial，并以定点形式进行解释。它显示了查找倒数的算法如何工作，以及如何处理有符号的除法和模数。
让我们考虑一下为什么 0.CCCCCCCC...（十六进制）或 0.110011001100... 二进制是 4/5。将二进制表示除以 4（右移 2 位），我们将得到 0.001100110011...，通过简单的检查可以将原来的值相加得到 0.111111111111...，它显然等于 1，同样的方式 0.9999999... 在十进制等于一。因此，我们知道 x + x/4 = 1，所以 5x/4 = 1，x=4/5。然后将其表示为十六进制的 CCCCCCCCCCCCD 以进行舍入（因为存在的最后一位之外的二进制数字将是 1）。
一般来说，乘法比除法快得多。因此，如果我们可以避免乘以倒数，我们可以显着加快除以常数
一个问题是我们不能精确地表示倒数（除非除法是 2 的幂，但在这种情况下，我们通常可以将除法转换为位移）。因此，为了确保正确答案，我们必须小心，倒数中的错误不会导致最终结果中的错误。
-3689348814741910323 是 0xCCCCCCCCCCCCCCCCCD，它是用 0.64 定点表示的刚好超过 4/5 的值。
当我们将 64 位整数乘以 0.64 定点数时，我们得到 64.64 结果。我们将值截断为 64 位整数（实际上将其舍入到零），然后执行进一步的移位，除以 4 并再次截断通过查看位级别，很明显我们可以将两个截断视为单个截断。
这显然给了我们至少一个除以 5 的近似值，但它是否给了我们一个正确四舍五入到零的准确答案？
为了得到准确的答案，误差需要足够小，以免将答案推到舍入边界上。
除以 5 的确切答案总是小数部分为 0、1/5、2/5、3/5 或 4/5。因此，相乘和移位结果中小于 1/5 的正误差永远不会将结果推到舍入边界上。
我们常量的误差是 (1/5) * 2-64。 i 的值小于 264，因此相乘后的误差小于 1/5。除以 4 后，误差小于 (1/5) * 2−2。
(1/5) * 2−2 < 1/5 所以答案总是等于进行精确除法并四舍五入到零。
不幸的是，这不适用于所有除数。
如果我们尝试将 4/7 表示为 0.64 定点数并从零四舍五入，我们最终会得到 (6/7) * 2-64 的错误。乘以略低于 264 的 i 值后，我们最终得到略低于 6/7 的误差，除以 4 后，我们最终得到略低于 1.5/7 的误差，大于 1/7。
因此，要正确实现除以 7，我们需要乘以 0.65 定点数。我们可以通过乘以定点数的低 64 位，然后加上原始数（这可能会溢出进位）然后通过进位循环来实现。
这是一个算法文档的链接，该算法生成我在 Visual Studio 中看到的值和代码（在大多数情况下），我假设它仍然在 GCC 中用于将变量整数除以常量整数。
http://gmplib.org/~tege/divcnst-pldi94.pdf
文章中，一个uword有N位，一个udword有2N位，n=分子=被除数，d=分母=除数，ℓ初始设置为ceil(log2(d))，shpre为预移位（在乘法前使用) = e = d 中尾随零位的数量，shpost 是移位后（在乘法之后使用），prec 是精度 = N - e = N - shpre。目标是使用预移位、乘法和后移位来优化 n/d 的计算。
向下滚动到图 6.2，它定义了如何生成 udword 乘数（最大大小为 N+1 位），但没有清楚地解释该过程。我将在下面解释这一点。
图 4.2 和图 6.2 显示了对于大多数除数，乘法器如何减少到 N 位或更少的乘法器。公式 4.5 解释了用于处理图 4.1 和 4.2 中的 N+1 位乘法器的公式是如何得出的。
在现代 X86 和其他处理器的情况下，乘法时间是固定的，因此预移位对这些处理器没有帮助，但它仍然有助于将乘法器从 N+1 位减少到 N 位。我不知道 GCC 或 Visual Studio 是否已经消除了 X86 目标的预移位。
回到图 6.2。只有当分母（除数）> 2^(N-1)（当ℓ == N => mlow = 2^(2N)）时，mlow 和 mhigh 的分子（除数）才可以大于 udword，在这种情况下n/d 的优化替换是比较（如果 n>=d，q = 1，否则 q = 0），因此不会生成乘数。 mlow 和 mhigh 的初始值将是 N+1 位，并且可以使用两个 udword/uword 除法来产生每个 N+1 位值（mlow 或 mhigh）。以 64 位模式下的 X86 为例：
; upper 8 bytes of dividend = 2^(ℓ) = (upper part of 2^(N+ℓ))
; lower 8 bytes of dividend for mlow = 0
; lower 8 bytes of dividend for mhigh = 2^(N+ℓ-prec) = 2^(ℓ+shpre) = 2^(ℓ+e)
dividend dq 2 dup(?) ;16 byte dividend
divisor dq 1 dup(?) ; 8 byte divisor
; ...
mov rcx,divisor
mov rdx,0
mov rax,dividend+8 ;upper 8 bytes of dividend
div rcx ;after div, rax == 1
mov rax,dividend ;lower 8 bytes of dividend
div rcx
mov rdx,1 ;rdx:rax = N+1 bit value = 65 bit value
您可以使用 GCC 进行测试。您已经了解了 j = i/5 的处理方式。看看 j = i/7 是如何处理的（应该是 N+1 位乘法器的情况）。
在大多数当前处理器上，乘法具有固定的时序，因此不需要预移位。对于 X86，最终结果是大多数除数的两个指令序列，以及像 7 这样的除数的五个指令序列（为了模拟 N+1 位乘法器，如等式 4.5 和 pdf 文件的图 4.2 所示）。示例 X86-64 代码：
; rbx = dividend, rax = 64 bit (or less) multiplier, rcx = post shift count
; two instruction sequence for most divisors:
mul rbx ;rdx = upper 64 bits of product
shr rdx,cl ;rdx = quotient
;
; five instruction sequence for divisors like 7
; to emulate 65 bit multiplier (rbx = lower 64 bits of multiplier)
mul rbx ;rdx = upper 64 bits of product
sub rbx,rdx ;rbx -= rdx
shr rbx,1 ;rbx >>= 1
add rdx,rbx ;rdx = upper 64 bits of corrected product
shr rdx,cl ;rdx = quotient
; ...
为了解释 5 指令序列，一个简单的 3 指令序列可能会溢出。让 u64() 表示高 64 位（商所需的全部）
mul rbx ;rdx = u64(dvnd*mplr)
add rdx,rbx ;rdx = u64(dvnd*(2^64 + mplr)), could overflow
shr rdx,cl
为了处理这种情况，cl = post_shift-1。 rax = 乘数 - 2^64，rbx = 股息。 u64() 是高 64 位。请注意，rax = rax<<1 - rax。商数是：
u64( ( rbx * (2^64 + rax) )>>(cl+1) )
u64( ( rbx * (2^64 + rax<<1 - rax) )>>(cl+1) )
u64( ( (rbx * 2^64) + (rbx * rax)<<1 - (rbx * rax) )>>(cl+1) )
u64( ( (rbx * 2^64) - (rbx * rax) + (rbx * rax)<<1 )>>(cl+1) )
u64( ( ((rbx * 2^64) - (rbx * rax))>>1) + (rbx*rax) )>>(cl ) )
mul rbx ; (rbx*rax)
sub rbx,rdx ; (rbx*2^64)-(rbx*rax)
shr rbx,1 ;( (rbx*2^64)-(rbx*rax))>>1
add rdx,rbx ;( ((rbx*2^64)-(rbx*rax))>>1)+(rbx*rax)
shr rdx,cl ;((((rbx*2^64)-(rbx*rax))>>1)+(rbx*rax))>>cl
我会从稍微不同的角度回答：因为是允许的。
C 和 C++ 是针对抽象机器定义的。编译器按照 as-if 规则将该程序从抽象机器转换为具体机器。
只要不改变抽象机指定的可观察行为，编译器就可以进行任何更改。没有合理的期望编译器会以最直接的方式转换你的代码（即使很多 C 程序员都这么认为）。通常，它这样做是因为编译器希望与直接方法相比优化性能（如其他答案中详细讨论的那样）。
如果在任何情况下编译器将正确的程序“优化”为具有不同可观察行为的东西，那就是编译器错误。
我们代码中的任何未定义行为（有符号整数溢出是一个经典示例）并且此合同无效。
关注公众号
不定期副业成功案例分享
-O3也是如此。编译器必须编写为所有可能的输入值提供正确结果的代码。这仅适用于具有-ffast-math的浮点数，并且 AFAIK 没有“危险”整数优化。 （启用优化后，编译器可能能够证明可能的值范围，例如，它可以使用仅适用于非负符号整数的东西。）-O0（但不是-Os），默认情况下也会启用模乘逆。其他编译器（如 clang）将 DIV 用于-O0处的非 2 次幂常量。相关：我想我在 my Collatz-conjecture hand-written asm answer 中包含了一段关于此的内容