这里是一个算法文档的链接,该算法生成的值和代码是我在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