c++ - 为什么在x86上除以3需要右移(以及其他奇数)？

Question

我具有以下C/C++函数:

unsigned div3(unsigned x) {
    return x / 3;
}

When compiled using clang 10的 -O3，结果为:

div3(unsigned int):
        mov     ecx, edi         # tmp = x
        mov     eax, 2863311531  # result = 3^-1
        imul    rax, rcx         # result *= tmp
        shr     rax, 33          # result >>= 33
        ret

我的理解是:除以3等于乘以乘积逆3-1 mod 232，即2863311531。
有些事情我还是不明白:

为什么我们需要完全使用ecx/rcx？我们不能直接将rax与edi相乘吗？

为什么我们要在64位模式下进行乘法运算？将eax和ecx相乘不是更快吗？

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

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

编辑1
对于那些不理解3-1 mod 232我的意思的人，我在这里谈论的是乘法逆。
例如:

// multiplying with inverse of 3:
15 * 2863311531      = 42949672965
42949672965 mod 2^32 = 5

// using fixed-point multiplication
15 * 2863311531      = 42949672965
42949672965 >> 33    = 5

// simply dividing by 3
15 / 3               = 5

因此，与42949672965乘以实际上等于除以3。我假设clang的优化实际上是基于模块化算法的，而实际上它是基于定点算法的。
编辑2
我现在已经知道，乘法逆仅可用于除法运算，而无余数。例如，将3-1乘以1等于3-1，而不是零。只有定点算法才具有正确的舍入。
不幸的是，即使在可能的情况下，clang也不使用模块化算术，在这种情况下，模块化算术只是一个 imul指令。以下函数具有与上面相同的编译输出。

unsigned div3(unsigned x) {
    __builtin_assume(x % 3 == 0);
    return x / 3;
}

(关于精确除法的定点乘法逆的规范问答，该逆适用于每种可能的输入: Why does GCC use multiplication by a strange number in implementing integer division?-不太重复，因为它仅涵盖数学，而不包括某些实现细节，如寄存器宽度和imul vs. mul。)

Answer 1

1. Can't we multiply rax with edi directly?

我们不能使用 imul rax, rdi，因为调用约定允许调用者在RDI的高位上留下垃圾。仅EDI部分包含该值。内联时这不是问题；编写32位寄存器确实会将零扩展到完整的64位寄存器，因此编译器通常不需要额外的指令即可对32位值进行零扩展。
(如果不能避免，最好使用 limitations on mov-elimination零扩展到另一个寄存器中)。
从字面上看，甚至没有问题，x86没有任何乘法指令对它们的输入之一进行零扩展以使您将32位和64位寄存器相乘。两个输入的宽度必须相同。

2. Why do we multiply in 64-bit mode?

(术语:所有这些代码都在64位模式下运行。您在问为什么64位操作数大小如此。)
您可以使用 mul edi将EAX与EDI相乘以在EDX:EAX上获得64位结果，但是在Intel CPU上 mul edi是3 uops，而在大多数现代x86-64 CPU上具有快速64位 imul。 (尽管 imul r64, r64在AMD Bulldozer系列以及某些低功耗CPU上速度较慢。) https://uops.info/和 https://agner.org/optimize/(指令表和Microarch PDF)
(有趣的事实: mul rdi实际上在Intel CPU上更便宜，只有2 oups。也许不必对整数乘法单元的输出进行额外的拆分，例如 mul edi将不得不拆分64位低半倍数乘法器输出分成EDX和EAX一半，但对于64x64 => 128位多的像素自然会发生这种情况。)
另外，您需要的部分在EDX中，因此您需要另一个 mov eax, edx来处理它。 (同样，因为我们正在寻找的是该函数的独立定义的代码，而不是在内联到调用方之后。)
GCC 8.3和更早版本的确使用32位 mul而不是64位 imul( https://godbolt.org/z/5qj7d5)。当Bulldozer系列和旧的Silvermont CPU更加相关时，对于 -mtune=generic来说并不疯狂，但是对于最近的GCC而言，这些CPU在过去更遥远，其通用调整选择反射(reflect)了这一点。不幸的是，GCC还浪费了 mov指令将EDI复制到EAX，使这种方式看起来更糟:/

# 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

mov eax, 0xAAAAAAAB/ mul edi只能是6 oups，但仍然比:

# 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位 0x00000000AAAAAAAB不能表示为32位符号扩展的立即数，因此 imul rax, rcx, 0xAAAAAAAB无法编码。这将意味着 0xFFFFFFFFAAAAAAAB。

3. Why are we using imul instead of mul? I thought modular arithmetic would be all unsigned.

它是未签名的。输入的符号仅影响结果的上半部分，但是 imul reg, reg不会产生结果的上半部分。只有operate形式的 mul和 imul是NxN => 2N的完全乘法，因此只有它们需要单独的有符号和无符号版本。
只有 imul才具有更快，更灵活的低半值形式。关于 imul reg, reg签署的唯一一件事是，它基于下半部分的签署溢出设置OF。仅拥有一个 mul r,r与 imul r,r唯一的区别是FLAGS输出是不值得花费更多的操作码和更多的晶体管的。
英特尔手册( https://www.felixcloutier.com/x86/imul)甚至指出了它可以用于无符号的事实。

4. What's up with the 33-bit rightshift at the end? I thought we can just drop the highest 32-bits.

不，如果以这种方式实现，则没有乘数常量可以为每个可能的输入x提供正确的正确答案。 “按原样”优化规则不允许近似，仅允许对程序使用的每个输入产生完全相同的可观察行为的实现。如果不知道 x的值范围而不是整个 unsigned的值范围，则编译器没有该选项。 ( -ffast-math仅适用于浮点；如果需要更快的整数数学近似值，请按如下所示手动进行编码):
请参阅 Why does GCC use multiplication by a strange number in implementing integer division?，以获取有关编译器用于通过编译时间常数进行精确除法的定点乘法逆方法的更多信息。
有关此示例在一般情况下不起作用的示例，请参见我对 Divide by 10 using bit shifts?答案的修改，其中建议

// 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);
}

当 div10(1073741829) = 107374183实际上是107374182时，它的第一个错误答案(如果从0向上循环)是 1073741829/10。(它应四舍五入，而不是像C整数除法那样四舍五入。)

从您的编辑中，我看到您实际上是在谈论使用乘法结果的下半部分，显然，该结果对于从UINT_MAX一直到精确的倍数都非常适用。
如您所说，当除法有余数时，例如截断为32位而不是 16 * 0xaaaaaaab时， 0xaaaaaab0 = 5。

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

是的，如果该数学方法可行，则编译器使用32位imul实现该方法是合法且最佳的。他们不寻求这种优化，因为这鲜为人知。如果值得在编译时间方面增加编译器代码甚至寻找优化，则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;
}