第41章 除以9

我们来看一个非常简单的函数:

int f(int a)
{
        return a/9;
};

41.1 x86

x86平台编译器的编译方法十分直白:

指令清单41.1 MSVC

_a$ = 8               ; size = 4
_f     PROC
     push   ebp
     mov    ebp, esp
     mov    eax, DWORD PTR _a$[ebp]
     cdq              ; sign extend EAX to EDX:EAX
     mov    ecx, 9
     idiv   ecx
     pop    ebp
     ret    0
_f   ENDP

IDIV指令是除法指令。它会从寄存器对EDX:EAX中提取被除数、从ECX寄存器中提取除数。计算结束以后,它把计算结果/商存储在EAX寄存器里,把余数存储在EDX寄存器。除法计算之后,商就位于EAX寄存器里,直接成为f()函数的返回值;因此没有其他值传递的操作。为了通知IDIV指令从EDX:EAX寄存器对中提取64位被除数,编译器在IDIV指令之前分派了CDQ指令。IDIV指令就会进行MOVSX那样的符号位处理和数据扩展处理。

启用编译器的优化选项之后,可得到下述程序:

指令清单41.2 采用MSVC 优化

_a$ = 8                       ; size = 4
_f     PROC
     mov     ecx, DWORD PTR _a$[esp-4]
     mov     eax, 954437177   ; 38e38e39H
     imul    ecx
     sar     edx, 1
     mov     eax, edx
     shr     eax, 31          ; 0000001fH
     add     eax, edx
     ret     0
_f     ENDP

编译器用乘法指令来变相实现除法运算。大家知道乘法会比除法运行快很多。使用这里介绍[1]的方法可以有效提高程序效率、节省时间开销。

在编译优化中,我们常常将其称为“强度减轻”的办法。

若使用GCC 4.4.1编译此程序,即使我们刻意不启用其优化选项,GCC的非优化编译结果也足以和MSVC的优化编译结果媲美。

指令清单41.3 不带优化的GCC 4.4.1

      public f
f     proc near

arg_0 = dword ptr  8

      push    ebp
      mov     ebp, esp
      mov     ecx, [ebp+arg_0]
      mov     edx, 954437177 ; 38E38E39h
      mov     eax, ecx
      imul    edx
      sar     edx, 1
      mov     eax, ecx
      sar     eax, 1Fh
      mov     ecx, edx
      sub     ecx, eax
      mov     eax, ecx
      pop     ebp
      retn
f     endp

41.2 ARM

ARM处理器和其他的RISC处理器一样,“纯洁”得不支持硬件级别的除法指令。此外,这种CPU还不难“直接”进行32位常量的乘法运算(32位opcode容纳不下32位常量)。因此,在进行除法运算时,编译器会混合加减法运算和位移运算、变相实现除法运算(详情请参阅第19章)。

本节引用参考书目Ltd94(第3.3节)的一个例子,介绍一个“32位数除以10”的例子,分别计算商和余数。

; takes argument in a1
; returns quotient in a1, remainder in a2
; cycles could be saved if only divide or remainder is required
    SUB    a2, a1, #10               ; keep (x-10) for later
    SUB    a1, a1, a1, lsr #2
    ADD    a1, a1, a1, lsr #4
    ADD    a1, a1, a1, lsr #8
    ADD    a1, a1, a1, lsr #16
    MOV    a1, a1, lsr #3
    ADD    a3, a1, a1, asl #2
    SUBS   a2, a2, a3, asl #1        ; calc (x-10) - (x/10)*10
    ADDPL  a1, a1, #1                ; fix-up quotient
    ADDMI  a2, a2, #10               ; fix-up remainder
    MOV    pc, lr
41.2.1 ARM模式下,采用Xcode 4.6.3(LLVM)优化
__text:00002C58 39 1E 08 E3 E3 18 43 E3 MOV      R1, 0x38E38E39
__text:00002C60 10 F1 50 E7             SMMUL    R0, R0, R1
__text:00002C64 C0 10 A0 E1             MOV      R1, R0,ASR#1
__text:00002C68 A0 0F 81 E0             ADD      R0, R1, R0,LSR#31
__text:00002C6C 1E FF 2F E1             BX       LR

这里的代码和采用优化算法时的MSVC与GCC基本相同。很明显,LLVM采用了相同的算法来处理常数。

细心的读者可能会问:既然ARM模式的单条指令不能把32位立即数赋值给寄存器,那么这个程序又是怎样做到单条MOV指令赋值的呢?实际上,原始指令并非是IDA显示的那种单条MOV指令。仔细观察您就会发现,那“条”指令占用了8个字节,而标准的ARM指令只有4个字节。原始的指令分两步进行32位赋值:首先用MOV指令将低16位(本例是常量0x8E39)复制到寄存器的低16位,再用MOVT指令把立即数的高16位复制到寄存器的高16位。IDA能够识别出这种指令组合,为了便于读者理解,把两条指令“排版”为一条“伪指令”。这8字节的MOV指令实际上是2条指令。

SMMUL是Signed Most Significant Word Multiply的简称。它是2个32位有符号数的乘法运算指令,会把64位结果的高32位保存在寄存器R0中,舍弃结果中的低32位。

MOV R1,R0,ASR#1是算术右移1位的运算指令。

而指令ADD R0,R1,R0,LSR#31的执行结果相当于将R0的值右移31位,并与R1的值相加,其结果保存在R0中。也就是:R0=R1+R0>>31。

在ARM模式的指令中没有单独的位移指令。不过,它可以在MOV、ADD、SUB以及RSB[2]指令中,使用“后缀”形式的参数调节符对第二个操作数进行位移运算。在使用位移调节符的时候,应当指定位移的确切位数。

ASR是算术右移Arithmetic Shift Right的简称。算术右移需要考虑符号位。

LSR是逻辑右移Logical Shift Right的简称。逻辑右移不考虑符号位。

41.2.2 Thumb-2模式下的Xcode 4.6.3优化(LLVM)
MOV             R1, 0x38E38E39
SMMUL.W         R0, R0, R1
ASRS            R1, R0, #1
ADD.W           R0, R1, R0,LSR#31
BX              LR

Thumb模式的指令集里有单独的位移运算指令。本例中的ASRS就是算术右移指令。

41.2.3 非优化的Xcode 4.6.3(LLVM) 以及Keil 6/2013

在没有启用优化选项的情况下,LLVM编译器不会生成上面那种混合运算指令,它会调用仿真库里的模拟运算函数__divsi3。

然而,无论是否启用优化选项,Keil编译器都只会调用库函数__aeabi_idivmod。

41.3 MIPS

出于某些原因,优化的GCC 4.4.5有除法指令。

指令清单41.4 优化的GCC 4.4.5(IDA)

f:
                li      $v0, 9
                bnez    $v0, loc_10
                div     $a0, $v0 ; branch delay slot
                break   0x1C00    ; "break 7" in assembly output and objdump

loc_10:
                mflo    $v0
                jr      $ra
                or      $at, $zero ; branch delay slot, NOP

本例出现了新指令BREAK。它是由编译器产生的异常处理指令,会在除数为零的情况下抛出错误信息。毕竟在正规的数学概念里,除数不可以是零。即使在启用优化选项的情况下,GCC也未能判断出本例“除数$V0用于不会是零”的切实情况,机械性地分派了异常处理的检测指令。BREAK指令只是一个在除数为零的情况下通知操作系统进行异常处理的指令。只要除数不是零,程序将会执行MFLO指令,把LO寄存器里的商复制到$V0寄存器。

这里顺便说明一下,乘法指令MUL会把积的高32位保存在HI寄存器中,把积的低32位保存在LO寄存器中。而除法指令DIV则会把商保存在LO寄存器,把余数在保存在HI寄存器中。

如果把源程序的除法指令修改为“a % 9”、计算余数,那么编译之后的程序就会用MFHI指令替换本例中的MFLO指令。

41.4 它是如何工作的

在引入2的n次方之后,除法运算可转换为乘法运算:

 \text{result}=\frac{\text{input}}{\text{divisor}}=\frac{\text{input}\cdot \frac{2^{n}}{\text{divisor}}}{2^{n}}=\frac{\text{input}\cdot M}{2^{n}}

这里的M是魔术因子(Magic coefficient)。其计算公式是

 M=\frac{2^{n}}{\text{divisor}}

最终,除法运算就转换成了

\text{result}=\frac{\text{input}\bullet M}{2^{n}}

“除以2的n次方”的运算可以直接通过右移操作实现。如果n小于32,那么中间运算的积的低n位(通常位于EAX或者RAX寄存器)就会被位移运算直接抹去;如果n大于或等于32,那么积的高半部分(通常位于EDX或者RDX寄存器)的数值都会受到影响。

可见,参数n的取值直接决定了转换运算的计算精度。

在进行有符号数的除法运算时,符号位也对计算精度及n的取值产生了显著影响。

下面这个例子将验证符号位的影响。

int f3_32_signed(int a)
{
         return a/3;
};

unsigned int f3_32_unsigned(unsigned int a)
{
         return a/3;
};

在无符号数的计算过程中,魔术因子是0xaaaaaaab。乘法的中间结果要除以2的33次方。

而在有符号数的计算过程中,魔术因子则是0x55555556,乘法的中间结果要除以2的32次方。虽然这里没有进行除法运算,但是根据前面的讨论可知:商的有效位取自于EDX寄存器。

请别忘记中间一步的乘法计算同样存在符号位的问题:积的高32位右移31位,将在EAX寄存器的最低数权位保留有符号数的符号位(正数为0,负数为1)。将符号位加入积的高32位值,可实现负数补码的“+1”修正

指令清单41.5 带优化的MSVC 2012

_f3_32_unsigned PROC
         mov      eax, -1431655765     ; aaaaaaabH
         mul      DWORD PTR _a$[esp-4] ; unsigned multiply
; EDX=(input*0xaaaaaaab)/2^32
         shr      edx, 1
; EDX=(input*0xaaaaaaab)/2^33
         mov      eax, edx
         ret      0
_f3_32_unsigned ENDP

_f3_32_signed PROC
         mov      eax, 1431655766              ; 55555556H
         imul     DWORD PTR _a$[esp-4] ; signed multiply
; take high part of product
; it is just the same as if to shift product by 32 bits right or to divide it by 2^32
         mov      eax, edx       ; EAX=EDX=(input*0x55555556)/2^32
         shr      eax, 31        ; 0000001fH
         add      eax, edx       ; add 1 if sign is negative
         ret      0
_f3_32_signed ENDP
41.4.1 更多的理论

其实,我们都知道,乘法和除法互为逆运算。因此下面的除法可以用乘法来代替。我们可以表述为

 \frac{x}{c}=x\frac{1}{c}

1/c可以称为乘法的逆运算,是c的倒数。可以用编译器来做提前运算。

但是这是为浮点计算用的,有没有整数呢?在模算术计算环境下,是可能会有的。CPU寄存器的长度是很规整的,要么32位要么64位,因此几乎所有的寄存器的算术操作都是对2的32次方或者2的64次方进行操作。

可以查阅War02的第10章第3节。

41.5 计算除数

41.5.1 变位系数#1

通常我们看到的代码可能就像这样:

        mov        eax, MAGICAL CONSTANT
        imul       input value
        sar        edx, SHIFTING COEFFICIENT ; signed division by 2xusing arithmetic shift right
        mov        eax, edx
        shr        eax, 31
        add        eax, edx

我们这里用大写字母M来代表32位的魔术因子,把变位系数记为C,把除数记为D

因此除数可以表示为

 D=\frac{2^{32+c}}{M}

指令清单41.6 优化的MSVC 2012代码

        mov        eax, 2021161081                      ; 78787879H
        imul       DWORD PTR _a$[esp-4]
        sar        edx, 3
        mov        eax, edx
        shr        eax, 31                              ; 0000001fH
        add        eax, edx

用公式可以表示为

 D=\frac{2^{32+3}}{2021161081}

这个数字超过了32位的表达范围。我们可使用Mathematica程序计算除数:

指令清单41.7 Wolfram Mathematica的计算结果

In[1]:=N[2^(32+3)/2021161081]
Out[1]:=17.

它可算出本例采用的除数是17。

在64位数据的除法运算中,计算除数的方法完全相同。只是不再采用2的32次方,而采用了2的64次方。

uint64_t f1234(uint64_t a)
{
         return a/1234;
};

指令清单41.8 64位下的优化MSVC2012

f1234    PROC
         mov     rax, 7653754429286296943              ; 6a37991a23aead6fH
         mul     rcx
         shr     rdx, 9
         mov     rax, rdx
         ret     0
f1234    ENDP

指令清单41.9 WolframMathematica的计算结果

In[1]:=N[2^(64+9)/16^^6a37991a23aead6f]
Out[1]:=1234.
41.5.2 变位系数#2

变位系数确实可能为零:

                mov      eax, 55555556h ; 1431655766
                imul     ecx
                mov      eax, edx
                shr      eax, 1Fh

这样,计算除数的方法就更简单一些了:

 D=\frac{2^{32}}{M}

就本例而言,除数的计算方法是:

 D=\frac{2^{32}}{1431655766}

再次使用Mathematica程序计算除数

指令清单41.10 Wolfram Mathematica的计算结果

In[1]:=N[2^32/16^^55555556]
Out[1]:=3.

最终求得除数为3。

41.6 练习题

请描述下述代码的功能。

指令清单41.11 采用MSVC 2010优化的代码

_a$ = 8
_f        PROC
          mov    ecx, DWORD PTR _a$[esp-4]
          mov    eax, -968154503 ; c64b2279H
          imul   ecx
          add    edx, ecx
          sar    edx, 9
          mov    eax, edx
          shr    eax, 31           ; 0000001fH
          add    eax, edx
          ret    0
_f        ENDP

指令清单41.12 ARM64位下采用GCC 4.9优化

f:
          mov    w1, 8825
          movk   w1, 0xc64b, lsl 16
          smull  x1, w0, w1
          lsr    x1, x1, 32
          add    w1, w0, w1
          asr    w1, w1, 9
          sub    w0, w1, w0, asr 31
          ret

答案请参见G.1.14。


[1] 可以看War02 pp.10-3中介绍的用乘法来代替除法的部分。

[2] 这些指令也称为“数据处理指令”。