为什么在x86上除以3需要右移和其他奇数? 编辑1 编辑2

问题描述

我具有以下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 2 32 相乘,即2863311531。

有些事情我还是不明白:

  1. 为什么我们需要完全使用ecx / rcx?我们不能直接将raxedi相乘吗?
  2. 为什么要在64位模式下进行乘法运算?将eaxecx相乘会更快吗?
  3. 我们为什么使用imul而不是mul?我以为模块化算术将是无符号的。
  4. 最后的33位右移是怎么回事?我以为我们可以删除最高的32位。

编辑1

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

// 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。)

解决方法

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

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

(如果无法避免的话,由于limitations on mov-elimination,最好零扩展到另一个寄存器中。)

从字面上更进一步地回答您的问题,不,x86没有任何乘法指令对它们的输入之一进行零扩展以使您可以将32位和64位寄存器相乘。两个输入的宽度必须相同。

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

(术语:所有这些代码都以64位模式运行。您要问为什么要使用64位 operand-size 。)

可以 mul edi EAX 与EDI相乘以得到EDX:EAX的64位结果,但是mul edi为3与大多数具有快速64位imul的现代x86-64 CPU相比,它具有Intel CPU的性能。 (尽管imul r64,r64在AMD Bulldozer系列和某些低功率CPU上速度较慢。)https://uops.info/https://agner.org/optimize/(指令表和Microarch PDF) (有趣的事实:mul rdi实际上是Intel CPU上的更便宜的,只有2 oups。也许与不必对整数乘法单元的输出进行额外的拆分有关,例如{{ 1}}必须将64位低半乘法器输出拆分为EDX和EAX一半,但这自然发生在64x64 => 128位mul上。)

您想要的零件也位于EDX中,因此您需要另一个mul edi来处理它。 (同样,因为我们正在查看的是该函数的独立定义的代码,而不是在内联到调用方之后。)

GCC 8.3和更早版本的 did 使用32位mov eax,edx而不是64位mulhttps://godbolt.org/z/5qj7d5)。当Bulldozer系列和旧的Silvermont CPU更加相关时,这对于imul来说并不疯狂,但是对于最近的GCC而言,这些CPU在过去更遥远,其通用调整选择反映了这一点。不幸的是,GCC还浪费了-mtune=generic指令,将EDI复制到EAX,使这种方式看起来更糟:/

mov

使用# 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时只有6 oups,但仍然比:

mul edi

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

  1. 为什么我们使用imul代替mul?我以为模块化算术将是无符号的。

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

只有imul具有更快,更灵活的仅对下半部分的格式。关于imul的唯一签名是,它根据下半部分的有符号溢出来设置OF。仅拥有一个imul reg,reg与FLAGS输出唯一的mul r,r是不值得花费更多的操作码和更多的晶体管的。

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

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

否,如果您以这种方式实现,则没有乘数常量可以为每个可能的输入imul r,r提供正确的正确答案。允许近似值,仅实现为程序使用的每个输入产生完全相同的可观察行为的实现。除了不知道x的整个范围之外,x的值范围不存在,编译器没有该选项。 ({unsigned仅适用于浮点;如果您想更快地近似整数数学,请按如下所示手动进行编码):

有关用于编译器按编译时间常数进行精确除法的定点乘法逆方法的更多信息,请参见Why does GCC use multiplication by a strange number in implementing integer division?

有关在一般情况下不能正常工作的示例,请参见我对Divide by 10 using bit shifts?上的答案的修改

-ffast-math

// 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); } 实际上是107374182时,它的第一个错误答案(如果从0向上循环)是div10(1073741829) = 107374183。(它应四舍五入,而不是像C整数除法那样取整为0。)>


从您的编辑中,我看到您实际上是在使用乘积结果的 low 一半,显然,该结果对于直到UINT_MAX的精确倍数都非常适用。

正如您所说,如果除法有余数,例如1073741829/10 = 16 * 0xaaaaaaab被截断为32位而不是0xaaaaaab0时。

5

是的,如果算术可行,那么对于编译器而言,使用32位imul实现该方法是合法且最佳的。他们不寻求这种优化,因为这鲜为人知。如果值得在编译时间方面增加编译器代码甚至寻找优化,则IDK值得一提,更不用说在开发人员时间中的编译器维护成本了。这不是运行时成本上的巨大差异,而且几乎不可能实现。不过很好。

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

但是,至少在div3_exact_only: imul eax,edi,0xAAAAAAAB # 1 uop,3c latency ret 之类的已知类型宽度中,您可以在源代码中完成此操作:

uint32_t
,

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

您必须考虑3^(-1) mod 3而不是0.3333333,其中0之前的.位于高32位,而3333位于较低的32位。 该定点运算可以很好地工作,但是结果显然移到了rax的上部,因此CPU必须在运算后再次将结果向下移。

为什么我们使用imul代替mul?我以为模块化算术都是无符号的。

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

a <= a * b

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

,

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

Why does GCC use multiplication by a strange number in implementing integer division?

它包含指向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是

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

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

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的结果来证明这一点。

相关问答

错误1:Request method ‘DELETE‘ not supported 错误还原:...
错误1:启动docker镜像时报错:Error response from daemon:...
错误1:private field ‘xxx‘ is never assigned 按Alt...
报错如下,通过源不能下载,最后警告pip需升级版本 Requirem...