为什么在 x86 上除以 3 需要右移(和其他奇怪的东西)?

Jan*_*tke 32 c++ assembly compilation x86-64 integer-division

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

unsigned div3(unsigned x) {
    return x / 3;
}
Run Code Online (Sandbox Code Playgroud)

使用 clang 10 at编译-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
Run Code Online (Sandbox Code Playgroud)

我所理解的是:除以 3 相当于乘以乘法逆 3 -1 mod 2 32,即 2863311531。

不过还是有些不明白的地方:

  1. 为什么我们需要使用ecx/rcx呢?我们不能乘raxedi直接?
  2. 为什么我们在 64 位模式下乘法?乘以eaxand不是更快ecx吗?
  3. 为什么我们使用imul而不是mul?我认为模算术都是无符号的。
  4. 最后的 33 位右移是怎么回事?我以为我们可以放弃最高的 32 位。

编辑 1

对于那些不明白我所说的 3 -1 mod 2 32是什么意思的人,我在这里谈论的是乘法逆。例如:

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
Run Code Online (Sandbox Code Playgroud)

所以乘以 42949672965 实际上相当于除以 3。我认为 clang 的优化是基于模算术的,而实际上它是基于定点算术的。

编辑 2

我现在意识到乘法逆只能用于没有余数的除法。例如,1 乘以 3 -1等于 3 -1,而不是零。只有定点算术具有正确的舍入。

不幸的是,clang 没有使用任何模算术,imul在这种情况下它只是一条指令,即使它可以。以下函数具有与上述相同的编译输出。

// 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
Run Code Online (Sandbox Code Playgroud)

(关于适用于每个可能输入的精确除法的定点乘法逆的规范问答:为什么 GCC 在实现整数除法时使用一个奇怪的数乘法? - 不完全重复,因为它只涵盖数学,而不是一些实现详细信息,例如寄存器宽度和 imul 与 mul。)

Pet*_*des 29

  1. 不能直接用 edi 乘以 rax 吗?

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

(由于mov-elimination限制,如果无法避免,则零扩展到不同的寄存器会更好)。

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

  1. 为什么我们在 64 位模式下乘法?

(术语:所有这些代码都在 64 位模式下运行。你问为什么 64 位操作数大小。)

可以 mul ediEAX与 EDI相乘以获得跨 EDX:EAX 拆分的 64 位结果,但mul edi在 Intel CPU 上为 3 uop,而大多数现代 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 uop。也许与不必对整数乘法单元的输出进行额外拆分有关,例如mul edi必须将 64 位低半乘法器输出拆分为 EDX 和 EAX 两半,但是对于 64x64 => 128 位 mul,这自然发生。)

你想要的部分也在 EDX 中,所以你需要另一个mov eax, edx来处理它。(同样,因为我们正在查看函数的独立定义的代码,而不是在内联到调用者之后。)

GCC 8.3 及更早版本确实使用 32 位mul而不是 64 位imulhttps://godbolt.org/z/5qj7d5)。-mtune=generic当推土机系列和旧的 Silvermont CPU 更相关时,这并不疯狂,但这些 CPU 过去更远用于最近的 GCC,其通用调整选择反映了这一点。不幸的是,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
Run Code Online (Sandbox Code Playgroud)

使用mov eax, 0xAAAAAAAB/只会是 6 uop mul 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
Run Code Online (Sandbox Code Playgroud)

不幸的是,64 位0x00000000AAAAAAAB不能表示为 32 位符号扩展立即数,因此imul rax, rcx, 0xAAAAAAAB不可编码。这将意味着0xFFFFFFFFAAAAAAAB.

  1. 为什么我们使用 imul 而不是 mul?我认为模算术都是无符号的。

它是未签名的。输入的符号性仅影响结果的高半部分,但imul reg, reg不会产生高半部分。只有一个操作数的形式mulimul充满乘法指令做的N×N => 2N,所以他们只需要单独的符号和无符号版本。

Onlyimul有更快更灵活的低半单形式。唯一有符号的imul reg, reg是它根据低半部分的有符号溢出设置 OF。仅仅为了与 FLAGS 输出mul r,r的唯一区别而花费更多的操作码和更多的晶体管是不值得的imul r,r

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

  1. 最后的 33 位右移是怎么回事?我以为我们可以放弃最高的 32 位。

不,x如果您以这种方式实现,则没有乘数常数可以为每个可能的输入提供准确的正确答案。 “as-if”优化规则不允许近似,只允许为程序使用的每个输入产生完全相同的可观察行为的实现。如果不知道x除全范围之外的值范围unsigned,编译器没有该选项。(-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);
}
Run Code Online (Sandbox Code Playgroud)

它的第一个错误的答案(如果你从0向上循环)是div10(1073741829) = 1073741831073741829/10实际是107374182.(这围捕,而不是朝着0像C整数除法应该。)


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

正如您所说,当除法有余数时它完全失败,例如16 * 0xaaaaaaab=0xaaaaaab0当截断为 32 位时,而不是5.

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

是的,如果这个数学计算成立,编译器用 32 位 imul 实现它是合法和最佳的。他们不寻找这种优化,因为它很少是一个已知的事实。IDK是否值得添加编译器代码甚至寻找优化,就编译时间而言,更不用说开发人员时间的编译器维护成本。这不是一个巨大的运行成本差异,而且它很少会是可能的。不过,这很好。

div3_exact_only:
    imul  eax, edi, 0xAAAAAAAB        # 1 uop, 3c latency
    ret
Run Code Online (Sandbox Code Playgroud)

但是,您可以在源代码中自己做一些事情,至少对于已知的类型宽度,例如uint32_t

uint32_t div3_exact_only(uint32_t x) {
    return x * 0xaaaaaaabU;
}
Run Code Online (Sandbox Code Playgroud)


Cos*_*nus 11

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

相反的3^(-1) mod 3,你要多想想0.3333333在那里0的前.位于高32位,在3333位于低32位。这个定点运算工作正常,但结果显然移到了 的上半部分rax,因此CPU必须在运算后再次将结果移下。

为什么我们使用 imul 而不是 mul?我认为模算术都是无符号的。

没有与MUL指令等效的IMUL指令。使用的IMUL变体需要两个寄存器:

a <= a * b
Run Code Online (Sandbox Code Playgroud)

没有任何MUL指令可以做到这一点。MUL指令更昂贵,因为它们将结果作为 128 位存储在两个寄存器中。当然,您可以使用旧指令,但这不会改变结果存储在两个寄存器中的事实。

  • 它不直接使用 rdi 的原因是,当它作为参数接收时,rdi 的高位不一定为 0。清除寄存器高 32 位的最佳方法是使用 32 位移动指令,该指令会自动清除它们。编译器可以使用“mov edi, edi”。但是 rcx 是可用的,而且事实证明,在某些情况下“mov ecx, edi”实际上是免费的(由寄存器重命名器处理),而“mov edi, edi”不一定是免费的。 (2认同)

rcg*_*ldr 8

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

为什么GCC在实现整数除法时使用乘以奇怪的数字?

它包含一个 pdf 文章的链接,解释了这一点(我的回答澄清了这篇 pdf 文章中没有很好解释的内容):

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

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

如果更改为,请查看生成的代码

unsigned div7(unsigned x) {
    return x / 7;
}
Run Code Online (Sandbox Code Playgroud)

所以为了解释这个过程,让 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
Run Code Online (Sandbox Code Playgroud)

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

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)
Run Code Online (Sandbox Code Playgroud)

所以乘数是 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
L = L-1 = 2
...                 visual studio generated code for div7, input is rcx
mov eax, 613566757
mul ecx
sub ecx, edx                   ; handle 2^32 bit
shr ecx, 1                     ; ...
lea eax, DWORD PTR [edx+ecx]   ; ...
shr eax, 2
Run Code Online (Sandbox Code Playgroud)

如果需要余数,则可以使用以下步骤:

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


gna*_*729 5

x/3 大约为 (x * (2^32/3)) / 2^32。因此,我们可以执行一次 32x32->64 位乘法,取较高的 32 位,并得到大约 x/3。

\n

存在一些错误,因为我们无法精确乘以 2^32/3,只能乘以该数字四舍五入为整数。我们使用 x/3 \xe2\x89\x88 (x * (2^33/3)) / 2^33 获得更高的精度。(我们不能使用 2^34/3,因为它 > 2^32)。事实证明,这足以在所有情况下准确地得到 x/3。如果输入是 3k 或 3k+2,您可以通过检查公式是否给出 k 的结果来证明这一点。

\n