SIMD 算法设计模式

SIMD 编程的核心不是”如何用 intrinsic”，而是”如何把问题变换成 SIMD 能高效处理的形状”。这篇文章总结了我在实践中反复遇到的几类设计模式，从数据布局到查表，从无分支选择到流压缩，每一个模式都对应着一类真实的性能优化场景。

一、SoA vs AoS——数据布局决定一切

在谈任何 SIMD 指令之前，必须先谈数据布局。这是我在 SIMD 优化中最深刻的一条教训：指令选择是战术问题，数据布局才是战略问题。

1.1 AoS 布局的困境

Array of Structures（AoS）是最自然的面向对象布局方式：

// AoS: Array of Structures
typedef struct {
    float x, y, z;
    float vx, vy, vz;
    float mass;
    int   id;
} Particle;

Particle particles[1024];

在内存中的排列：

[x0 y0 z0 vx0 vy0 vz0 mass0 id0] [x1 y1 z1 vx1 vy1 vz1 mass1 id1] ...

如果我们只需要更新所有粒子的 x 坐标，SIMD 需要从每个结构体中”挑出” x 字段。这意味着：

加载一个缓存行（64 字节），但只用其中 4 字节（x 字段）。
要么用 gather 指令按索引加载（慢），要么手动 load + shuffle 重组（繁琐）。
缓存利用率极低——带宽被大量无用字段占据。

1.2 SoA 布局的优势

Structure of Arrays（SoA）将每个字段拆成独立数组：

// SoA: Structure of Arrays
typedef struct {
    float x[1024];
    float y[1024];
    float z[1024];
    float vx[1024];
    float vy[1024];
    float vz[1024];
    float mass[1024];
    int   id[1024];
} ParticleSystem;

ParticleSystem ps;

内存中：

[x0 x1 x2 x3 x4 x5 x6 x7 ...] [y0 y1 y2 y3 ...]  [z0 z1 z2 z3 ...]

更新所有 x 坐标？直接对 ps.x 数组做连续的 SIMD 加载和运算：

// SoA 布局下的向量化更新——连续内存，完美 SIMD
void update_positions_soa(ParticleSystem *ps, float dt, int n) {
    for (int i = 0; i < n; i += 8) {
        __m256 x  = _mm256_load_ps(&ps->x[i]);
        __m256 vx = _mm256_load_ps(&ps->vx[i]);
        __m256 dt_vec = _mm256_set1_ps(dt);
        x = _mm256_fmadd_ps(vx, dt_vec, x);   // x += vx * dt
        _mm256_store_ps(&ps->x[i], x);
    }
}

对比 AoS 布局下的等价操作：

// AoS 布局下——散乱的内存访问，难以向量化
void update_positions_aos(Particle *particles, float dt, int n) {
    for (int i = 0; i < n; i++) {
        particles[i].x += particles[i].vx * dt;
    }
}

编译器也许能自动向量化 AoS 版本，但生成的代码通常需要额外的 gather/scatter 或交叉存取指令，效率远不如 SoA 版本。

1.3 AoSoA——折中方案

纯 SoA 的缺点是失去了空间局部性——如果算法需要同时访问同一粒子的多个字段，SoA 会在多个数组之间跳跃。AoSoA（Array of Structure of Arrays）是折中方案：

// AoSoA: 每个块内是 SoA，块的大小等于 SIMD 宽度
#define BLOCK 8  // AVX2: 256-bit / 32-bit = 8 floats
typedef struct {
    float x[BLOCK];
    float y[BLOCK];
    float z[BLOCK];
    float vx[BLOCK];
    float vy[BLOCK];
    float vz[BLOCK];
} ParticleBlock;

ParticleBlock blocks[1024 / BLOCK];

这样，一个块内所有粒子的 x 是连续的（SIMD 友好），同时同一粒子的 x 和 y 在相邻内存中（空间局部性好）。

1.4 Data-Oriented Design 的核心思想

SoA 布局是 Data-Oriented Design（DOD）的核心实践之一。DOD 的哲学很简单：

按访问模式组织数据，而不是按”对象”的逻辑关系。
热/冷分离：频繁访问的字段放一起，偶尔访问的分开。
为硬件而设计：缓存行是 64 字节，SIMD 寄存器是 256/512 位——数据布局应该匹配这些物理约束。

这不是一种”优化技巧”，而是一种设计哲学的转变。游戏引擎（Unity DOTS / Unreal Mass Framework）和高性能计算领域已经广泛采用这种思路。

个人观点：我见过太多团队在 AoS 布局上拼命调指令、换算法，收效甚微。如果你的数据布局对 SIMD 不友好，那么在指令层面做的任何努力都是在弥补结构性的缺陷。先改布局，再谈指令。

二、Mask + Blend 替代 if-else

分支是 SIMD 的天敌。向量寄存器里的 8 个元素可能各自需要走不同的分支——这在标量代码中用 if-else 很自然，但在 SIMD 中必须换一种思路。

2.1 基本模式：compare → mask → blend

核心思想是同时计算两个分支的结果，然后用 mask 选择：

// 标量版本：条件赋值
for (int i = 0; i < n; i++) {
    if (a[i] > threshold) {
        result[i] = a[i] * 2;     // true 分支
    } else {
        result[i] = a[i] + 100;   // false 分支
    }
}

SIMD 版本：

// SIMD 版本：无分支条件赋值
void conditional_assign_avx2(const int32_t *a, int32_t *result,
                             int32_t threshold, int n) {
    __m256i thresh_vec = _mm256_set1_epi32(threshold);
    __m256i two_vec    = _mm256_set1_epi32(2);
    __m256i hundred    = _mm256_set1_epi32(100);

    for (int i = 0; i < n; i += 8) {
        __m256i va = _mm256_load_si256((__m256i*)&a[i]);

        // 1. Compare: 生成 mask（全1 或 全0）
        __m256i mask = _mm256_cmpgt_epi32(va, thresh_vec);

        // 2. 计算两个分支的结果
        __m256i true_val  = _mm256_mullo_epi32(va, two_vec);   // a * 2
        __m256i false_val = _mm256_add_epi32(va, hundred);     // a + 100

        // 3. Blend: 按 mask 选择
        // mask 为全1的位置取 true_val，全0 的位置取 false_val
        __m256i res = _mm256_blendv_epi8(false_val, true_val, mask);

        _mm256_store_si256((__m256i*)&result[i], res);
    }
}

2.2 _mm256_blendv_epi8 的工作原理

_mm256_blendv_epi8(a, b, mask) 的语义：

对于每个字节 i:
    if (mask[i] 的最高位 == 1)
        result[i] = b[i]
    else
        result[i] = a[i]

注意它是按字节粒度的，且只看最高位。_mm256_cmpgt_epi32 生成的 mask 每个 int32 元素要么是 0xFFFFFFFF（全1），要么是 0x00000000（全0），所以各字节的最高位一致，可以安全用于 blendv。

2.3 浮点数的条件选择

对浮点数，可以直接用 _mm256_blendv_ps：

// 浮点 clamp: result = max(min_val, min(x, max_val))
// 等价于: if (x < min_val) result = min_val;
//         else if (x > max_val) result = max_val;
//         else result = x;
void clamp_avx2(const float *x, float *result,
                float min_val, float max_val, int n) {
    __m256 vmin = _mm256_set1_ps(min_val);
    __m256 vmax = _mm256_set1_ps(max_val);

    for (int i = 0; i < n; i += 8) {
        __m256 vx = _mm256_load_ps(&x[i]);

        // 直接用 min/max 指令——这本身就是无分支的
        __m256 res = _mm256_max_ps(vmin, _mm256_min_ps(vx, vmax));

        _mm256_store_ps(&result[i], res);
    }
}

实际上，_mm256_min_ps 和 _mm256_max_ps 内部就是 compare + blend 的硬件实现。很多看起来需要条件分支的操作，仔细分析后都可以用 min/max/abs 等指令直接替代。

2.4 多分支的处理

当条件不止两个时，逐层叠加 mask：

// 三段分类: x < 0 → -1, x == 0 → 0, x > 0 → 1 (signum)
__m256i signum_avx2(__m256i x) {
    __m256i zero = _mm256_setzero_si256();
    __m256i pos_mask = _mm256_cmpgt_epi32(x, zero);    // x > 0
    __m256i neg_mask = _mm256_cmpgt_epi32(zero, x);    // x < 0

    __m256i ones     = _mm256_set1_epi32(1);
    __m256i neg_ones = _mm256_set1_epi32(-1);

    // pos_mask 处取 1，否则取 0
    __m256i result = _mm256_and_si256(pos_mask, ones);
    // neg_mask 处取 -1，OR 上去
    result = _mm256_or_si256(result, _mm256_and_si256(neg_mask, neg_ones));

    return result;
}

这里用 AND + OR 代替了 blend——当选择值是常量时，AND mask 比 blendv 更高效（1 cycle vs 2 cycles on many microarchitectures）。

个人观点：mask + blend 模式最大的思维转变在于——你必须接受”浪费性地计算两个分支”。在标量思维中这很反直觉，但在 SIMD 世界里，“一条指令算 8 个没用的结果”远比”8 次分支预测”要快。

三、pshufb——16 路并行 4-bit 查表

pshufb（Packed Shuffle Bytes）可能是 SSSE3 以来最强大的单条指令。它的本质是一个 16 路并行的 4-bit 索引查表操作。

3.1 _mm_shuffle_epi8 的工作原理

// __m128i _mm_shuffle_epi8(__m128i table, __m128i index)
//
// 对于每个字节 i (0..15):
//   if (index[i] 的最高位 == 1)
//       result[i] = 0            // 最高位为 1 → 输出 0
//   else
//       result[i] = table[index[i] & 0x0F]  // 低 4 位作为索引

关键点： - table 是一个 16 字节的查找表。 - index 的每个字节的低 4 位是查找索引（0-15），最高位是”置零”标志。 - 16 路查找同时完成，延迟仅 1 个时钟周期。

3.2 应用一：SIMD popcount

popcount（计算一个整数中 1 的个数）是 pshufb 最经典的应用。思路：将每个字节拆成高 4 位和低 4 位，分别查表得到 1 的个数，再相加。

#include <immintrin.h>

// 标量 popcount（参考实现）
int popcount_scalar(const uint8_t *data, int n) {
    int count = 0;
    for (int i = 0; i < n; i++) {
        count += __builtin_popcount(data[i]);
    }
    return count;
}

// SSE3 pshufb popcount
// 每次处理 16 字节，每个字节查表得到其 popcount
__m128i popcount_16bytes(__m128i v) {
    // 查找表：nibble_count[i] = popcount(i)，i = 0..15
    const __m128i lut = _mm_setr_epi8(
        0, 1, 1, 2, 1, 2, 2, 3,
        1, 2, 2, 3, 2, 3, 3, 4
    );
    const __m128i low_mask = _mm_set1_epi8(0x0F);

    __m128i lo = _mm_and_si128(v, low_mask);           // 低 4 位
    __m128i hi = _mm_and_si128(_mm_srli_epi16(v, 4), low_mask);  // 高 4 位

    __m128i cnt_lo = _mm_shuffle_epi8(lut, lo);  // 查表: popcount(低4位)
    __m128i cnt_hi = _mm_shuffle_epi8(lut, hi);  // 查表: popcount(高4位)

    return _mm_add_epi8(cnt_lo, cnt_hi);  // 每字节的 popcount
}

// 统计整个数组的 1 的总数
uint64_t popcount_array_sse(const uint8_t *data, int n) {
    __m128i total = _mm_setzero_si128();
    __m128i zero  = _mm_setzero_si128();

    int i = 0;
    for (; i + 15 < n; ) {
        __m128i acc = _mm_setzero_si128();
        // 内层循环：累加最多 255 次，避免 uint8 溢出
        int inner_end = (i + 255 * 16 < n) ? i + 255 * 16 : n;
        for (; i + 15 < inner_end; i += 16) {
            __m128i v = _mm_loadu_si128((__m128i*)&data[i]);
            acc = _mm_add_epi8(acc, popcount_16bytes(v));
        }
        // 将 uint8 累加器扩展到更宽的类型
        total = _mm_add_epi64(total, _mm_sad_epu8(acc, zero));
    }

    // 提取结果
    uint64_t result = (uint64_t)_mm_extract_epi64(total, 0)
                    + (uint64_t)_mm_extract_epi64(total, 1);

    // 处理剩余元素
    for (; i < n; i++) {
        result += __builtin_popcount(data[i]);
    }
    return result;
}

这里有一个精妙的细节：_mm_sad_epu8(acc, zero) 将 16 个 uint8 值求和，结果放在两个 uint64 字段中。这是 SIMD 中”水平求和”的常见技巧。

3.3 应用二：Base64 编码的 SIMD 加速

Base64 编码的核心是将 6-bit 值映射到 ASCII 字符。标量版本用查找表：

// 标量 Base64 编码
static const char base64_table[] =
    "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/";

void base64_encode_scalar(const uint8_t *src, char *dst, int len) {
    for (int i = 0; i < len; i += 3) {
        uint32_t val = (src[i] << 16) | (src[i+1] << 8) | src[i+2];
        dst[0] = base64_table[(val >> 18) & 0x3F];
        dst[1] = base64_table[(val >> 12) & 0x3F];
        dst[2] = base64_table[(val >>  6) & 0x3F];
        dst[3] = base64_table[(val      ) & 0x3F];
        dst += 4;
    }
}

SIMD 版本不能直接用 pshufb（因为查找表有 64 个条目，超过 16 字节），但可以用分段查表 + 偏移的方式：

// SIMD Base64 编码的核心：将 6-bit 值转换为 ASCII
// 分为 5 个区间：
//   0-25  → 'A'-'Z' (加 65)
//   26-51 → 'a'-'z' (加 71)
//   52-61 → '0'-'9' (减 4)
//   62    → '+'     (减 19)
//   63    → '/'     (减 16)
__m256i base64_encode_6bit_to_ascii_avx2(__m256i input) {
    // 用 compare + mask 判断每个 6-bit 值落在哪个区间
    __m256i result;

    // 构建偏移量：根据输入范围选择不同的偏移
    __m256i range_0_25  = _mm256_cmpgt_epi8(_mm256_set1_epi8(26), input);
    __m256i range_26_51 = _mm256_andnot_si256(range_0_25,
        _mm256_cmpgt_epi8(_mm256_set1_epi8(52), input));

    // 选择偏移值
    __m256i offset = _mm256_set1_epi8(71);  // 默认: a-z 范围
    offset = _mm256_blendv_epi8(offset, _mm256_set1_epi8(65), range_0_25);
    // 对 52-63 范围需要更精细的处理...

    result = _mm256_add_epi8(input, offset);
    return result;
}

实际的 SIMD Base64 实现（如 Wojciech Mula 的方案）还涉及 3 字节到 4 字节的重排（12 个输入字节变成 16 个 6-bit 索引），这部分也大量使用 pshufb 进行字节重排。

3.4 pshufb 的 AVX2 注意事项

_mm256_shuffle_epi8 在 AVX2 中是两个独立的 128-bit 操作，不会跨 lane 查找：

// AVX2 pshufb: 高 128-bit lane 和低 128-bit lane 各自独立
// index[0..15]  → 从 table[0..15]  查找
// index[16..31] → 从 table[16..31] 查找
// 不能用 index[0] 去查 table[20]！

这是 AVX2 最常见的坑之一。如果需要 32 字节的查找表，必须将同一个 16 字节表复制到两个 lane 中。

四、SIMD 前缀和与 Stream Compaction

4.1 标量前缀和

前缀和（prefix sum / scan）是并行计算中最基础的原语之一：

// 标量前缀和
void prefix_sum_scalar(const int32_t *input, int32_t *output, int n) {
    output[0] = input[0];
    for (int i = 1; i < n; i++) {
        output[i] = output[i-1] + input[i];
    }
}

这是一个典型的串行依赖——output[i] 依赖 output[i-1]。怎么向量化？

4.2 SIMD 前缀和：log(n) 步算法

在一个 SIMD 寄存器内（比如 AVX2 的 8 个 int32 元素），前缀和可以用 log2(8) = 3 步完成：

// AVX2 寄存器内前缀和（8 个 int32 元素）
__m256i prefix_sum_epi32(__m256i v) {
    // 步骤 1: 偏移 1，加
    // [a, b, c, d, e, f, g, h]
    // + [0, a, b, c, d, e, f, g]  (左移 1 个元素)
    // = [a, a+b, b+c, c+d, d+e, e+f, f+g, g+h]
    __m256i shifted;

    // 注意：AVX2 没有跨 lane 的字节移位
    // 需要用 permute + blend 来实现

    // Step 1: shift by 1 element (4 bytes)
    // _mm256_slli_si256 只在 128-bit lane 内移位
    shifted = _mm256_slli_si256(v, 4);
    // 还需要把 lane 0 的最高元素搬到 lane 1 的最低位
    // 用 permute2x128 + blend 实现跨 lane 传递
    __m256i cross = _mm256_permute2x128_si256(v, v, 0x08);
    cross = _mm256_srli_si256(cross, 12);  // 取出 lane 0 的最后一个元素
    shifted = _mm256_or_si256(shifted, cross);
    v = _mm256_add_epi32(v, shifted);

    // Step 2: shift by 2 elements (8 bytes)
    shifted = _mm256_slli_si256(v, 8);
    cross = _mm256_permute2x128_si256(v, v, 0x08);
    cross = _mm256_srli_si256(cross, 8);
    shifted = _mm256_or_si256(shifted, cross);
    v = _mm256_add_epi32(v, shifted);

    // Step 3: shift by 4 elements (跨 lane)
    shifted = _mm256_permute2x128_si256(_mm256_setzero_si256(), v, 0x02);
    v = _mm256_add_epi32(v, shifted);

    return v;
}

注意这里的跨 lane 操作——AVX2 的 256-bit 寄存器被分成两个 128-bit lane，大部分 shuffle/shift 操作不能跨 lane。这使得 SIMD 前缀和比理论上复杂得多。

4.3 大数组的前缀和

对于大数组，策略是分块处理：

每个 SIMD 寄存器内做局部前缀和。
提取每块的总和。
对块总和做前缀和（可以递归或标量处理）。
将块前缀和广播回每个块，加到局部结果上。

// 大数组前缀和的分块策略（伪代码）
void prefix_sum_large(const int32_t *input, int32_t *output, int n) {
    const int BLOCK = 8;  // AVX2 每组 8 个 int32
    int num_blocks = n / BLOCK;

    // 阶段 1: 各块内前缀和，同时记录块总和
    int32_t block_sums[num_blocks];
    for (int b = 0; b < num_blocks; b++) {
        __m256i v = _mm256_load_si256((__m256i*)&input[b * BLOCK]);
        __m256i ps = prefix_sum_epi32(v);
        _mm256_store_si256((__m256i*)&output[b * BLOCK], ps);
        // 提取最后一个元素作为块总和
        block_sums[b] = _mm256_extract_epi32(ps, 7);
    }

    // 阶段 2: 对块总和做前缀和（块数通常较少，标量即可）
    for (int b = 1; b < num_blocks; b++) {
        block_sums[b] += block_sums[b - 1];
    }

    // 阶段 3: 将块前缀和广播加回
    for (int b = 1; b < num_blocks; b++) {
        __m256i add_val = _mm256_set1_epi32(block_sums[b - 1]);
        __m256i ps = _mm256_load_si256((__m256i*)&output[b * BLOCK]);
        ps = _mm256_add_epi32(ps, add_val);
        _mm256_store_si256((__m256i*)&output[b * BLOCK], ps);
    }
}

4.4 Stream Compaction

Stream compaction 是另一个核心原语：给定一组数据和一个 mask，只保留 mask 为 true 的元素，连续输出。

标量版本很简单：

// 标量 stream compaction
int compact_scalar(const int32_t *input, int32_t *output,
                   const int *mask, int n) {
    int j = 0;
    for (int i = 0; i < n; i++) {
        if (mask[i]) {
            output[j++] = input[i];
        }
    }
    return j;
}

SIMD 版本需要一个查找表，将 mask 的 bit 模式映射到 permute 索引：

// AVX2 stream compaction: 使用 permutevar 重排
// 预计算的查找表：对于每个 8-bit mask，存储对应的 permute 索引
// 例如 mask = 0b10110001 → 保留位置 0, 4, 5, 7
// → permute 索引 = [0, 4, 5, 7, ?, ?, ?, ?]
// 这个表有 256 个条目，每个 32 字节

// 简化版本：对 4 个 int32 做 compaction（SSE）
int compact_sse(const int32_t *input, int32_t *output, __m128i mask) {
    int bitmask = _mm_movemask_ps(_mm_castsi128_ps(mask));

    // 预计算的 shuffle 表（4-bit mask → 4 个索引）
    static const int32_t shuffle_table[16][4] = {
        {0,0,0,0}, {0,0,0,0}, {1,0,0,0}, {0,1,0,0},  // 0000-0011
        {2,0,0,0}, {0,2,0,0}, {1,2,0,0}, {0,1,2,0},  // 0100-0111
        {3,0,0,0}, {0,3,0,0}, {1,3,0,0}, {0,1,3,0},  // 1000-1011
        {2,3,0,0}, {0,2,3,0}, {1,2,3,0}, {0,1,2,3},  // 1100-1111
    };
    static const int count_table[16] = {
        0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
    };

    __m128i shuf = _mm_loadu_si128((__m128i*)shuffle_table[bitmask]);
    __m128i data = _mm_loadu_si128((__m128i*)input);
    __m128i result = _mm_permutevar_ps(_mm_castsi128_ps(data),shuf);
    _mm_storeu_si128((__m128i*)output, _mm_castps_si128(result));

    return count_table[bitmask];
}

在 AVX-512 中，stream compaction 有硬件原生支持——_mm512_mask_compress_epi32，一条指令完成。这是 AVX-512 最令人兴奋的能力之一。

五、Gather/Scatter 与跨 Lane 操作

5.1 AVX2 Gather

AVX2 引入了 gather 指令，可以按索引从内存加载数据：

// _mm256_i32gather_epi32: 按索引加载 8 个 int32
// 等价于: for (int i = 0; i < 8; i++) result[i] = base[index[i]]
void gather_example(const int32_t *table, const int32_t *indices,
                    int32_t *output, int n) {
    for (int i = 0; i < n; i += 8) {
        __m256i idx = _mm256_load_si256((__m256i*)&indices[i]);
        __m256i val = _mm256_i32gather_epi32(table, idx, 4);
        _mm256_store_si256((__m256i*)&output[i], val);
    }
}

5.2 AVX2 Gather 的性能现实

AVX2 的 gather 指令在微架构层面并不是一条”真正的”向量加载——它通常被分解为多个标量加载操作。在 Haswell/Broadwell 上，vgatherdps（8x float gather）的吞吐量约为每 9-12 个周期一次，而 8 次独立标量加载可能也差不多。

实测对比：

// 手动 gather vs 硬件 gather 性能对比
__m256i manual_gather(const int32_t *base, __m256i indices) {
    int32_t idx[8];
    _mm256_storeu_si256((__m256i*)idx, indices);
    return _mm256_setr_epi32(
        base[idx[0]], base[idx[1]], base[idx[2]], base[idx[3]],
        base[idx[4]], base[idx[5]], base[idx[6]], base[idx[7]]);
}

微架构	硬件 gather (cycles/op)	手动 gather (cycles/op)	备注
Haswell	~12	~10-14	差别不大
Skylake	~7	~10-14	硬件 gather 开始有优势
Ice Lake	~5	~10-14	硬件 gather 明显更快
Zen 3	~6	~10-14	AMD 也优化了 gather

结论：在 Skylake 及更新的处理器上，硬件 gather 通常值得使用。在 Haswell 上，手动 gather 可能更好。

5.3 跨 Lane 操作的代价

AVX2 的 256-bit 寄存器在硬件上是两个 128-bit lane 的组合。大部分 shuffle 操作只在 lane 内进行。跨 lane 的操作有额外代价：

// Lane 内操作（快，1 cycle）:
_mm256_shuffle_epi32(v, imm)    // 每个 lane 内 32-bit shuffle
_mm256_shufflelo_epi16(v, imm)  // 低 lane 内 16-bit shuffle
_mm256_unpacklo_epi32(a, b)     // lane 内交叉排列

// 跨 Lane 操作（慢，3 cycles on Haswell）:
_mm256_permutevar8x32_epi32(v, idx)  // 全 256-bit 32-bit permute
_mm256_permute2x128_si256(a, b, imm) // 交换/选择 128-bit lane
_mm256_permute4x64_epi64(v, imm)     // 全 256-bit 64-bit permute

跨 lane 操作在不同微架构上的延迟：

操作	Haswell	Skylake	Zen 2
`vperm2i128`	3	3	1
`vpermd`	3	3	3
`vpermq`	3	3	2

设计原则：尽量在 lane 内完成操作，只在必要时才跨 lane。如果算法需要频繁跨 lane，考虑是否可以改变数据布局来避免。

5.4 Scatter 的缺失与替代

AVX2 没有 scatter 指令。如果需要按索引写入，只能用标量：

// AVX2 没有 scatter，必须手动实现
void manual_scatter(int32_t *base, __m256i indices, __m256i values) {
    int32_t idx[8], val[8];
    _mm256_storeu_si256((__m256i*)idx, indices);
    _mm256_storeu_si256((__m256i*)val, values);

    for (int i = 0; i < 8; i++) {
        base[idx[i]] = val[i];
    }
}

AVX-512 才引入了 scatter 指令：_mm512_i32scatter_epi32。

六、AVX-512 新能力

AVX-512 不仅仅是”更宽的 AVX2”。它引入了几个根本性的新概念。

6.1 Mask Register（k0-k7）

AVX-512 有 8 个专用的 mask 寄存器 k0-k7，每一位对应一个向量元素。几乎每条 AVX-512 指令都可以带一个 mask 操作数：

// AVX-512 masking: 两种模式

// 1. Zero masking: mask 为 0 的元素变为 0
__m512i result_z = _mm512_maskz_add_epi32(mask, a, b);
// result[i] = (mask & (1 << i)) ? a[i] + b[i] : 0

// 2. Merge masking: mask 为 0 的元素保留原值
__m512i result_m = _mm512_mask_add_epi32(src, mask, a, b);
// result[i] = (mask & (1 << i)) ? a[i] + b[i] : src[i]

这比 AVX2 的 blend 方式优雅得多——不需要单独的 compare + blend 步骤，mask 直接嵌入运算指令。

// AVX2 的条件加法（3 条指令）
__m256i mask = _mm256_cmpgt_epi32(a, threshold);
__m256i sum  = _mm256_add_epi32(a, b);
__m256i result = _mm256_blendv_epi8(a, sum, mask);  // mask ? sum : a

// AVX-512 的条件加法（2 条指令，且 mask 是寄存器）
__mmask16 mask = _mm512_cmpgt_epi32_mask(a, threshold);
__m512i result = _mm512_mask_add_epi32(a, mask, a, b);

6.2 Compress / Expand 指令

这是 AVX-512 最具革命性的指令之一——硬件级 stream compaction：

// _mm512_mask_compress_epi32: 将 mask 为 1 的元素压缩到连续位置
// 输入: [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p]
// mask: 1 0 1 1 0 0 1 0 0 1 0 0 0 1 0 1
// 输出: [a, c, d, g, j, n, p, ?, ?, ?, ?, ?, ?, ?, ?, ?]

void compress_example(const int32_t *input, int32_t *output,
                      __mmask16 mask) {
    __m512i data = _mm512_load_si512(input);
    _mm512_mask_compressstoreu_epi32(output, mask, data);
    // output 中连续存放所有 mask 为 1 的元素
}

// Expand 是反操作：将连续数据展开到 mask 为 1 的位置
void expand_example(const int32_t *input, int32_t *output,
                    __mmask16 mask) {
    __m512i data = _mm512_loadu_si512(input);
    __m512i result = _mm512_mask_expand_epi32(
        _mm512_setzero_si512(), mask, data);
    _mm512_store_si512(output, result);
}

有了 compress 指令，很多之前需要复杂查找表的操作变得极其简洁：

// 过滤数组中的正数——AVX-512 版本
int filter_positive_avx512(const int32_t *input, int32_t *output, int n) {
    int j = 0;
    __m512i zero = _mm512_setzero_si512();

    for (int i = 0; i < n; i += 16) {
        __m512i data = _mm512_loadu_si512(&input[i]);
        __mmask16 mask = _mm512_cmpgt_epi32_mask(data, zero);

        _mm512_mask_compressstoreu_epi32(&output[j], mask, data);
        j += _mm_popcnt_u32(mask);  // 计算保留了多少个元素
    }
    return j;
}

6.3 Conflict Detection

vpconflictd / vpconflictq 检测向量元素中是否有重复值——这对于直方图等”散乱写入”操作非常重要：

// 构建直方图时的冲突检测
void histogram_avx512(const int32_t *data, int32_t *hist, int n) {
    for (int i = 0; i < n; i += 16) {
        __m512i indices = _mm512_loadu_si512(&data[i]);

        // conflict[j] 的 bit k 为 1 表示 indices[k] == indices[j] 且 k < j
        __m512i conflict = _mm512_conflict_epi32(indices);

        // 无冲突的元素可以直接 gather-add-scatter
        __mmask16 no_conflict = _mm512_cmpeq_epi32_mask(
            conflict, _mm512_setzero_si512());

        __m512i old_vals = _mm512_mask_i32gather_epi32(
            _mm512_setzero_si512(), no_conflict, indices, hist, 4);
        __m512i new_vals = _mm512_add_epi32(old_vals, _mm512_set1_epi32(1));
        _mm512_mask_i32scatter_epi32(hist, no_conflict, indices, new_vals, 4);

        // 有冲突的元素需多轮处理，每轮取出一批不重复元素
        __mmask16 remaining = ~no_conflict & 0xFFFF;
        while (remaining) {
            int idx = __builtin_ctz(remaining);
            remaining &= (remaining - 1);
            hist[data[i + idx]]++;
        }
    }
}

个人观点：AVX-512 的 compress/expand 和 conflict detection 解决了 SIMD 编程中两个最棘手的问题——不规则输出和散乱写入。但 AVX-512 的部署现状仍然令人头疼：频率降低（thermal throttling）、不是所有 CPU 都支持、AMD 到 Zen 4 才完整支持。在 2025 年，我的建议是：如果你能控制部署环境（服务器端），积极使用 AVX-512；如果是客户端软件，还是老老实实写 AVX2 + 回退路径。

七、ARM NEON vs x86 SSE/AVX 的跨平台策略

7.1 NEON 与 SSE/AVX 的核心差异

Apple Silicon 的崛起让跨平台 SIMD 变成了现实需求。NEON 和 SSE/AVX 有几个关键差异：

特性	x86 SSE/AVX	ARM NEON
寄存器宽度	128/256/512 bit	128 bit (SVE: 128-2048)
寄存器数量	16 (SSE) / 32 (AVX-512)	32
整数乘法	受限（无 32x32→64）	丰富
饱和算术	有（paddusb 等）	有且更全面
条件执行	blend/mask	类似（BSL 指令）
FMA	AVX2+	有
pshufb 等价	SSSE3+	vtbl/vtbx
Gather/Scatter	AVX2+ / AVX-512	SVE2

NEON 的 vtbl / vtbx 是 pshufb 的等价物，但支持更大的查找表（最多 4 个寄存器 = 64 字节）。

7.2 手动跨平台抽象

最原始的方式是用 #ifdef 包裹不同平台的 intrinsic，但维护成本随代码量增长而爆炸：

#if defined(__x86_64__)
    #include <immintrin.h>
    typedef __m128i vec128i;
    #define vec_add_i32(a, b)     _mm_add_epi32(a, b)
    #define vec_shuffle_bytes(t,i) _mm_shuffle_epi8(t, i)
#elif defined(__aarch64__)
    #include <arm_neon.h>
    typedef int32x4_t vec128i;
    #define vec_add_i32(a, b)     vaddq_s32(a, b)
    #define vec_shuffle_bytes(t,i) vqtbl1q_u8(t, i)
#endif

7.3 Google Highway

Highway 是目前最推荐的跨平台 SIMD 库。它的核心设计理念：

// Highway 示例：跨平台 SIMD 代码
#include "hwy/highway.h"
namespace hn = hwy::HWY_NAMESPACE;

void AddArrays(const float* a, const float* b, float* result, size_t n) {
    const hn::ScalableTag<float> d;  // 自动选择最宽的可用 SIMD
    const size_t lanes = hn::Lanes(d);

    size_t i = 0;
    for (; i + lanes <= n; i += lanes) {
        auto va = hn::Load(d, a + i);
        auto vb = hn::Load(d, b + i);
        hn::Store(hn::Add(va, vb), d, result + i);
    }
    // 处理剩余元素
    for (; i < n; i++) {
        result[i] = a[i] + b[i];
    }
}

Highway 的优势： - 编译时分发：针对每个目标 ISA 生成最优代码。 - 运行时分发：同一个二进制文件可以在不同 CPU 上运行最优路径。 - 覆盖全面：支持 x86（SSE2-AVX-512）、ARM（NEON/SVE）、WASM SIMD、RISC-V V。 - 接口一致：不需要 #ifdef 地狱。

7.4 SIMDe——头文件级兼容层

SIMDe（SIMD Everywhere）的思路不同——它用头文件将 x86 intrinsic 翻译为其他平台的等价操作：

// 使用 SIMDe: 写 x86 intrinsic，跨平台编译
#include <simde/x86/avx2.h>

void process(const float *input, float *output, int n) {
    for (int i = 0; i < n; i += 8) {
        simde__m256 v = simde_mm256_loadu_ps(&input[i]);
        v = simde_mm256_mul_ps(v, simde_mm256_set1_ps(2.0f));
        simde_mm256_storeu_ps(&output[i], v);
    }
}
// 在 ARM 上编译时，simde_mm256_* 会被翻译为 NEON 指令

SIMDe 的适用场景是已有大量 x86 intrinsic 代码，想以最小改动支持 ARM。但翻译出来的代码未必能完美匹配 ARM 指令特性。

7.5 std::experimental::simd

C++ 标准正在努力将 SIMD 纳入标准库（P1928）。GCC 和 Clang 都有实验性支持：

#include <experimental/simd>
namespace stdx = std::experimental;

void add_arrays(const float *a, const float *b, float *c, int n) {
    using V = stdx::native_simd<float>;
    constexpr int lanes = V::size();

    int i = 0;
    for (; i + lanes <= n; i += lanes) {
        V va(a + i, stdx::element_aligned);
        V vb(b + i, stdx::element_aligned);
        V vc = va + vb;
        vc.copy_to(c + i, stdx::element_aligned);
    }
    for (; i < n; i++) c[i] = a[i] + b[i];
}

个人观点：在 2025 年中期的现实选择中，我推荐 Google Highway。它成熟、活跃维护、覆盖全面。SIMDe 适合迁移遗留代码。std::experimental::simd 还需要几年才能稳定。如果项目只需要支持 x86，直接写 intrinsic 也完全可以——代码更透明，性能更可预测。

八、完整示例：各模式的标量 / SSE / AVX2 / AVX-512 对比

8.1 条件替换：将负数替换为零（ReLU）

#include <immintrin.h>
#include <stdint.h>
#include <string.h>

// ===== 标量版本 =====
void relu_scalar(float *data, int n) {
    for (int i = 0; i < n; i++) {
        if (data[i] < 0.0f) data[i] = 0.0f;
    }
}

// ===== SSE 版本 =====
void relu_sse(float *data, int n) {
    __m128 zero = _mm_setzero_ps();
    int i = 0;
    for (; i + 3 < n; i += 4) {
        __m128 v = _mm_loadu_ps(&data[i]);
        v = _mm_max_ps(v, zero);
        _mm_storeu_ps(&data[i], v);
    }
    for (; i < n; i++) {
        if (data[i] < 0.0f) data[i] = 0.0f;
    }
}

// ===== AVX2 版本 =====
void relu_avx2(float *data, int n) {
    __m256 zero = _mm256_setzero_ps();
    int i = 0;
    for (; i + 7 < n; i += 8) {
        __m256 v = _mm256_loadu_ps(&data[i]);
        v = _mm256_max_ps(v, zero);
        _mm256_storeu_ps(&data[i], v);
    }
    for (; i < n; i++) {
        if (data[i] < 0.0f) data[i] = 0.0f;
    }
}

// ===== AVX-512 版本 =====
void relu_avx512(float *data, int n) {
    __m512 zero = _mm512_setzero_ps();
    int i = 0;
    for (; i + 15 < n; i += 16) {
        __m512 v = _mm512_loadu_ps(&data[i]);
        v = _mm512_max_ps(v, zero);
        _mm512_storeu_ps(&data[i], v);
    }
    for (; i < n; i++) {
        if (data[i] < 0.0f) data[i] = 0.0f;
    }
}

8.2 查表：字节分类（判断 ASCII 字符类型）

// 将 ASCII 字符分类: 数字→1, 大写字母→2, 小写字母→3, 其他→0
// 利用 pshufb 的 4-bit 查表能力

// ===== 标量版本 =====
void classify_ascii_scalar(const uint8_t *input, uint8_t *output, int n) {
    for (int i = 0; i < n; i++) {
        uint8_t c = input[i];
        if (c >= '0' && c <= '9')      output[i] = 1;
        else if (c >= 'A' && c <= 'Z') output[i] = 2;
        else if (c >= 'a' && c <= 'z') output[i] = 3;
        else                           output[i] = 0;
    }
}

// ===== SSE 版本（用减法 + 无符号比较技巧）=====
// 核心思路：sub + saturating compare 实现范围判断
// c in [lo, hi] 等价于 (uint8_t)(c - lo) <= (hi - lo)
// 对每个范围生成 mask，然后逐层 blend
// 完整实现见第三节的 pshufb 字符分类方案

8.3 水平求和

// 将向量内所有元素求和——SIMD 常见的收尾操作

// ===== SSE: 4 个 float 求和 =====
float hsum_sse(__m128 v) {
    __m128 shuf = _mm_movehdup_ps(v);    // [v1, v1, v3, v3]
    __m128 sums = _mm_add_ps(v, shuf);   // [v0+v1, ?, v2+v3, ?]
    shuf = _mm_movehl_ps(shuf, sums);    // [v2+v3, ?, ?, ?]
    sums = _mm_add_ss(sums, shuf);       // [v0+v1+v2+v3, ?, ?, ?]
    return _mm_cvtss_f32(sums);
}

// ===== AVX2: 8 个 float 求和 =====
float hsum_avx2(__m256 v) {
    // 先将高 128 加到低 128
    __m128 hi = _mm256_extractf128_ps(v, 1);
    __m128 lo = _mm256_castps256_ps128(v);
    __m128 sum128 = _mm_add_ps(lo, hi);
    return hsum_sse(sum128);
}

// ===== AVX-512: 16 个 float 求和 =====
float hsum_avx512(__m512 v) {
    return _mm512_reduce_add_ps(v);  // AVX-512 有专用的 reduce 指令
}

九、工程实践中的陷阱

9.1 常见问题速查表

问题	症状	原因	解决方案
段错误	对齐加载崩溃	数据未对齐到 16/32/64 字节	用 `_mm_loadu` 或对齐分配
性能不升反降	SIMD 比标量慢	频繁的 pack/unpack 开销	改数据布局为 SoA
结果不正确	高位元素垃圾值	尾部未填充的数据	用 maskload 或填充到对齐长度
编译器优化消失	-O2 下标量也很快	编译器自动向量化了	检查 `-fno-tree-vectorize` 对比
AVX 频率降低	AVX-512 吞吐不如预期	Thermal throttling	profile 实测，别只看理论值
跨 lane 瓶颈	pshufb 结果错误	AVX2 pshufb 不跨 128-bit lane	表复制到两个 lane
Gather 无加速	gather 不比标量快	老微架构的 gather 实现差	Haswell 上用手动 gather
Denormal 性能	运算突然变慢	产生了非规格化浮点数	设置 FTZ/DAZ 标志
精度问题	结果与标量版本不同	SIMD 的运算顺序不同	FMA 改变了结合性
编译器报错	intrinsic 未定义	缺少 ISA 编译标志	`-mavx2 -mfma` 等

9.2 对齐分配

// 对齐内存分配——必须用对应的 free 释放

// C11 标准方式
#include <stdlib.h>
float *data = aligned_alloc(32, n * sizeof(float));  // 32 字节对齐
// free(data);  // 标准 free 可以释放

// POSIX 方式
#include <stdlib.h>
float *data;
posix_memalign((void**)&data, 32, n * sizeof(float));
// free(data);

// Intel 方式
#include <immintrin.h>
float *data = (float*)_mm_malloc(n * sizeof(float), 32);
// _mm_free(data);  // 必须用 _mm_free

// Windows 方式
float *data = (float*)_aligned_malloc(n * sizeof(float), 32);
// _aligned_free(data);  // 必须用 _aligned_free

9.3 处理尾部元素

当数组长度不是 SIMD 宽度的倍数时，有三种常见策略：

// 方法 1: 标量回退（最安全）
int i = 0;
for (; i + 7 < n; i += 8) { /* AVX2 向量处理 */ }
for (; i < n; i++) { /* 标量处理剩余 */ }

// 方法 2: Masked load/store（AVX2）
if (i < n) {
    int remaining = n - i;
    __m256i mask = _mm256_cmpgt_epi32(
        _mm256_set1_epi32(remaining),
        _mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7));
    __m256i v = _mm256_maskload_epi32((int*)&data[i], mask);
    // 处理 v...
    _mm256_maskstore_epi32((int*)&data[i], mask, v);
}

// 方法 3: 重叠处理（最高效但需要 n >= SIMD_WIDTH）
// 对幂等操作，最后一组从 n-8 开始，与前一组重叠
if (n >= 8) {
    __m256 v = _mm256_loadu_ps(&data[n - 8]);
    v = _mm256_max_ps(v, _mm256_setzero_ps());
    _mm256_storeu_ps(&data[n - 8], v);
}

9.4 编译器标志和运行时检测

// 编译时：指定目标 ISA
// gcc/clang: -msse4.2 -mavx2 -mfma -mavx512f -mavx512bw
// MSVC: /arch:AVX2 或 /arch:AVX512

// 运行时检测 CPU 特性
#include <cpuid.h>

int has_avx2(void) {
    unsigned int eax, ebx, ecx, edx;
    __cpuid_count(7, 0, eax, ebx, ecx, edx);
    return (ebx >> 5) & 1;  // bit 5 of EBX = AVX2
}

int has_avx512f(void) {
    unsigned int eax, ebx, ecx, edx;
    __cpuid_count(7, 0, eax, ebx, ecx, edx);
    return (ebx >> 16) & 1;  // bit 16 of EBX = AVX-512F
}

// 函数指针分发
typedef void (*relu_fn)(float *, int);

relu_fn get_relu(void) {
    if (has_avx512f()) return relu_avx512;
    if (has_avx2())    return relu_avx2;
    return relu_sse;
}

9.5 Denormal 处理

// 非规格化浮点数会导致某些微架构性能骤降
// 设置 FTZ (Flush to Zero) 和 DAZ (Denormals Are Zero)
#include <pmmintrin.h>
void disable_denormals(void) {
    _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
    _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);
}

十、性能测量与验证方法论

10.1 基准测试的正确姿势

#include <time.h>
#include <stdio.h>

static inline uint64_t rdtsc(void) {
    unsigned int lo, hi;
    __asm__ volatile ("rdtsc" : "=a" (lo), "=d" (hi));
    return ((uint64_t)hi << 32) | lo;
}

void benchmark(const char *name, void (*fn)(float*, int),
               float *data, int n, int iterations) {
    for (int i = 0; i < 10; i++) fn(data, n);  // 预热

    struct timespec start, end;
    clock_gettime(CLOCK_MONOTONIC, &start);
    uint64_t cycles_start = rdtsc();

    for (int i = 0; i < iterations; i++) fn(data, n);

    uint64_t cycles_end = rdtsc();
    clock_gettime(CLOCK_MONOTONIC, &end);

    double elapsed = (end.tv_sec - start.tv_sec)
                   + (end.tv_nsec - start.tv_nsec) / 1e9;
    double cycles_per_elem = (double)(cycles_end - cycles_start)
                           / iterations / n;
    printf("%-20s %8.2f ms  %6.2f cycles/elem\n",
           name, elapsed * 1000.0 / iterations, cycles_per_elem);
}

10.2 验证正确性

个人观点：SIMD 代码一定要先写一个简单的标量参考实现，然后逐步向量化。每向量化一步就和标量版本对比。我见过太多人直接写全向量化版本，然后花三天 debug 一个 shuffle 索引错误。

10.3 使用 perf 分析

# SIMD 指令使用情况
perf stat -e instructions,cycles,fp_arith_inst_retired.256b_packed_single ./my_program

# 缓存命中率（数据布局优化相关）
perf stat -e cache-references,cache-misses,L1-dcache-load-misses ./my_program

十一、设计模式总结与选择指南

11.1 模式速查

模式	适用场景	关键指令	性能提升
SoA 布局	数据密集型批量处理	load/store	2-8x（减少内存浪费）
Mask+Blend	条件分支	cmpgt + blendv	1.5-4x（消除分支预测失败）
pshufb 查表	4-bit 范围映射	shuffle_epi8	4-16x（16 路并行查表）
前缀和	累积求和 / scan 操作	shift + add	2-4x
Stream Compaction	过滤/压缩数据	compress (AVX-512)	3-8x
Gather	间接寻址 / 查大表	i32gather	1-3x（依赖微架构）

11.2 决策流程

面对一个想要 SIMD 优化的算法，我建议按以下顺序思考：

数据布局能改成 SoA 吗？ 如果能，这通常是收益最大的单一改动。
有条件分支吗？ 用 mask + blend 替代。如果只有 min/max/abs，直接用对应指令。
有小查找表（16 字节以内）吗？ 用 pshufb。
有规约操作（求和/最大值等）吗？ 先向量化主循环，最后水平规约。
有数据过滤/压缩需求吗？ AVX-512 用 compress，AVX2 用查找表 + permute。
有散乱访问吗？ 先尝试改布局消除散乱访问；如果不行，用 gather。
需要跨平台吗？ 用 Highway 或 SIMDe。

11.3 何时不要用 SIMD

SIMD 不是万能药。以下情况手写 SIMD 可能不值得：

编译器已经做得很好：简单的循环（如数组加法），GCC/Clang 的自动向量化可能已经生成了接近最优的代码。先检查 -O2 -ftree-vectorize 的输出。
数据量太小：SIMD 有启动开销（寄存器设置等），如果数组只有几十个元素，收益可能被开销抵消。
算法内在串行：某些递推关系的依赖链无法打破，SIMD 帮不上忙。
分支密集且不规则：如果 mask 中只有 1/8 是 true，你在”浪费地”计算 7/8 的分支。

个人观点：SIMD 优化的最大敌人不是技术难度，而是工程复杂度。一段 SIMD 代码的维护成本是等价标量代码的 3-5 倍。在动手之前，先确认这段代码是真正的性能热点，而且自动向量化确实不够好。

十二、参考资料

Intel Intrinsics Guide：https://www.intel.com/content/www/us/en/docs/intrinsics-guide/
Agner Fog, “Instruction tables”：https://www.agner.org/optimize/instruction_tables.pdf
Wojciech Mula, “SIMD-friendly algorithms for substring searching”：http://0x80.pl/articles/simd-strfind.html
Daniel Lemire, “Stream VByte: Faster Byte-Oriented Integer Compression”
Google Highway 项目：https://github.com/google/highway
SIMDe 项目：https://github.com/simd-everywhere/simde
Geoff Langdale & Daniel Lemire, “Parsing Gigabytes of JSON per Second”（simdjson）
Wojciech Mula, “Base64 encoding and decoding with SIMD”：http://0x80.pl/notesen/2016-01-12-sse-base64-encoding.html
Intel Architecture Instruction Set Extensions Programming Reference
ARM NEON Programmer’s Guide：https://developer.arm.com/documentation/den0018/a

算法系列导航：上一篇：无分支编程 | 下一篇：随机化算法

相关阅读：SIMD 字符串搜索 | 向量化哈希

目录