在SIMD功能问世之前,x86兼容的处理器早就集成了FPU。自SSE2开始,SIMD拓展指令集提供了单指令多数据浮点计算的指令。这类指令支持的数据格式同样是IEEE 754。
随着硬件的发展,最近推出的x86/x86-64编译器越来越多地使用SIMD指令,都不再怎么分配FPU指令了。
这不得不说是一种好消息。毕竟SIMD的指令更为简单。
本章将基于第17章的源程序进行讲解。
#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 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编译出的程序与之相同,本文不再介绍。
在使用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所示。
图27.1 OllyDbg: MOVSD指令把变量a传递给XMM1
图27.2 OllyDbg:DIVSD进行除法运算、把商存储于XMM1
图27.3 OllyDbg:MULSD进行乘法运算、把积存储于XMM0
图27.4 OllyDbg:ADDSD计算XMM0和XMM1的和
从图27.5可以看到:虽然OllyDbg把单个XMM寄存器视为一对double数据的寄存器,但是它只使用了寄存器的低半部分。很明显,OllyDbg根据SSE2指令的后缀“-SD”判断出了数据类型,并据此进行了数据整理。实际上,OllyDbg还能够根据SSE指令进行判断,把XMM寄存器的数据整理为4个float浮点数、或者是16个字节。
图27.5 OllyDbg:FLD把函数返回值传递给ST(0)
#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)。
#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.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
使用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
这个程序用栈向函数传递变量a和b,并且把函数返回值存储在ST(0)寄存器里。除此之外它和64位的程序没有区别。
使用OllyDbg调试这个程序,则可以观察COMISD指令比较数值、设置/清零CF/PF标识位的具体过程,如图27.6所示。
图27.6 OllyDbg:COMISD指令调整CF和PF标识位
本书在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寄存器中。
现在,我们再次利用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型的数据。
在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。