3. 状态向量模拟器的核心算法

状态向量模拟器把 n 个量子比特的量子态表示为一个长度为 2ⁿ 的复数数组 α，其中 α[i] 对应计算基态 |i⟩ 的振幅。单比特门作用在比特 k 上时，只把下标中第 k 位为 0 和 1 的振幅两两配对，做一次 2×2 复数矩阵乘法。这就是“全振幅更新”的原始形态：一次门操作要遍历整个 2ⁿ 数组，但每次只触碰两个元素。理解这一点的关键，是看清“配对规则”和“内存步长”之间的关系。

设比特编号从 0 开始，且约定比特 k 对应下标二进制表示中的第 k 位（最低位为比特 0）。那么对任意下标 i，若其第 k 位为 0，则与它配对的伙伴下标是 i | (1<<k)，两者相差 stride = 2^k。当 k = 0 时，配对元素相邻，步长为 1；当 k = n−1 时，配对元素相隔 2ⁿ⁻¹，跨越半个数组。这个 stride 决定了访存是否连续，也决定了缓存命中率的高低。

直接按上述公式实现，最朴素的循环是外层遍历 i 从 0 到 2ⁿ−1，内层判断第 k 位。这样会做 2ⁿ 次判断，且当 stride 很大时，α[i] 与 α[j] 落在不同的缓存行甚至不同的内存页，每次访问都可能触发缓存缺失。更常见的做法是三层循环：外层按块（block）划分，中层遍历块内偏移，内层遍历 stride，使同一块内的访问尽量连续。下面给出一个以“块—偏移—步长”组织的单比特门更新骨架，它把 stride 作为最内层循环的步进单位，从而让相邻迭代访问相邻地址：

// 单比特门作用在比特 k 上，state 为长度 2^n 的复数数组
void apply_1q(complex* state, int n, int k, complex U[2][2]) {
    int stride = 1 << k;          // 配对元素的距离
    int block  = stride << 1;     // 一个配对块的长度
    for (int base = 0; base < (1 << n); base += block) {
        for (int off = 0; off < stride; ++off) {
            int i = base + off;       // 第 k 位为 0
            int j = i + stride;       // 第 k 位为 1
            complex a = state[i], b = state[j];
            state[i] = U[0][0]*a + U[0][1]*b;
            state[j] = U[1][0]*a + U[1][1]*b;
        }
    }
}

这段代码的访存模式是：内层 off 递增时，i 和 j 都连续递增，因此对 state[i] 和 state[j] 的访问都是顺序的，缓存行利用率高。当 k 较小时，block 很小，外层 base 频繁跳转，但每次跳转的跨度也小；当 k 较大时，block 很大，内层循环很长，顺序访问的优势更明显。真正的问题出现在多比特门和门序列上：如果每个门都独立遍历整个数组，那么 n 个比特上的门序列会产生 O(m·2ⁿ) 的访存总量，其中 m 是门数。门融合（gate fusion）正是为了减少这个总量。

门融合的基本思想是：把作用在相同或相邻比特上、且彼此之间没有其他比特门干扰的门，合并成一个更大的等效矩阵，再一次性作用到状态向量上。例如两个单比特门 U 和 V 依次作用在同一个比特 k 上，可以合并为 W = V·U，只做一次遍历。更一般地，若一组门只涉及 t 个比特，且这 t 个比特在门序列中连续出现，就可以把这组门合并成一个 2^t×2^t 的矩阵，然后按 t 比特门的方式更新。融合的收益来自两方面：一是减少遍历次数，二是把多次小步长访问合并成一次大步长访问，改善访存局部性。但融合也有代价：2^t×2^t 矩阵的存储和计算量随 t 指数增长，因此 t 通常限制在 4 到 8 之间。

下面这张图把“全振幅更新”的配对关系、stride 随比特位置的变化，以及门融合前后的遍历次数对比放在同一张图里，便于对照理解。

从图中可以看到，stride 越大，配对元素在内存中越分散，缓存行利用率越低。实际实现中常用的优化手段有三类。第一类是循环重排：把 stride 放到最内层，让相邻迭代访问相邻地址，正如前面代码所示。第二类是分块（tiling）：把状态向量切成若干块，使每块能放入 L1 或 L2 缓存，然后在块内完成尽可能多的门操作，减少对主存的往返。第三类是门融合：把多个门合并成一次遍历，从根本上减少访存次数。三者往往组合使用：先做门融合降低遍历次数，再用分块和循环重排改善每次遍历的局部性。

门融合的实现需要解决两个问题：哪些门可以融合，以及融合后的矩阵如何作用到状态向量。对于第一个问题，常用的判据是：若两个门的作用比特集合相同，或一个门的作用比特集合包含在另一个门的作用比特集合内，且它们之间没有其他门作用于这些比特，则可以融合。例如 CNOT(q0,q1) 和 X(q0) 可以融合，因为 X(q0) 只作用于 q0，而 CNOT 作用于 q0 和 q1，q0 在 CNOT 的作用集合内。融合后的矩阵维度是 2^t，其中 t 是涉及比特的总数。对于第二个问题，融合后的 t 比特门更新可以看作单比特门更新的推广：把 t 个比特对应的下标位提取出来，形成一个 2^t 维的子空间，然后在这个子空间上做矩阵乘法。实现时通常用位运算把下标映射到子空间索引，再用查表或直接计算完成更新。

下面给出一个融合后 t 比特门更新的简化实现思路，假设融合涉及比特集合 bits = {b0, b1, ..., b(t-1)}，且这些比特在下标中的位置已知。核心是把下标 i 中这些比特的值提取出来作为子空间索引，其余比特保持不变：

// 融合后 t 比特门更新（简化示意）
// bits: 涉及的比特位置数组，t: 比特数，M: 2^t × 2^t 矩阵
void apply_tq(complex* state, int n, int* bits, int t, complex* M) {
    int dim = 1 << t;
    int total = 1 << n;
    // 预计算每个子空间索引对应的基下标（其余比特为 0）
    int* base = new int[dim];
    for (int s = 0; s < dim; ++s) {
        int idx = 0;
        for (int b = 0; b < t; ++b)
            if (s & (1 << b)) idx |= (1 << bits[b]);
        base[s] = idx;
    }
    // 遍历所有“其余比特”的组合
    int mask = 0;
    for (int b = 0; b < t; ++b) mask |= (1 << bits[b]);
    for (int rest = 0; rest < total; ++rest) {
        if (rest & mask) continue;          // 只处理其余比特为 0 的基
        complex tmp[256];                    // t ≤ 8 时 dim ≤ 256
        for (int s = 0; s < dim; ++s)
            tmp[s] = state[rest | base[s]];
        for (int s = 0; s < dim; ++s) {
            complex acc = 0;
            for (int r = 0; r < dim; ++r)
                acc += M[s*dim + r] * tmp[r];
            state[rest | base[s]] = acc;
        }
    }
    delete[] base;
}

这段代码的访存模式是：外层 rest 遍历所有“其余比特为 0”的基，内层 s 遍历子空间。由于 base[s] 是预先算好的，内层对 state 的访问是“按 base 数组跳转”的，跳转模式取决于 bits 的分布。如果 bits 是连续的（例如 q0、q1、q2），base[s] 就是连续的 0 到 7，访问连续；如果 bits 分散（例如 q0、q3、q5），base[s] 的跳转就会变得不规则，缓存效率下降。因此，门融合时不仅要考虑融合哪些门，还要考虑融合后涉及的比特是否连续。实践中，调度器会优先融合作用在连续比特上的门，或者通过比特重映射（qubit reordering）把频繁交互的比特放到相邻位置，从而改善访存局部性。

综合来看，状态向量模拟器的核心算法可以概括为一条优化链：全振幅更新定义了“每次门操作必须遍历整个状态向量”这一基本约束；门融合通过合并门来减少遍历次数；内存访问模式优化（循环重排、分块、比特重映射）则让每次遍历的访存尽可能连续。三者共同决定了模拟器在 n 较大时的实际性能。理解 stride 与缓存行的关系，是判断一个实现是否高效的关键切入点：当 k 接近 n−1 时，单比特门更新会退化为跨半个数组的随机访问，此时必须依靠分块或重映射来补救；而当融合后的比特集合连续时，t 比特门更新反而能获得比多次单比特门更好的局部性。

在实际工程中，可以用一个简单的指标来评估访存效率：每次门操作实际读取的缓存行数与理论最少缓存行数之比。理论最少缓存行数约为 2ⁿ·16 字节 / 64 字节 = 2ⁿ/4（假设每个复数 16 字节，缓存行 64 字节）。如果实际读取行数远高于此，说明 stride 过大或访问不连续，需要引入分块或重映射。这个指标可以在开发阶段用性能计数器（如 cache-misses）粗略观测，作为优化是否有效的判断依据。