我一直在阅读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

但是看看生成除法。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 ?它是如何产生这个神奇的数字的,为什么一切都有效?


当前回答

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

一个问题是我们不能精确地表示倒数(除非除法是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。i的值小于264,因此相乘后的误差小于1/5。除4后误差小于(1/5)* 2−2。

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


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

如果我们试图将4/7表示为一个0.64的定点数,舍入为零,我们最终会得到(6/7)* 2-64的误差。在乘以一个略低于264的i值后,我们最终得到的误差略低于6/7,在除以4后,我们最终得到的误差略低于1.5/7,大于1/7。

因此,为了正确地实现除7,我们需要乘以一个0.65的定点数。我们可以通过乘以定点数的下64位来实现这一点,然后加上原始数(这可能会溢出到进位),然后做一个进位旋转。

其他回答

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

一个问题是我们不能精确地表示倒数(除非除法是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。i的值小于264,因此相乘后的误差小于1/5。除4后误差小于(1/5)* 2−2。

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


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

如果我们试图将4/7表示为一个0.64的定点数,舍入为零,我们最终会得到(6/7)* 2-64的误差。在乘以一个略低于264的i值后,我们最终得到的误差略低于6/7,在除以4后,我们最终得到的误差略低于1.5/7,大于1/7。

因此,为了正确地实现除7,我们需要乘以一个0.65的定点数。我们可以通过乘以定点数的下64位来实现这一点,然后加上原始数(这可能会溢出到进位),然后做一个进位旋转。

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

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

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

Let's consider for a minute why 0.CCCCCCCC... (hex) or 0.110011001100... binary is 4/5. Divide the binary representation by 4 (shift right 2 places), and we'll get 0.001100110011... which by trivial inspection can be added the original to get 0.111111111111..., which is obviously equal to 1, the same way 0.9999999... in decimal is equal to one. Therefore, we know that x + x/4 = 1, so 5x/4 = 1, x=4/5. This is then represented as CCCCCCCCCCCCD in hex for rounding (as the binary digit beyond the last one present would be a 1).

我会从一个稍微不同的角度回答:因为它是允许这样做的。

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

The compiler is allowed to make ANY changes as long as it doesn't change the observable behaviour as specified by the abstract machine. There is no reasonable expectation that the compiler will transform your code in the most straightforward way possible (even when a lot of C programmer assume that). Usually, it does this because the compiler wants to optimize the performance compared to the straightforward approach (as discussed in the other answers at length). If under any circumstances the compiler "optimizes" a correct program to something that has a different observable behaviour, that is a compiler bug. Any undefined behaviour in our code (signed integer overflow is a classical example) and this contract is void.

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

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

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

因此,只有在编译时不知道除数时,才会在输出中看到div或idiv。

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

这里是一个算法文档的链接,该算法生成的值和代码是我在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))时才能大于一个udword,在这种情况下,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