一、强度削弱
之前在偶尔看gcc对于除以一个常数的表达式生成的汇编代码中，发现一条除法表达式生成的汇编指令非常多，这些指令中没有乘法操作，比较明显的特征就是进行了一个大整数的乘法，之后是移位啊、减法啊什么的操作，虽然不知道是什么意思，但是感觉很厉害的样子。
后台就觉得这应该是一个优化，搜索了一下，看到gcc中的相关实现是基于一篇文章的描述实现，事实上，gcc中这部分代码的作者就是该文章的一个作者，结合gcc的代码、生成的汇编指令及这篇文章的内容惊人的一致，所以要理解这部分代码，这篇文章是基础。
理论上的东西在计算机工程中可能会收到各种限制，例如说理论上从不用考虑整数分解的复杂度，不用考虑浮点数的精度，整数的截断、指令中操作数长度等工程性的问题，而这些则是计算机工程中需要考虑的问题。
二、基本的想法
基本的想法就是将除法转换为整数乘法，即首先将计算结果放大，而不重要的部分将会受到计算机指令及操作数的影响而被自动丢弃，这个“不重要部分”需要满足的基本条件就是足够小，永远到底有多远，足够小到底有多小。直观上可以这么想，对于正常的整数除法(n/d)*x=y，每当x增加d，此时的y增加1。我们可以设计一个满足这个特性的式子n*m来代替n/d，n*m也满足对于n*m*x，每当x增加d时，n*m*x增加1+delta，也就是说n*m也增加1并且多了一点点，这一点点需要足够的小，小到满足文章中给出的基础表达式：
2^(N+l) <= m*d <= 2^(N + l) + 2^l
这里出现的N是为了便于使用大部分体系结构中都存在的两个整数相乘，结果存放在两个整数寄存器中，例如intel体系下两个int相乘的结果放在eax和edx中，使用N这样就可以直接从寄存器中获得这个移位前的值。上面的表达式，对于任意一个-2^N到2^N之间的数，m*d右移N+l位，它和1的差距都小于1/2^N，从而保证误差在1之内，利用整数的向下取整即可获得正确的结果。
还有的限制是：
所有的操作数及中间结果不能溢出一个int的大小，这一点是一个必须小心谨慎保证的一个变量。
这里最为复杂的表达式是当m>=2^m时，其计算方法为
t1=MULUH(m-2^N,n)
q=SRL(t1+SRL(n-t1,1),shpost-1)
其实它的第二个式子就是
SRL(t1+(n-t1)/2,shpost-1)=SRL((n+t1)/2,shpost-1)=SRL((n+t1),shpost)
由于t1是m*n结果中由小于2^N部分和m乘法之后的高int内容，加上第二个式子中计算的m本身大于2^N的部分，就是m*n中高int中的内容，最后再右移post，即最开始表达式中的l，总共移动了N+l个bit，满足文章最开始的定理限制。
t1+(n-t1)/2
这个需要解释下：
由于
t1=MULUH(m-2^N,n) <= MULUH(2^N,n) < n
所以n-t1>0，进而由于
t1+(n-t1)/2=(n+t1)/2，由于2<2^N，所有(n+t1)/2能够保存在一个Nbit的int中。
再说明下m为什么最多占用N+1个bit。在choose_multiplier函数中
由于 m*d<2^(N+post)(1+2^(-prec))所以
m<2^(N+post)(1+2^(-prec))/d
由于2^(l-1) <d,所以
m<2^(N+post)(1+2^(-prec))/d<2^(N+post)(1+2^(-prec))/2^(l-1)<2^(N+post)(1)/2^(1-l)=2^(N+post+1-l)
加上N+post<= l+prec，所以
m<2^(prec+1)
当prec为无符号时，此时prec为32，m最大为33bit。
三、有符号处理
基本和这个类似，但是有负数的算数右移和正数右移不同，
tsecer@harry :echo $((-5>>2)) $((5>>2))
-2 1
这一点可以这么理解，内码表示统一看作负数则这些数是按照
0 1 0x7FFFFFFF 0x80000000 …… 0xFFFFFFFE(-2) 0xFFFFFFFF(-1)
当移位时始终是向左取整，转换为负数之后就是向负无穷偏移了。为了这种情况，需要专门为这种移位进行一个处理，就是负数移位需要加上1，正数移位则不需要。
四、测试代码例子
tsecer@harry :cat div2mul.c
unsigned int uiDPos(unsigned int x)
{
return x / 21;
}
int iDPos(int x)
{
return x / 21;
}
int uiDNeg(unsigned int x)
{
return x / -21;
}
int iDNeg(int x)
{
return x / -21;
}
tsecer@harry :gcc div2mul.c -S
tsecer@harry :cat div2mul.s
.file "div2mul.c"
.text
.globl uiDPos
.type uiDPos, @function
uiDPos:
.LFB2:
pushq %rbp
.LCFI0:
movq %rsp, %rbp
.LCFI1:
movl %edi, -4(%rbp)
movl -4(%rbp), %eax
movl %eax, -12(%rbp)
movl $-2045222521, -16(%rbp)
movl -16(%rbp), %eax
mull -12(%rbp)
movl %edx, %ecx
movl -12(%rbp), %eax
subl %ecx, %eax
shrl %eax
leal (%rcx,%rax), %eax
shrl $4, %eax
leave
ret
.LFE2:
.size uiDPos, .-uiDPos
.globl iDPos
.type iDPos, @function
iDPos:
.LFB3:
pushq %rbp
.LCFI2:
movq %rsp, %rbp
.LCFI3:
movl %edi, -4(%rbp)
movl -4(%rbp), %ecx
movl $818089009, -12(%rbp)
movl -12(%rbp), %eax
imull %ecx
sarl $2, %edx
movl %ecx, %eax
sarl $31, %eax
movl %edx, %ecx
subl %eax, %ecx
movl %ecx, %eax
leave
ret
.LFE3:
.size iDPos, .-iDPos
.globl uiDNeg
.type uiDNeg, @function
uiDNeg:
.LFB4:
pushq %rbp
.LCFI4:
movq %rsp, %rbp
.LCFI5:
movl %edi, -4(%rbp)
movl -4(%rbp), %eax
cmpl $-21, %eax
setae %al
movzbl %al, %eax
leave
ret
.LFE4:
.size uiDNeg, .-uiDNeg
.globl iDNeg
.type iDNeg, @function
iDNeg:
.LFB5:
pushq %rbp
.LCFI6:
movq %rsp, %rbp
.LCFI7:
movl %edi, -4(%rbp)
movl -4(%rbp), %ecx
movl $818089009, -12(%rbp)
movl -12(%rbp), %eax
imull %ecx
sarl $2, %edx
movl %ecx, %eax
sarl $31, %eax
subl %edx, %eax
leave
ret
.LFE5:
.size iDNeg, .-iDNeg
.section .eh_frame,"a",@progbits
.Lframe1:
.long .LECIE1-.LSCIE1
.LSCIE1:
.long 0x0
.byte 0x1
.string "zR"
.uleb128 0x1
.sleb128 -8
.byte 0x10
.uleb128 0x1
.byte 0x3
.byte 0xc
.uleb128 0x7
.uleb128 0x8
.byte 0x90
.uleb128 0x1
.align 8
.LECIE1:
.LSFDE1:
.long .LEFDE1-.LASFDE1
.LASFDE1:
.long .LASFDE1-.Lframe1
.long .LFB2
.long .LFE2-.LFB2
.uleb128 0x0
.byte 0x4
.long .LCFI0-.LFB2
.byte 0xe
.uleb128 0x10
.byte 0x86
.uleb128 0x2
.byte 0x4
.long .LCFI1-.LCFI0
.byte 0xd
.uleb128 0x6
.align 8
.LEFDE1:
.LSFDE3:
.long .LEFDE3-.LASFDE3
.LASFDE3:
.long .LASFDE3-.Lframe1
.long .LFB3
.long .LFE3-.LFB3
.uleb128 0x0
.byte 0x4
.long .LCFI2-.LFB3
.byte 0xe
.uleb128 0x10
.byte 0x86
.uleb128 0x2
.byte 0x4
.long .LCFI3-.LCFI2
.byte 0xd
.uleb128 0x6
.align 8
.LEFDE3:
.LSFDE5:
.long .LEFDE5-.LASFDE5
.LASFDE5:
.long .LASFDE5-.Lframe1
.long .LFB4
.long .LFE4-.LFB4
.uleb128 0x0
.byte 0x4
.long .LCFI4-.LFB4
.byte 0xe
.uleb128 0x10
.byte 0x86
.uleb128 0x2
.byte 0x4
.long .LCFI5-.LCFI4
.byte 0xd
.uleb128 0x6
.align 8
.LEFDE5:
.LSFDE7:
.long .LEFDE7-.LASFDE7
.LASFDE7:
.long .LASFDE7-.Lframe1
.long .LFB5
.long .LFE5-.LFB5
.uleb128 0x0
.byte 0x4
.long .LCFI6-.LFB5
.byte 0xe
.uleb128 0x10
.byte 0x86
.uleb128 0x2
.byte 0x4
.long .LCFI7-.LCFI6
.byte 0xd
.uleb128 0x6
.align 8
.LEFDE7:
.ident "GCC: (GNU) 4.1.2 20070115 (prerelease) (SUSE Linux)"
.section .note.GNU-stack,"",@progbits
tsecer@harry :
五、gcc代码
/* Choose a minimal N + 1 bit approximation to 1/D that can be used to
replace division by D, and put the least significant N bits of the result
in *MULTIPLIER_PTR and return the most significant bit.
The width of operations is N (should be <= HOST_BITS_PER_WIDE_INT), the
needed precision is in PRECISION (should be <= N).
PRECISION should be as small as possible so this function can choose
multiplier more freely.
The rounded-up logarithm of D is placed in *lgup_ptr. A shift count that
is to be used for a final right shift is placed in *POST_SHIFT_PTR.
Using this function, x/D will be equal to (x * m) >> (*POST_SHIFT_PTR),
where m is the full HOST_BITS_PER_WIDE_INT + 1 bit multiplier. */
static
unsigned HOST_WIDE_INT
choose_multiplier (unsigned HOST_WIDE_INT d, int n, int precision,
rtx *multiplier_ptr, int *post_shift_ptr, int *lgup_ptr)
{
HOST_WIDE_INT mhigh_hi, mlow_hi;
unsigned HOST_WIDE_INT mhigh_lo, mlow_lo;
int lgup, post_shift;
int pow, pow2;
unsigned HOST_WIDE_INT nl, dummy1;
HOST_WIDE_INT nh, dummy2;
/* lgup = ceil(log2(divisor)); */
lgup = ceil_log2 (d);
gcc_assert (lgup <= n);
pow = n + lgup;
pow2 = n + lgup - precision;
/* We could handle this with some effort, but this case is much
better handled directly with a scc insn, so rely on caller using
that. */
gcc_assert (pow != 2 * HOST_BITS_PER_WIDE_INT);
/* mlow = 2^(N + lgup)/d */
if (pow >= HOST_BITS_PER_WIDE_INT)
{
nh = (HOST_WIDE_INT) 1 << (pow - HOST_BITS_PER_WIDE_INT);
nl = 0;
}
else
{
nh = 0;
nl = (unsigned HOST_WIDE_INT) 1 << pow;
}
div_and_round_double (TRUNC_DIV_EXPR, 1, nl, nh, d, (HOST_WIDE_INT) 0,
&mlow_lo, &mlow_hi, &dummy1, &dummy2);
/* mhigh = (2^(N + lgup) + 2^N + lgup - precision)/d */
if (pow2 >= HOST_BITS_PER_WIDE_INT)
nh |= (HOST_WIDE_INT) 1 << (pow2 - HOST_BITS_PER_WIDE_INT);
else
nl |= (unsigned HOST_WIDE_INT) 1 << pow2;
div_and_round_double (TRUNC_DIV_EXPR, 1, nl, nh, d, (HOST_WIDE_INT) 0,
&mhigh_lo, &mhigh_hi, &dummy1, &dummy2);
gcc_assert (!mhigh_hi || nh - d < d);
gcc_assert (mhigh_hi <= 1 && mlow_hi <= 1);
/* Assert that mlow < mhigh. */
gcc_assert (mlow_hi < mhigh_hi
|| (mlow_hi == mhigh_hi && mlow_lo < mhigh_lo));
/* If precision == N, then mlow, mhigh exceed 2^N
(but they do not exceed 2^(N+1)). */
/* Reduce to lowest terms. */
for (post_shift = lgup; post_shift > 0; post_shift--)
{
unsigned HOST_WIDE_INT ml_lo = (mlow_hi << (HOST_BITS_PER_WIDE_INT - 1)) | (mlow_lo >> 1);
unsigned HOST_WIDE_INT mh_lo = (mhigh_hi << (HOST_BITS_PER_WIDE_INT - 1)) | (mhigh_lo >> 1);
if (ml_lo >= mh_lo)
break;
mlow_hi = 0;
mlow_lo = ml_lo;
mhigh_hi = 0;
mhigh_lo = mh_lo;
}
*post_shift_ptr = post_shift;
*lgup_ptr = lgup;
if (n < HOST_BITS_PER_WIDE_INT)
{
unsigned HOST_WIDE_INT mask = ((unsigned HOST_WIDE_INT) 1 << n) - 1;
*multiplier_ptr = GEN_INT (mhigh_lo & mask);
return mhigh_lo >= mask;
}
else
{
*multiplier_ptr = GEN_INT (mhigh_lo);
return mhigh_hi;
}
}