最近我一直在琢磨一件事：量子计算模拟器到底难不难写？如果抛开 Qiskit、Cirq 这类现成框架，完全从零开始，用纯 Python 搭建一个能跑量子线路、能测量、能复现经典实验的模拟器，代码量会不会大得吓人？带着这个疑问，我把 DREAMVFIA 这个开源项目整理成了完整的实现笔记。这篇是上篇，先把最核心的框架、数据结构和基本算子讲透——量子态怎么表示、量子门怎么作用、测量怎么坍缩、贝尔态怎么跑通。整个实现只用 NumPy 一个第三方库，核心代码不足三百行，非常适合两类读者：一类是想搞懂量子计算底层数学逻辑的 Python 开发者，另一类是正在学量子信息、想找一个能动手改的模拟器教学项目的学生。看完这篇，你会对“量子比特就是一个复向量”、“量子门就是矩阵乘法”、“测量就是概率采样”这三句话有切身的体感。

1. 项目定位与整体设计思路

1.1 为什么非要自己写一个模拟器

市面上的量子模拟器其实不少。Qiskit Aer 能模拟上百个量子比特，Cirq 背靠 Google 的量子生态，ProjectQ 性能也很强。那为什么还要自己做？我的原因很简单，也很“笨”：不自己实现一遍，量子计算的很多概念永远是黑箱。

用 Qiskit 跑一个贝尔态，三行代码出结果，但你不会知道 1000 次测量里为什么恰好一半是 00、一半是 11。这个“为什么”不是靠文档能消化的，必须亲手把量子态向量写出来、把 H 门乘进去、把 CNOT 门作用上去，亲眼看到振幅是怎么从[1, 0, 0, 0]变成[0.707, 0, 0, 0.707]的，才算真正理解。

DREAMVFIA 的定位就是教学型模拟器，它刻意不追求性能和功能堆砌，而是追求可读性和可扩展性。每一个矩阵、每一步位运算都没有隐藏魔法，读代码的过程就是复习量子计算基础的过程。对于给学生上课的场景，这种“能拆开看内部结构”的代码比任何商业框架都有价值。

1.2 分层架构与模块规划

写模拟器之前，我先把物理模型映射成软件模块。量子计算的运行过程可以分成四件事：准备量子态、施加量子门、测量、统计结果。于是 DREAMVFIA 的核心也分成四层：

这样的分层有一个明显的好处：每一层都可以独立测试。我写代码的时候是严格按层推进的，先完成量子态模块并跑通单元测试，再写量子门，再写线路调度，最后做测量。任何一个环节出错，都能快速定位到具体的层，而不是在几百行代码里大海捞针。

1.3 技术选型：Python 与 NumPy 的组合逻辑

选 Python 不是因为性能，而是因为它能把线性代数写得像数学公式。量子计算本质上是复向量空间上的线性变换，NumPy 提供的complex128复数类型、向量化运算和广播机制，几乎是为此量身定制的。

选 NumPy 而非纯 Python 列表也是同理。模拟器的核心操作是矩阵乘法和振幅遍历，如果用嵌套列表手动实现，不仅代码冗长，还容易在“列表拷贝”这种细节上踩坑。NumPy 的ndarray切片和拷贝语义清晰，配合dtype=np.complex128，一个振幅向量的内存占用和读写效率都远优于原生列表。

当然，代价是有的。Python 的循环逐索引处理 2 的 n 次方个振幅时，性能天花板很低。这个问题的应对策略我在第 5 节专门讲，这里先卖个关子。总而言之，教学项目的第一目标是“让人看懂”，在这个前提下，Python + NumPy 就是最优解。

2. 量子计算模拟的核心原理

2.1 量子比特：从比特到状态向量

经典比特只有 0 和 1 两个状态，量子比特则多了一个“叠加”的可能性。数学上，一个量子比特的状态用一个二维复向量表示：

|ψ⟩ = α|0⟩ + β|1⟩

其中 α 和 β 都是复数，分别叫“|0⟩ 态振幅”和“|1⟩ 态振幅”。它们满足归一化条件|α|² + |β|² = 1，这个条件的物理含义是：对量子比特进行测量，得到 0 的概率是|α|²，得到 1 的概率是|β|²，总概率必须等于 1。

你可以把量子比特想象成一根箭头在二维复平面里的方向向量，只不过这根箭头的“长度”固定为 1，但“方向”可以连续变化。经典比特只能指向两个固定方向（0 或 1），而量子比特可以指向圆周上的任意一个位置。正是这种“连续自由度”，让 n 个量子比特能编码 2 的 n 次方个复数振幅。

n 个量子比特的联合状态是一个 2ⁿ 维的复向量：

|ψ⟩ = Σ c_i |i⟩，其中 i 从 0 到 2ⁿ-1

每个下标 i 对应一个计算基态，比如三个量子比特时下标 5（二进制 101）对应的基态是|101⟩。整个模拟器最核心的数据结构，就是存放这 2ⁿ 个复数振幅的一维数组。

2.2 量子门：酉矩阵与矩阵乘法

量子门是对量子态施加的变换，数学上就是一个酉矩阵（酉矩阵满足 U†U = I，即共轭转置乘自身等于单位阵）。酉矩阵保证变换前后向量的长度不变，也就是归一化条件不会因为施加量子门而被破坏。

单比特量子门是一个 2×2 酉矩阵，作用在量子比特上就是一次矩阵乘法。以最经典的 Hadamard 门（简称 H 门）为例：

H = 1/√2 [[1, 1], [1, -1]]

H 门的作用是把|0⟩变成(|0⟩ + |1⟩)/√2，也就是将确定态变成等概率叠加态。这个门的矩阵写法是：

H|0⟩ = 1/√2 [[1, 1], [1, -1]] [1, 0]ᵀ = [1/√2, 1/√2]ᵀ

看到没有，本质上就是高中数学里的矩阵乘向量。Pauli-X 门（就是量子版 NOT 门）的矩阵是[[0, 1], [1, 0]]，作用在|0⟩上得到|1⟩，作用在|1⟩上得到|0⟩，和经典取反一模一样。

两比特门中最重要的是 CNOT 门（受控非门）。它有两个输入：控制比特和目标比特。如果控制比特是 1，就对目标比特取反；如果控制比特是 0，目标比特保持不变。CNOT 是一个 4×4 矩阵，具体形式在实现一节里给出。

这里有一个初学者最容易忽略的点：量子门必须保证可逆。经典与非门会丢弃信息，但量子门不会，因为酉矩阵一定可逆。这也是为什么模拟器的所有门操作都只是数值上的线性变换，不存在“if 分支丢弃数据”这种写法（测量除外）。

2.3 测量：Born 规则与波函数坍缩

测量是量子计算里最特殊的一环，也是模拟器里唯一有“随机性”和“破坏性”的操作。对一个量子比特测量，得到结果 k 的概率由振幅的模平方决定：

P(k) = |c_k|²

这个规则叫 Born 规则，它把抽象的复数振幅和现实世界里的观测概率联系了起来。模拟器的测量模块要做两件事：

第一，按概率随机采样。用np.random.random()生成一个 [0, 1) 的随机数，再计算累积概率，落在哪个区间就返回哪个测量结果。

第二，坍缩量子态。测量一旦发生，量子态不再是叠加态，而是“坍缩”到与测量结果对应的基态上。在数学上，这相当于把所有不满足测量结果的振幅置为零，然后重新归一化。

测量是不可逆的破坏性操作，这就是为什么模拟器里测量之后量子态会“变掉”。如果多次测量同一个量子比特，第一次测量会决定后续所有结果——这一点在模拟器里天然成立，因为坍缩后其他振幅已经归零了。

3. DREAMVFIA 核心代码实现

3.1 量子态的数据结构

量子态是整个模拟器的心脏。我选择把 n 个量子比特的状态存成一个长度为 2ⁿ 的一维复数数组，下标 i 的二进制表示对应计算基态|i⟩的振幅。初始化时全部振幅置 0，仅下标 0 处置 1，表示所有量子比特都处在|0...0⟩态。

import numpy as np class QuantumState: def __init__(self, n_qubits): self.n_qubits = n_qubits self.dim = 1 << n_qubits # 2^n self.amp = np.zeros(self.dim, dtype=np.complex128) self.amp[0] = 1.0 + 0.0j # |00...0> def probability(self, index): return abs(self.amp[index]) ** 2

dtype=np.complex128是必须写清楚的。如果不指定，NumPy 默认用float64，复数振幅会被截断成实数，整个模拟器直接废掉。这里还有一个隐含约定：下标 0 对应所有比特为 0 的状态，这也是绝大多数量子模拟器的默认约定。

3.2 量子门定义与线性代数基础

我在gates.py里把常用量子门定义成函数，返回 NumPy 矩阵。这样做的可读性比直接写一坨带数字的数组好得多，后续如果要加参数化门（比如相位门），只需要在这个函数里加一个角度参数。

def hadamard(): return np.array([[1, 1], [1, -1]]) / np.sqrt(2) def pauli_x(): return np.array([[0, 1], [1, 0]]) def pauli_z(): return np.array([[1, 0], [0, -1]]) def cnot(): return np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])

注意 CNOT 矩阵的行列顺序：输入态按|00⟩, |01⟩, |10⟩, |11⟩排列，矩阵把|10⟩映射到|11⟩，把|11⟩映射到|10⟩，这就是“控制比特为 1 时目标比特取反”的矩阵表达。

这些矩阵都是酉矩阵。写完之后建议用np.allclose(U @ U.conj().T, np.eye(len(U)))自检一下，这是我在开发中坚持做的第一层验证。

3.3 单比特门的施加与比特序

有了量子态和量子门，最关键的问题来了：怎么把一个 2×2 的矩阵“作用到”一个 2ⁿ 维向量中的某一个量子比特上？

我的做法是逐对索引处理。假设目标比特是 t，那么任意下标 i 都有一个“孪生兄弟”j = i ^ (1 << t)，它们的二进制位只在第 t 位不同。施加单比特门 U 时，对每一对 (i, j) 做局部矩阵乘法：

def apply_single_qubit_gate(state, gate, target): new_amp = state.amp.copy() u00, u01, u10, u11 = gate[0, 0], gate[0, 1], gate[1, 0], gate[1, 1] for i in range(state.dim): j = i ^ (1 << target) if i < j: # 每对只处理一次 a_i, a_j = state.amp[i], state.amp[j] new_amp[i] = u00 * a_i + u01 * a_j new_amp[j] = u10 * a_i + u11 * a_j state.amp = new_amp

这里有一个非常重要的约定：比特序是“小端”的，即 qubit 0 对应二进制的最低位。比如 2 个量子比特时下标 1（二进制 01）表示 qubit 0 为 1、qubit 1 为 0。

为什么i < j就能保证每对只处理一次？因为当第 t 位是 0 时，j = i + 2^t > i；当第 t 位是 1 时，j = i - 2^t < i。所以循环只在“第 t 位为 0 的那个下标”处完成这对振幅的更新，不会重复计算，也不会漏掉任何一对。

3.4 CNOT 门的施加逻辑

CNOT 是两比特门，作用于控制比特 c 和目标比特 t。当控制位为 1 时，交换目标位为 0 和 1 的振幅；控制位为 0 时不做任何事。

def apply_cnot(state, control, target): new_amp = state.amp.copy() for i in range(state.dim): if ((i >> control) & 1) == 1: # 控制位为 1 j = i ^ (1 << target) if i < j: new_amp[i] = state.amp[j] new_amp[j] = state.amp[i] state.amp = new_amp

同样用i < j避免把同一对振幅交换两次。这里再强调一次：交换操作必须基于原始state.amp读取，写入到new_amp，否则会在同一轮循环里把已经交换过的值又读出来，造成“二次交换”，结果完全错误。

3.5 测量与坍缩的实现

测量的实现分三步：计算目标比特为 1 的总概率、按概率随机采样、坍缩并归一化。

def measure(state, qubit_index): prob_one = 0.0 for i in range(state.dim): if (i >> qubit_index) & 1: prob_one += abs(state.amp[i]) ** 2 outcome = 1 if np.random.random() < prob_one else 0 mask = 1 << qubit_index if outcome == 1: for i in range(state.dim): if (i & mask) == 0: state.amp[i] = 0.0 else: for i in range(state.dim): if (i & mask) != 0: state.amp[i] = 0.0 norm = np.sqrt(np.sum(np.abs(state.amp) ** 2)) state.amp /= norm return outcome

坍缩后必须要做归一化。因为置零之后向量的模长不再等于 1，如果不除以新的模长，后续所有的概率计算都会失真。state.amp /= norm这一步是在模拟“测量后量子态重新成为合法量子态”的物理过程。

4. 实操验证：跑通你的第一个量子线路

4.1 环境准备与项目结构

代码只依赖 Python 3.9+ 和 NumPy。建议用虚拟环境隔离依赖：

python -m venv venv source venv/bin/activate pip install numpy pytest

项目结构我按模块职责拆成下面这样，每个文件只干一件事，测试文件也一一对应：

dreamvfia/ ├── dreamvfia/ │ ├── __init__.py │ ├── state.py │ ├── gates.py │ ├── circuit.py │ └── measure.py ├── examples/ │ └── bell_state.py └── tests/ ├── test_state.py ├── test_gates.py └── test_bell_state.py

circuit.py我封装了一个简单的线路类，把量子态、量子门和测量串起来：

from .state import QuantumState from . import gates from .gates import hadamard, pauli_x from .measure import measure as measure_qubit class QuantumCircuit: def __init__(self, n_qubits): self.n_qubits = n_qubits self.state = QuantumState(n_qubits) self.instructions = [] def h(self, target): self.instructions.append(("h", target)) apply_single_qubit_gate(self.state, hadamard(), target) def x(self, target): self.instructions.append(("x", target)) apply_single_qubit_gate(self.state, pauli_x(), target) def cnot(self, control, target): self.instructions.append(("cnot", control, target)) apply_cnot(self.state, control, target) def measure(self, target): return measure_qubit(self.state, target)

线路类把每一次操作记录到instructions列表里，这一步对调试验证很有用。跑完一个实验，可以打印这个列表，确认门是按预期顺序被施加的。

4.2 构建贝尔态并测量

贝尔态（Bell State）是量子纠缠里最经典的例子。线路只有两步：先把 qubit 0 用 H 门变成叠加态，再用 CNOT 门让 qubit 1 与 qubit 0 纠缠：

from dreamvfia import QuantumCircuit qc = QuantumCircuit(2) qc.h(0) qc.cnot(control=0, target=1)

运行之后直接看一眼振幅向量：

print(qc.state.amp.round(4)) # [0.7071+0.j 0.+0.j 0.+0.j 0.7071+0.j]

这个输出信息量很大。下标 0 和下标 3 的振幅非零且都等于1/√2，说明量子态是(|00⟩ + |11⟩)/√2——两个量子比特处于纠缠态。下标 1 和下标 2 的振幅为 0，意味着无论怎么测量，都不可能得到 01 或 10 的结果。

再对两个量子比特做 1000 次测量，统计结果分布：

from collections import Counter results = Counter() shots = 1000 for _ in range(shots): qc = QuantumCircuit(2) qc.h(0) qc.cnot(0, 1) b0 = qc.measure(0) b1 = qc.measure(1) results[(b0, b1)] += 1 print(results) # 典型输出: Counter({(0, 0): 498, (1, 1): 502})

完美复现了贝尔态的理论预言：只可能得到 00 或 11，各约 50% 概率，永远不可能出现 01 或 10。这背后就是量子纠缠的数学本质——测量 qubit 0 会瞬间决定 qubit 1 的状态，因为它们的振幅被 CNOT 门耦合在了一起。

4.3 单元测试与自检方法

我一个人写模拟器时最大的感受是：量子计算太容易“看起来对”了。很多错误不会让程序崩溃，只会让概率统计悄悄偏离一点点。所以自检必须做成自动化测试。下面是我在项目里保留的几组核心测试：

def test_hadamard_on_zero(): qc = QuantumCircuit(1) qc.h(0) p0 = abs(qc.state.amp[0]) ** 2 p1 = abs(qc.state.amp[1]) ** 2 assert abs(p0 - 0.5) < 1e-12 assert abs(p1 - 0.5) < 1e-12 def test_x_gate_flips_state(): qc = QuantumCircuit(1) qc.x(0) assert np.allclose(qc.state.amp, [0, 1]) def test_cnot_on_control_one(): qc = QuantumCircuit(2) qc.x(0) # qubit0 置为 |1⟩ qc.cnot(0, 1) assert np.allclose(qc.state.amp, [0, 0, 0, 1]) # |11⟩ def test_h_twice_is_identity(): qc = QuantumCircuit(1) qc.h(0) qc.h(0) assert np.allclose(qc.state.amp, [1, 0])

这几个测试覆盖了绝大部分底层逻辑。尤其是H门连续作用两次等于单位变换这一点，如果单比特门的比特序搞错了，这个测试会立刻报错。我在开发过程中就是靠这组测试，把 index 操作和矩阵方向上的各种小错误一网打尽的。

5. 开发中踩过的坑与排查经验

5.1 复数精度与归一化

第一个坑是我自己踩得最深的：state.amp.copy()拷贝出来的数组，类型是complex128，但如果在初始化时写了np.zeros(self.dim, dtype=float)，后面所有复数振幅都会被截断成实部，H 门算出来的1/√2倒是没问题，但相位门这类依赖虚部的门会全军覆没。

还有一个隐藏坑：概率计算必须用abs(x) ** 2，而不是x ** 2。复数(a+bj)的平方是(a+bj)² = a² - b² + 2abj，得到的结果还是一个复数，取出来当概率用简直是在随机数生成器上跳舞。正确的算法是(a+bj) * (a-bj)，也就是模平方，abs()函数做的就是这件事。

5.2 比特序混乱

比特序是模拟器开发里最容易出问题、也最难排查的问题。同一个振幅数组，如果 qubit 0 对应最低位，下标 1 表示|01⟩；如果对应最高位，下标 1 表示|10⟩。两种约定本身没有对错，但代码内部必须全程一致。

我建议在项目 README 的第一行就写明“qubit i 对应二进制第 i 位（最低位为 qubit 0）”，并且在所有变更振幅的代码里都用位运算(i >> target) & 1提取比特状态，不要靠手算数值。检查比特序问题的一个实用技巧：构造一个|1000⟩态（下标 8），分别对 qubit 0 和 qubit 3 施加 X 门，看结果下标是 9 还是 1——这一步就能快速确认你的比特序约定。

5.3 动态范围与数值稳定性

量子线路深了之后，振幅会变成非常小的浮点数，比如 30 层线路后振幅可能只有 1e-15 量级。这时候如果某个振幅应该是 0 但实际是 1e-16，测量概率计算中abs()之后几乎无影响，但归一化时可能会引入微小误差。

我的处理方法是：在关键测试里统一用np.allclose(actual, expected, atol=1e-12)而不是精确相等，给数值误差留出余地。同时在测量坍缩之后立刻重新归一化，避免误差在多次测量中累积。

5.4 性能瓶颈在哪里

最后聊一下性能。这个实现每个门都是遍历 2ⁿ 个振幅，时间复杂度 O(2ⁿ)，12 量子比特以下体验流畅，20 比特以上开始卡顿。内存方面，complex128每个振幅占 16 字节，可以估算一下：

到 30 量子比特，个人电脑基本就到极限了。教学场景里 10~15 比特完全够用，这也是 DREAMVFIA 现阶段的目标。如果真需要模拟更大的规模，优化方向是把循环改成 NumPy 向量化操作，或者对稀疏状态用字典存储非零振幅，再进一步可以上 GPU 或分布式计算——这些都是下篇会展开的内容。

5.5 问题速查表

6. 扩展方向与下篇预告

6.1 从“能跑”到“好用”的三步

当前实现能跑通贝尔态，但距离“好用”还有三步要走。第一步是补齐量子门库，把 Toffoli 门、SWAP 门、相位门 S/T、参数化旋转门 Rx/Ry/Rz 都加进去，这会覆盖大部分量子算法的教学需求。第二步是支持部分测量和复位操作，很多算法（比如量子隐形传态）中间需要测量某个比特并据此决定后续操作，这要求线路模型支持经典比特和控制流。第三步是测量统计的可视化，用 Matplotlib 画概率直方图，教学演示时直观很多。

6.2 性能优化路线

如果要往更大规模走，我有几条明确的优化路线。最直接的是把单比特门施加从 for 循环改成向量化操作，用数组切片和矩阵广播一次性更新所有振幅对，实测能提速几十倍。其次是引入稀疏振幅表示，很多算法（比如 GHZ 态制备）的中间态只有少部分非零振幅，用字典或稀疏矩阵存储能大幅扩展规模。最后是并行化，不同振幅对之间的更新完全独立，天然适合多线程或 GPU 加速，这也是 DREAMVFIA 下篇要展示的重头戏。

6.3 下篇能学到的内容

说了这么多，这篇文章其实只覆盖了“静态量子线路”——量子态按固定的门序列演化，最后一次性测量。下篇我们会做一个质的飞跃：实现带经典控制流的动态线路，跑通量子隐形传态和 Deutsch-Jozsa 算法，再把性能优化做一个完整的基准测试对比。到那时候，你会看到同一个模拟器从 12 比特模拟到 20 比特以上，靠的完全是向量化技巧和数据结构设计，而不是换语言。

我从这个项目里最大的体会是：模拟量子计算机这件事，门槛远没有想象中高，但细节的坑远比想象中多。一个看似简单的小问题——比如比特序约定——就能让你调试一晚上。好在只要把每个环节拆开、弄懂、测试，一切都变得清晰可控。希望这篇上篇笔记能帮你迈过“量子计算很神秘”这道心理门槛，如果你自己也跟着写了一遍，可以在评论区聊聊，尤其是你踩过的那些坑，大概率我当年也踩过。