除法运算

有符号div,无符号idiv

标 题: 【原创】第一章:1.9、除法与取模运算的识别与优化原理
作 者: AOnePass
时 间: 2010-09-02,23:59:39

链 接: http://bbs.pediy.com/showthread.php?t=119744

在计算机的0和1的世界里，最复杂的整数运算莫过于除法的运算了，对于这一节内容的演练断断续续的花费了笔者近一周的时间，

笔者一直比较头痛的问题是怎样让初学者来理解这复杂而微妙的优化过程，但是经过多次演练后，最终还是打算从最基本的数学概念开始，

笔者个人认为这是能快速理解除法优化的唯一捷径。

    这是笔者有史以来花费心思最多的一篇文章，希望各位读者能珍惜笔者的心得结晶。

除法与倒数相乘

何为倒数相乘？很简单，编译器世界中倒数相乘的中心思想其实就是用乘法来代替除法运算。

它的原理很简单，就是将被除数乘以除数的倒数，其公式为x/y = x*(1/y)，我们拿10/2作为例子，我可以得出以下推论：

由 公式x/y = x*(1/y) 可得 10/2 = 10*(1/2) = 10*0.5

编译器也正是由这个公式才得以将除法转换为除法，但是编译器为什么要这样做呢？原因同样很简单，因为乘法的运算速度会比除法快4倍左右，在现有的Intel指令集中，就属除法指令最慢，因此将其优化为乘法的理由显得很充分，但是这么做似乎有不妥之处，让我们看看这是为什么……

int _tmain(int argc, _TCHAR* argv[])
{
    int   a = 5/2;
    float b = 5.0/2.0;

  printf("    5/2 = %d
5.0/2.0 = %1.1f",a,b);
  return 0;
}
-Export--------------------------------
    5/2 = 2
5.0/2.0 = 2.5

    

很明显，我们没能充分考虑到浮点类型，一般情况下，在C语言中1除以任何数其结果皆为0，这个问题显得比较严重，我们要怎样才能解决它呢？

倒数相乘与定点运算的配合

为了使除法的倒数相乘优化成为可能，编译器使用了定点运算方案来表示小数。

    那么何为定点运算，定点运算有什么特点呢？

    我们都知道一般情况下我们的小数都是用浮点类型来表示的，有一位会记录当前的小数点位于那里，当然还有其他的一些转换规则，这些都不是本文所关心的，我们只需知道浮点类型的小数可以位于任意一位，也就是说“小数点是浮动的”。

    而定点运算根据字面意思来理解就是“小数点是固定的”，这种小数的定点表示法有很多优点，首当其冲的就是效率上的提高。当然，作为代价，同样也必须承受随之而来的精度上的丢失。

    那么，这个定点表示法又是怎样运作的呢？它怎么就能保存小数信息呢？这部分内容很难寻找，经过大量近似资料的启示与笔者的试验，最终才证实其原理其实很简单，这首先还要从定点表示法的小数点位置与精度说起。

    首先，编译器一般都将小数点定位在第一位，因为要表示一个数的倒数，那么用小数表示的话必然是一个小于1的数，例如8的倒数与12345678的倒数分别为0.125与~0.000000081，因此其整数部分恒为0。

    其次，就是精度问题，例如我们用一个4bit大小的数来表述小数的话，那么它的精度就是0.0625，其表示的整数每加1或减1，其表示的小数就随之增加或减少0.0625，如下所示：

确定精度：
┏━━┳━━━━━┳━━━━━━┳━━━━━━┳━━━━━━━┓
┃Bit   ┃1         ┃2           ┃3           ┃4             ┃
┣━━   ╋━━━━━╋━━━━━━╋━━━━━━╋━━━━━━━┫
┃精度┃2^-1 [0.5]┃2^-2 [0.25] ┃2^-3 [0.125]┃2^-4 [0.0625] ┃
┗━━┻━━━━━┻━━━━━━┻━━━━━━┻━━━━━━━┛

例子：
 0001 = 1 = 1*0.0625 = 0.0625
 0010 = 2 = 2*0.0625 = 0.125
 0101 = 5 = 5*0.0625 = 0.3125

    而对于32位的编译器来说，只不过是将上面的例子的位数扩充了一下，因此精度也就随之提高，我们直接看例子：

精度：2^-32 = 0.00000000023283064365386962890625

例子：
┏━━┳━━━━━━┳━━━━━━━━━━┳━━━━━━━━━┓
┃除数┃倒数(小数)  ┃定点运算 : 小数/精度┃舍去小数点的16进制┃
┣━━╋━━━━━━╋━━━━━━━━━━╋━━━━━━━━━┫
┃3   ┃0.3333333333┃1431655765.190167756┃0x55555555        ┃
┣━━╋━━━━━━╋━━━━━━━━━━╋━━━━━━━━━┫
┃5   ┃0.2         ┃858993459.2000000000┃0x33333333        ┃
┣━━╋━━━━━━╋━━━━━━━━━━╋━━━━━━━━━┫
┃11  ┃0.0909090909┃390451572.3245912064┃0x1745D174        ┃
┣━━╋━━━━━━╋━━━━━━━━━━╋━━━━━━━━━┫
┃59  ┃0.0169491525┃72796055.68241664000┃0x0456C797        ┃
┣━━╋━━━━━━╋━━━━━━━━━━╋━━━━━━━━━┫
┃99  ┃0.0101010101┃43383508.03606568960┃0x0295FAD4        ┃
┗━━┻━━━━━━┻━━━━━━━━━━┻━━━━━━━━━┛

    有了上面的基础，我们就可正式进入除法优化的神秘领地了……

除法运算的识别与优化原理

    大多数情况下除法是最好识别的，因为其特征很明显。另外，在上一节《乘法的识别与优化原理》笔者详细的讲解了编译器利用位移指令(sar、shr)的优化情况，

很明显的，除法必然也会存在这种优化，但是鉴于其基础知识已经讲解且原理比较简单，因此除法的位移优化在本小节中将一带而过，

我们将重点放在倒数相乘相关的优化上。

    按照惯例，请读者阅读以下源代码，并猜测编译器会将其进行怎样的优化：

int _tmain(int argc, _TCHAR* argv[])
{
    printf("03=%d",argc/3);  //-*
    printf("05=%d",argc/5);  // |-> 注意，这些除数的倒数在上一段“1.9.2节”中有提及
    printf("11=%d",argc/11); // |
    printf("59=%d",argc/59); //-*
    printf("04=%d",argc/4);
    printf("64=%d",argc/64);
    return 0;
}

    下面我们先看看它的Debug版部分反汇编代码：
    (再次提醒一下，Debug的反汇编代码并不是完整的，与逻辑无关的部分代码都已经被笔者删除了，如有疑问请查看第1.3节)

.text:004113A0     push ebp
.text:004113A1     mov ebp, esp
.text:004113A3     sub esp, 0C0h
......             ......
.text:004113BE     mov eax, [ebp+argc]
.text:004113C1     cdq                                 ; 此指令的意思是将edx扩展为eax的高位，将其扩展为64位
.text:004113C2     mov ecx, 3                          ; ecx = 3(除数)
.text:004113C7     idiv ecx                            ; edx:eax = eax/ecx = argc/3
.text:004113CB     push eax
.text:004113CC     push offset Format                  ; "03=%d"
.text:004113D1     call ds:__imp__printf
.text:004113D7     add esp, 8
......             ......
.text:004113E1     mov eax, [ebp+argc]
.text:004113E4     cdq
.text:004113E5     mov ecx, 5
.text:004113EA     idiv ecx
.text:004113EE     push eax
.text:004113EF     push offset a05D                    ; "05=%d"
.text:004113F4     call ds:__imp__printf
.text:004113FA     add esp, 8
......             ......
......             ......
.text:0041144A     mov eax, [ebp+argc]
.text:0041144D     cdq                                 ; 扩展为64位 edx:eax
.text:0041144E     and edx, 3                          ; 选取edx的低2位(0x3 = 00000011)
.text:00411451     add eax, edx                        ; 当被除数小于除数时，将结果清零【注1】
.text:00411451                                         ; 
.text:00411451                                         ;  -= 注1 =-
.text:00411451                                         ; (以下为扩展知识，非必须掌握)
.text:00411451                                         ;     这个操作其实其最终目的就是控制特殊情况下的计算结果，
.text:00411451                                         ; 《黑客反汇编揭秘》对edx的解释是“数值符号位”，这是错误
.text:00411451                                         ;的！在此提醒各位读者注意……
.text:00411451                                         ;     我们知道，其实这种除数为正的除法计算无非就四种情况：
.text:00411451                                         ; 1、被除数>除数
.text:00411451                                         ; 2、被除数<除数
.text:00411451                                         ; 3、负的被除数>除数
.text:00411451                                         ; 4、负的被除数<除数
.text:00411451                                         ;     根据整型除法的规则，任何被除数如果除以比自己大的除数
.text:00411451                                         ; 其结果皆为0(例如1/2=0或1/9999=0)。
.text:00411451                                         ;     但是作为编译器来讲，必须要考虑到以上所有情况，按照我
.text:00411451                                         ; 们正常人的思路，既然后面必须使用sar指令，那么我们只需要保
.text:00411451                                         ; 证这种特殊情况下的除法计算在右移后为0即可。
.text:00411451                                         ;     这其实很简单，只需要给被除数加上一个比除数小一的数即
.text:00411451                                         ; 可，原理如下所示：
.text:00411451                                         ; -16/4 = (-16+3)>>2 = -13>>2 = -4
.text:00411451                                         ; -5/4  = (-5+3) >>2 = -2 >>2 = -1
.text:00411451                                         ; -4/4  = (-4+3) >>2 = -1 >>2 = -1
.text:00411451                                         ; -3/4  = (-3+3) >>2 =  0 >>2 =  0
.text:00411451                                         ; -2/4  = (-2+3) >>2 =  1 >>2 =  0
.text:00411451                                         ;     通过上面的几个例子相信各位都能明白之中的原理了，但是
.text:00411451                                         ; 我们并没有将被除数为正数的情况考虑进去。而通过编译器生成
.text:00411451                                         ; 的汇编代码则给了我们很好的提示，就是利用and指令在高位取
.text:00411451                                         ; 值，这也就是为什么要将被除数扩展为64位的根本原因。
.text:00411451                                         ;     通过cdq将eax扩展为64位后，cpu将以edx:eax的方式来表示
.text:00411451                                         ; 数据，那么edx肯定就是高位了，而我们的编译器编译出来的代码
.text:00411451                                         ; 始终都是32位的，因此edx寄存器理论上是根本用不上的。
.text:00411451                                         ;     但是有一种情况下例外，就是被除数为负数时，edx寄存器里
.text:00411451                                         ; 面的值是0xFFFFFFFF(二进制全1)，反之则二进制全0。
.text:00411451                                         ;     利用这个特性编译器使用了与(and)运算来巧妙地达到了目的，
.text:00411451                                         ; 这样一来如果被除数为负数，那么将任何数与edx进行and运算都得
.text:00411451                                         ; 那个数本身，反之则恒为0。
.text:00411451                                         ;     总结起来原理虽非常简单，但也不得不感叹人类智慧的伟大…
.text:00411451
.text:00411451
.text:00411453     sar eax, 2                          ; 带符号右移2位，相当于除以4
.text:00411458     push eax
.text:00411459     push offset a04D                    ; "04=%d"
.text:0041145E     call ds:__imp__printf
.text:00411464     add esp, 8
......             ......
.text:0041146E     mov eax, [ebp+argc]
.text:00411471     cdq
.text:00411472     and edx, 0Fh
.text:00411475     add eax, edx
.text:00411477     sar eax, 4
.text:0041147C     push eax
.text:0041147D     push offset a16D                    ; "16=%d"
.text:00411482     call ds:__imp__printf
.text:00411488     add esp, 8
......             ......
......             ......
.text:004114B6     xor eax, eax
......             ......
.text:004114C8     mov esp, ebp
.text:004114CA     pop ebp
.text:004114CB     retn

    通过上面例子可知，虽然Debug模式下未开启任何优化，不过在除数为2的次方的情况下，编译器仍对其做了不小的优化，那么Release版就更是可想而知了：

.text:00401000     push esi
.text:00401001     mov esi, [esp+4+argc]
.text:00401005     mov eax, 55555556h                  ; 注意这里，比我们在“1.9.2节”中算出来的值大1
.text:0040100A     imul esi                            ; argc*(1/3)，结果保存在edx:eax组成的64位寄存器里
.text:0040100C     mov eax, edx                        ; 去掉小数部分【注2】
.text:0040100C                                         ;
.text:0040100C                                         ; -= 注2 =-
.text:0040100C                                         ;     记得我们在上面学习定点计算时讲过一句话“编译器一
.text:0040100C                                         ; 般都将小数点定位在第一位，因为要表示一个数的倒数，那
.text:0040100C                                         ; 么用小数表示的话必然是一个小于1的数”。
.text:0040100C                                         ;     因此，当imul指令变向将其扩展为64位后，其小数点的
.text:0040100C                                         ; 位置是不会变得，所以始终保证小数点都在一个固定的位置，
.text:0040100C                                         ; 在这里也就是eax里的最高位上。
.text:0040100C                                         ;     由此可知，如果想舍去小数，那么只需要清除eax里的值
.text:0040100C                                         ; 即可。
.text:0040100C                                         ;     当然，纵观汇编指令的上下文可知，这一步操作还有另
.text:0040100C                                         ; 外的一个目的--备份，将edx中的值备份一份，以便于简单高
.text:0040100C                                         ; 效的取出其符号位。
.text:0040100C                                         ;
.text:0040100C
.text:0040100E     shr eax, 1Fh                        ; 取符号位
.text:00401011     push edi
.text:00401012     mov edi, ds:__imp__printf
.text:00401018     add eax, edx                        ; 将符号位加到结果上【注3】
.text:00401018                                         ; 
.text:00401018                                         ; -= 注3 =-
.text:00401018                                         ;     对于这步操作应该会让很多读者感到困惑，为什么要加
.text:00401018                                         ; 上符号位呢？其实很简单，就是为了遵守整数除法的规则。
.text:00401018                                         ; 让我们先看看下面几个小例子：
.text:00401018                                         ; ┏━━┳━━━━━━━━━━━━┳━━━━━━━━━┓
.text:00401018                                         ; ┃算式┃浮点结果[16进制值(64位)]┃整型结果[16进制值]┃
.text:00401018                                         ; ┣━━╋━━━━━━━━━━━━╋━━━━━━━━━┫
.text:00401018                                         ; ┃ 8/3┃ 2.67[00000002 AAAAAAB0]┃ 2[0000 0002]     ┃
.text:00401018                                         ; ┣━━╋━━━━━━━━━━━━╋━━━━━━━━━┫
.text:00401018                                         ; ┃-8/3┃-2.67[FFFFFFFD 55555550]┃-2[FFFF FFFE]     ┃
.text:00401018                                         ; ┗━━┻━━━━━━━━━━━━┻━━━━━━━━━┛
.text:00401018                                         ;     通过上面的例子我们不难发现，其实这个加法只是为了兼
.text:00401018                                         ; 容负数，因为负数的编码规则与正数不同，仅此而已。
.text:00401018                                         ; 
.text:00401018
.text:0040101A     push eax
.text:0040101B     push offset Format                  ; "03=%d"
.text:00401020     call edi ; __imp__printf
.text:00401022     mov eax, 66666667h                  ; 注意这里，比我们在“1.9.2节”中算出来的值大了1倍多
.text:00401027     imul esi
.text:00401029     sar edx, 1                          ; 将整数部分右移1位(既除2)，恢复到原来的数值。
.text:0040102B     mov ecx, edx
.text:0040102D     shr ecx, 1Fh
.text:00401030     add ecx, edx
.text:00401032     push ecx
.text:00401033     push offset a05D                    ; "05=%d"
.text:00401038     call edi ; __imp__printf
.text:0040103A     mov eax, 2E8BA2E9h
.text:0040103F     imul esi
.text:00401041     sar edx, 1
.text:00401043     mov eax, edx
.text:00401045     shr eax, 1Fh
.text:00401048     add eax, edx
.text:0040104A     push eax
.text:0040104B     push offset a11D                    ; "11=%d"
.text:00401050     call edi ; __imp__printf
.text:00401052     mov eax, 22B63CBFh                  ; 比我们在“1.9.2节”中算出来的值大了8倍多
.text:00401057     imul esi
.text:00401059     sar edx, 3                          ; 将整数部分右移3位(既除8)，恢复到原来的数值。
.text:0040105C     mov ecx, edx
.text:0040105E     shr ecx, 1Fh
.text:00401061     add ecx, edx
.text:00401063     push ecx
.text:00401064     push offset a59D                    ; "59=%d"
.text:00401069     call edi ; __imp__printf
.text:0040106B     mov eax, esi
.text:0040106D     cdq                                 ;-*
.text:0040106E     and edx, 3                          ; |
.text:00401071     add eax, edx                        ; |-> 与在Debug版中的优化一模一样，为了节省篇幅后面就省略掉了
.text:00401073     sar eax, 2                          ; |
.text:00401076     push eax                            ;-*
.text:00401077     push offset a04D                    ; "04=%d"
.text:0040107C     call edi ; __imp__printf
......             ......
.text:0040108A     push offset a16D                    ; "16=%d"
......             ......
.text:0040109D     push offset a64D                    ; "64=%d"
......             ......
.text:004010A4     add esp, 38h
.text:004010A7     pop edi
.text:004010A8     xor eax, eax
.text:004010AA     pop esi
.text:004010AB     retn

    后面除数为2的次方的优化没什么可说的，与Debug版里一模一样，但是上面的除法优化可能是有些读者很困惑，这究竟是为什么呢？我们先看看下面这张表：

┏━━┳━━━━━━┳━━━━━━━━━━┳━━━━━━━━━╥━━━━━━━━┳━━━━━━┳━━━━━━┓
┃除数┃倒数(小数)  ┃定点运算 : 小数/精度┃舍去小数点的16进制║反汇编中对应的值┃(1/除数)×2 ┃(1/除数)×8 ┃
┣━━╋━━━━━━╋━━━━━━━━━━╋━━━━━━━━━╫━━━━━━━━╋━━━━━━╋━━━━━━┫
┃3   ┃0.3333333333┃1431655765.190167756┃0x55555555        ║0x55555556      ┃-           ┃-           ┃
┣━━╋━━━━━━╋━━━━━━━━━━╋━━━━━━━━━╫━━━━━━━━╋━━━━━━╋━━━━━━┫
┃5   ┃0.2         ┃858993459.2000000000┃0x33333333        ║0x66666667      ┃0x66666666  ┃-           ┃
┣━━╋━━━━━━╋━━━━━━━━━━╋━━━━━━━━━╫━━━━━━━━╋━━━━━━╋━━━━━━┫
┃11  ┃0.0909090909┃390451572.3245912064┃0x1745D174        ║0x2E8BA2E9      ┃0x2E8BA2E8  ┃-           ┃
┣━━╋━━━━━━╋━━━━━━━━━━╋━━━━━━━━━╫━━━━━━━━╋━━━━━━╋━━━━━━┫
┃59  ┃0.0169491525┃72796055.68241664000┃0x0456C797        ║0x22B63CBF      ┃-           ┃0x22B63CBE  ┃
┗━━┻━━━━━━┻━━━━━━━━━━┻━━━━━━━━━╨━━━━━━━━┻━━━━━━┻━━━━━━┛

    由上表可知，不管我们反汇编中体现出来的值最终比实际计算出来的值大几倍，而且都是2的次方，例如11的倒数对应的是大2倍，而59对应的是大8倍。

    其次反汇编中体现出来的值不管比原值大几倍或是大致没变，都要比原值大1。

    着看起来似乎有些怪异，我们先想想后面这种情况，想想为什么非要大个1呢？它为什么不大2或100呢？这其实不难理解，很明显的与精度有关，我们拿除数为59的定点运算结果72796055.68241664000来说，在将其转为16进制后，它的小数会直接被舍去，而按照四舍五
入的原则其实是应该加1的，但是实际并没有这样做。

    编译器在进行优化时考虑到了这点，为了弥补后面丢失的小数精度，所以在最终生成代码时就统一都给加了一个1。

    第二种加1的问题解决了，思维活跃的读者这时应该可以猜到第一种情况也是精度问题了。事实也正像读者们所想象的那样，将原定点表示法计算出来的值乘以一个2的倍数，就是出于精度的考虑。

    让我们共同回忆一下定点表示法的特点：

(1) 计算更加简便快速
(2) 有精度问题

    很明显的，当我们试图用定点表示法表示越小的小数时，其精度问题就越明显，就拿一个4位的定点小数来讲，其最小精度为0.0625，也就是说如果我们想表示0.05的话是不行的，那么除了增加位数以外，还有没有什么其他方法呢？答案是肯定的！

    做法很简单，我们只需将0.05乘以5，那么我们就可以用0x4来表示它了，而且精度一丁点也没丢失，再次使用的时候只需再将其除以5即可。

    不过令人沮丧的，上面的思路虽然原理是对的，但是实际情况要复杂一些，因为编译器最终是要将除法优化掉的，所以我们必须保证在数据还原时不能出现除法，这很明显要用到位移了，而位移的特点是只能转换除数为2的次方的除法，因此这就需要我们保证我们在变幻小数时将其增加的倍数也必须是2的次方。

    我们拿1876523938/876523938这个例子来讲，这次我们要将其乘以4以期减少误差，其运算逻辑如下：

参考数据：

    此算式的浮点结果    = 2.140870
    此算式的整除结果    = 2
    1876523938          = 0x6FD97BA2
    1/876523938定点表示 = 0x00000004
    4/876523938定点表示 = 0x00000013

以普通方式计算：

    1876523938/876523938 = 1876523938 * (1/876523938) = 0x6FD97BA2 * 0x00000004 = 0x00000001 BF65EE88
    0x00000001 E0C00000 取整 = 0x00000001 = 1
    结果为1
                 
以编译器优化后方式计算：

    1876523938/876523938 = 1876523938 * (4/876523938) = 0x6FD97BA2 * 0x00000013 = 0x00000008 4D242D06
    0x00000008 4D242D06 取整 = 0x00000008
    0x00000008 >> 2 = 0x00000002 = 2
    结果为2

    由此可见，这种计算前增加倒数值(称之为‘倒数向上圆整’)后，在将结果减去对应倍数的值(称之为‘商向下圆整’)的方法可以比较完美的解决精度问题。

    到此，我们有关于除法逆向的学习就可以告一段落了，下面我们就趁热打铁，再看看另一个与除法相关的知识点。