为什么除以 3 需要在 x86 上右移（和其他奇怪之处）？

最后的 33 位右移是怎么回事？我以为我们可以去掉最高的 32 位。

相反，您必须更多地考虑 before 在上 32 位中的位置和 the 在下 32 位的位置。 这个定点运算工作正常，但结果明显偏移到上半部分，因此 CPU 必须在运算后再次将结果向下移。3^(-1) mod 30.33333330.3333rax

为什么我们使用 imul 而不是 mul？我以为模块化算术都是无符号的。

没有等同于该指令的指令。使用的变体采用两个寄存器：MULIMULIMUL

a <= a * b

没有这样做的指令。 指令更昂贵，因为它们将结果存储为两个寄存器中的 128 位。 当然，您可以使用传统指令，但这并不能改变结果存储在两个寄存器中的事实。MULMUL

1. 我们不能直接将 rax 与 edi 相乘吗？

我们不能，因为调用约定允许调用方在 RDI 的高位中留下垃圾;只有 EDI 部分包含该值。内联时这不是问题;写入 32 位寄存器会隐式地对整个 64 位寄存器进行零扩展，因此编译器通常不需要额外的指令来对 32 位值进行零扩展。imul rax, rdi

（如果无法避免，则将零扩展到不同的寄存器会更好，因为对移动消除的限制）。

从字面上理解你的问题，不，x86 没有任何乘法指令可以零扩展它们的一个输入，让你乘以 32 位和 64 位寄存器。两个输入的宽度必须相同。

2. 为什么我们要以 64 位模式乘法？

（术语：所有这些代码都以 64 位模式运行。你问为什么是 64 位操作数大小。

您可以将 EAX 与 EDI 相乘，以获得跨 EDX：EAX 的 64 位结果，但在 Intel CPU 上是 3 uops，而大多数现代 x86-64 CPU 具有快速的 64 位。（虽然在 AMD Bulldozer 系列和一些低功耗 CPU 上速度较慢。https://uops.info/ 和 https://agner.org/optimize/（说明书和微架构 PDF） （有趣的事实：实际上在英特尔 CPU 上更便宜，只有 2 个 uops。也许与不必对整数乘法单元的输出进行额外拆分有关，例如必须将 64 位低半乘法器输出拆分为 EDX 和 EAX 两半，但这自然会发生在 64x64 => 128 位 mul 中。mul edimul ediimulimul r64, r64mul rdimul edi

此外，您想要的部分在 EDX 中，因此您需要另一个来处理它。（同样，因为我们查看的是函数的独立定义的代码，而不是在内联到调用方之后。mov eax, edx

GCC 8.3 及更早版本确实使用了 32 位而不是 64 位 （https://godbolt.org/z/5qj7d5）。当 Bulldozer 系列和旧的 Silvermont CPU 更相关时，这并不疯狂，但这些 CPU 对于最近的 GCC 来说已经更久远了，其通用的调整选择反映了这一点。不幸的是，GCC 还浪费了将 EDI 复制到 EAX 的指令，使这种方式看起来更糟：/mulimul-mtune=genericmov

# gcc8.3 -O3  (default -mtune=generic)
div3(unsigned int):
        mov     eax, edi                 # 1 uop, stupid wasted instruction
        mov     edx, -1431655765         # 1 uop  (same 32-bit constant, just printed differently)
        mul     edx                      # 3 uops on Sandybridge-family
        mov     eax, edx                 # 1 uop
        shr     eax                      # 1 uop
        ret
                                  # total of 7 uops on SnB-family

只有 6 个 uops 带 / ，但仍然比：mov eax, 0xAAAAAAABmul edi

# gcc9.3 -O3  (default -mtune=generic)
div3(unsigned int):
        mov     eax, edi                # 1 uop
        mov     edi, 2863311531         # 1 uop
        imul    rax, rdi                # 1 uop
        shr     rax, 33                 # 1 uop
        ret
                      # total 4 uops, not counting ret

遗憾的是，64 位不能表示为 32 位符号扩展的 immediate，因此不可编码。这意味着.0x00000000AAAAAAABimul rax, rcx, 0xAAAAAAAB0xFFFFFFFFAAAAAAAB

3. 为什么我们使用 imul 而不是 mul？我以为模块化算术都是无符号的。

它是未签名的。输入的符号仅影响结果的高半部分，但不会产生高半部分。只有 和 的单操作数形式是 NxN => 2N 的全乘法，因此只有它们需要单独的有符号和无符号版本。imul reg, regmulimul

只有更快、更灵活的低半形式。唯一有符号的问题是它根据低半部分的有符号溢出设置 OF。花费更多的操作码和更多的晶体管只是为了拥有与FLAGS输出的唯一区别是不值得的。imulimul reg, regmul r,rimul r,r

英特尔的手册（https://www.felixcloutier.com/x86/imul）甚至指出，它可以用于未签名的事实。

4. 最后的 33 位右移是怎么回事？我以为我们可以去掉最高的 32 位。

不，如果你以这种方式实现它，就没有乘数常数可以为每个可能的输入 x 给出确切正确的答案。“假设”优化规则不允许近似值，只允许为程序使用的每个输入生成完全相同的可观察行为的实现。如果不知道 full range 以外的值范围，编译器就没有该选项。（仅适用于浮点;如果您想要更快的整数数学近似值，请手动编码它们，如下所示）：xunsigned-ffast-math

请参阅为什么 GCC 在实现整数除法时使用乘法乘以奇数？有关编译器用于按编译时常量进行精确除法的定点乘法反方法的详细信息。

有关这在一般情况下不起作用的示例，请参阅我对使用位移除 10 的答案的编辑？

// Warning: INEXACT FOR LARGE INPUTS
// this fast approximation can just use the high half,
// so on 32-bit machines it avoids one shift instruction vs. exact division
int32_t div10(int32_t dividend)
{
    int64_t invDivisor = 0x1999999A;
    return (int32_t) ((invDivisor * dividend) >> 32);
}

它的第一个错误答案（如果你从 0 向上循环）是 when 实际上是107374182。（它向上舍入，而不是像 C 整数除法那样向 0 方向发展。div10(1073741829) = 1073741831073741829/10

---

从您的编辑中，我看到您实际上是在谈论使用乘法结果的下半部分，这显然非常适合精确倍数一直到 UINT_MAX。

正如你所说，当除法有余数时，它完全失败了，例如 = 截断为 32 位时，而不是 .16 * 0xaaaaaaab0xaaaaaab05

unsigned div3_exact_only(unsigned x) {
    __builtin_assume(x % 3 == 0);  // or an equivalent with if() __builtin_unreachable()
    return x / 3;
}

是的，如果该数学方法成功，编译器使用 32 位 imol 实现该数学将是合法且最佳的。他们不寻找这种优化，因为它很少是已知的事实。IDK 是否值得添加编译器代码来寻找优化，就编译时间而言，更不用说开发人员时间的编译器维护成本了。这在运行时成本上没有太大的差异，而且很少可能。不过，这很好。

div3_exact_only:
    imul  eax, edi, 0xAAAAAAAB        # 1 uop, 3c latency
    ret

但是，您可以在源代码中自己完成此操作，至少对于已知的类型宽度，例如：uint32_t

uint32_t div3_exact_only(uint32_t x) {
    return x * 0xaaaaaaabU;
}

如果你看看我对上一个问题的回答：

为什么 GCC 在实现整数除法时使用奇数乘法？

它包含一个指向 pdf 文章的链接，该文章对此进行了解释（我的答案澄清了这篇 pdf 文章中没有很好地解释的内容）：

https://gmplib.org/~tege/divcnst-pldi94.pdf

请注意，某些除数（例如 7）需要额外的精度，乘法器通常需要 33 位，乘积通常需要 65 位，但这可以通过单独处理 2^32 位和 3 个附加指令来避免，如我之前的回答和下面的答案所示。

如果更改为

unsigned div7(unsigned x) {
    return x / 7;
}

因此，为了解释该过程，让 L = ceil（log2（divisor））。对于上面的问题，L = ceil（log2（3）） == 2。右移位计数最初为 32+L = 34。

为了生成具有足够位数的乘法器，需要生成两个潜在的乘法器：mhi 将是要使用的乘法器，移位计数将为 32+L。

mhi = (2^(32+L) + 2^(L))/3 = 5726623062
mlo = (2^(32+L)        )/3 = 5726623061

然后检查是否可以减少所需的位数：

while((L > 0) && ((mhi>>1) > (mlo>>1))){
    mhi = mhi>>1;
    mlo = mlo>>1;
    L   = L-1;
}
if(mhi >= 2^32){
    mhi = mhi-2^32
    L   = L-1;
    ; use 3 additional instructions for missing 2^32 bit
}
... mhi>>1 = 5726623062>>1 = 2863311531
... mlo>>1 = 5726623061>>1 = 2863311530  (mhi>>1) > (mlo>>1)
... mhi    = mhi>>1 = 2863311531
... mlo    = mhi>>1 = 2863311530
... L = L-1 = 1
... the next loop exits since now (mhi>>1) == (mlo>>1)

所以乘数是 mhi = 2863311531 和移位计数 = 32+L = 33。

在现代 X86 上，乘法和移位指令是恒定时间，因此没有必要将乘法器 （mhi） 减少到小于 32 位，因此将上面的 while（...） 更改为 if（...）。

在 7 的情况下，循环在第一次迭代时退出，并且需要 3 条额外的指令来处理 2^32 位，因此 mhi 为 <= 32 位：

L = ceil(log2(7)) = 3
mhi = (2^(32+L) + 2^(L))/7 = 4908534053
mhi = mhi-2^32 = 613566757

设 ecx = dividend，简单的方法可能会在 add 上溢出：

mov eax, 613566757             ; eax = mhi
mul ecx                        ; edx:eax = ecx*mhi
add edx, ecx                   ; edx:eax = ecx*(mhi + 2^32), potential overflow
shr edx, 3

为避免潜在的溢出，请注意 eax = eax*2 - eax：

(ecx*eax)           =   (ecx*eax)<<1)             -(ecx*eax)
(ecx*(eax+2^32))    =   (ecx*eax)<<1)+  (ecx*2^32)-(ecx*eax)
(ecx*(eax+2^32))>>3 =  ((ecx*eax)<<1)+  (ecx*2^32)-(ecx*eax)     )>>3
                    = (((ecx*eax)   )+(((ecx*2^32)-(ecx*eax))>>1))>>2

所以实际的代码，使用 u32（） 表示上 32 位：

...                 visual studio generated code for div7, dividend is ecx
mov eax, 613566757
mul ecx                        ; edx = u32( (ecx*eax) )
sub ecx, edx                   ; ecx = u32(            ((ecx*2^32)-(ecx*eax))        )
shr ecx, 1                     ; ecx = u32(           (((ecx*2^32)-(ecx*eax))>>1)    )
lea eax, DWORD PTR [edx+ecx]   ; eax = u32( (ecx*eax)+(((ecx*2^32)-(ecx*eax))>>1)    )
shr eax, 2                     ; eax = u32(((ecx*eax)+(((ecx*2^32)-(ecx*eax))>>1))>>2)

---

如果需要余数，则可以使用以下步骤：

mhi and L are generated based on divisor during compile time
...
quotient  = (x*mhi)>>(32+L)
product   = quotient*divisor
remainder = x - product

x/3 约为 （x * （2^32/3）） / 2^32。因此，我们可以执行单个 32x32->64 位乘法，取更高的 32 位，得到大约 x/3。

有一些错误，因为我们不能精确地乘以 2^32/3，只能乘以这个四舍五入为整数的数字。我们使用 x/3 ≈ （x * （2^33/3）） / 2^33 获得更高的精度。（我们不能使用 2^34/3，因为那是 2^32 >）。事实证明，这足以在所有情况下得到 x/3。如果输入为 3k 或 3k+2，则通过检查公式是否给出 k 的结果来证明这一点。