别被量子劝退：手写一个 Deutsch-Jozsa，看看它到底快在哪

量子计算的书翻开第一页，看到狄拉克符号和希尔伯特空间，我就合上了。不是我不想学，是那种「从数学公理出发推导一切」的讲法，实在让人头疼。

但后来我发现，Deutsch-Jozsa 算法这东西，用纯经典的 Python 代码就能把核心逻辑跑通，而且跑完之后你会立刻明白量子加速到底快在哪——不需要先啃完一整本线性代数。

这篇文章我们就手写一个 Deutsch-Jozsa 模拟器，用经典计算机模拟量子计算的过程，然后亲眼看看它凭什么比经典算法快。

问题本身：一个看似无聊但信息量极大的场景

先来看我们要解决什么问题。

假设你有一个黑盒函数 f(x)，输入是一个 n 位的二进制串（比如 "1010"），输出是 0 或 1。这个函数被承诺只可能是两种情况之一：

恒定函数：对任意输入，永远返回 0，或者永远返回 1；
平衡函数：对于一半的输入返回 0，对另一半返回 1。

你的任务是判断这个黑盒到底是恒定的还是平衡的。

经典算法怎么做？最坏情况下，你需要查询这个黑盒 2^(n-1) + 1 次。为什么？假设 n=3，总共有 8 种输入。你运气极差，前 4 次查询都返回了 0，这时候你仍然无法下结论——万一下一次返回 1 呢？只有当你查到第 5 次，如果还是 0，你才能断定它是恒定的（因为平衡函数最多只能有一半返回 0）。

Deutsch-Jozsa 算法只需要查询 1 次。不管 n 多大，一次就够了。

这听起来像作弊，但它是真的。我们把它拆开看。

量子比特的经典模拟：用数组表示叠加态

在真实量子计算机里，一个量子比特可以同时处于 0 和 1 的叠加态。在经典计算机里模拟它，最直接的办法是用向量和矩阵。

单个量子比特的状态用一个二维复数向量表示：

import numpy as np

# |0⟩ 状态
state_0 = np.array([1, 0], dtype=complex)

# |1⟩ 状态
state_1 = np.array([0, 1], dtype=complex)

多个量子比特用张量积组合。两个量子比特就有 4 种基态：|00⟩、|01⟩、|10⟩、|11⟩，对应一个 4 维向量。

# 两个量子比特都初始化为 |0⟩
state_2q = np.kron(state_0, state_0)  # [1, 0, 0, 0]

量子门是酉矩阵，作用在状态向量上就是矩阵乘法。Hadamard 门是最关键的一个，它把确定的基态变成等概率叠加：

H = np.array([[1, 1],
              [1, -1]]) / np.sqrt(2)

# 对第一个量子比特施加 H 门，第二个保持不变
H_on_first = np.kron(H, np.eye(2))
state_after_H = H_on_first @ state_2q

这就是经典模拟量子计算的基本思路：状态是向量，运算是矩阵乘法，测量是概率采样。计算量随量子比特数指数增长（n 个量子比特需要 2^n 维向量），但模拟小规模完全可行。

Deutsch-Jozsa 的核心电路：就三步

Deutsch 算法的原始版本只处理 1 个输入比特，Jozsa 把它推广到了 n 个比特。电路结构如下：

|0⟩^⊗n  ── H^⊗n ── U_f ── H^⊗n ── 测量
|1⟩      ── H ────── U_f ──────────── 忽略

我们用 n=3 来具体实现。初始化分两部分：前 n 个量子比特全部置为 |0⟩，最后一个辅助量子比特置为 |1⟩。

def initialize_state(n):
    """初始化 n 个工作比特 + 1 个辅助比特，辅助比特从 |1⟩ 开始"""
    # 工作比特全是 |0⟩
    work_state = np.array([1, 0], dtype=complex)
    for _ in range(n - 1):
        work_state = np.kron(work_state, np.array([1, 0], dtype=complex))
    
    # 辅助比特 |1⟩
    ancilla = np.array([0, 1], dtype=complex)
    
    return np.kron(work_state, ancilla)

第一步，对所有量子比特施加 Hadamard 门。

def apply_hadamard_all(state, num_qubits):
    """对所有量子比特同时加 H 门"""
    H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
    H_all = H
    for _ in range(num_qubits - 1):
        H_all = np.kron(H_all, H)
    return H_all @ state

施加完 H 门后，量子态变成了所有可能输入的等权重叠加。对于 n=3 且带辅助比特的情况，工作比特部分包含了从 |000⟩ 到 |111⟩ 共 8 种基态，每种振幅相等。

第二步是 oracle，也就是黑盒函数 U_f。在量子电路里，U_f 的作用是：把输入 |x⟩|y⟩ 变成 |x⟩|y ⊕ f(x)⟩。这个操作必须用酉矩阵表示。

def make_oracle_matrix(n, f):
    """根据经典函数 f 构建 oracle 酉矩阵"""
    dim = 2 ** (n + 1)  # n 个工作比特 + 1 个辅助比特
    U_f = np.zeros((dim, dim), dtype=complex)
    
    for x in range(2 ** n):
        for y in range(2):  # 辅助比特可以是 0 或 1
            # 输入基态的索引：x 在高位，y 在低位
            input_idx = (x << 1) | y
            # 输出：y 变成 y XOR f(x)
            output_y = y ^ f(x)
            output_idx = (x << 1) | output_y
            U_f[output_idx, input_idx] = 1
    
    return U_f

这里的关键细节：oracle 矩阵的每一列只有一个 1，其余全是 0，这保证它是酉矩阵（实际上是置换矩阵）。f(x) 可以是任意满足恒定或平衡条件的函数——我们可以在构造时自由选择。

第三步，对工作比特再次施加 H 门（辅助比特不动）。

def apply_hadamard_work(state, n):
    """只对前 n 个工作比特加 H 门"""
    H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
    H_work = H
    for _ in range(n - 1):
        H_work = np.kron(H_work, H)
    # 辅助比特保持恒等
    H_full = np.kron(H_work, np.eye(2))
    return H_full @ state

最后测量前 n 个工作比特。如果全部测到 0，说明函数是恒定的；只要有一个测到 1，就是平衡的。

def measure_work_qubits(state, n):
    """测量前 n 个工作比特，返回测量结果的概率分布"""
    dim = 2 ** n
    probs = np.zeros(dim)
    
    for i in range(dim):
        # 提取工作比特为 i、辅助比特为 0 和 1 的概率
        prob = abs(state[i * 2]) ** 2 + abs(state[i * 2 + 1]) ** 2
        probs[i] = prob
    
    return probs

拼起来跑一遍，看看到底发生了什么

def deutsch_jozsa(n, f):
    """模拟 Deutsch-Jozsa 算法，n 是工作比特数，f 是经典函数"""
    state = initialize_state(n)
    state = apply_hadamard_all(state, n + 1)  # 所有比特加 H
    U_f = make_oracle_matrix(n, f)
    state = U_f @ state
    state = apply_hadamard_work(state, n)
    probs = measure_work_qubits(state, n)
    return probs

现在构造两个测试函数：

# 恒定函数：永远返回 0
def constant_zero(x):
    return 0

# 平衡函数：输入 x 的二进制表示中 1 的个数是奇数则返回 1
def balanced_parity(x):
    return bin(x).count('1') % 2

跑一下 n=3 的情况：

n = 3
probs_constant = deutsch_jozsa(n, constant_zero)
probs_balanced = deutsch_jozsa(n, balanced_parity)

print("恒定函数测量结果：", probs_constant)
print("平衡函数测量结果：", probs_balanced)

输出会是：

恒定函数测量结果： [1. 0. 0. 0. 0. 0. 0. 0.]
平衡函数测量结果： [0. 0. 0. 0. 0. 0. 0. 1.]

对于恒定函数，100% 概率测到 |000⟩；对于平衡函数，0% 概率测到 |000⟩，具体哪个态被测量到取决于函数，但关键是 |000⟩ 的概率一定是 0。

这就意味着：你只需要看测量结果是不是全 0，就能一次判定函数类型。

为什么一次就够了？干涉是关键

经典算法慢的原因是：每次查询只能获得一个 f(x) 的值，信息是局部的。

量子算法里，第一次 H 门把所有输入叠加起来，oracle 在同一个查询中计算了所有 f(x) 的值，但这些值被编码在了量子态的相位里，而不是直接可读的振幅里。第二次 H 门做了一件事：让这些相位信息通过量子干涉转化为可测量的振幅。

具体来说，对于恒定函数，所有路径的相位变化一致，干涉是相长的，振幅全部集中在 |00...0⟩ 上。对于平衡函数，一半路径相位翻转，干涉是相消的，|00...0⟩ 的振幅被精确抵消到零。

这就是量子加速的本质：不是并行计算了所有可能，而是通过干涉提取了全局性质。并行计算本身不够——如果你并行算出了 2^n 个结果，测量时只能随机得到一个，跟经典算法没区别。干涉让你能提取「所有结果的某种整体模式」，这才是量子的杀手锏。

那这个算法到底有什么用？

说实话，Deutsch-Jozsa 本身没什么实际用途。谁会在现实中遇到一个被承诺「要么恒定要么平衡」的黑盒？

但它是量子算法设计的教科书级示范，展示了几个核心技巧：

相位反冲：oracle 通过辅助比特把 f(x) 的值编码到输入态的相位里，而不是直接修改输入态的值；
干涉提取全局信息：用第二次 H 门把相位差转化为振幅差，让测量结果反映函数的整体性质；
单次查询获取指数级信息：虽然不能读出所有 f(x)，但足以回答一个全局性问题。

这些技巧在后续真正有用的量子算法（比如 Shor 算法、Grover 算法）里反复出现。Deutsch-Jozsa 就是那个最小的、能完整展示这些技巧的工作模型。

把上面这段代码跑一遍，改改 n 的值，换几个不同的平衡函数试试，你会发现结果永远成立。这种「亲手验证一个反直觉结论」的体验，比看十页数学推导都管用。

常见问题

为什么经典算法最坏要查 2^(n-1)+1 次，而不是 2^n 次？

因为平衡函数的定义是一半输出 0、一半输出 1。最坏情况下，你前 2^(n-1) 次查询全返回了同一种结果（比如全是 0），此时你仍然不能确定——如果是平衡函数，正好有一半是 0，你运气差全查到了这一半。但第 2^(n-1)+1 次查询如果是 0，那就超过一半了，只能是恒定函数。所以最坏情况就是正好踩到那个边界上。

模拟代码里 oracle 矩阵的维度为什么是 2^(n+1)？

n 个工作比特加上 1 个辅助比特，总共 n+1 个量子比特。每个量子比特有 2 个基态，所以整个系统的希尔伯特空间维度是 2^(n+1)。oracle 矩阵必须在这个空间上作用，并且是酉矩阵。代码里对每一对 (x, y) 计算输出索引，就是在这个 2^(n+1) 维空间里做置换。

n 很大的时候这段代码还能跑吗？

跑不了。n=10 时状态向量就有 2048 维，矩阵是 2048×2048，内存还勉强撑得住。n=20 时维度超过 200 万，经典计算机直接爆炸。这就是为什么经典模拟量子计算只能用于小规模验证——但小规模足够你理解算法原理了。