SIMD

SIMD 为什么快、快到哪里为止；tiling 怎样把瓶颈从内存拉回 CPU；哪些循环能向量化、哪些不能、不能的时候怎么改。

适用范围：C++，Linux x86-64 + GCC/Clang，以 AVX2 和 AVX-512 为主线。微架构参数（延迟、端口数）以 Intel Skylake 为例，属于实现细节，随 CPU 型号变化。

实测状态：本文全部代码、耗时和流量数字均未实测。数字按公开的微架构参数估算，用于说明量级。

#	章节	主题
一	指令集与寄存器	SSE / AVX2 / AVX-512 / NEON / SVE，load 与 gather
二	SIMD 的加速来源	指令数、分支数、上限在哪
三	tiling：让数据留在 cache 里	cache tiling、register tiling
四	能用 SIMD 的场景	归约、掩码、字节查找、过滤压缩
五	用不了或不划算的场景	依赖链、指针追逐、冲突写入、变长编码
六	数据布局	AoS / SoA / AoSoA
七	实现细节	对齐、尾部、多累加器、gather、NT store、浮点
八	三种写法	自动向量化、intrinsics、可移植库
九	部署陷阱	运行时分发、ODR、vzeroupper、降频
十	验证方法	编译器报告、反汇编、perf
十一	面试速答卡	高频问题浓缩答案

一、指令集与寄存器

SIMD（Single Instruction, Multiple Data）用一条指令同时处理一个向量寄存器里的多个元素。每个元素占据的位置叫一个 lane。

指令集	向量宽度	寄存器	float / double 个数	说明
SSE2	128 位	xmm0–xmm15	4 / 2	x86-64 基线，所有 64 位 CPU 都有
AVX2 + FMA	256 位	ymm0–ymm15	8 / 4	`-march=x86-64-v3`，Haswell（2013）起
AVX-512	512 位	zmm0–zmm31，k0–k7	16 / 8	`-march=x86-64-v4`，有独立的掩码寄存器
NEON	128 位	v0–v31	4 / 2	AArch64 基线
SVE / SVE2	128–2048 位	z0–z31，p0–p15	由硬件决定	向量长度在运行时确定

AVX-512 相对 AVX2 的差别不只是宽度：

比较指令直接产出 k 掩码寄存器，每个 lane 占 1 位
几乎所有指令都支持掩码（mask 保留原值，maskz 清零）
新增 compress / expand、冲突检测（AVX-512CD）、64 位整数乘法（AVX-512DQ）

连续读取与 gather

向量寄存器从内存取数有两种方式：

方式	取到的元素	标量等价写法	指令
load	从一个地址开始的连续 N 个	`v[i] = p[i]`	`vmovups`、`vmovdqu`
gather	下标向量指定的 N 个任意位置	`v[i] = base[idx[i]]`	`vgatherdps`、`vpgatherdd`

load：一次访存，取连续 4 个

  p    → | a0 | a1 | a2 | a3 | a4 | a5 | ...
           └─────────────────┘
  向量  = | a0 | a1 | a2 | a3 |

gather：idx = (0, 5, 2, 9)，每个 lane 各访存一次

  base → | a0 | a1 | a2 | a3 | a4 | a5 | a6 | a7 | a8 | a9 | ...
           ↓         ↓              ↓                   ↓
  向量  = | a0 | a5 | a2 | a9 |        按 idx 的顺序排列

#include <immintrin.h>
#include <cstddef>
#include <cstdint>

// out[i] = table[idx[i]]
void lookup(const float* table, const int32_t* idx, float* out, size_t n) {
  size_t i = 0;
  for (; i + 8 <= n; i += 8) {
    __m256i vi = _mm256_loadu_si256((const __m256i*)(idx + i));  // 8 个下标
    __m256  v  = _mm256_i32gather_ps(table, vi, 4);              // 地址 = table + idx * 4
    _mm256_storeu_ps(out + i, v);
  }
  for (; i < n; ++i) out[i] = table[idx[i]];
}
// g++ -O3 -mavx2

_mm256_i32gather_ps 的第三个参数是下标的缩放系数，取 1、2、4 或 8，这里是 sizeof(float)。

scatter 是反方向的写入：base[idx[i]] = v[i]。AVX2 只有 gather，scatter 需要 AVX-512F。

gather 的 N 个元素仍要 N 次访存，代价见第七节。

二、SIMD 的加速来源

每个元素的指令数和分支数下降

void saxpy(float* __restrict y, const float* __restrict x, float a, size_t n) {
  for (size_t i = 0; i < n; ++i) y[i] += a * x[i];
}
// g++ -O3 -mavx2 -mfma

向量化后的循环体形态如下（典型输出形态，未实测）：

.L3:
    vmovups      ymm1, [rsi+rax]          ; 读 8 个 x
    vfmadd213ps  ymm1, ymm0, [rdi+rax]    ; ymm1 = a*x + y
    vmovups      [rdi+rax], ymm1          ; 写 8 个 y
    add          rax, 32
    cmp          rax, rdx
    jne          .L3

	标量	AVX2	AVX-512
每轮处理元素数	1	8	16
每轮指令数	约 6	约 6	约 6
每个元素的循环分支数	1	1/8	1/16

循环体的指令条数相同，每条指令处理的元素数是 8 倍或 16 倍。

数据相关的分支变成掩码

if (cond) x = a; else x = b; 在向量里写成「比较得掩码，按掩码选择」。循环体内没有跳转，分支预测失败的代价（每次约 15~20 周期）消失。第四节例 2给出完整代码。

上限：内存带宽

SIMD 只加快计算。数据不在 cache 里时，循环的速度由内存带宽决定，向量再宽也要等数据到达。

循环	瓶颈	向量宽度加倍的收益
数据在 L1/L2 内	CPU 执行端口	接近 2 倍
数据在 L3 内	介于两者之间	小于 2 倍
数据超出 L3	内存带宽	接近 0

把第三行的循环变成第一行，是 tiling 的工作。

三、tiling：让数据留在 cache 里

loop tiling（循环分块，也叫 cache blocking）改变循环的遍历顺序：把一小块数据用完再换下一块。计算量不变，数据从内存搬进 cache 的次数下降。

下面用矩阵乘法 C = A × B 说明，N = 2048，元素为 double，单个矩阵 2048 × 2048 × 8 = 32 MB。

版本 1：ijk

for (int i = 0; i < N; ++i)
  for (int j = 0; j < N; ++j) {
    double s = 0;
    for (int k = 0; k < N; ++k)
      s += A[i*N + k] * B[k*N + j];   // B 按列走，步长 N*8 = 16 KB
    C[i*N + j] = s;
  }

B[k*N + j] 在内存里不连续，相邻两次访问相距 N 个元素，向量化只能用 gather
s += 是浮点归约，不加 -ffast-math 时编译器不能改变加法顺序

版本 2：ikj

for (int i = 0; i < N; ++i)
  for (int k = 0; k < N; ++k) {
    double a = A[i*N + k];
    for (int j = 0; j < N; ++j)
      C[i*N + j] += a * B[k*N + j];   // C、B 都连续，无归约
  }

内层循环能向量化，但每个 i 都要扫完整个 B：

项目	估算
纯计算	2048³ ≈ 8.6e9 次 FMA；AVX2 每条 4 个 double、每周期 2 条，3 GHz 下约 0.36 s
内存流量	每个 `i` 扫一遍 B（32 MB），共 2048 × 32 MB ≈ 65 GB
内存下限	单线程带宽按 10~~20 GB/s 计，约 3~~6 s

内存耗时约为计算耗时的 10 倍，换成 AVX-512 几乎没有收益。

版本 3：cache tiling

#include <algorithm>
constexpr int T = 128;
for (int ii = 0; ii < N; ii += T)
 for (int kk = 0; kk < N; kk += T)
  for (int jj = 0; jj < N; jj += T)
   for (int i = ii; i < std::min(ii + T, N); ++i)
    for (int k = kk; k < std::min(kk + T, N); ++k) {
      double a = A[i*N + k];
      for (int j = jj; j < std::min(jj + T, N); ++j)
        C[i*N + j] += a * B[k*N + j];
    }

版本 2：固定 i，扫过整个 B            版本 3：固定一个 tile，128 个 i 复用它

        j →                                  j →
   ┌────────────────────┐              ┌────┬────┬────┬────┐
 k │ → → → → → → → → →  │            k │####│    │    │    │  #### = 128×128×8
 ↓ │ → → → → → → → → →  │            ↓ ├────┼────┼────┼────┤        = 128 KB
   │    ...2048 行...    │              │    │    │    │    │        放得进 L2
   │ → → → → → → → → →  │              ├────┼────┼────┼────┤
   └────────────────────┘              │    │    │    │    │
   32 MB，超出 L2 和 L3                  └────┴────┴────┴────┘

版本 2 里，扫到 B 的后半部分时前半部分已被挤出 cache，下一个 i 又从 B[0] 开始，数据要重新从内存搬。B 的每一行被从内存搬入 2048 次。

版本 3 里，B 的一个 tile 被第 1 个 i 搬进 L2，后面 127 个 i 直接从 L2 读。B 的每个 tile 只在外层 ii 换块时重新搬入，共 2048 / 128 = 16 次。

	B 每份数据从内存搬入的次数	B 的总流量	A + B + C 总流量
版本 2	2048	约 65 GB	约 65 GB
版本 3	16	约 0.5 GB	约 1.5 GB

tile 大小的选择：A、B、C 三个 tile 合计 3 × 128 KB = 384 KB，小于 L2（1~2 MB，实现相关）。T 太小则循环开销占比上升，T 太大则 tile 放不进 L2。

register tiling

版本 3 的内层循环里，每条 FMA 仍然配 2 次 load（C 和 B）和 1 次 store（C）。把 C 的一个 4×8 小块放进 8 个 ymm 寄存器，跨整个 k 循环不落内存：

#include <immintrin.h>
// C[0..4)[0..8) += A[0..4)[0..K) * B[0..K)[0..8)
// lda / ldb / ldc 是各矩阵一行的元素个数
void kernel_4x8(const double* A, int lda, const double* B, int ldb,
                double* C, int ldc, int K) {
  __m256d c[4][2];
  for (int r = 0; r < 4; ++r) {
    c[r][0] = _mm256_loadu_pd(C + r*ldc);
    c[r][1] = _mm256_loadu_pd(C + r*ldc + 4);
  }
  for (int k = 0; k < K; ++k) {
    __m256d b0 = _mm256_loadu_pd(B + k*ldb);
    __m256d b1 = _mm256_loadu_pd(B + k*ldb + 4);
    for (int r = 0; r < 4; ++r) {
      __m256d a = _mm256_broadcast_sd(A + r*lda + k);
      c[r][0] = _mm256_fmadd_pd(a, b0, c[r][0]);
      c[r][1] = _mm256_fmadd_pd(a, b1, c[r][1]);
    }
  }
  for (int r = 0; r < 4; ++r) {
    _mm256_storeu_pd(C + r*ldc,     c[r][0]);
    _mm256_storeu_pd(C + r*ldc + 4, c[r][1]);
  }
}
// g++ -O3 -mavx2 -mfma

	版本 3 内层（每 8 条 FMA）	`kernel_4x8`（每 8 条 FMA）
load	16	6（2 次 B + 4 次 A）
store	8	0
独立依赖链	无链，受 store 端口限制	8 条

累加器取 8 个的依据：Skylake 上 FMA 延迟 4 周期、每周期可执行 2 条。让两个端口都不空转，需要 4 × 2 = 8 条互相独立的依赖链。

OpenBLAS 和 BLIS 是这两层的组合：外层按 L2/L3 大小分块，最内层是按寄存器数量定大小的 micro-kernel。

tiling 的其他收益

收益	原因
TLB miss 减少	一个 tile 只跨几十个 4 KB 页
冲突 miss 减少	N 为 2 的幂时各行映射到同一组 cache set；把 tile 拷到连续缓冲区（packing）可消除
多线程可扩展	每个线程的 tile 在各自的 L2 里，不争抢共享的内存带宽

四、能用 SIMD 的场景

循环能向量化需要同时满足五个条件：

条件	不满足时的表现
数据连续、类型相同	只能用 gather，收益小
迭代之间独立，或只有可结合的归约	依赖链限制速度
分支能转成掩码	不同元素要走不同的代码路径
指针不重叠	编译器生成运行时检查，或放弃
运算有对应的向量指令	退回标量

例 1：逐元素运算与归约

double dot(const double* p, const double* q, size_t n) {
  double s = 0;
  for (size_t i = 0; i < n; ++i) s += p[i] * q[i];
  return s;
}

只用 -O3 时，浮点加法必须保持 ((s + x0) + x1) + ... 的顺序，整个循环是一条串行依赖链。让它向量化有三种办法：

办法	命令或写法	影响范围
放宽浮点语义	`-ffast-math` 或 `-fassociative-math -fno-signed-zeros -fno-trapping-math`	整个翻译单元
标注这一个循环	`#pragma omp simd reduction(+:s)`，编译加 `-fopenmp-simd`	单个循环
手写多个累加器	见第七节	单个函数

整数归约满足结合律，-O3 直接向量化。

例 2：分支转掩码

// 只累加买单的数量
int64_t sum_buy(const int32_t* qty, const uint8_t* side, size_t n) {
  int64_t s = 0;
  for (size_t i = 0; i < n; ++i)
    if (side[i] == 'B') s += qty[i];
  return s;
}

买卖方向随机时，这个 if 约 50% 预测失败。

标量的无分支写法：

int64_t sum_buy_branchless(const int32_t* qty, const uint8_t* side, size_t n) {
  int64_t s = 0;
  for (size_t i = 0; i < n; ++i) {
    int32_t mask = -(int32_t)(side[i] == 'B');  // 买单 0xFFFFFFFF，卖单 0
    s += qty[i] & mask;                         // 卖单贡献 0
  }
  return s;
}

比较结果是 0 或 1，取负后成为全 0 或全 1，再与 qty 按位与。每个元素执行同样的指令。

AVX2 写法（一次 8 个元素）：

#include <immintrin.h>
#include <cstdint>
#include <cstddef>

int64_t sum_buy_avx2(const int32_t* qty, const uint8_t* side, size_t n) {
  const __m128i vB = _mm_set1_epi8('B');
  __m256i acc = _mm256_setzero_si256();                 // 4 个 int64 累加器
  size_t i = 0;
  for (; i + 8 <= n; i += 8) {
    __m128i s8  = _mm_loadl_epi64((const __m128i*)(side + i)); // 读 8 个字节
    __m128i m8  = _mm_cmpeq_epi8(s8, vB);                      // 相等 0xFF，否则 0x00
    __m256i m32 = _mm256_cvtepi8_epi32(m8);                    // 符号扩展为 32 位掩码
    __m256i q   = _mm256_loadu_si256((const __m256i*)(qty + i));
    __m256i sel = _mm256_and_si256(q, m32);                    // 卖单清零

    // int32 扩成 int64 再累加，避免溢出
    __m256i lo = _mm256_cvtepi32_epi64(_mm256_castsi256_si128(sel));
    __m256i hi = _mm256_cvtepi32_epi64(_mm256_extracti128_si256(sel, 1));
    acc = _mm256_add_epi64(acc, _mm256_add_epi64(lo, hi));
  }
  __m128i s2 = _mm_add_epi64(_mm256_castsi256_si128(acc),
                             _mm256_extracti128_si256(acc, 1));
  int64_t sum = _mm_cvtsi128_si64(s2) + _mm_extract_epi64(s2, 1);

  for (; i < n; ++i)                                           // 尾部不足 8 个
    sum += qty[i] & -(int32_t)(side[i] == 'B');
  return sum;
}
// g++ -O3 -mavx2

一组数据走一遍：

步骤	lane 0	1	2	3	4	5	6	7
`side`	B	S	B	B	S	S	B	S
`m8`	FF	00	FF	FF	00	00	FF	00
`m32`	FFFFFFFF	0	FFFFFFFF	FFFFFFFF	0	0	FFFFFFFF	0
`qty`	10	20	30	40	50	60	70	80
`sel`	10	0	30	40	0	0	70	0

lo = (10, 0, 30, 40)，hi = (0, 0, 70, 0)，acc 增加 (10, 0, 100, 40)，合并得 150。

掩码扩展必须用符号扩展：0xFF 作为有符号字节是 -1，扩展到 32 位是 0xFFFFFFFF。零扩展（_mm256_cvtepu8_epi32）得到 0x000000FF，只保留 qty 的低 8 位。

AVX-512 写法的循环体（片段，需要 AVX-512F + BW + VL）：

__m128i   s8 = _mm_loadu_si128((const __m128i*)(side + i));   // 16 个字节
__mmask16 m  = _mm_cmpeq_epi8_mask(s8, vB);                   // 每个 lane 1 位
__m512i   q  = _mm512_maskz_loadu_epi32(m, qty + i);          // 读取并清零卖单
acc = _mm512_add_epi64(acc, _mm512_cvtepi32_epi64(_mm512_castsi512_si256(q)));
acc = _mm512_add_epi64(acc, _mm512_cvtepi32_epi64(_mm512_extracti64x4_epi64(q, 1)));

比较直接产出 __mmask16，不需要扩展掩码宽度。

原始的 if 写法在 -O3 -mavx2 下通常会被编译器自动转换（if-conversion），用 -fopt-info-vec-optimized 确认。

例 3：字节查找

在 FIX 报文里找分隔符 SOH（0x01）：

#include <immintrin.h>
#include <bit>
#include <cstddef>
#include <cstdint>

// 返回 [p, p+n) 中第一个 0x01 的下标，没有则返回 n
size_t find_soh(const char* p, size_t n) {
  const __m256i needle = _mm256_set1_epi8(0x01);
  size_t i = 0;
  for (; i + 32 <= n; i += 32) {
    __m256i  v = _mm256_loadu_si256((const __m256i*)(p + i));
    uint32_t m = _mm256_movemask_epi8(_mm256_cmpeq_epi8(v, needle));
    if (m) return i + std::countr_zero(m);
  }
  for (; i < n; ++i) if (p[i] == 0x01) return i;
  return n;
}
// g++ -O3 -std=c++20 -mavx2

步骤	指令	作用
比较	`vpcmpeqb`	32 个字节各自与 `0x01` 比较
压缩	`vpmovmskb`	每个字节的最高位拼成 32 位整数
定位	`tzcnt`	最低的 1 所在位就是第一个命中的下标

每 32 字节只有 1 次比较和 1 个分支。同一模式的应用：

glibc 的 memchr、strlen
simdjson 定位引号和结构字符
abseil flat_hash_map（SwissTable）一次比对一组 16 个控制字节

例 4：过滤并压缩

列式扫描：把满足条件的行号收集起来。

#include <immintrin.h>
#include <bit>
#include <cstddef>
#include <cstdint>

// 把 price > th 的下标写入 out，返回个数；out 的容量至少为 n + 16
size_t filter_gt(const float* price, size_t n, float th, uint32_t* out) {
  const __m512 vth = _mm512_set1_ps(th);
  __m512i idx = _mm512_setr_epi32(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15);
  size_t k = 0, i = 0;
  for (; i + 16 <= n; i += 16) {
    __mmask16 m = _mm512_cmp_ps_mask(_mm512_loadu_ps(price + i), vth, _CMP_GT_OQ);
    _mm512_storeu_si512(out + k, _mm512_maskz_compress_epi32(m, idx));
    k  += std::popcount((unsigned)m);
    idx = _mm512_add_epi32(idx, _mm512_set1_epi32(16));
  }
  for (; i < n; ++i) if (price[i] > th) out[k++] = (uint32_t)i;
  return k;
}
// g++ -O3 -std=c++20 -mavx512f

vpcompressd 把掩码为 1 的 lane 向低位挤紧
每次整条写出 16 个元素，只有前 popcount(m) 个有效，下一次写入从 out + k 覆盖无效部分
先压缩到寄存器再写出：vpcompressd 以内存为目的操作数的形式在 AMD Zen 4 上由微码实现，速度很慢（实现相关）

AVX2 没有 compress 指令，做法是用 movemask 的结果查一张 256 项的置换表，再用 vpermd 重排。

其他能用的场景

场景	手法
小查找表（hex、base64、半字节 popcount）	`pshufb` 当作 16 项查表
排序	排序网络、双调归并、AVX-512 快排（x86-simd-sort）
校验与加密	CRC32 用 `pclmulqdq`，AES-NI，SHA-NI
图像、音频、3D 数学	像素混合、FIR 滤波、4×4 矩阵乘
数学函数	glibc libmvec 提供 `sin`、`cos`、`exp`、`log`、`pow` 的向量版，需要 `-ffast-math`

五、用不了或不划算的场景

判断流程

flowchart TD
    A[热点循环] --> B{数据连续?}
    B -- 否 --> B1[改数据布局: AoS 转 SoA]
    B -- 是 --> C{本轮依赖上一轮的结果?}
    C -- 是 --> C1{是可结合的归约?}
    C1 -- 是 --> C2[多累加器向量化]
    C1 -- 否 --> C3[换方向: 多条独立序列并排算]
    C -- 否 --> D{下一个地址依赖本次读到的值?}
    D -- 是 --> D1[不适合 SIMD, 优化 cache 局部性]
    D -- 否 --> E{分支能转成掩码?}
    E -- 否 --> E1[拆成两趟: SIMD 找位置, 标量处理]
    E -- 是 --> F{数据在 cache 内?}
    F -- 否 --> F1[先做 tiling 或合并多趟扫描]
    F -- 是 --> G[向量化]

例 1：循环携带依赖

单条序列的指数移动平均（EMA）：

for (size_t i = 1; i < n; ++i)
  ema[i] = alpha * x[i] + (1 - alpha) * ema[i-1];

ema[i] 依赖 ema[i-1]。关键路径上每个元素至少一条 FMA 的延迟（4 周期），n 个元素至少 4n 周期，与向量宽度无关。

SIMD 有两个使用方向：

方向	含义	对 EMA
沿序列	一个向量装同一条序列的相邻 8 个元素	不可行，相邻元素有依赖
跨序列	一个向量装 8 条序列在同一时刻的元素	可行，序列之间独立

跨序列的写法，8 个品种各占一个 lane：

#include <immintrin.h>
#include <cstddef>

// x 按时间步存放：x[t*8 + s]，s 为品种编号；out 长度为 8
void ema8(const float* x, size_t T, float alpha, float* out) {
  const __m256 va = _mm256_set1_ps(alpha);
  const __m256 vb = _mm256_set1_ps(1 - alpha);
  __m256 e = _mm256_loadu_ps(x);
  for (size_t t = 1; t < T; ++t) {
    __m256 ax = _mm256_mul_ps(va, _mm256_loadu_ps(x + t*8)); // 不在依赖链上
    e = _mm256_fmadd_ps(vb, e, ax);                          // 依赖链只有这一条 FMA
  }
  _mm256_storeu_ps(out, e);
}
// g++ -O3 -mavx2 -mfma

每个时间步仍是 4 周期，但一次产出 8 个品种的结果。alpha * x[t] 不依赖 e，放在依赖链之外提前算好。

例 2：指针追逐

for (Node* p = head; p; p = p->next) sum += p->val;

下一个节点的地址要等本次 load 完成才知道。链表、std::map、跳表都属于这一类。瓶颈是每次 cache miss 的延迟，计算量本身很小。gather 指令需要事先知道全部地址，用不上。

例 3：间接写入有冲突

for (size_t i = 0; i < n; ++i) ++hist[key[i]];

同一个向量里有两个相同的 key 时：

步骤	lane 0（key=5）	lane 1（key=5）
gather	读到 `hist[5]` = 3	读到 `hist[5]` = 3
加 1	4	4
scatter	写 `hist[5]` = 4	写 `hist[5]` = 4

正确结果是 5，得到 4。AVX-512CD 的 vpconflictd 能检测同一向量内的重复下标，检测加上分情况处理的开销通常使它不比标量快。

哈希表批量探测 table[hash(key[i])] 同理：地址随机，时间花在 cache miss 上。

例 4：变长编码

varint、UTF-8、protobuf 中，第 i+1 个元素的起点取决于第 i 个元素的长度，元素位置本身构成依赖链。

用 SIMD 需要改算法或改格式：

方案	做法
两趟处理（simdjson）	第一趟用 SIMD 找出所有结构字符的位置，第二趟按位置逐个解析
长度与数据分离（Stream VByte）	所有元素的长度字段集中存放，读 1 个字节即可知道 4 个整数的长度，再用 `pshufb` 一次解出

例 5：指针别名

void add(float* a, const float* b, size_t n) {
  for (size_t i = 0; i < n; ++i) a[i] += b[i];
}

调用方传入 add(p + 1, p, n) 时，a[i] 的写入会改变后面要读的 b[i+1]，向量版和标量版结果不同。编译器无法排除这种调用，只能在循环前插入重叠检查，并生成向量和标量两份循环（loop versioning）。

void add(float* __restrict a, const float* __restrict b, size_t n);

__restrict 向编译器保证两者不重叠，检查和标量副本都可以去掉。

其他不适用的场景

场景	原因	处理办法
循环里有非内联调用、虚函数、锁、I/O	没有向量版本，且有副作用	把调用移出循环；数学函数用 libmvec
整数除法 `a[i] / b[i]`	x86 没有向量整数除法指令	除数为常量时改乘法加移位；或转成 double
64 位整数乘法	AVX2 没有，`vpmullq` 需要 AVX-512DQ	用 32 位乘法拼
AoS 布局	同一字段在内存里不连续	改成 SoA，见第六节
元素很少（如 5 档盘口）	装载和尾部处理的开销大于收益	保持标量
大数组上的简单扫描	已受内存带宽限制	tiling，或把多趟扫描合并成一趟
提前退出且边界未知	多读的部分可能跨进未映射的页	按向量宽度对齐读取；GCC 14 起自动向量化支持部分形式（实现相关）
严格浮点语义	归约顺序、NaN 传播不能改变	对单个循环放宽，或手写

六、数据布局

AoS 与 SoA

// AoS：Array of Structures
struct Order {
  double  price;     // 8
  int32_t qty;       // 4
  uint8_t side;      // 1
  // 3 字节填充
  uint64_t id;       // 8
};                   // sizeof = 24
std::vector<Order> orders;

AoS 内存：| price qty side pad id | price qty side pad id | price ...
           ↑                       ↑
           相邻两个 price 相距 24 字节，读 4 个 price 要用 gather

// SoA：Structure of Arrays
struct Orders {
  std::vector<double>   price;
  std::vector<int32_t>  qty;
  std::vector<uint8_t>  side;
  std::vector<uint64_t> id;
};

SoA 内存：price: | p0 p1 p2 p3 p4 p5 ... |   一条 vmovupd 读 4 个
          qty:   | q0 q1 q2 q3 q4 q5 ... |

	AoS	SoA
只扫一个字段	每个 cache line（64 字节）含不到 3 个 price	每个 cache line 含 8 个 price
向量化	需要 gather	连续 load
访问单个对象的全部字段	1 次 cache miss	每个字段 1 次
增删元素	操作 1 个数组	操作多个数组，要保持同步

按列扫描的场景用 SoA，按对象随机访问的场景用 AoS。

AoSoA

struct OrderBlock {            // 每块 8 个订单
  double   price[8];
  int32_t  qty[8];
  uint8_t  side[8];
  uint64_t id[8];
};
std::vector<OrderBlock> blocks;

块内是 SoA，可以向量化；块与块之间是 AoS，一个订单的全部字段落在相邻的几个 cache line 内。块大小取向量宽度的整数倍。

七、实现细节

对齐

指令	对地址的要求	不满足时
`_mm256_load_ps` / `_mm512_load_ps`	32 / 64 字节对齐	触发 #GP，进程收到 SIGSEGV
`_mm256_loadu_ps` / `_mm512_loadu_ps`	无	正常执行

数据已对齐时，loadu 与 load 速度相同（Nehalem 之后的 Intel CPU，实现相关）。开销来自跨越边界的访问：

情况	代价
访问落在一个 cache line 内	无额外代价
跨 cache line	需要访问两个 line，延迟增加
跨 4 KB 页	需要两次地址翻译，代价更高

一个 zmm 寄存器恰好 64 字节。地址不是 64 的倍数时，每一次 512 位访问都跨 cache line。AVX-512 代码的数据按 64 字节对齐：

alignas(64) float buf[1024];
float* p = static_cast<float*>(std::aligned_alloc(64, n * sizeof(float)));  // n*sizeof(float) 需为 64 的倍数

尾部处理

n 不是向量宽度的整数倍时，剩余元素有三种处理方式：

方式	做法	限制
标量尾循环	剩余元素逐个处理	无；n 很小时尾循环占比高
掩码读写	`_mm512_maskz_loadu_ps(m, p)`，被屏蔽的 lane 不访问内存，不会触发缺页异常	AVX2 的 `vmaskmov` 较慢（实现相关）
重叠最后一个向量	对 `[n-8, n)` 再处理一次	只适用于幂等的操作（逐元素映射到独立输出、max、min），不适用于求和

归约要用多个累加器

单个累加器的向量归约仍是一条依赖链：每轮迭代要等上一轮的加法完成。Skylake 上 FMA 延迟 4 周期、每周期可执行 2 条，单累加器只用到了 1/8 的执行能力。

#include <immintrin.h>
#include <cstddef>

double dot_avx2(const double* p, const double* q, size_t n) {
  __m256d s0 = _mm256_setzero_pd(), s1 = s0, s2 = s0, s3 = s0;
  size_t i = 0;
  for (; i + 16 <= n; i += 16) {
    s0 = _mm256_fmadd_pd(_mm256_loadu_pd(p + i),      _mm256_loadu_pd(q + i),      s0);
    s1 = _mm256_fmadd_pd(_mm256_loadu_pd(p + i + 4),  _mm256_loadu_pd(q + i + 4),  s1);
    s2 = _mm256_fmadd_pd(_mm256_loadu_pd(p + i + 8),  _mm256_loadu_pd(q + i + 8),  s2);
    s3 = _mm256_fmadd_pd(_mm256_loadu_pd(p + i + 12), _mm256_loadu_pd(q + i + 12), s3);
  }
  __m256d s  = _mm256_add_pd(_mm256_add_pd(s0, s1), _mm256_add_pd(s2, s3));
  __m128d lo = _mm_add_pd(_mm256_castpd256_pd128(s), _mm256_extractf128_pd(s, 1));
  double  r  = _mm_cvtsd_f64(_mm_add_sd(lo, _mm_unpackhi_pd(lo, lo)));
  for (; i < n; ++i) r += p[i] * q[i];
  return r;
}
// g++ -O3 -mavx2 -mfma

累加器个数	每 4 周期完成的 FMA 数	端口利用率
1	1	1/8
4（上面的代码）	4	1/2
8	8	1

打满两个 FMA 端口需要 8 个累加器。这与 kernel_4x8 取 8 个累加器是同一个依据。

gather 与 scatter

gather 和 scatter 的定义见第一节。gather 省掉的是下标提取和结果拼装的指令，不减少内存访问次数：8 个元素仍是 8 次 load。

数据在 cache 内时，gather 相对标量读取的收益很小
数据不在 cache 内时，时间花在 cache miss 上，gather 无收益
受 Downfall（CVE-2022-40982）影响的 Intel CPU（第 6 到第 11 代酷睿）打上微码补丁后，gather 明显变慢（实现相关）

能把数据排成连续的，就不用 gather。

non-temporal store

普通 store 写一个不在 cache 里的 line 时，CPU 先把这个 line 从内存读进来（Read For Ownership），再修改。整条 line 都要被覆盖时，这次读取是浪费。

#include <immintrin.h>
#include <cstddef>

// dst 必须 32 字节对齐，n 为 8 的倍数
void fill_nt(float* dst, float v, size_t n) {
  const __m256 x = _mm256_set1_ps(v);
  for (size_t i = 0; i < n; i += 8)
    _mm256_stream_ps(dst + i, x);   // 不经过 cache，不触发 RFO
  _mm_sfence();                     // 保证后续 store 的顺序
}
// g++ -O3 -mavx

适用	不适用
写入量超过 cache 容量，写完后短期内不再读	写完马上要读：数据不在 cache 里，读取会 miss

denormal

绝对值小于 FLT_MIN 的非零浮点数是 denormal。运算遇到 denormal 时 CPU 转入微码处理，单次运算的耗时上升到上百周期（实现相关）。IIR 滤波、EMA 这类衰减到接近 0 的序列容易触发。

#include <xmmintrin.h>
#include <pmmintrin.h>
_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);           // 结果为 denormal 时置 0
_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);   // 输入为 denormal 时当作 0

这两个标志位在 MXCSR 寄存器里，每个线程各有一份，新线程要单独设置。

浮点结果一致性

向量化的归约改变了加法顺序，结果的末几位与标量版不同；AVX2 路径（4 个 lane）和 AVX-512 路径（8 个 lane）的结果也互不相同。

场景	影响
回测与实盘跑在不同 CPU 上	同一份数据算出的信号末位不同，阈值判断可能翻转
单元测试用 `==` 比较浮点	换编译选项或换机器后失败
多副本状态机要求逐位一致	各副本必须走同一条代码路径

要求逐位一致时，固定向量宽度和累加器个数，不使用运行时分发；或者用定点整数。

八、三种写法

写法	优点	缺点
自动向量化	代码不变，可移植	是否成功取决于编译器版本和选项，容易在重构后悄悄失效
intrinsics	指令选择完全可控	绑定指令集，每个宽度写一份
可移植库（`std::simd`、Highway、xsimd）	一份代码对应多个指令集	特殊指令（compress、`pshufb` 查表）的覆盖不完整

帮助自动向量化

手段	作用
`__restrict`	排除别名
`#pragma omp simd`（`-fopenmp-simd`）	声明迭代之间无依赖
`#pragma omp simd reduction(+:s)`	允许重排这一个归约
循环计数用 `size_t` 或 `ptrdiff_t`	避免 32 位无符号下标回绕导致分析失败
循环体内不调用非内联函数	调用会阻止向量化
`-fopt-info-vec-missed`	打印失败原因

可移植库

libstdc++ 自 GCC 11 起提供 <experimental/simd>（Parallelism TS 2）；C++26 把 SIMD 类型纳入标准库，各实现的支持进度不同（实现相关）。

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

float sum(const float* p, size_t n) {
  using V = stdx::native_simd<float>;      // 宽度由 -march 决定
  V acc = 0.0f;
  size_t i = 0;
  for (; i + V::size() <= n; i += V::size())
    acc += V(p + i, stdx::element_aligned);
  float s = stdx::reduce(acc);
  for (; i < n; ++i) s += p[i];
  return s;
}
// g++ -O3 -std=c++17 -march=x86-64-v3

九、部署陷阱

运行时分发

用 -mavx2 编译的二进制在不支持 AVX2 的 CPU 上执行到 AVX2 指令时收到 SIGILL。要让同一个二进制适配多种 CPU，按函数指定指令集，运行时选择：

#include <cstddef>

__attribute__((target("avx2")))
size_t find_soh_avx2(const char* p, size_t n);     // 函数体内可用 AVX2 intrinsics
size_t find_soh_scalar(const char* p, size_t n);

size_t find_soh(const char* p, size_t n) {
  static const bool has_avx2 = __builtin_cpu_supports("avx2");
  return has_avx2 ? find_soh_avx2(p, n) : find_soh_scalar(p, n);
}
// g++ -O3（不加 -mavx2）

GCC 的 __attribute__((target_clones("avx2", "default"))) 自动生成多个版本，并在加载时通过 ifunc 选择。

内联函数与不同的编译选项

a.cpp（-mavx2 编译）   ──┐
                         ├── 都包含头文件里的 inline float norm(const Vec&)
b.cpp（不加 -mavx2）   ──┘

两个目标文件里各有一份 norm，链接器只保留一份。保留的是 AVX2 版本时，b.cpp 里本应只在老 CPU 上走的标量路径也会调用到它，触发 SIGILL。

处理办法：用到特定指令集的代码放在独立的 .cpp 里，不包含会被其他翻译单元共用的内联函数；或者给每个指令集版本套一层不同的命名空间。

SSE 与 AVX 混用

ymm 寄存器的高 128 位非零时执行旧式 SSE 编码的指令，会产生额外开销：Haswell 及更早的 CPU 是状态切换（约 70 周期），Skylake 是每条 SSE 指令附带虚假依赖（实现相关）。

来源	处理
自己的代码	整个翻译单元用 `-mavx` 以上编译，标量浮点指令也会用 VEX 编码
调用未用 AVX 编译的库	调用前执行 `vzeroupper`；GCC/Clang 在用到 ymm 的函数返回前自动插入

AVX-512 降频

CPU	行为（实现相关）
Intel Skylake-X / Cascade Lake	执行重型 512 位指令时核心降频，恢复需要毫秒级时间
Intel Ice Lake 及之后	降幅很小
AMD Zen 4	512 位指令拆成两个 256 位执行，不降频
AMD Zen 5	桌面和服务器型号有完整的 512 位数据通路

在降频明显的 CPU 上，一小段 AVX-512 代码会拖慢它前后的全部标量代码。延迟敏感的路径上用 -mprefer-vector-width=256 把自动向量化限制在 256 位。

Amdahl 定律

整体加速比 = 1 / ((1 - p) + p / s)      p：被加速部分占总时间的比例，s：该部分的加速比

p	s	整体加速比
60%	8	2.1
90%	8	4.7
60%	16	2.3

向量化之前先用 perf record 确认热点占比。p 为 60% 时，把宽度从 8 加到 16 只把整体加速比从 2.1 提到 2.3。

十、验证方法

g++ -O3 -march=x86-64-v3 -fopt-info-vec-optimized -fopt-info-vec-missed a.cpp

打印每个循环是否向量化，以及失败的原因。

objdump -d --no-show-raw-insn a.out | grep -c ymm

统计用到 ymm 寄存器的指令条数，确认生成的是 256 位代码。

perf stat -e cycles,instructions,branch-misses,L1-dcache-load-misses,LLC-load-misses ./a.out

指标	向量化生效时的变化	tiling 生效时的变化
`instructions`	明显下降	基本不变
`branch-misses`	下降（分支转掩码时）	基本不变
`LLC-load-misses`	基本不变	明显下降
IPC（instructions / cycles）	受内存限制时很低	上升

llvm-mca -mcpu=skylake loop.s

静态分析一段汇编的端口压力和依赖链瓶颈，用于判断累加器个数是否足够。

十一、面试速答卡

SIMD 为什么快 一条指令处理 8 或 16 个元素，每个元素摊到的指令数和分支数下降；数据相关的分支转成掩码后不再有预测失败。

SIMD 的上限在哪 内存带宽。数据不在 cache 里时，向量再宽也要等数据到达。

tiling 为什么能加速 计算量不变，改变遍历顺序，让一小块数据在 cache 里被反复使用后再换下一块。矩阵乘法 N=2048、T=128 时，内存流量从约 65 GB 降到约 1.5 GB。

register tiling 是什么 把输出的一个小块放在向量寄存器里跨循环累加，省掉每次迭代的 load 和 store。累加器个数 = FMA 延迟 × 每周期可执行条数，Skylake 上是 4 × 2 = 8。

什么样的循环能向量化 数据连续、迭代独立或只有可结合的归约、分支能转掩码、指针不重叠、运算有向量指令。

什么样的循环不能 三类：依赖链（EMA、前缀和、指针追逐）；地址不可预知或有冲突（直方图、哈希探测、变长编码）；语义限制（别名、浮点顺序、副作用）。

不能向量化时怎么改 依赖链：换方向，多条独立序列并排算。AoS：改 SoA。变长编码：拆成两趟，或把长度字段与数据分离。别名：加 __restrict。

浮点求和为什么不自动向量化 浮点加法不满足结合律，重排会改变结果。用 -ffast-math、#pragma omp simd reduction 或手写多累加器。

load 和 loadu 的区别 load 要求对齐，否则触发异常；loadu 无要求。数据已对齐时两者速度相同，开销来自跨 cache line 和跨页。

AVX-512 一定比 AVX2 快吗 不一定。受内存带宽限制的循环没有收益；Skylake-X 上会降频；向量化部分占比不高时受 Amdahl 定律限制。

向量化之后结果和标量不一样，是 bug 吗 归约的加法顺序变了，末几位不同是预期行为。要求逐位一致时固定代码路径，或用定点整数。