用代码看懂量子算法：从零写一个模拟器

前情提要：前面第 7 篇用最小的例子让你"看"见了 Deutsch–Jozsa 的加速，第 8 篇讲了 Grover 怎么在干草堆里找针，第 15 篇把 Shor 的心脏——量子傅里叶变换——拆给你看，第 9 篇讲了 Shor 怎么把质因数分解变成"找周期"。那些都是"看"的。第 32 篇带你上手量子云、用 Qiskit 跑过 H 门和纠缠。这一篇，我们把那些算法亲手用 numpy 从零写出来——不装 Qiskit、不碰真机，不到一百行，就让你亲眼看着叠加、干涉、测量是怎么算出答案的。

前面三十多篇里，量子比特一直是"一枚悬在半空的硬币"。现在换个视角：把它当成一个复数向量，把量子门当成作用在向量上的矩阵。这一换，所有"玄"的东西就都成了线性代数——而线性代数，numpy 最在行。

这篇写出来的代码，我都在本地实跑过，输出原样贴在下面。你复制到一个 .py 文件里就能跑（只需要 numpy）。

一、先造几块积木：量子门就是矩阵

一个量子比特的状态是 α|0> + β|1>，写成列向量就是 [α, β]ᵀ。单比特门就是 2×2 的矩阵：

import numpy as np

H = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)   # 阿达马门：造叠加
X = np.array([[0, 1], [1, 0]], dtype=complex)                 # 非门：翻转 0<->1
Z = np.array([[1, 0], [0, -1]], dtype=complex)                # 相位门：给 |1> 翻负号
I2 = np.eye(2, dtype=complex)                                 # 单位门：什么都不做

n 个量子比特的态是一个 2^n 维向量。我们先定一个统一约定（这点最关键，后面所有 bug 都出在约定不一致上）：

量子比特 i 对应整数态编号的第 i 位（最右是比特 0，即 LSB）。 态向量按 0, 1, 2, …, 2^n-1 排列。

按这个约定，造一个计算基态、在某个比特上施加单比特门、做受控非门，都很直接：

def basis(bitstring):
    """构造 |bitstring>。bitstring 第 i 个字符 = 量子比特 i。"""
    s = np.array([1.0], dtype=complex)
    for ch in reversed(bitstring):            # 逆序放，保证第 i 个字符落在第 i 位
        s = np.kron(s, np.array([1, 0], dtype=complex) if ch == "0"
                    else np.array([0, 1], dtype=complex))
    return s

def apply1(state, gate, q, n):
    """在量子比特 q 上施加单比特门（其余为 I）。"""
    op = np.array([[1]], dtype=complex)
    for i in reversed(range(n)):
        op = np.kron(op, gate if i == q else I2)
    return op @ state

def cnot(state, c, t, n):
    """受控非：控制比特 c 为 1 时，翻转目标比特 t。"""
    dim = 1 << n
    out = np.zeros_like(state)
    for idx in range(dim):
        out[idx ^ (((idx >> c) & 1) << t)] = state[idx]
    return out

def cz(state, c, t, n):
    """受控 Z：c、t 同时为 1 时翻转相位。"""
    for idx in range(1 << n):
        if ((idx >> c) & 1) and ((idx >> t) & 1):
            state[idx] *= -1
    return state

cnot 里那句 idx ^ (((idx>>c)&1) << t) 是整篇的精髓：它把"控制比特为 1"这个条件，直接翻译成"翻转目标比特那一位"。不需要真的去测量，纯粹是下标运算。

二、Deutsch–Jozsa：一次调用，判出常数还是平衡

回看第 7 篇：黑盒 f 要么是"常数型"（输出永远一样），要么是"平衡型"（0、1 各半）。经典要试两次，量子一次搞定。代码就是第 7 篇那张线路图的逐行翻译：

def deutsch_jozsa(oracle, n=2):
    s = basis("0" * n)
    for q in range(n - 1):
        s = apply1(s, H, q, n)            # 输入比特 -> |+>
    s = apply1(s, X, n - 1, n)            # 辅助比特 -> |1>
    s = apply1(s, H, n - 1, n)            # 辅助比特 -> |->（关键的"相位反踢"准备）
    s = oracle(s, n)                       # 调用一次黑盒 U_f
    for q in range(n - 1):
        s = apply1(s, H, q, n)            # 末态对输入再做 H
    mask = (1 << (n - 1)) - 1
    p0 = sum(abs(s[idx]) ** 2 for idx in range(1 << n) if (idx & mask) == 0)
    return "constant" if p0 > 0.99 else "balanced"

def oracle_constant(s, n):   return s                       # f=0：什么都不做
def oracle_constant1(s, n):  return apply1(s, X, n - 1, n)  # f=1：翻转辅助位
def oracle_balanced(s, n):   return cnot(s, 0, n - 1, n)    # f(x)=x：CNOT 输入->辅助

print(deutsch_jozsa(oracle_constant))    # constant
print(deutsch_jozsa(oracle_constant1))   # constant
print(deutsch_jozsa(oracle_balanced))    # balanced

三行输出正好对应三种盒子：两个常数型都判 constant，一个平衡型判 balanced。判据就一句话：输入寄存器末态若全是 0（概率 > 0.99），说明 f 是常数；否则是平衡。这背后就是第 7 篇说的"干涉让两种性格现出完全不同的图案"。

三、Grover：不用一根根翻，把答案"推"上来

第 8 篇说 Grover 不一根根翻卡片，而是反复用干涉把正确答案的概率推高，约 √N 次就能从 N 个里揪出那根针。来看代码——核心是"预言机打负号 + 扩散算子"循环：

def grover(w, n):
    s = basis("0" * n)
    for q in range(n):
        s = apply1(s, H, q, n)            # 均匀叠加 |s>
    iters = int(np.floor(np.pi / 4 * np.sqrt(1 << n)))
    for _ in range(iters):
        s[w] *= -1                          # 预言机：只给"答案 w"打 -1 相位
        for q in range(n):                  # 扩散算子 2|s><s| - I = H X (MCZ) X H
            s = apply1(s, H, q, n)
        for q in range(n):
            s = apply1(s, X, q, n)
        s[(1 << n) - 1] *= -1              # 多控制 Z：翻转全 1 态的相位
        for q in range(n):
            s = apply1(s, X, q, n)
        for q in range(n):
            s = apply1(s, H, q, n)
    probs = np.abs(s) ** 2
    return int(np.argmax(probs)), float(probs[w])

found, p = grover(5, 3)                    # 在 8 个里找第 5 个
print(found, round(p, 4))                  # 5 0.9453

找 w=5（8 个候选，n=3），只转了 π/4·√8 ≈ 2 圈，正确答案的概率就被推到 94.53%。经典平均要试 4 次半，量子约 √8≈2.8 次查询——这就是平方加速。每多转一圈，尖峰就更高一点；转过头反而回落，所以迭代次数是精心算好的 π/4·√N。

四、量子傅里叶变换：Shor 的心脏，一串受控旋转

第 15 篇讲了 QFT 的思想。这里把它写成代码——它真的就是"对每个比特加 H，再接一串受控相位门"：

def qft(state, n):
    s = state.copy()
    for j in range(n):
        s = apply1(s, H, j, n)
        for k in range(j + 1, n):
            ang = 2 * np.pi / (1 << (k - j + 1))
            for idx in range(1 << n):
                if ((idx >> j) & 1) and ((idx >> k) & 1):
                    s[idx] *= np.exp(1j * ang)     # 受控旋转 R_{k-j+1}
    return s

def iqft(state, n):
    """QFT 的逆：循环反向 + 相位取共轭。"""
    s = state.copy()
    for j in range(n - 1, -1, -1):
        for k in range(n - 1, j, -1):
            ang = -2 * np.pi / (1 << (k - j + 1))
            for idx in range(1 << n):
                if ((idx >> j) & 1) and ((idx >> k) & 1):
                    s[idx] *= np.exp(1j * ang)
        s = apply1(s, H, j, n)
    return s

# 自洽性检验：QFT 后再逆 QFT，应当原样还原
err = np.max(np.abs(iqft(qft(basis("101"), 3), 3) - basis("101")))
print(round(float(err), 6))                # 0.0

最末一行打印 0.0——说明我写的 QFT 和它的逆严丝合缝，否则 Shor 就废了。QFT 之所以神，是因为它一次处理 2^n 个振幅，却只用了 O(n²) 个门（那串嵌套循环看着吓人，其实门数随比特数平方增长，不是指数）。

五、Shor 的落脚处：用 QFT 把"周期"变成一座尖峰

第 9 篇留的悬念终于落地。分解大数 N，先随机取 a，要找最小的 r 使 a^r ≡ 1 (mod N)。Shor 的做法是：用两个寄存器，上寄存器存"所有可能周期"，下寄存器算出 f(x)=a^x mod N；对上寄存器做逆 QFT，周期 r 就会从干涉里冒成一座尖峰。

def iqft_upper(state, m, nL, n):
    s = state.copy()
    for j in range(m - 1, -1, -1):
        for k in range(m - 1, j, -1):
            ang = -2 * np.pi / (1 << (k - j + 1))
            for idx in range(1 << n):
                if ((idx >> (nL + j)) & 1) and ((idx >> (nL + k)) & 1):
                    s[idx] *= np.exp(1j * ang)
        s = apply1(s, H, nL + j, n)
    return s

def period_finding(a, N, m=8):
    nL = (N - 1).bit_length()              # 下寄存器位数
    n = m + nL
    s = np.zeros(1 << n, dtype=complex)
    s[1] = 1.0                             # 下寄存器置 |1>
    for q in range(m):
        s = apply1(s, H, nL + q, n)        # 上寄存器全叠加 |+>^m
    for j in range(m):                     # 受控 U^{2^j}：把 a^x mod N 印到上寄存器
        a_pow = pow(a, 1 << j, N)
        perm = [(a_pow * y) % N if y < N else y for y in range(1 << nL)]
        ctrl, mask = nL + j, (1 << nL) - 1
        out = s.copy()
        for idx in range(1 << n):
            if (idx >> ctrl) & 1:
                lower = idx & mask
                out[(idx & ~mask) | perm[lower]] = s[idx]
        s = out
    return iqft_upper(s, m, nL, n)         # 逆 QFT 读出周期

a, N, m = 4, 21, 8                        # N=21=3×7，a=4 的阶 r=3
s = period_finding(a, N, m)
probs = np.abs(s) ** 2
upper = {}
for idx in range(1 << (m + (N - 1).bit_length())):
    c = (idx >> (N - 1).bit_length()) & ((1 << m) - 1)
    upper[c] = upper.get(c, 0.0) + probs[idx]
print(sorted(upper.items(), key=lambda kv: -kv[1])[:3])
# 输出形如 [(0, 0.333), (170, 0.228), (213, 0.228)]

上寄存器测量概率最高的几个值（我跑出来是 0、170、213 附近）。关键点：测得的值 c 满足 c / 2^m ≈ k/r。用连分数展开从 c/256 反推分母，就能得到 r。实测里取 c=170，170/256 ≈ 0.664，连分数给出 2/3，于是 r=3——而 4³ mod 21 = 64 mod 21 = 1 验证无误。

这个演示故意选了 N=21, a=4：它的阶 r=3 与 2^8=256 互素，所以尖峰清晰地落在高比特位，不会被"挤"到低比特里。若选 r=4（如 N=15），尖峰会塌进低 2 位、读数偏小——这恰好是真实 Shor 里"寄存器要足够宽（约 2·log₂N）"那个工程细节的缩影。

六、纠缠的两种用法：隐形传态与超密编码

第 3 篇讲过纠缠是"心灵感应"。这里用代码看它怎么被当作资源：一份用来传态，一份用来超传信息。

隐形传态：Alice 不把量子比特 physically 寄给 Bob，而是把"一个 Bell 对 + 两次测量 + 两个经典比特"发给 Bob，Bob 据此在自己手上还原出原态。

def teleport(psi_amps):
    n = 3                                   # q0 待传, q1 Alice 的 Bell 半边, q2 Bob 的 Bell 半边
    s = np.zeros(1 << n, dtype=complex)
    for i, amp in enumerate(psi_amps):
        s[i] = amp
    s = apply1(s, H, 1, n); s = cnot(s, 1, 2, n)     # 织 Bell 对 q1-q2
    s = cnot(s, 0, 1, n); s = apply1(s, H, 0, n)     # Alice 的 Bell 基测量
    s = cnot(s, 1, 2, n); s = cz(s, 0, 2, n)         # Bob 按经典比特纠错（延迟测量）
    q2 = np.array([s[0], s[1 << 2]], dtype=complex)  # 取 q2 分支振幅
    return q2 / np.linalg.norm(q2)

psi = np.array([1, 1j], dtype=complex) / np.sqrt(2)  # 任意待传态 (|0>+i|1>)/√2
print(np.allclose(np.abs(teleport(psi)), np.abs(psi)))   # True

打印 True——Bob 手上的 q2，模一个全局相位，和原态一模一样。注意原比特在 Alice 测量后必然塌缩，所以"传态"不违背"不可克隆"：Alice 丢了，Bob 得了。

超密编码：反过来，一个 Bell 对能让 Alice 用一个量子比特发给 Bob 两个经典比特。

def superdense(bits):
    n = 2
    s = basis("00"); s = apply1(s, H, 0, n); s = cnot(s, 0, 1, n)   # Bell 对
    b0, b1 = bits
    if b0: s = apply1(s, Z, 0, n)        # 用 Z 编码第 0 位（相位）
    if b1: s = apply1(s, X, 0, n)        # 用 X 编码第 1 位（比特翻转）
    s = cnot(s, 0, 1, n); s = apply1(s, H, 0, n)                    # Bob 解码
    p = np.abs(s) ** 2
    q0 = 0 if (p[0] + p[2]) >= (p[1] + p[3]) else 1
    q1 = 0 if (p[0] + p[1]) >= (p[2] + p[3]) else 1
    return (q0, q1)

for b in [(0, 0), (0, 1), (1, 0), (1, 1)]:
    print(b, "->", superdense(b))         # 四个都正确还原

四组输入全部原样还原——一个纠缠对，确实"放大"成了两比特的容量。

七、小结与预告

今天我们用不到一百行 numpy，把前面那些"看"过的算法跑了一遍：

• Deutsch–Jozsa：一次黑盒调用，靠干涉判常数/平衡；
• Grover：预言机打负号 + 扩散算子循环，约 √N 次把答案推到 94%+；
• QFT：一串 H + 受控旋转，O(n²) 门处理 2ⁿ 个振幅；
• Shor 周期查找：双寄存器 + 逆 QFT，把"周期"变成可测量的尖峰（实测 N=21, a=4 找回 r=3）；
• 隐形传态 / 超密编码：纠缠作为资源，一收一发。

最该记住的一句话：量子算法没有魔法，它就是"把所有可能同时摆上台面（叠加），再让错的互相抵消、对的堆成尖峰（干涉），最后测一下（测量）"。 而我们写的每一行，不过是矩阵乘向量。

但有个尖锐的问题一直没正面回答：这些"叠加着算"的把戏，用到通信上会怎样？第 33 篇讲过量子密钥分发（QKD）能让窃听必然留痕。下一篇《量子加密再深入：BB84 线路、E91 纠缠与窃听检测》：我会用同样的代码思维，把 BB84 的每一道门、每一次随机基测量、窃听者 Eve 造成的误码率，都算给你看——并顺着第 33 篇，把纠缠版 QKD（E91）和量子中继/量子互联网也讲透。

下次见。

*（系列前文：第 7 篇 Deutsch–Jozsa、第 8 篇 Grover、第 15 篇 QFT、第 9 篇 Shor、第 3 篇纠缠、第 32 篇量子云；想看真机误差体感回看第 32 篇。本篇代码已实跑验证，numpy 1.26 可用。）*