分支预测、restrict 与 SIMD

周三上午,你在国内某头部私募的衍生品自营桌做期权 Greeks 重估。沪深300 ETF(510300.SH)期权链上挂着 480 张合约,每张合约要算一次 Delta、Gamma、Vega、Theta,每个 Greek 都是一次 100 万路径的 Monte Carlo——480 × 4 × 1M = 19.2 亿条路径。上一课你已经用 monotonic_buffer_resource 把 simulate_one_path_arena 的每路径分配压到零,绝对时间从盘后 17:00 算到 17:42。组长今天追问:能不能压到 17:10 之前?——因为风控会在 17:15 截稿,而下午两点行情大动后他想再跑一遍敏感度。算法没动、布局已优、分配器已换,剩下的钱藏在哪里?

藏在三层机器级低效里:​​分支预测​​(branch prediction)误预测、​​指针别名​​(pointer aliasing)阻断自动向量化、​​标量内循环​​没用上 256 位 YMM 寄存器。本课依次拆开:先把可预测的分支保留、不可预测的分支重塑;再用 __restrict__ 解开编译器的手脚;最后写一段 AVX2 内核,并把它接回上一课的 Monte Carlo 定价器做收官。

排序数组比未排序数组快:分支预测的真实代价

经典演示:1M 个 [0, 256) 范围的随机整数,数一下大于 128 的有几个。同样的代码、同样的数据、同样的结果——​​排序后比未排序快 3-5 倍​​。

std::vector<int> make_data(std::size_t n, std::uint32_t seed = 42) {
    std::mt19937 rng(seed);
    std::uniform_int_distribution<int> dist(0, 255);
    std::vector<int> v(n);
    for (std::size_t i = 0; i < n; ++i) v[i] = dist(rng);
    return v;
}

std::int64_t count_above_128(const std::vector<int>& v) {
    std::int64_t c = 0;
    for (int x : v) {
        if (x > 128) {        // unpredictable when v is unsorted; predictable when sorted
            ++c;
        }
    }
    return c;
}

static void BM_CountAbove_Unsorted(benchmark::State& state) {
    auto v = make_data(1'000'000);
    for (auto _ : state) {
        benchmark::DoNotOptimize(count_above_128(v));
    }
}
BENCHMARK(BM_CountAbove_Unsorted);

static void BM_CountAbove_Sorted(benchmark::State& state) {
    auto v = make_data(1'000'000);
    std::sort(v.begin(), v.end());
    for (auto _ : state) {
        benchmark::DoNotOptimize(count_above_128(v));
    }
}
BENCHMARK(BM_CountAbove_Sorted);

机制:Ice Lake 一代 CPU 的分支预测器在​​可预测​​模式下命中率 > 95%,在​​随机​​模式下不到 60%;每次误预测要刷掉 15-20 个流水线 cycle。排序后的数组,前半段全 ≤128、后半段全 > 128,中间只有一次跳变——预测器几轮迭代就锁定模式;未排序时它每次都在赌币,赌输一半。perf stat -e branch-misses 把这件事直接量出来:未排序版本的误预测次数是排序版本的 20-30 倍。

实操结论:​​热循环里不可预测的分支极贵​​。三条出路是 (a) 用无分支算术改写,(b) 通过排序 / 分区 / 分桶让分支变可预测,(c) 当确实可预测时把方向告诉编译器。

[[likely]] / [[unlikely]]:把已知的偏向告诉编译器

C++20 标准化了 [[likely]] / [[unlikely]] 属性;本模块目标是 C++17,在 GCC / Clang 上用 __builtin_expect 等价回退。同样需要可移植的 __restrict__ 包装,顺手给出。

// Portable hint macros (C++17 fallback; C++20 has [[likely]] / [[unlikely]])
#if defined(__GNUC__) || defined(__clang__)
  #define LIKELY(x)   __builtin_expect(!!(x), 1)
  #define UNLIKELY(x) __builtin_expect(!!(x), 0)
  #define RESTRICT    __restrict__
#else
  #define LIKELY(x)   (x)
  #define UNLIKELY(x) (x)
  #define RESTRICT
#endif

// Hot-path example: in a tick handler, the no-event-this-tick branch dominates.
void on_tick(const Tick& t) {
    if (UNLIKELY(t.is_trade_event())) {
        process_trade(t);
    } else {
        update_quote_state(t);
    }
}

// Aliasing-aware AXPY: developer promises x and y do not overlap.
void axpy(double* RESTRICT y, const double* RESTRICT x,
          double a, std::size_t n) {
    for (std::size_t i = 0; i < n; ++i) y[i] += a * x[i];
}

提示干两件事:一,引导编译器把​​可能成立​​的分支编译为顺延执行(fall-through, 不付跳转代价、对热代码缓存友好),不可能的分支编译为跳出到冷区的 forward jump;二,在 CPU 预测器还没学会一条新分支时(典型场景:刚冷启动、刚 perf 重置),提供静态退路。

​重要警告:错的提示比无提示更糟​​。提示一错,编译器把冷路径放到顺延位置、热路径变成跳转;CPU 静态预测器也跟着读反。把 LIKELY / UNLIKELY 留给剖析器已经证实比例 ≥ 95/5 的分支,其他地方一概不写。在 on_tick 上,真实 tick 流里"行情更新而非成交"占 99%+,所以 UNLIKELY(t.is_trade_event()) 是被业务事实加持的。

restrict 与自动向量化:解开编译器的手脚

接着看 axpy——BLAS-1 的 y = a*x + y,量化代码里到处是它。​​不​​写 __restrict__ 时,把它放进 godbolt.org(g++ -O3 -mavx2),编译器吐出的是​​带运行时重叠检查的标量循环​​:用 mulsd / addsd 一次处理一个 double。原因是 C++ 允许指针别名——axpy(p, p, 2.0, 100) 在语言上完全合法,此时写 y[i] 会改下一轮 x[i+1] 读到的值,形成一条循环携带依赖(loop-carried dependency),向量化的前提就此崩塌。编译器只能假设最坏,生成两份代码(标量保底 + 向量化重叠检查)或者干脆全程标量。

加上 __restrict__ 后,你向编译器​​承诺​ x 与 y 不重叠。编译器随即吐出向量化版本:把 a 用 _mm256_set1_pd 广播到一个 __m256d,每轮 load 4 个 x 的 double、_mm256_mul_pd 乘 a、_mm256_add_pd 加到当前 4 个 y 上、_mm256_storeu_pd 写回。内循环迭代次数从 n 变成 n/4,吞吐近 4 倍。__restrict__ 是 GCC / Clang 扩展(C99 标准里 restrict 是关键字,C++ 标准未引入),用 #define RESTRICT __restrict__ 包装就跨编译器可用,MSVC 上对应 __declspec(restrict)(详见 cppreference)。

合同是单向且严酷的:如果你撒谎、x 与 y 确实重叠,​​行为未定义​​——可能数值出错、可能段错误、可能在 release 下偶发崩溃。__restrict__ 是"我保证"型保证,不是"如果"型检查。

AVX2 内联函数:一周期处理 4 个 double

x86_64 从 Sandy Bridge(2011, AVX)到 Haswell(2013, AVX2)就有了 256 位 YMM 寄存器,一根能装 4 个 double 或 8 个 float,AVX2 指令一周期完成一根寄存器上的并行运算。<immintrin.h> 把这些指令暴露为返回 __m256d / __m256 / __m256i 的 C 函数。需要默记的就六个:

• _mm256_loadu_pd(ptr) —— 从 ptr 起载入连续 4 个 double 入 YMM;u 表示 unaligned,Sandy Bridge 之后基本零代价,优先用它。
• _mm256_storeu_pd(ptr, v) —— 反向写回。
• _mm256_add_pd(a, b) —— 通道(lane)逐对相加。
• _mm256_mul_pd(a, b) —— 通道逐对相乘。
• _mm256_setzero_pd() —— 4 通道全置 0.0。
• _mm256_set1_pd(a) —— 广播 scalar a 到 4 通道。

水平归约(horizontal reduction, 把 4 通道加成一个 scalar)的成语写法:取高 128 位、加到低 128 位、再水平加一次、提取标量。一次写好封进 hsum_pd,后续忘了细节也无妨。

工作例:tick 流的单遍均值方差核

量化策略每天要对实时 tick 算 streaming mean / streaming variance / Sharpe——给一段 double,一遍循环算 sum 与 sumsq,均值是 sum / n、方差是 sumsq / n − mean²。先写标量参考,再写手向量化版,两者数值要可比。

MeanVar mean_var_scalar(const double* x, std::size_t n) {
    double sum = 0.0;
    double sumsq = 0.0;
    for (std::size_t i = 0; i < n; ++i) {
        sum   += x[i];
        sumsq += x[i] * x[i];
    }
    double mean = sum / static_cast<double>(n);
    double var  = sumsq / static_cast<double>(n) - mean * mean;
    return {mean, var};
}

#include <immintrin.h>

struct MeanVar { double mean; double var; };

static inline double hsum_pd(__m256d v) {
    // sum the 4 lanes of v into a single double
    __m128d lo = _mm256_castpd256_pd128(v);
    __m128d hi = _mm256_extractf128_pd(v, 1);
    __m128d s  = _mm_add_pd(lo, hi);
    __m128d sh = _mm_shuffle_pd(s, s, 1);
    return _mm_cvtsd_f64(_mm_add_sd(s, sh));
}

MeanVar mean_var_avx2(const double* RESTRICT x, std::size_t n) {
    __m256d vsum   = _mm256_setzero_pd();
    __m256d vsumsq = _mm256_setzero_pd();
    std::size_t i = 0;
    for (; i + 4 <= n; i += 4) {
        __m256d v = _mm256_loadu_pd(x + i);
        vsum   = _mm256_add_pd(vsum,   v);
        vsumsq = _mm256_add_pd(vsumsq, _mm256_mul_pd(v, v));
    }
    double sum   = hsum_pd(vsum);
    double sumsq = hsum_pd(vsumsq);
    for (; i < n; ++i) {                 // scalar tail for n not divisible by 4
        sum   += x[i];
        sumsq += x[i] * x[i];
    }
    double mean = sum / static_cast<double>(n);
    double var  = sumsq / static_cast<double>(n) - mean * mean;
    return {mean, var};
}

1M 个 double(初始化 x[i] = 100.0 + 0.01 * (i % 200) - 1.0;,与上一课 SoA 工作例同源)上跑 Google Benchmark:标量约 X ms、-O3 -mavx2 自动向量化约 X/2 ms、手写 AVX2 约 X/3 ms。手写版赢自动向量化的两点常出在尾部处理(n 不整除 4 时的 scalar tail 调度)与双累加器(vsum / vsumsq 间无依赖,流水线更深)。

数值层面要点两条:​​(1)​ scalar 与手写 AVX2 的 mean 应当​​按字节相同​​——加法顺序差别只发生在 SIMD 累加器之间,而 4 通道求和顺序与 scalar 不同,所以 var 可能差 1-2 ULP;​​(2)​ 自动向量化版本因为编译器自由重排(浮点加法不满足结合律),mean 与 var 都可能差 1-2 ULP。这不是 bug——是浮点数学的事实。要 bit-exact 复现的生产代码用 pairwise summation 或 Kahan summation,代价标量 2x、但同样可向量化,具体留给数值方法模块(2.2.x)。

编译与剖析

把 axpy 拷进 godbolt 看汇编;把 capstone 编译并以 taskset 锁核跑 perf stat,数据才稳定。

# Inspect axpy emit on godbolt.org or locally:
g++ -std=c++17 -O3 -mavx2 -S -masm=intel axpy.cpp -o axpy.s
less axpy.s   # look for `vfmadd` / `vmulpd` / `vaddpd` lines (vectorised)
              # vs scalar `mulsd` / `addsd` (not vectorised)

# Build the capstone Monte Carlo binary:
g++ -std=c++17 -O3 -mavx2 -march=haswell -g mc_capstone.cpp -lbenchmark -lpthread -o mc_capstone

# Run it pinned to one core so the measurement is stable:
taskset -c 2 perf stat -e cycles,instructions,cache-misses,branch-misses ./mc_capstone

-march=haswell 是覆盖最广的安全选择——AVX2 全面支持,且兼容到 2013 年以后的 Intel 与 2017 年以后的 AMD。已知部署机型是 Ice Lake-SP 时可用 -march=icelake-server,但本课不依赖。

AVX-512(__m512d, _mm512_*)在 Skylake-SP 时代曾因频率回退(clock throttling)劝退延迟敏感场景,Ice Lake-SP 之后基本修好,但国内 quant 桌出于跨机型稳定与 AMD EPYC 兼容(Zen 4 / 2022 才广泛支持 AVX-512),仍以 AVX2 为主流。Apple Silicon 与 AWS Graviton 用 ARM NEON <arm_neon.h>(float64x2_t 是 128 位、两通道),迁移机械但不在本课范围。FMA3 指令 _mm256_fmadd_pd 一周期内做一次"乘加 + 一次舍入",对 AXPY 是天然下一步,要求 -mfma;本课内核为与 scalar 字节一致不用它。

收官:把 L1 → L2 → L3 → L4 全链应用在 Monte Carlo 定价器上

回到开头的 480 合约 × 1M 路径问题。L3 的 simulate_one_path_arena 把每条路径的临时数组从堆搬上栈;L2 的 SoA 把 1M 条路径的终值压成一根连续 std::pmr::vector<double>;L4 的 mean_var_avx2 把折现回报的求和从标量变成 AVX2 内核(mean_var_avx2 一次过求 sum 与 sumsq,后者顺手给出 Monte Carlo 标准误差的样本方差)。

具体改造:外层 for (i = 0; i < N_PATHS; ++i) 不再边算边累加 std::max(ST − K, 0.0),而是把每条路径的终值存进一个 std::pmr::vector<double> path_ends(N_PATHS, &arena);,出完所有路径再调一次 mean_var_avx2,得到 mean_payoff 与 var_payoff,折现:price = exp(-r * T) * mean_payoff,标准误差 se = exp(-r * T) * sqrt(var_payoff / N_PATHS)。

沪深300 ETF(510300.SH)European call 参数 S = 4.20, K = 4.30, r = 0.024, sigma = 0.18, T = 0.5,N_paths = 1_000_000, N_steps = 252,seed 42(沿用上一课以便对账)。在 Intel Xeon Gold 6342 / 6354(Ice Lake-SP)单核锁定下,L4 capstone 相对 3.4.2-时代基线应该 5-10x 提速,价格落在基线的 ±1 标准误差(约 0.4 / sqrt(N_paths) × S × σ × sqrt(T) ≈ 0.4 × 1e-3 × 4.20 × 0.18 × 0.707 ≈ 2e-4 量级)内。perf stat 上对照基线,期望看到 branch-misses 降低、instructions per cycle 上升、cache-misses 因 SoA 终值数组的顺序访问保持低位。

正确性是承重项:​​优化改的是性能,不是数值​​。如果价格偏出标准误差,先回退最后一步、再二分定位——通常是手写 SIMD 的尾部循环边界写歪,或者 SoA 终值数组的填入索引错位。

操作规则与下一课

本课的层级:​​先 perf stat 量出可预测度,再 godbolt 看 emit、用 __restrict__ 解开自动向量化,实在不够再手写 intrinsic​​。CFFEX 股指期货策略的 tick-to-trade 关键路径常落在 5 μs 以内,信号 + 风控 + 决策合起来 < 1 μs——这预算下哪怕一段 SIMD 内核省 200 ns 都值得。但本课 capstone 对应的不是 tick-to-trade,是​​非实时、CPU 密集的批处理​​:Greeks 重估、VaR 重估、波动率曲面拟合——延迟容忍、吞吐至上,吞吐翻倍意味着研究迭代翻倍。

到这里,本模块在​​单核单线程​​层面收官:L1 教你看见瓶颈、L2 教你把数据摆好、L3 教你把分配器换掉、L4 教你把指令换掉。3.4.4 并发与网络要把同一条 capstone 在多核上再放大一个量级——线程池、std::execution::par_unseq、SPSC 队列、内核旁路网络(DPDK / io_uring),都是 3.4.4 的事;3.4.5 交易系统则把这条优化过的内核嵌入完整策略框架。先把单核吞吐天花板钉死,再谈并行。

练习