GCC为什么在整数除法实现中使用一个奇怪的数字来进行乘法运算？

Question

GCC为什么在整数除法实现中使用一个奇怪的数字来进行乘法运算？

269

我一直在阅读有关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文件，它并不包含任何除法操作！相反，它使用位移和魔法数进行某种黑魔法。这里是一个计算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呢？它是如何生成这个神奇的数字的，为什么一切都能正常工作？

- qiubit

33

gcc优化常数除法，尝试用2、3、4、5、6、7、8进行除法运算，你很可能会看到每种情况下的代码都非常不同。 - Jabberwocky

33

注意：魔数“-3689348814741910323”在使用“uint64_t”类型时转换为“CCCCCCCCCCCCCCCD”，大约相当于(2^64)*4/5。 - chux - Reinstate Monica

36

编译器不会因为关闭优化而故意生成低效的代码。它将进行一个不涉及代码重排或变量消除的微小"优化"，例如单个源代码语句将被翻译为该操作最高效的代码。编译器优化考虑的是周围的代码而不仅仅是单条语句。 - Clifford

25

阅读这篇精彩的文章：劳动分工 - Jester

12

有些编译器实际上会生成效率低下的代码，因为禁用了优化。特别是为了方便调试（如在单独的代码行上设置断点），它们会这样做。事实上，GCC相当不寻常，因为它没有真正的“无优化”模式，因为许多优化都是固定开启的。你可以通过GCC来看到这一点。另一方面，Clang和MSVC会在“-O0”时发出“div”指令。 - Cody Gray

显示剩余17条评论

5个回答

142

将一个数除以5等价于将它乘以1/5，同样也等价于将其乘以4/5并右移两位。所涉及的值是十六进制下的CCCCCCCCCCCCCCCD，如果放在十六进制小数点后面，它就是4/5的二进制表示（即四分之四五的二进制是0.110011001100循环的——请参见下文了解原因）。我想你可以从这里开始！你可能想看看定点算术（注意最后会舍入到整数）。

至于为什么，乘法比除法快，当除数固定时，这是一条更快的路线。倒数乘法教程详细介绍了如何运作，解释了定点数的术语，展示了如何找到倒数的算法，以及如何处理有符号除法和模数。

让我们考虑一下为什么0.CCCCCCCC...（十六进制）或0.110011001100...（二进制）是四分之四五。将二进制表示除以4（向右移动2位），我们得到0.001100110011...，显然可以通过简单的检查与原始数相加得到0.111111111111...，它显然等于1，就像十进制中的0.9999999...等于一样。因此，我们知道x + x/4 = 1，所以5x/4 = 1，x=4/5。然后将其表示为CCCCCCCCCCCCD（舍入为最后一个二进制位之后的数字是1）。

- abligh

2

@user2357112 可以随意发表自己的答案，但我不同意。您可以将乘法视为64.0位乘以0.64位，得到128位定点答案，其中最低的64位被丢弃，然后除以4（正如我在第一段中指出的那样）。您可能能够提出另一种模数算术答案，同样可以解释位移，但我相当确定这个解释是有效的。 - abligh

8

这个值实际上是“CCCCCCCCCCCCCCCD”。最后的D很重要，确保当结果被截断时，精确除法会得出正确的答案。 - plugwash

7

没关系。我没有看到他们正在获取128位乘法结果的上64位；这不是大多数语言可以做到的，所以一开始我没有意识到这种情况。如果明确提到如何取128位结果的上64位相当于乘以一个定点数并向下舍入，那么这个答案会更好。（此外，最好解释为什么必须是4/5而不是1/5，以及为什么我们必须将4/5向上舍入而不是向下。） - user2357112

2

据我所知，您需要计算出需要多大的误差才能使除以5向上跨越舍入边界，然后将其与您计算中的最坏情况误差进行比较。可以假设gcc开发人员已经这样做，并得出结论它总是会给出正确的结果。 - plugwash

4

实际上，您可能只需要检查5个可能的最高输入值，如果这些值正确舍入，那么其他所有值也应该正确。 - plugwash

显示剩余5条评论

65

通常乘法比除法更快。因此，如果我们可以用倒数取代除法，就可以显著加速一个常数的除法。

问题在于，我们无法精确表示倒数（除非除数是二的幂，但在这种情况下，我们通常只需将除法转换为位移）。因此，为了确保正确的答案，我们必须小心，以确保我们的倒数中的误差不会导致最终结果出错。

-3689348814741910323 是 0xCCCCCCCCCCCCCCCD，在 0.64 固定点表示中略大于4/5。

当我们将一个64位整数乘以一个0.64固定点数时，我们得到一个64.64结果。我们将值截断为64位整数（有效地四舍五入至零），然后执行进一步的移位，将其除以四并再次截断。通过查看位级别，很明显我们可以将这两个截断视为单个截断。

这显然为我们提供了至少近似于除以5的结果，但它是否给出了向零正确舍入的精确答案？

要获得精确答案，误差需要足够小，以不会推动答案超过舍入边界。

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

我们常数的误差为（1/5）* 2^-64。i的值小于2^64，因此在乘法后的误差小于1/5。除以4后，误差小于(1/5)*2^-2。

(1/5)*2^-2 < 1/5，因此答案总是等于进行精确除法并向零舍入。

不幸的是，这并不适用于所有除数。

如果我们尝试使用向零舍入方式将4/7表示为0.64固定点数，那么误差将为(6/7) * 2^-64。在乘以略小于2⁶⁴的i值之后，我们得到的误差只有不到6/7，在除以四之后，我们得到的误差仅略小于1.5/7，大于1/7。

因此，要正确实现除以7，我们需要乘以一个0.65固定点数。我们可以通过乘以固定点数的低64位，然后加上原始数字（这可能会溢出到进位位），然后进行带进位的旋转来实现该操作。

- plugwash

9

这个答案将模数乘法逆元从“看起来比我想花时间学习的数学更复杂”转换成了一些有意义的东西。因为易于理解，所以加1分。我从未需要做其他事情，只需使用编译器生成的常量，因此我只是浏览其他解释这个数学概念的文章。 - Peter Cordes

2

我在代码中完全看不到任何与模运算有关的内容。不知道其他评论者从哪里得出这个结论。 - plugwash

4

这里的取模是指对2的n次幂取模，就像寄存器中的所有整数运算一样。详见：https://zh.wikipedia.org/wiki/%E6%A8%A1%E5%8F%8D%E7%B4%A0%E5%8F%96%E5%80%BC#%E5%BA%94%E7%94%A8 - Peter Cordes

5

据我所知，模乘逆元用于精确除法，但并不适用于一般的除法。 - harold

5

@PeterCordes 乘以定点数倒数？我不知道大家都怎么称呼它，但我可能会这样称呼它，因为这个名称相当描述性。 - harold

显示剩余5条评论

14

这里有一个算法文档的链接，它生成了我在Visual Studio（大多数情况下）看到的值和代码，并且我假设它仍然用于GCC对变量整数除以常量整数的操作。

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

在文章中，uword有N位，udword有2N位，n = 分子 = 被除数，d= 分母 = 除数，ℓ最初设置为ceil(log2(d))，shpre是预移位（用于乘法之前）= e = d中的尾随零位数，shpost是后移位（用于乘法之后），prec是精度 = 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。当分母（除数）> 2^(N-1)（当ℓ == N => mlow = 2^(2N)）时，mlow和mhigh的分子（被除数）可能大于udword，此时n/d的优化替换是一个比较操作（如果n>=d，则q = 1，否则q = 0），因此不生成乘法器。mlow和mhigh的初始值将是N+1位，每个N+1位值可以使用两个udword/uword除法产生（mlow或mhigh）。以64位模式下的X86为例：

; 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测试这个。你已经看到了i/5的处理方式。再来看看i/7的处理方式（应该是N+1位乘法器的情况）。

在大多数现代处理器上，乘法器都有固定的时序，因此不需要预先进行移位。对于X86，大多数除数的最终结果是一个两条指令的序列，而像7这样的除数则需要五条指令的序列（以模拟pdf文件中的等式4.5和图4.2所示的N+1位乘法器）。以下是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 = multiplier - 2^64，rbx = dividend。其中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

- rcgldr

1

那篇论文描述了在gcc中实现它，因此我认为可以安全地假设相同的算法仍在使用。 - Peter Cordes

1

那篇1994年的论文描述了在gcc中实现它，所以gcc有时间更新其算法。以防其他人没有时间检查一下URL中的94是什么意思。 - Ed Grimm

0

我会从稍微不同的角度回答：因为允许这样做。

C和C++是针对抽象机器定义的。编译器按照as-if规则将程序转换为具体机器。

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

- dmeister

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Sneftel · Accepted Answer

整数除法是现代处理器上执行的最慢的算术运算之一，延迟高达数十个周期，吞吐量差。(对于x86，请参见Agner Fog's instruction tables and microarch guide)。

如果您提前知道除数，可以通过将其替换为一组其他操作（乘法、加法和移位）来避免除法，这些操作具有等效的效果。即使需要多个操作，通常仍然比整数除法本身快得多。

用这种方式实现C的/运算符而不是使用涉及div的多条指令序列只是GCC默认使用常量除法的方法。它不需要跨操作进行优化，甚至对调试也没有任何影响。（使用-Os进行小型代码大小时，GCC会使用div。）使用乘法逆元而不是除法就像使用lea而不是mul和add。

作为结果，如果除数在编译时未知，则输出中只会看到div或idiv。

有关编译器如何生成这些序列的信息以及让您自己生成它们的代码（除非您正在使用一个愚蠢的编译器，否则几乎肯定是不必要的），请参见libdivide。