《高级语言内的单指令多数据流计算(SIMD)》整理

侯思松整理出了一个非常好的文章——《高级语言内的单指令多数据流计算(SIMD)》。我这里对其内容做进一步整理，使得其能够被更容易地接受。

我们知道传统计算机的指令集架构主要完成的是基本的算术逻辑计算、条件分支、I/O访问等基础功能。但是目前，在很多领域光靠这些基础指令还不足以完成高密集计算，因此DSP、流处理器等各种适合高性能计算的处理器诞生。而在90年代，Intel出了MMX（多媒体指令集扩展）用于加速应用多媒体领域的高密集计算。随后，x86处理器架构中又诞生了SSE、AMD自带的3D NOW!（+），SSE2、SSE3、SSSE3（仅Intel支持）、SSE4.1、SSE4.2、SSE4A（仅AMD支持）、AVX等更高级的SIMD指令集扩展。而SIMD的意思是“单指令多数据”，也就是说一条指令可以让多条数据独立而并行地获得处理。比如我们可以同时让4对整型数据同时做加法计算，而只需要1个指令周期甚至更少。

除了x86体系之外，其它传统处理器架构也引入了SIMD指令集扩展，比如Power架构的AltiVec、ARMv7架构中所引入的NEON技术、MIPS架构所引入的MIPS-3D等等。

然而还有很多中低端的处理器并没有专门的SIMD指令集，这个时候我们可以采取某些特定的算法，以比较小的代价来模拟SIMD所能达到的功效。下面，我们将以32位处理器为例，介绍一些算法。

情景一：将一组8位整型数据做逻辑右移一位。一般做法是：

extern void naive1(void);
void naive1(void)
{
    // Initialize
    __attribute__((aligned(4))) unsigned char buffer[4] = { 255, 128, 11, 6 };
    
    // Calculate
    for(int i = 0; i < 4; i++)
        buffer[i] >>= 1;
    
    // Output results
    printf("Results: ");
    for(int i = 0; i < 4; i++)
        printf("0x%.2X  ", buffer[i]);
    puts("");
}

下面提供上面算法的仿SIMD算法：

extern void simd_test1(void);
void simd_test1(void)
{
    // Initialize
    unsigned buffer = 0x060b80ff;
    
    // Calculate
    buffer &= 0xfefefefe;
    buffer >>= 1;
    
    // Output results
    printf("The result is: 0x%.8X\n", buffer);
}

我们可以看到，原始版本需要以单字节访存4次，然后做移位操作；而仿SIMD的优化版本仅仅做一次4字节的访存操作，然后做两步算术逻辑计算就完事了，非常精炼。

当然，从这个例子我们还能再扩展到左移的情况。下面呢，我们将扩展算术右移一位的情况，这个是原文中没有的，先看原始算法：

extern void naive1(void);
void naive1(void)
{
    // Initialize
    __attribute__((aligned(4))) char buffer[4] = { -1, -128, 11, 6 };
    
    // Calculate
    for(int i = 0; i < 4; i++)
        buffer[i] >>= 1;
    
    // Output results
    printf("Results: ");
    for(int i = 0; i < 4; i++)
        printf("0x%.2X  ", buffer[i]);
    puts("");
}

再看看仿SIMD的优化版：

extern void simd_test1(void);
void simd_test1(void)
{
    // Initialize
    unsigned buffer = 0x060b80ff;
    
    // Calculate
    unsigned mask = buffer & 0x80808080;
    buffer &= 0xfefefefe;
    buffer >>= 1;
    buffer |= mask;
    
    // Output results
    printf("The result is: 0x%.8X\n", buffer);
}

由于算术右移需要考虑符号位。因此，我们需要用一个mask变量来存放原始数据的符号值，然后将最终移位结果“或”上 符号位。这样就比逻辑右移要多2步操作，不过即便如此，其开销比访存4次仍然小不少。

情景二：

对一个无符号8位整型做一次补操作。即，255 - x（x为一个无符号8位整型数据）。原始算法为：

extern void naive2(void);
void naive2(void)
{
    // Initialize
    __attribute__((aligned(4))) unsigned char buffer[4] = { 255, 128, 11, 6 };
    
    // Calculate
    for(int i = 0; i < 4; i++)
        buffer[i] = 255 - buffer[i];
    
    // Output results
    printf("Results: ");
    for(int i = 0; i < 4; i++)
        printf("0x%.2X  ", buffer[i]);
    puts("");
}

优化版本为：

extern void simd_test2(void);
void simd_test2(void)
{
    // Initialize
    unsigned buffer = 0x060b80ff;
    
    // Calculate
    buffer = ~buffer;
    
    // Output results
    printf("The result is: 0x%.8X\n", buffer);
}

这个算法相对来说还是比较容易理解的。因为这操作本身就是取0～255整数范围内的补。 

情景三：

计算两个无符号整型的算术平均数。这里有两种情况，一种是当碰到两数的和是奇数的话，向下取整；另一种是向上取整。

我们先看向下取整的情况，即(x + y) / 2：

原始算法为：

extern void naive3(void);
void naive3(void)
{
    // Initialize
    __attribute__((aligned(4))) unsigned char buffer1[4] = { 255, 128, 11, 33 };
    __attribute__((aligned(4))) unsigned char buffer2[4] = { 100, 129, 19, 55 };
    __attribute__((aligned(4))) unsigned char dstBuffer[4];
    
    // Calculate
    for(int i = 0; i < 4; i++)
        dstBuffer[i] = ((unsigned)buffer1[i] + (unsigned)buffer2[i]) / 2;
    
    // Output results
    printf("Results: ");
    for(int i = 0; i < 4; i++)
        printf("0x%.2X  ", dstBuffer[i]);
    puts("");
}

仿SIMD的优化算法为：

extern void simd_test3(void);
void simd_test3(void)
{
    // Initialize
    unsigned buffer1 = 0x210b80ff;
    unsigned buffer2 = 0x37138164;
    unsigned dstBuffer;
    
    // Calculate
    dstBuffer = (buffer1 & buffer2) + (((buffer1 ^ buffer2) & 0xfefefefe) >> 1);
    
    // Output results
    printf("The result is: 0x%.8X\n", dstBuffer);
}

这个算法源自FFmpeg，算法本身非常巧妙。

然后，我们看一下向上取整的算法，即：sum = (x + y); result = sum / 2;if(sum % 1 != 0)result += 1;

原始算法为：

extern void naive3(void);
void naive3(void)
{
    // Initialize
    __attribute__((aligned(4))) unsigned char buffer1[4] = { 255, 128, 11, 33 };
    __attribute__((aligned(4))) unsigned char buffer2[4] = { 100, 129, 19, 55 };
    __attribute__((aligned(4))) unsigned char dstBuffer[4];
    
    // Calculate
    for(int i = 0; i < 4; i++)
    {
        unsigned sum = (unsigned)buffer1[i] + (unsigned)buffer2[i];
        unsigned char result = sum / 2;
        dstBuffer[i] = (sum & 1) == 0? result : result + 1;
    }
    
    // Output results
    printf("Results: ");
    for(int i = 0; i < 4; i++)
        printf("0x%.2X  ", dstBuffer[i]);
    puts("");
}

仿SIMD的优化算法为：

extern void simd_test3(void);
void simd_test3(void)
{
    // Initialize
    unsigned buffer1 = 0x210b80ff;
    unsigned buffer2 = 0x37138164;
    unsigned dstBuffer;
    
    // Calculate
    dstBuffer = (buffer1 | buffer2) - (((buffer1 ^ buffer2) & 0xfefefefe) >> 1);
    
    // Output results
    printf("The result is: 0x%.8X\n", dstBuffer);
}

情景四：

将两个字节按比例做混合，即，dst=(a * (255 - s) + b * s) / 255：