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

我一直在阅读div和mul汇编操作，我决定用C写一个简单的程序来看看它们的实际操作:

文件division.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 ?它是如何产生这个神奇的数字的，为什么一切都有效?

32050 次浏览
小开

整数除法是在现代处理器上可以执行的最慢的算术运算之一，延迟可达数十个周期，吞吐量很差。(对于x86，请参见Agner Fog的指令表和微拱指南)。

如果您事先知道除法，那么可以通过将除法替换为一组具有等效效果的其他操作(乘法、加法和移位)来避免除法。即使需要数个运算，它通常仍然比整数除法本身快得多。

以这种方式实现C /操作符，而不是使用包含div的多指令序列，这只是GCC按常量除法的默认方式。它不需要跨操作进行优化，甚至在调试时也不需要更改任何内容。(不过，使用-Os来实现较小的代码大小确实会让GCC使用div。)使用乘法逆而不是除法，就像使用lea而不是mul和add

因此，如果在编译时除数未知，则只能在输出中看到div或idiv。

有关编译器如何生成这些序列的信息，以及让你自己生成它们的代码(几乎肯定是不必要的，除非你使用的是脑残编译器)，请参阅libdivide。

小开

除以5就等于乘以1/5，也就是乘以4/5然后右移2位。有关的值是十六进制的CCCCCCCCCCCCCCCD，这是4/5的二进制表示，如果放在十六进制点之后(即四分之五的二进制是0.110011001100循环-原因见下文)。我想你可以接手了!你可能想要检查定点算法(尽管注意它在结尾被舍入为整数)。

至于为什么，乘法比除法快，当除数固定时，这是一条更快的路径。

关于它如何工作的详细描述，请参阅倒数乘法教程，并从定点的角度进行解释。它展示了如何寻找倒数的算法，以及如何处理有符号除法和取模。

让我们考虑一下为什么0.CCCCCCCC...(十六进制)或0.110011001100...二进制是4/5。将二进制表示除4(右移2位)，我们将得到0.001100110011...，通过简单的检查，可以将原始的0.111111111111...相加，得到0.111111111111...，显然等于1，同样的方式，十进制的0.9999999...等于1。因此，我们知道x + x/4 = 1，所以5x/4 = 1， x=4/5。然后用十六进制表示为CCCCCCCCCCCCD进行四舍五入(因为在最后一位之后的二进制数字将是1)。

小开

一般来说，乘法比除法快得多。所以如果我们可以不用乘以倒数来代替我们可以用一个常数来大大加快除法的速度

一个问题是我们不能精确地表示倒数(除非除法是2的幂，但在这种情况下，我们通常可以将除法转换为位移位)。因此，为了确保正确的答案，我们必须注意，倒数中的错误不会导致最终结果的错误。

-3689348814741910323是0xcccccccccccccd，这是一个用0.64固定点表示的略高于4/5的值。

当我们用一个64位整数乘以一个0.64的定点数时，我们得到一个64.64的结果。我们将值截断为一个64位整数(实际上是四舍五入到零)，然后执行进一步的移位，除以4，再次截断。通过查看位级，很明显，我们可以将两次截断视为单个截断。

显然，这至少给了我们一个除以5的近似值，但它能给我们一个精确的趋近于0的答案吗?

为了得到一个精确的答案，误差需要足够小，不会把答案推到四舍五入的边界。

除5的准确答案总是有一个小数部分，0,1 / 5,2 / 5,3 /5或4/5。因此，在相乘和移位结果中小于1/5的正误差将永远不会使结果超过舍入边界。

常量的误差是(1/5)* 2-64年。我的值小于264，因此相乘后的误差小于1/5。除4后，误差小于(1/5)* 2-; 2。

(1/5) * 2-; 2 <1/5所以答案总是等于做一个精确的除法并四舍五入到0。

不幸的是，这并不适用于所有因子。

如果我们试图将4/7表示为一个0.64固定点数，舍入为0，我们最终会得到(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是pre-shift(用于乘前)= e = d中后面的零位数，shpost是post-shift(用于乘后)，prec是precision = 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。mlow和mhigh的分子(被除数)只有在分母(除数)>2^(N-1)(当ℓ== N =>mlow = 2^(2N))，在这种情况下，n/d的优化替换是一个比较(如果n>=d, q = 1，否则q = 0)，因此不会产生乘数。mlow和mhigh的初始值将是N+1位，并且可以使用两个udword/uword分割来生成每个N+1位值(mlow或mhigh)。以X86 64位模式为例:

; 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++是针对抽象机器定义的。编译器按照似规则将该程序从抽象机转换为具体机。

  • 编译器可以做任何改变，只要它不改变抽象机器指定的可观察行为。没有合理的期望编译器会以最直接的方式转换你的代码(即使很多C程序员都这么认为)。通常，它这样做是因为编译器想要优化与直接方法相比的性能(在其他答案中详细讨论)。
  • 如果在任何情况下，编译器将一个正确的程序“优化”为具有不同的可观察行为，那就是编译器错误。
  • 在我们的代码中任何未定义的行为(带符号整数溢出是一个经典的例子)，这个合同是无效的。

提交答案