第25章 SIMD

SIMD意为“单指令流多数据流”,其全称是“Single Instruction, Multiple Data”。顾名思义,这类指令可以一次处理多个数据。在x86的CPU里,SIMD子系统和FPU都由专用模块实现。

不久之后,x86 CPU通过MMX指令率先整合了SIMD的运算功能。支持这种技术的CPU里都有8个64位寄存器,即MM0~MM7。每个MMX寄存器都是8字节寄存器,可以容纳2个32位数据,或者4个16位数据。使用SIMD指令进行操作数的计算时,它可以把8个8位数据分为4组数据同时运算。

平面图像可视为一种由二维数组构成的数据结构。在美工人员调整图像亮度的时候,图像编辑程序得对每个像素的亮度系数进行加减法的运算。简单地说,图像的每个像素都有灰度系数,而且灰度系统是8位数据(1个字节)。因而,每执行一个SIMD指令就能够同时调整8个像素的灰度(即亮度)。为了满足这种需要,处理器厂商后来专门推出了基于SIMD技术的饱和度调整指令。这种饱和度调整指令甚至有越界保护功能,能够避免亮度调整时可能会产生的因子上溢(overflow)和下溢(underflow)等问题。

在刚刚推出MMX技术的时候,SIMD借用了FPU的寄存器(别名)。这有利于FPU或MMX的混合运算。有人说Intel这样处理是为了节约晶体管,但是其实际原因更为简单:老式操作系统并不会利用新增的CPU寄存器,还是会把数据保存到FPU上。所以说,只有同时满足“支持MMX技术的CPU”+“老式操作系统”+“利用MMX指令集的程序”这3个条件,才能充分体现SIMD的优势。

SSE技术是SIMD的扩展技术,它把SIMD寄存器扩展为128位寄存器。支持SSE技术的CPU已经有了单独的SIMD通用寄存器,不再复用FPU的寄存器。

AVX是另外一种SIMD技术,它的通用寄存器都是256位寄存器。

现在转入正题,看看使用SIMD的具体应用。

SIMD技术当然可用于内存复制(memcpy)和内存比较(memcmp)等用途。

此外它还可用于DES的运算。DES(Data Encryption Standard)是分组对称密码算法。它采用了64位的分组和56位的密钥,将64位的输入经过一系列变换得到64位的输出。如果在电路的与、或、非门和导线实现DES模块,电路规模将会是非常庞大。基于并行分组密码算法的Bitslice DES[1]应运而生。这种算法可由单指令并行处理技术/SIMD实现。我们已经知道,x86平台的unsigned int型数据占用32位空间。因此,在进行unsigned int型数据的64位数据和56位密钥的演算时,可以以把64位中间运算结果和运算密钥都分为2个32位数据,再进行分步的并行处理。

在Oracle存储密码和hash的算法就是DES。为了进行有关演示,您还可以研究一下暴力破解Oracle RDBMS密码和hash的小程序。这个程序利用了bitslice DES算法,针对SSE2和AVX指令集进行了优化。对其稍加修改,则可用于128位/256位消息模块的并行加密运算。如需查看这个程序的源代码,请参阅http://conus.info/utils/ops_SIMD/

25.1 矢量化

矢量化[2]泛指将多个标量数组计算为(转换成)一个矢量数组的技术。进行这种计算时,循环体从输入数组中取值,进行某种运算后生成最终数组。这种算法只对数组中的单个元素进行一次运算。“(并行)矢量化技术”是矢量化处理的并行计算技术。

矢量化并不是最近才出现的技术。本书的作者在1998年的Cray Y-MP超级计算机系列产品中就见到过“Cray Y-MP EL”的羽量级lite版本。有兴趣的读者可以去他们在线超级计算机博物馆去看看:http://www.cray-cyber.org/

举例来说,下述循环就采用了矢量化技术:

for (i = 0; i < 1024; i++)
{
     C[i] = A[i]*B[i];
}

上述代码从数组A和数组B取值,把这两个数组元素进行乘法运算,并把求得的积存储在数组C里。

如果输入数组中的每个元素都是32位int型数据,那么就可以把A的四个元素放在128位的MMX寄存器中,再把数组B的四个元素放在另一个MMX寄存器里,然后通过PMULLD(Multiply Packed Signed Dword Integers and Store Low Result)和PMULHW(Multiply Packed Signed Integers and Store High Result)指令进行并行乘法运算。这可以一次获得4个64位的积。

这样一来,循环体的执行次数不再是1024,而是1024/4。换句话说,程序的性能将提高4倍。

某些编译器具有简单的矢量化自动优化功能。Intel C++编译器就有这样的智能优化功能,详情请参见http://www.intel.com/intelpress/sum_vmmx.htm

25.1.1 用于加法计算

本节将使用下述程序作为编译对象。

int f (int sz, int *ar1, int *ar2, int *ar3)
{
         for (int i=0; i<sz; i++)
                  ar3[i]=ar1[i]+ar2[i];
         return 0; 
};

Intel C++

我们通过下述指令用win32版本的Intel C++编译器编译上述代码:

icl intel.cpp /QaxSSE2 /Faintel.asm /Ox

利用IDA打开上面生成的可执行文件,可看到:

; int __cdecl f(int, int *, int *, int *)
                   public ?f@@YAHHPAH00@Z
?f@@YAHHPAH00@Z proc near

var_10 = dword ptr -10h
sz     = dword ptr  4
ar1    = dword ptr  8
ar2    = dword ptr  0Ch
ar3    = dword ptr  10h

       push    edi
       push    esi
       push    ebx
       push    esi
       mov     edx, [esp+10h+sz]
       test    edx, edx
       jle     loc_15B
       mov     eax, [esp+10h+ar3]
       cmp     edx, 6
       jle     loc_143
       cmp     eax, [esp+10h+ar2]
       jbe     short loc_36
       mov     esi, [esp+10h+ar2]
       sub     esi, eax
       lea     ecx, ds:0[edx*4]
       neg     esi
       cmp     ecx, esi
       jbe     short loc_55

loc_36: ; CODE XREF: f(int,int *,int *,int *)+21
       cmp     eax, [esp+10h+ar2]
       jnb     loc_143
       mov     esi, [esp+10h+ar2]
       sub     esi, eax
       lea     ecx, ds:0[edx*4]
       cmp     esi, ecx
       jb      loc_143

loc_55: ; CODE XREF: f(int,int *,int *,int *)+34
       cmp     eax, [esp+10h+ar1]
       jbe     short loc_67
       mov     esi, [esp+10h+ar1]
       sub     esi, eax
       neg     esi
       cmp     ecx, esi
       jbe     short loc_7F

loc_67: ; CODE XREF: f(int,int *,int *,int *)+59
       cmp     eax, [esp+10h+ar1]
       jnb     loc_143
       mov     esi, [esp+10h+ar1]
       sub     esi, eax
       cmp     esi, ecx
       jb      loc_143

loc_7F: ; CODE XREF: f(int,int *,int *,int *)+65
       mov     edi, eax        ;edi = ar3
       and     edi, 0Fh        ; is ar3 16-byte aligned?
       jz      short loc_9A    ; yes
       test    edi, 3
       jnz     loc_162
       neg     edi
       add     edi, 10h
       shr     edi, 2

loc_9A: ; CODE XREF: f(int,int *,int *,int *)+84
       lea     ecx, [edi+4]
       cmp     edx, ecx
       jl      loc_162
       mov     ecx, edx
       sub     ecx, edi
       and     ecx, 3
       neg     ecx
       add     ecx, edx
       test    edi, edi
       jbe     short loc_D6
       mov     ebx, [esp+10h+ar2]
       mov     [esp+10h+var_10], ecx
       mov     ecx, [esp+10h+ar1]
       xor     esi, esi

loc_C1: ; CODE XREF: f(int,int *,int *,int *)+CD
       mov     edx, [ecx+esi*4]
       add     edx, [ebx+esi*4]
       mov     [eax+esi*4], edx
       inc     esi
       cmp     esi, edi
       jb      short loc_C1
       mov     ecx, [esp+10h+var_10]
       mov     edx, [esp+10h+sz]

loc_D6: ; CODE XREF: f(int,int *,int *,int *)+B2
       mov     esi, [esp+10h+ar2]
       lea     esi, [esi+edi*4] ; is ar2+i*4 16-byte aligned?
       test    esi, 0Fh
       jz      short loc_109   ; yes!
       mov     ebx, [esp+10h+ar1]
       mov     esi, [esp+10h+ar2]

loc_ED: ; CODE XREF: f(int,int *,int *,int *)+105
       movdqu  xmm1, xmmword ptr [ebx+edi*4]
       movdqu  xmm0, xmmword ptr [esi+edi*4] ; ar2+i*4 is not 16-byte aligned, so load it to XMM0
       paddd   xmm1, xmm0
       movdqa  xmmword ptr [eax+edi*4], xmm1 add edi, 4
       cmp     edi, ecx
       jb      short loc_ED
       jmp     short loc_127

loc_109: ; CODE XREF: f(int,int *,int *,int *)+E3
       mov     ebx, [esp+10h+ar1]
       mov     esi, [esp+10h+ar2]

loc_111: ; CODE XREF: f(int,int *,int *,int *)+125
       movdqu  xmm0, xmmword ptr [ebx+edi*4]
       paddd   xmm0, xmmword ptr [esi+edi*4]
       movdqa  xmmword ptr [eax+edi*4], xmm0
       add     edi, 4
       cmp     edi, ecx
       jb      short loc_111

loc_127: ; CODE XREF: f(int,int *,int *,int *)+107
         ; f(int,int *,int *,int *)+164
       cmp     ecx, edx
       jnb     short loc_15B
       mov     esi, [esp+10h+ar1]
       mov     edi, [esp+10h+ar2]

loc_133: ; CODE XREF: f(int,int *,int *,int *)+13F
       mov     ebx, [esi+ecx*4]
       add     ebx, [edi+ecx*4]
       mov     [eax+ecx*4], ebx
       inc     ecx
       cmp     ecx, edx
       jb      short loc_133
       jmp     short loc_15B

loc_143: ; CODE XREF: f(int,int *,int *,int *)+17
         ; f(int,int *,int *,int *)+3A ...
       mov     esi, [esp+10h+ar1]
       mov     edi, [esp+10h+ar2]
       xor     ecx, ecx

loc_14D: ; CODE XREF: f(int,int *,int *,int *)+159
       mov     ebx, [esi+ecx*4]
       add     ebx, [edi+ecx*4]
       mov     [eax+ecx*4], ebx
       inc     ecx
       cmp     ecx, edx
       jb      short loc_14D

loc_15B: ; CODE XREF: f(int,int *,int *,int *)+A
         ; f(int,int *,int *,int *)+129 ...
       xor     eax, eax
       pop     ecx
       pop     ebx
       pop     esi
       pop     edi
       retn

loc_162: ; CODE XREF: f(int,int *,int *,int *)+8C
         ; f(int,int *,int *,int *)+9F
       xor     ecx, ecx
       jmp     short loc_127
?f@@YAHHPAH00@Z endp

其中,与SSE2有关的指令有:

所以上述SSE2指令的运行条件有2个:需要进行加法处理的操作数至少有4对;指针ar3的地址已经向16字节边界对齐。

而且,如果ar2的地址也向16字节边界对齐,则会执行下述指令:

movdqu  xmm0, xmmword ptr [ebx+edi*4] ; ar1+i*4
paddd   xmm0, xmmword ptr [esi+edi*4] ; ar2+i*4
movdqa  xmmword ptr [eax+edi*4], xmm0 ; ar3+i*4

另外,MOVDQU指令会把ar2的值加载到XMM0寄存器。虽然这条指令不要求指针地址对齐,但是性能略微慢些:

movdqu  xmm1, xmmword ptr [ebx+edi*4] ; ar1+i*4
movdqu  xmm0, xmmword ptr [esi+edi*4] ; ar2+i*4 is not 16-byte aligned, so load it to xmm0
paddd   xmm1, xmm0
movdqa  xmmword ptr [eax+edi*4], xmm1 ; ar3+i*4

其他代码没有涉及与SSE2有关的指令。

GCC

在指定-O3选项并且启用SSE2支持(-msse2选项)的情况下,GCC编译器也能够进行简单的矢量化智能处理。[4]

我们使用GCC 4.4.1编译上述程序,可得到:

; f(int, int *, int *, int *)
               public _Z1fiPiS_S_
_Z1fiPiS_S_ proc near

var_18      = dword ptr -18h
var_14      = dword ptr -14h
var_10      = dword ptr -10h
arg_0       = dword ptr  8
arg_4       = dword ptr  0Ch
arg_8       = dword ptr  10h
arg_C       = dword ptr  14h

            push    ebp
            mov     ebp, esp
            push    edi
            push    esi
            push    ebx
            sub     esp, 0Ch
            mov     ecx, [ebp+arg_0]
            mov     esi, [ebp+arg_4]
            mov     edi, [ebp+arg_8]
            mov     ebx, [ebp+arg_C]
            test    ecx, ecx
            jle     short loc_80484D8
            cmp     ecx, 6
            lea     eax, [ebx+10h]
            ja      short loc_80484E8

loc_80484C1: ; CODE XREF: f(int,int *,int *,int *)+4B
             ; f(int,int *,int *,int *)+61 ...
            xor     eax, eax
            nop
            lea     esi, [esi+0]

loc_80484C8: ; CODE XREF: f(int,int *,int *,int *)+36
            mov     edx, [edi+eax*4]
            add     edx, [esi+eax*4]
            mov     [ebx+eax*4], edx
            add     eax, 1
            cmp     eax, ecx
            jnz     short loc_80484C8

loc_80484D8: ; CODE XREF: f(int,int *,int *,int *)+17
             ; f(int,int *,int *,int *)+A5
            add     esp, 0Ch
            xor     eax, eax
            pop     ebx
            pop     esi
            pop     edi
            pop     ebp
            retn

            align8

loc_80484E8: ; CODE XREF: f(int,int *,int *,int *)+1F
            test    bl, 0Fh
            jnz     short loc_80484C1
            lea     edx, [esi+10h]
            cmp     ebx, edx
            jbe     loc_8048578

loc_80484F8: ; CODE XREF: f(int,int *,int *,int *)+E0
            lea     edx, [edi+10h]
            cmp     ebx, edx
            ja      short loc_8048503
            cmp     edi, eax
            jbe     short loc_80484C1

loc_8048503: ; CODE XREF: f(int,int *,int *,int *)+5D
            mov     eax, ecx
            shr     eax, 2
            mov     [ebp+var_14], eax
            shl     eax, 2
            test    eax, eax
            mov     [ebp+var_10], eax
            jz      short loc_8048547
            mov     [ebp+var_18], ecx
            mov     ecx, [ebp+var_14]
            xor     eax, eax
            xor     edx, edx
            nop

loc_8048520: ; CODE XREF: f(int,int *,int *,int *)+9B
            movdqu  xmm1, xmmword ptr [edi+eax]
            movdqu  xmm0, xmmword ptr [esi+eax]
            add     edx, 1
            paddd   xmm0, xmm1
            movdqa  xmmword ptr [ebx+eax], xmm0
            add     eax, 10h
            cmp     edx, ecx
            jb      short loc_8048520
            mov     ecx, [ebp+var_18]
            mov     eax, [ebp+var_10]
            cmp     ecx, eax
            jz      short loc_80484D8

loc_8048547: ; CODE XREF: f(int,int *,int *,int *)+73
            lea     edx, ds:0[eax*4]
            add     esi, edx
            add     edi, edx
            add     ebx, edx
            lea     esi, [esi+0]

loc_8048558: ; CODE XREF: f(int,int *,int *,int *)+CC
            mov     edx, [edi]
            add     eax, 1
            add     edi, 4
            add     edx, [esi]
            add     esi, 4
            mov     [ebx], edx
            add     ebx, 4
            cmp     ecx, eax
            jg      short loc_8048558
            add     esp, 0Ch
            xor     eax, eax
            pop     ebx
            pop     esi
            pop     edi
            pop     ebp
            retn

loc_8048578: ; CODE XREF: f(int,int *,int *,int *)+52
            cmp     eax, esi
            jnb     loc_80484C1
            jmp     loc_80484F8
_Z1fiPiS_S_ endp

GCC生成的代码和Intel C++十分相似,只是严谨性略差。

25.1.2 用于内存复制

我们回顾一下本书14.2节的memcpy()函数:

#include <stdio.h>

void my_memcpy (unsigned char* dst, unsigned char* src, size_t cnt)
{
     size_t i;
     for (i=0; i<cnt; i++)
          dst[i]=src[i];
};

用GCC 4.9.1进行优化编译,可得到如下所示的代码。

指令清单25.1 Optimizing GCC 4.9.1 x64

my_memcpy:
; RDI = destination address
; RSI = source address
; RDX = size of block
        test    rdx, rdx
        je      .L41
        lea     rax, [rdi+16]
        cmp     rsi, rax
        lea     rax, [rsi+16]
        setae   cl
        cmp     rdi, rax
        setae   al
        or      cl, al
        je      .L13
        cmp     rdx, 22
        jbe     .L13
        mov     rcx, rsi
        push    rbp
        push    rbx
        neg     rcx
        and     ecx, 15
        cmp     rcx, rdx
        cmova   rcx, rdx
        xor     eax, eax
        test    rcx, rcx
        je      .L4
        movzx   eax, BYTE PTR [rsi]
        cmp     rcx, 1
        mov     BYTE PTR [rdi], al
        je      .L15
        movzx   eax, BYTE PTR [rsi+1]
        cmp     rcx, 2
        mov     BYTE PTR [rdi+1], al
        je      .L16
        movzx   eax, BYTE PTR [rsi+2]
        cmp     rcx, 3
        mov     BYTE PTR [rdi+2], al
        je      .L17
        movzx   eax, BYTE PTR [rsi+3]
        cmp     rcx, 4
        mov     BYTE PTR [rdi+3], al
        je      .L18
        movzx   eax, BYTE PTR [rsi+4]
        cmp     rcx, 5
        mov     BYTE PTR [rdi+4], al
        je      .L19
        movzx   eax, BYTE PTR [rsi+5]
        cmp     rcx, 6
        mov     BYTE PTR [rdi+5], al
        je      .L20
        movzx   eax, BYTE PTR [rsi+6]
        cmp     rcx, 7
        mov     BYTE PTR [rdi+6], al 
        je      .L21
        movzx   eax, BYTE PTR [rsi7+]
        cmp     rcx, 8
        mov     BYTE PTR [rdi+7], al 
        je      .L22
        movzx   eax, BYTE PTR [rsi+8]
        cmp     rcx, 9
        mov     BYTE PTR [rdi+8], al 
        je      .L23
        movzx   eax, BYTE PTR [rsi+9]
        cmp     rcx, 10
        mov     BYTE PTR [rdi+9], al 
        je      .L24
        movzx   eax, BYTE PTR [rsi+10]
        cmp     rcx, 11
        mov     BYTE PTR [rdi+10], al 
        je      .L25
        movzx   eax, BYTE PTR [rsi+11]
        cmp     rcx, 12
        mov     BYTE PTR [rdi+11], al 
        je      .L26
        movzx   eax, BYTE PTR [rsi+12]
        cmp     rcx, 13
        mov     BYTE PTR [rdi+12], al 
        je      .L27
        movzx   eax, BYTE PTR [rsi+13]
        cmp     rcx, 15
        mov     BYTE PTR [rdi+13], al 
        jne     .L28
        movzx   eax, BYTE PTR [rsi+14]
        mov     BYTE PTR [rdi+14], al
        mov     eax, 15
.L4:
        mov     r10, rdx
        lea     r9, [rdx-1]
        sub     r10, rcx
        lea     r8, [r10-16]
        sub     r9, rcx
        shr     r8, 4
        add     r8, 1
        mov     r11, r8
        sal     r11, 4
        cmp     r9,14
        jbe     .L6
        lea     rbp, [rsi+rcx]
        xor     r9d, r9d
        add     rcx, rdi
        xor     ebx, ebx
.L7:
        movdqa  xmm0, XMMWORD PTR [rbp+0+r9]
        add     rbx, 1
        movups  XMMWORD PTR [rcx+r9], xmm0
        add     r9, 16
        cmp     rbx, r8
        jb      .L7
        add     rax, r11
        cmp     r10, r11
        je      .L1
.L6:
        movzx   ecx, BYTE PTR [rsi+rax]
        mov     BYTE PTR [rdi+rax], cl
        lea     rcx, [rax+1]
        cmp     rdx, rcx
        jbe     .L1
        movzx   ecx, BYTE PTR [rsi+1+rax]
        mov     BYTE PTR [rdi+1+rax], cl
        lea     rcx, [rax+2]
        cmp     rdx, rcx
        jbe     .L1
        movzx   ecx, BYTE PTR [rsi+2+rax]
        mov     BYTE PTR [rdi+2+rax], cl
        lea     rcx, [rax+3]
        cmp     rdx, rcx
        jbe     .L1
        movzx   ecx, BYTE PTR [rsi+3+rax]
        mov     BYTE PTR [rdi+3+rax], cl
        lea     rcx, [rax+4]
        cmp     rdx, rcx
        jbe     .L1
        movzx   ecx, BYTE PTR [rsi+4+rax]
        mov     BYTE PTR [rdi+4+rax], cl
        lea     rcx, [rax+5]
        cmp     rdx, rcx
        jbe     .L1
        movzx   ecx, BYTE PTR [rsi+5+rax]
        mov     BYTE PTR [rdi+5+rax], cl
        lea     rcx, [rax+6]
        cmp     rdx, rcx
        jbe     .L1
        movzx   ecx, BYTE PTR [rsi+6+rax]
        mov     BYTE PTR [rdi+6+rax], cl
        lea     rcx, [rax+7]
        cmp     rdx, rcx
        jbe     .L1
        movzx   ecx, BYTE PTR [rsi+7+rax]
        mov     BYTE PTR [rdi+7+rax], cl
        lea     rcx, [rax+8]
        cmp     rdx, rcx
        jbe     .L1
        movzx   ecx, BYTE PTR [rsi+8+rax]
        mov     BYTE PTR [rdi+8+rax], cl
        lea     rcx, [rax+9]
        cmp     rdx, rcx
        jbe     .L1
        movzx   ecx, BYTE PTR [rsi+9+rax]
        mov     BYTE PTR [rdi+9+rax], cl
        lea     rcx, [rax+10]
        cmp     rdx, rcx
        jbe     .L1
        movzx   ecx, BYTE PTR [rsi+10+rax]
        mov     BYTE PTR [rdi+10+rax], cl
        lea     rcx, [rax+11]
        cmp     rdx, rcx
        jbe     .L1
        movzx   ecx, BYTE PTR [rsi+11+rax]
        mov     BYTE PTR [rdi+11+rax], cl
        lea     rcx, [rax+12]
        cmp     rdx, rcx
        jbe     .L1
        movzx   ecx, BYTE PTR [rsi+12+rax]
        mov     BYTE PTR [rdi+12+rax], cl
        lea     rcx, [rax+13]
        cmp     rdx, rcx
        jbe     .L1
        movzx   ecx, BYTE PTR [rsi+13+rax]
        mov     BYTE PTR [rdi+13+rax], cl
        lea     rcx, [rax+14]
        cmp     rdx, rcx
        jbe     .L1
        movzx   edx, BYTE PTR [rsi+14+rax]
        mov     BYTE PTR [rdi+14+rax], dl
.L1:
        pop     rbx
        pop     rbp
.L41:
        rep ret
.L13:
        xor     eax, eax
.L3:
        movzx   ecx, BYTE PTR [rsi+rax]
        mov     BYTE PTR [rdi+rax], cl
        add     rax, 1
        cmp     rax, rdx
        jne     .L3
        rep ret
.L28:
        mov     eax, 14
        jmp     .L4
.L15:
        mov     eax, 1
        jmp     .L4
.L16:
        mov     eax, 2
        jmp     .L4
.L17:
        mov     eax, 3
        jmp     .L4
.L18:
        mov     eax, 4
        jmp     .L4
.L19:
        mov     eax, 5
        jmp     .L4
.L20:
        mov     eax, 6
        jmp     .L4
.L21:
        mov     eax, 7
        jmp     .L4
.L22:
        mov     eax, 8
        jmp     .L4
.L23:
        mov     eax, 9
        jmp     .L4
.L24:
        mov     eax, 10
        jmp     .L4
.L25:
        mov     eax, 11
        jmp     .L4
.L26:
        mov     eax, 12
        jmp     .L4
.L27:
        mov     eax, 13
        jmp     .L4

25.2 SIMD实现strlen()

根据msdn(http://msdn.microsoft.com/en-us/library/y0dh78ez(VS.80).aspx

),在C/C++源程序中插入特定的宏即可调用SIMD指令。在MSVC编译器的库文件之中,intrin.h文件就使用了这种宏。

我们可以把字符串进行分组处理,使用SIMD指令实现strlen()函数。这种算法比常规实现算法快2~2.5倍。它每次把16个字符加载到1个XMM寄存器里,然后检测字符串的结束标志—数值为零的结束符。[5]

size_t strlen_sse2(const char *str)
{
    register size_t len = 0;
    const char *s=str;
    bool str_is_aligned=(((unsigned int)str)&0xFFFFFFF0) == (unsigned int)str;

    if (str_is_aligned==false)
         return strlen (str);

    __m128i xmm0 = _mm_setzero_si128();
    __m128i xmm1;
    int mask = 0;

    for (;;) 
    {
         xmm1 = _mm_load_si128((__m128i *)s);
         xmm1 = _mm_cmpeq_epi8(xmm1, xmm0);
         if ((mask = _mm_movemask_epi8(xmm1)) != 0)
         {
             unsigned long pos;
             _BitScanForward(&pos, mask);
             len += (size_t)pos;

             break; 
         }
         s += sizeof(__m128i);
         len += sizeof(__m128i);
    };

    return len; 
}

使用MSVC 2010(启用/Ox选项)编译上述程序,可得到如下所示的代码。

指令清单25.2 Optimizing MSVC 2010

_pos$75552 = -4             ; size = 4
_str$ = 8                   ; size = 4
?strlen_sse2@@YAIPBD@Z PROC ; strlen_sse2

     pushebp
     mov      ebp, esp
     and      esp, -16      ; fffffff0H
     mov      eax, DWORD PTR _str$[ebp]
     sub      esp, 12       ; 0000000cH
     push     esi
     mov      esi, eax
     and      esi, -16      ; fffffff0H
     xor      edx, edx
     mov      ecx, eax
     cmp      esi, eax
     je       SHORT $LN4@strlen_sse
     lea      edx, DWORD PTR [eax+1]
     npad     3 ; align next label
$LL11@strlen_sse:
     mov      cl, BYTE PTR [eax]
     inc      eax
     test     cl, cl
     jne      SHORT $LL11@strlen_sse
     sub      eax, edx
     pop      esi
     mov      esp, ebp
     pop      ebp
     ret      0
$LN4@strlen_sse:
     movdqa   xmm1, XMMWORD PTR [eax]
     pxor     xmm0, xmm0
     pcmpeqb  xmm1, xmm0
     pmovmskb eax, xmm1
     test     eax, eax
     jne      SHORT $LN9@strlen_sse
$LL3@strlen_sse:
     movdqa   xmm1, XMMWORD PTR [ecx+16]
     add      ecx, 16                ; 00000010H
     pcmpeqb  xmm1, xmm0
     add      edx, 16                ; 00000010H
     pmovmskb eax, xmm1
     test     eax, eax
     je       SHORT $LL3@strlen_sse
$LN9@strlen_sse:
     bsf      eax, eax
     mov      ecx, eax
     mov      DWORD PTR _pos$75552[esp+16], eax
     lea      eax, DWORD PTR [ecx+edx]
     pop      esi
     mov      esp, ebp
     pop      ebp
     ret      0
?strlen_sse2@@YAIPBD@Z ENDP          ; strlen_sse2

程序首先检查str指针是否已经向16字节边界对齐。如果这个地址没有对齐,就调用常规的strlen()函数计算字符串长度。

然后通过MOVDQA指令把16字节数据加载到XMM1寄存器。

部分读者可能会问:为什么不直接使用不要求指针对齐的MOVDQU指令?

这种说法似乎有道理:如果指针已经向16字节边界对齐了,我们可以使用MOVDQA指令;否则,就使用速度慢些的MOVDQU指令。

但是本文意在强调:

在Windows NT系列操作系统,以及其他一些操作系统里,内存的最小分页单位是4KB(4096字节)。虽然每个win32程序表面上都有(虚拟)独立的4GB内存可以使用,但是实际上部分地址空间受到物理内存和分页方法的限制。如果程序试图访问不存在的内存块,将会引发错误。这是虚拟内存的工作方式决定的,详情请参见http://en.wikipedia.org/wiki/Page_(computer_memory)

所以,如果函数一次加载16字节数据,可能会跨越虚拟内存的内存块边界。例如说,操作系统可能在地址为0x00c0000的物理内存处给这个程序分配了8192(0x2000)字节的内存块。那么这个程序的内存就在0x008c0000与0x008c1fff之间。

在这个内存块之后,即位于0x008c2000地址以后的物理内存都不可被程序访问。换而言之,操作系统没有给这个程序分配那块内存。如果程序访问这块不可用的内存地址就将引发错误。

我们再来研究一下这个程序。程序可能在内存块的末端存储5个字符。这种分配方法是合情合理的。

0x008c1ff8
0x008c1ff9
0x008c1ffa
0x008c1ffb
0x008c1ffc
0x008c1ffd
0x008c1ffe
0x008c1fff
’h’
’e’
’l’
’l’
’o’
’\x00’
random noise
random noise

当采用常规算法的strlen()函数在处理这个字符串时,首先把指针指向“hello”的第一个字节,即地址0x008c1ff8。strlen()函数会逐字节地读取字符串。在遇到0x008c1ffd的0字节的时候,函数就结束了。

假如我们编写的strlen()函数不考虑字符串的地址对齐问题,无论字符串指针地址是否已经向16字节边界对齐,每次都读取16字节。那么它就可能要从0x008c1ff8一直读取到0x008c2008,这将引发错误。当然,我们应当尽力避免这种错误。

所以仅在字符串地址向16字节边界对齐的情况下,才可以使用基于SIMD技术的算法。另外,操作系统在进行内存分页时也是向16字节边界对齐的,所以在字符串首地址同样向16字节边界对齐时,我们的函数不会访问未分配的内存空间。

下面继续介绍我们的函数。

_mm_setzero_si128():清空XMM0寄存器的指令,相当于pxor xmm0, xmm0。

_mm_load_si128():MOVDQA的宏。它从指定地址里加载16字节数据到XMM1寄存器。

_mm_cmpeq_epi8():PCMPEQB的宏。它以字节为单位,比较两个XMM专用寄存器的值。如果某个字节相等,该字节返回值为0xff,否则返回值为0。举例来说,运行pcmpeqdb xmm1, xmm0指令的情况如下表所示。

运行前 XMM1 11223344556677880000000000000000
运行前 XMM0 11ab3444007877881111111111111111
运行后 XMM1 ff0000ff0000ffff0000000000000000

本例先用“pxor xmm0, xmm0”将xmm0寄存器清零,然后把从字符串摘取出的16字节与XMM0寄存器的16个0相比较。

下一条指令是_mm_movemask_epi8(),即PMOVVMSKB指令。这个指令通常与PCMPEQB配合使用:

pmovmskb eax, xmm1

在XMM1的第一个字节的最高位是1的情况下,设置EAX寄存器的最高位为1。也就是说,如果XMM1寄存器的值是0xff,那么EAX寄存器的最高位还会被赋值为1。

它还会在XMM1的第二个字节是0xff的情况下,设置EAX寄存器的次高位为1。

总而言之,这个指令此处的作用是“判断XMM1的每个字节是否是0xFF”。如果XMM1的第n个字节是0xFF,则设置EAX的第n位为1,否则EAX的第n位为0。

另外,不要忘记我们的算法可能会遇到的实际问题:指针指向的字符串可能包含有效值、结束符和脏数据。例如:

15
14
13
12
11
10
9~3
2
1~0
“h”
“e”
“l”
“l”
“o”
0
脏数据
0
脏数据

一次读取16字节到XMM寄存器,并将它与XMM0里的零进行比较,可能会得到这样的计算结果(从左到右是LSB到MSB的顺序):

XMM1: 0000ff00000000000000ff0000000000

也就是说它找到了2个以上的零。这是情理之中的结果。

在这种情况下PMOVMSKB指令将设置EAX为(二进制):0010000000100000b。

很明显,我们的函数必须以第一个结束符为准,忽略掉脏数据里的结束符。

下一条BSF(Bit Scan Forward)指令将保留操作符二进制值里的第一个1,将其余位清零,并把结果存储到第一个操作符里。

在EAX=0010000000100000b的情况下,执行“bsf eax, eax”,那么EAX将将保留第5位的1(从0开始数)。

MSVC里这条指令对应的宏是_BitScanForward。

代码的其他部分就容易理解了。如果找到了0字节,就把它的位置编号累加到计数器里,从而获取字符串长度。

应当注意到MSVC编译器优化掉了循环体的部分结构。

Intel Core I7处理器推出了SSE 4.2。它提供的新的指令大大简化了字符串的处理。有兴趣的读者可参见:http://www.strchr.com/strcmp_and_strlen_using_sse_4.2


[1] 请参见http://www.darkside.com.au/bitslice/。

[2] Vectorization,又称“向量化”。请参见http://en.wikipedia.org/wiki/Vectorization_ (computer_science)。

[3] 有关地址对齐的详细信息,请参见:http://en.wikipedia.org/wiki/Data_structure_ alignment。

[4] https://gcc.gnu.org/projects/tree-ssa/vectorization.html。

[5] 这种算法借鉴了网上的资料http://www.strchr.com/sse2_optimised_strlen。