注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Tsecer的回音岛

Tsecer的博客

 
 
 

日志

 
 

gcc如何将常量除法转换为乘法及移位  

2014-09-30 21:58:28|  分类: Gcc源代码分析 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
一、强度削弱
之前在偶尔看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;
    }
}
  评论这张
 
阅读(683)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017