第27章 SIMD与浮点数的并行运算

在SIMD功能问世之前,x86兼容的处理器早就集成了FPU。自SSE2开始,SIMD拓展指令集提供了单指令多数据浮点计算的指令。这类指令支持的数据格式同样是IEEE 754。

随着硬件的发展,最近推出的x86/x86-64编译器越来越多地使用SIMD指令,都不再怎么分配FPU指令了。

这不得不说是一种好消息。毕竟SIMD的指令更为简单。

本章将基于第17章的源程序进行讲解。

27.1 样板程序

#include <stdio.h>

double f (double a, double b)
{
          return a/3.14 + b*4.1;
};

int main() 
{
          printf ("%f\n", f(1.2, 3.4));
};
27.1.1 x64

指令清单27.1 Optimizing MSVC 2012 x64

__real@4010666666666666 DQ 04010666666666666r   ; 4.1
__real@40091eb851eb851f DQ 040091eb851eb851fr   ; 3.14

a$ = 8
b$ = 16
f        PROC
         divsd   xmm0, QWORD PTR __real@40091eb851eb851f
         mulsd   xmm1, QWORD PTR __real@4010666666666666
         addsd   xmm0, xmm1
         ret     0
f        ENDP

程序使用XMM0~XMM3传递参数中的前几个浮点数,其余参数都通过栈传递。有关浮点数的参数传递规范,请参见MSDN:https://msdn.microsoft.com/en-us/library/zthk2dkh.aspx

变量a存储在XMM0寄存器里,变量b存储在XMM1寄存器里。XMM系列寄存器都是128位寄存器,而双精度double型浮点数都是64位数据。所以这个程序只使用了寄存器的低半部分。

DIVSD是SSE指令,它的全称是“Divide Scalar Double-Precision Floating-Point Values”,即标量双精度浮点除法。DIVSD指令以第一个操作数中的低64位中的双精度浮点数做被除数,以第二操作数的低64位中的双精度浮点数做除数进行除法运算,商将保存到目标地址(第一操作数)的低64位中,目标地址的高64位保存不变。

在上述程序中可以看到,编译器以IEEE 754格式封装浮点数常量。

MULSD和ADDSD分别是乘法和加法的运算指令。

上述函数返回的数据类型是double,由XMM0寄存器回传。

由MSVC编译器进行非优化编译,可得如下所示的代码。

指令清单27.2 MSVC 2012 x64

__real@4010666666666666 DQ 04010666666666666r   ; 4.1 
__real@40091eb851eb851f DQ 040091eb851eb851fr   ; 3.14

a$ = 8
b$ = 16

f        PROC
         movsdx  QWORD PTR [rsp+16], xmm1
         movsdx  QWORD PTR [rsp+8], xmm0
         movsdx  xmm0, QWORD PTR a$[rsp]
         divsd   xmm0, QWORD PTR __real@40091eb851eb851f
         movsdx  xmm1, QWORD PTR b$[rsp]
         mulsd   xmm1, QWORD PTR __real@4010666666666666
         addsd   xmm0, xmm1
         ret     0
f        ENDP

上述程序还可以进行优化。请注意:传入的参数被保存到了“阴影空间”(可参见本书8.2.1节)。但是只有寄存器的低半部分数据——即double型64位数据会被存在阴影空间里。

GCC编译出的程序与之相同,本文不再介绍。

27.1.2 x86

在使用MSVC 2012编译x86平台的可执行程序时,编译器会分配SSE2指令。

指令清单27.3 Non-optimizing MSVC 2012 x86

tv70 = -8       ;size=8
_a$ = 8         ;size=8
_b$ = 16        ;size=8
_f      PROC
        push    ebp
        mov     ebp, esp
        sub     esp, 8
        movsd   xmm0, QWORD PTR _a$[ebp]
        divsd   xmm0, QWORD PTR __real@40091eb851eb851f
        movsd   xmm1, QWORD PTR _b$[ebp]
        mulsd   xmm1, QWORD PTR __real@4010666666666666
        addsd   xmm0, xmm1
        movsd   QWORD PTR tv70[ebp], xmm0
        fld     QWORD PTR tv70[ebp]
        mov     esp, ebp
        pop     ebp
        ret     0
_f      ENDP

指令清单27.4 Optimizing MSVC 2012 x86

tv67 = 8        ; size = 8
_a$ = 8         ; size = 8
_b$ = 16        ; size = 8
_f      PROC
        movsd   xmm1, QWORD PTR _a$[esp-4]
        divsd   xmm1, QWORD PTR __real@40091eb851eb851f
        movsd   xmm0, QWORD PTR _b$[esp-4]
        mulsd   xmm0, QWORD PTR __real@4010666666666666
        addsd   xmm1, xmm0
        movsd   QWORD PTR tv67[esp-4], xmm1
        fld     QWORD PTR tv67[esp-4]
        ret     0
_f      ENDP

这种x86程序与前面介绍的64位程序十分相似。它们的区别主要体现在参数的调用规范上:

① x86程序会像第17章的FPU例子那样使用栈来传递浮点数,而不是像64位程序那样使用XMM寄存器传递浮点型参数。

② x86程序的浮点型数据的运算结果保存在ST(0)寄存器里——这个值是通过局部变量tv、从XMM寄存器复制到ST(0)寄存器的。

现在,我们使用OllyDbg调试前面那个优化编译的32位程序,如图27.1到图27.4所示。

..\tu\2701.tif{}

图27.1 OllyDbg: MOVSD指令把变量a传递给XMM1

..\tu\2702.tif{}

图27.2 OllyDbg:DIVSD进行除法运算、把商存储于XMM1

..\tu\2703.tif{}

图27.3 OllyDbg:MULSD进行乘法运算、把积存储于XMM0

..\tu\2704.tif{}

图27.4 OllyDbg:ADDSD计算XMM0和XMM1的和

从图27.5可以看到:虽然OllyDbg把单个XMM寄存器视为一对double数据的寄存器,但是它只使用了寄存器的低半部分。很明显,OllyDbg根据SSE2指令的后缀“-SD”判断出了数据类型,并据此进行了数据整理。实际上,OllyDbg还能够根据SSE指令进行判断,把XMM寄存器的数据整理为4个float浮点数、或者是16个字节。

..\tu\2705.tif{}

图27.5 OllyDbg:FLD把函数返回值传递给ST(0)

27.2 传递浮点型参数

#include <math.h>
#include <stdio.h>

int main () 
{
        printf ("32.01 ^ 1.54 = %lf\n", pow (32.01,1.54));
        return 0;
}

在函数间负责传递浮点数的寄存器是XMM0~XMM3寄存器。

指令清单27.5 Optimizing MSVC 2012 x64

$SG1354 DB       '32.01 ^ 1.54 = %lf', 0aH, 00H

__real@40400147ae147ae1 DQ 040400147ae147ae1r    ; 32.01
__real@3ff8a3d70a3d70a4 DQ 03ff8a3d70a3d70a4r    ; 1.54

Main    PROC
        sub      rsp, 40                                ; 00000028H
        movsdx   xmm1, QWORD PTR __real@3ff8a3d70a3d70a4
        movsdx   xmm0, QWORD PTR __real@40400147ae147ae1
        call     pow
        lea      rcx, OFFSET FLAT:$SG1354
        movaps   xmm1, xmm0
        movd     rdx, xmm1
        call     printf
        xor      eax, eax
        add       rsp, 40                               ; 00000028H
        ret      0
main    ENDP

无论是Intel指令白皮书[1]还是AMD指令白皮书[2],都没有记载MOVSDX指令。它就是MOVSD指令。这就是说x86指令集的指令可能对应着多个名称(详情请参阅附录A.6.2)。而此处出现了这个指令,说明微软的研发人员希望用指令名称来标明操作符的数据类型,避免产生混淆。它是把浮点值传递给XMM寄存器的低半部分的指令。

pow()函数从XMM0和XMM1中读取参数,再把返回结果存储在XMM0寄存器。函数返回值又刻意通过RDX寄存器传递给printf()函数,不过这是为什么?坦白地讲,我不知道个中缘由。或许这是printf()的特性——可受理不同类型参数——造成的。

指令清单27.6 Optimizing GCC 4.4.6 x64

.LC2:
         .string "32.01 ^ 1.54 = %lf\n"
main:
         sub     rsp, 8
         movsd   xmm1, QWORD PTR .LC0[rip]
         movsd   xmm0, QWORD PTR .LC1[rip]
         call    pow
         ; result is now in XMM0
         mov     edi, OFFSET FLAT:.LC2
         mov     eax, 1 ; number of vector registers passed
         call    printf
         xor     eax, eax
         add     rsp, 8
         ret

.LC0:
         .long   171798692
         .long   1073259479
.LC1:
         .long   2920577761
         .long   1077936455

GCC的编译手法更易懂一些。它使用XMM0寄存器向printf()函数传递浮点型参数。此外,在向printf()函数传递参数的时候,程序将EAX寄存器的值设置为1。这就相当于告诉函数“有1个参数在矢量寄存器里”。这种特例完全遵循了有关规范(请参阅Mit13)。

27.3 浮点数之间的比较

#include <stdio.h>

double d_max (double a, double b)
{
         if (a>b)
                  return a;
         return b; 
};

int main()
{
         printf ("%f\n", d_max (1.2, 3.4));
         printf ("%f\n", d_max (5.6, -4));
};
27.3.1 x64

指令清单27.7 Optimizing MSVC 2012 x64

a$ = 8
b$ = 16
d_max   PROC
        comisd      xmm0, xmm1
        ja          SHORT $LN2@d_max
        movaps      xmm0, xmm1
$LN2@d_max:
        fatret      0
d_max   ENDP

MSVC的优化编译结果简明易懂。

其中,COMISD的全称是“Compare Scalar Ordered Double-Precision Floating-Point Values and Set EFLAGS”。它是比较双精度标量并设置标志位的指令。

MSVC的非优化编译结果同样不难阅读,只是程序的效率低了一些。

指令清单27.8 MSVC 2012 x64

a$ = 8
b$ = 16
d_max   PROC
        movsdx   QWORD PTR [rsp+16], xmm1
        movsdx   QWORD PTR [rsp+8], xmm0
        movsdx   xmm0, QWORD PTR a$[rsp]
        comisd   xmm0, QWORD PTR b$[rsp]
        jbe      SHORT $LN1@d_max
        movsdx   xmm0, QWORD PTR a$[rsp]
        jmp      SHORT $LN2@d_max
$LN1@d_max:
        movsdx   xmm0, QWORD PTR b$[rsp]
$LN2@d_max:
        fatret   0
d_max   ENDP

然而GCC 4.4.6优化得更为彻底。它直接使用MAXSD指令(Return Maximum Scalar Double-Precision Floating-Point Value),一步就可完成任务。

指令清单27.9 Optimizing GCC 4.4.6 x64

d_max:
          maxsd   xmm0, xmm1
          ret
27.3.2 x86

使用MSVC 2012对上述代码进行优化编译,可得到如下所示的代码。

指令清单27.10 Optimizing MSVC 2012 x86

_a$ = 8          ; size = 8
_b$ = 16         ; size = 8
_d_max  PROC
        movsd    xmm0, QWORD PTR _a$[esp-4]
        comisd   xmm0, QWORD PTR _b$[esp-4]
        jbe      SHORT $LN1@d_max
        fld      QWORD PTR _a$[esp-4]
        ret      0
$LN1@d_max:
        fld      QWORD PTR _b$[esp-4]
        ret      0
_d_max  ENDP

这个程序用栈向函数传递变量ab,并且把函数返回值存储在ST(0)寄存器里。除此之外它和64位的程序没有区别。

使用OllyDbg调试这个程序,则可以观察COMISD指令比较数值、设置/清零CF/PF标识位的具体过程,如图27.6所示。

..\tu\2706.tif{}

图27.6 OllyDbg:COMISD指令调整CF和PF标识位

27.4 机器精度

本书在22.2.2节就介绍过机器精度的具体计算算法。在这一节里,我们把它编译为如下所示的64位应用程序。

指令清单27.11 Optimizing MSVC 2012 x64

v$ = 8
calculate_machine_epsilon PROC
         movsdx   QWORD PTR v$[rsp], xmm0
         movaps   xmm1, xmm0
         inc      QWORD PTR v$[rsp]
         movsdx   xmm0, QWORD PTR v$[rsp]
         subsd    xmm0, xmm1
         ret      0
calculate_machine_epsilon ENDP

由于INC指令无法对128位的XMM寄存器里的值进行操作,程序首先把寄存器的值传递到了内存中的某个地址。

然而这也不是必须如此。如果在此处直接使用 ADDSD(Add Scalar Double-Precision Floating-Point Values)指令,就能够直接对XMM寄存器的低64位进行加法运算,同时保持寄存器的高64位不变。或许,我们不得不说MSVC 2012还有改善的余地。

我们再看这个程序。在完成递增运算之后,程序把返回值再次回传给了XMM寄存器,然后对XMM寄存器进行减法(比较)操作。SUBSD的全称是“Substract Scalar Double-Precision Floating-Point Values”。它只对128位XMM寄存器的低64位进行操作。最后,减法运算的返回值保存在XMM0寄存器中。

27.5 伪随机数生成程序(续)

现在,我们再次利用MSVC 2012编译指令清单22.1的源程序。这个版本的MSVC编译器能够直接分配FPU的SIMD指令。

指令清单27.12 Optimizing MSVC 2012

__real@3f800000 DD 03f800000r              ; 1

tv128 = -4
_tmp$ = -4
?float_rand@@YAMXZ PROC
         push   ecx
         call   ?my_rand@@YAIXZ
; EAX=pseudorandom value
         and    eax, 8388607               ; 007fffffH
         or     eax, 1065353216            ; 3f800000H
; EAX=pseudorandom value & 0x007fffff | 0x3f800000
; store it into local stack:
         mov    DWORD PTR _tmp$[esp+4], eax
; reload it as float point number:
         movss  xmm0, DWORD PTR _tmp$[esp+4]
; subtract 1.0:
         subss  xmm0, DWORD PTR __real@3f800000
; move value to ST0 by placing it in temporary variable...
         movss  DWORD PTR tv128[esp+4], xmm0
; ... and reloading it into ST0:
         fld    DWORD PTR tv128[esp+4]
         pop    ecx
         ret    0
?float_rand@@YAMXZ ENDP

编译器大量使用了带有“-SS”后缀的指令。这个后缀是“Scalar Single”的缩写。“Scalar(标量)”表明寄存器里只有一个值,“Single”则表明该值是单精度float型的数据。

27.6 总结

在XMM寄存器里存储IEEE 754格式格式浮点数的时候,本章的所有程序都只用到了128位寄存器空间的低64位。

这是由操作指令决定的。本章中的指令都带有-SD(Scalar Double-Precision)后缀。因此,所有指令操作的数据类型都是IEEE 754格式的双精度浮点数。这种数据只会占用XMM寄存器的低64位。

SIMD指令集比FPU的指令更易理解。过去的那些FPU指令距离人类语言较远,而且它们还往往涉及FPU寄存器的栈结构,远不如SIMD指令直观。

如果把本章源程序的数据类型都从double改为float,那么操作浮点数的指令都将改为-SS后缀(Scalar Single-Precision)。您将看到MOVSS/COMISS/ADDSS等指令。

“Scalar”就是“标量”的意思,它表示寄存器里只含有一个值。如果一个寄存器里存有多个值,那么对这些寄存器进行操作的指令的名字里都会带有“packed”字样。

还需要强调的是,SSE2指令的操作数都是IEEE 754格式的64位数据,而FPU内部则使用80位的数据格式存储数据。所以FPU的误差较小、精度较高。


[1] 请参阅Int13。

[2] 请参阅AMD13a。