量子计算机上的多项式时间大数分解算法

本文基于 Peter W. Shor 1995 年论文《Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer》整理。原稿篇幅较长，涵盖量子计算模型、可逆逻辑、量子傅里叶变换、大数分解与离散对数五大部分。本文聚焦于大数分解算法的核心思路与数学推导，对量子计算模型与可逆计算部分做精简介绍，力求以现代标准记号呈现这一里程碑式算法的完整脉络。

1. 引言：从经典困难性到量子突破

1.1 计算复杂性与 Church 论题

计算科学的奠基性成果之一，是 Church、Turing 与 Post 在 1936 年分别独立建立的可计算性与不可计算性的区分。其核心是 Church 论题：一切计算装置都可以被图灵机模拟。

随着实用计算机的发展，人们意识到"可计算"与"不可计算"的区分过于粗糙。计算机科学家更关心的是：一个函数能以多高的效率被计算。普遍接受的折中方案是以多项式时间与超多项式时间来划分高效与低效。将问题规模（输入位数）记为 (n)，若算法的步骤数关于 (n) 至多为多项式，则称该算法是多项式时间的。

定量 Church 论题（Strong Church's Thesis）断言：任何物理计算装置都可以被图灵机以多项式步数模拟。然而，量子力学的引入对这一论题提出了挑战。

1.2 量子计算的早期发展

Benioff（1980）首先研究了计算与量子力学的交互，证明了可逆的幺正演化足以实现图灵机的计算能力。Feynman（1982）进一步提出量子力学可能在计算上本质上难以被经典模拟。Deutsch（1985）则首次明确提出：使用量子力学原理的计算机是否能比经典计算机更高效地计算？

此后，Deutsch–Jozsa（1992）与 Simon（1994）分别给出了量子计算机相对经典计算机具有指数级加速的问题，但这些问题都属于"预言机问题"（oracle problem），尚未触及实际重要的计算问题。

Shor 的突破在于：大整数分解和离散对数这两个被广泛认为经典困难的问题，可以在量子计算机上以多项式时间求解。这两个问题正是 RSA 等公钥密码体系的安全性基础。

2. 量子计算模型

2.1 量子比特与叠加态

考虑一个由 (n) 个二态分量组成的系统。经典物理中，完整描述该系统状态只需 (n) 个比特；而在量子物理中，需要 (2^n - 1) 个复数。

具体地，系统的状态处于 (2^n) 维复向量空间中。对每个经典状态（如 (|011\cdots 0\rangle)），都有对应的基向量。系统的一般状态是这些基向量的叠加态：

$$
|\psi\rangle = \sum_{i=0}^{2^n - 1} a_i |S_i\rangle,
$$

其中振幅 (a_i) 为复数，满足归一化条件 (\sum_i |a_i|^2 = 1)。对系统进行测量时，观察到基态 (|S_i\rangle) 的概率为 (|a_i|^2)，而测量会将系统投影到该基态上。

2.2 量子门与幺正变换

量子计算中的状态变换必须是幺正变换（unitary transformation），即满足 (U^\dagger U = I) 的矩阵。计算仅允许局部幺正变换——即作用在固定数量比特上的变换。

通用量子门集由所有一位比特门加上一个二位门（受控非门，CNOT）构成。CNOT 门的作用是：当且仅当第一位为 1 时，对第二位取反。

例如，以下是一个二位量子门的真值表：

$$
|00\rangle \to |00\rangle, \qquad |01\rangle \to |01\rangle,
$$

$$
|10\rangle \to \frac{1}{\sqrt{2}}(|10\rangle + |11\rangle), \qquad |11\rangle \to \frac{1}{\sqrt{2}}(|10\rangle - |11\rangle).
$$

其对应的矩阵为

$$
U = \begin{pmatrix}
1 & 0 & 0 & 0 \
0 & 1 & 0 & 0 \
0 & 0 & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \
0 & 0 & \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}}
\end{pmatrix}.
$$

2.3 量子干涉

量子计算的一个关键特性是干涉。假设系统处于叠加态

$$
\frac{1}{\sqrt{2}}|10\rangle - \frac{1}{\sqrt{2}}|11\rangle,
$$

对上述门 (U) 作用后，结果为

$$
\frac{1}{2}(|10\rangle + |11\rangle) - \frac{1}{2}(|10\rangle - |11\rangle) = |11\rangle.
$$

原本 (|10\rangle) 和 (|11\rangle) 各自都有被观测到的概率，但在叠加态下，(|10\rangle) 的概率振幅完全相消，观测结果确定为 (|11\rangle)。这就是量子干涉的力量——通过精心设计算法，使"正确答案"的振幅相长叠加，"错误答案"的振幅相消。

2.4 量子门阵列

量子门阵列（quantum gate array）由一系列量子门通过逻辑"导线"连接而成。输入经过一系列量子门后，在最后一步测量所有比特的值，得到输出。

与经典门阵列类似，量子门阵列也需要满足一致性条件：门阵列的设计必须能在经典多项式时间内生成，且门矩阵中的元素必须是可计算数。

量子图灵机与量子门阵列若允许小概率错误，则在多项式时间内可计算相同的函数类，记为 BQP（Bounded-error Quantum Polynomial time）。

3. 可逆计算与模幂运算

3.1 可逆计算的必要性

量子计算中的幺正变换天然要求完全可逆：知道输出状态就能唯一确定输入状态。这与经典计算机通过耗散能量实现不可逆计算截然不同。

幸运的是，Lecerf（1963）与 Bennett（1973）已经证明：任何确定性计算都可以被转化为可逆计算，只需在时间上付出常数倍的代价。

3.2 Toffoli 门与 Fredkin 门

可逆计算的通用门有两个经典代表：

Toffoli 门（受控受控非门）：当且仅当前两位都为 1 时，对第三位取反。若将第三位初始化为 1，则输出为前两位输入的 NAND。由于 NAND 是经典通用门，Toffoli 门是可逆计算的通用门。

Fredkin 门（受控交换门）：当第一位为 0 时交换后两位，为 1 时保持不变。通过适当设置辅助位，可以实现 AND 与 NOT 运算，因此也是通用门。

3.3 Bennett 方法：清除中间结果

可逆计算的一个关键问题是：计算过程中产生的中间"记录"（RECORD）必须被清除，否则会影响量子干涉。

Bennett 的方法如下：

  1. 复制输入 (x)。
  2. 保留一份输入，用可逆门计算 (F(x))，将中间结果存入 RECORD 寄存器。
  3. 将 (F(x)) 复制到预设为零的输出寄存器。
  4. 逆向运行第二步的计算，清除 RECORD 和 OUTPUT 寄存器。

最终效果是将 (x) 映射为 (F(x))，同时清除了所有中间垃圾。

若还需清除输入 (x)（即用 (F(x)) 替换 (x)），则还需要 (F^{-1}) 的多项式时间算法。方法是：先正向计算得到 ((x, F(x)))，再可逆地计算 ((F(x), F^{-1}(F(x))) = (F(x), x))，然后逆向运行得到 (F(x))。

3.4 可逆模幂运算

模幂运算是量子分解算法的瓶颈子程序：给定 (n, x, a)，计算 (x^a \bmod n)。

经典上，通过反复平方法，计算 (x^{2^i} \bmod n)（(i < l)，其中 (l) 为位数），然后将 (a) 的二进制展开中为 1 的对应项相乘。这需要 (O(l^3)) 时间（使用初等乘法）或 (O(l^2 \log l \log\log l)) 时间（使用 Schönhage–Strassen 算法）。

将上述过程可逆化后，空间复杂度可控制在 (O(l))。具体伪代码如下：

power := 1
for i = 0 to l - 1
    if (a_i == 1) then
        power := power * x^{2^i} (mod n)
    endif
endfor

其中 (x^{2^i} \bmod n) 可以预先经典计算并嵌入门阵列结构中。输入 (a) 保持不变，输出为 (x^a \bmod n)。

这个门阵列可以进一步分解为模乘子程序，而模乘又可以通过模加的反复累加来实现。清除输入 (b) 时，利用 (\gcd(c, n) = 1) 时存在模逆 (c^{-1} \bmod n) 的事实，通过逆向模减将 (b) 归零。

4. 量子傅里叶变换

4.1 定义

量子傅里叶变换是 Shor 算法的核心。对于 (q = 2^l)，变换 (A_q) 将基态 (|a\rangle) 映射为

$$
|a\rangle \mapsto \frac{1}{\sqrt{q}} \sum_{c=0}^{q-1} \exp!\left(\frac{2\pi i a c}{q}\right) |c\rangle.
$$

即幺正矩阵的第 ((a, c)) 元素为 (\frac{1}{\sqrt{q}} \exp(2\pi i a c / q))。

4.2 量子线路实现

当 (q = 2^l) 时，量子傅里叶变换可以用两类门高效实现：

(R_j) 门：作用在第 (j) 个比特上，

$$
R_j = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \ 1 & -1 \end{pmatrix}.
$$

(S_{j,k}) 门（(j < k)）：作用在第 (j) 和第 (k) 个比特上，

$$
S_{j,k} = \begin{pmatrix}
1 & 0 & 0 & 0 \
0 & 1 & 0 & 0 \
0 & 0 & 1 & 0 \
0 & 0 & 0 & e^{i\theta_{k-j}}
\end{pmatrix},
$$

其中 (\theta_{k-j} = \pi / 2^{k-j})。

对 (l) 个比特，按以下顺序施加门（从左到右）：

$$
R_{l-1}; S_{l-2,l-1}; R_{l-2}; S_{l-3,l-1}; S_{l-3,l-2}; R_{l-3}; \cdots; R_1; S_{0,l-1}; \cdots; S_{0,1}; R_0.
$$

即从 (R_{l-1}) 到 (R_0) 逆序施加 (R_j)，在 (R_{j+1}) 与 (R_j) 之间施加所有 (S_{j,k})（(k > j)）。总共需要 (l(l-1)/2 + l = l(l+1)/2) 个门——关于输入位数是多项式量级。

4.3 正确性验证

从 (|a\rangle = |a_{l-1}\cdots a_0\rangle) 到 (|b\rangle = |b_{l-1}\cdots b_0\rangle) 的振幅中，各 (R_j) 门的因子 (\frac{1}{\sqrt{2}}) 相乘得到总因子 (\frac{1}{\sqrt{q}})。相位部分为

$$
\sum_{0 \le j \le k < l} \frac{\pi , a_j , b_k}{2^{k-j}}.
$$

设 (c) 为 (b) 的比特翻转（bit-reversal），经过变量替换后，上式化为

$$
\sum_{j,k=0}^{l-1} \frac{2\pi}{2^l} \cdot 2^j a_j \cdot 2^k c_k = \frac{2\pi a c}{q},
$$

恰好是量子傅里叶变换所需的相位。

4.4 近似傅里叶变换

当 (k - j) 很大时，(S_{j,k}) 门施加的相位因子 (e^{i\pi/2^{k-j}}) 非常接近 1。Coppersmith（1994）证明：可以忽略这些微小相位因子，定义一个近似傅里叶变换，其精度仍足以用于大数分解。这大幅减少了所需的量子门数量。

5. 大数分解算法

5.1 从分解到求阶

直接分解 (n) 的问题可以归约为求阶问题：给定 (x)（(\gcd(x, n) = 1)），求最小正整数 (r) 使得

$$
x^r \equiv 1 \pmod{n}.
$$

归约方法如下：随机选取 (x \bmod n)，用量子算法求其阶 (r)，然后计算 (\gcd(x^{r/2} - 1, n))。

原理：由于 (x^r \equiv 1 \pmod{n})，有

$$
(x^{r/2} - 1)(x^{r/2} + 1) \equiv 0 \pmod{n}.
$$

因此 (\gcd(x^{r/2} - 1, n)) 通常是 (n) 的非平凡因子。算法失败的条件是 (r) 为奇数或 (x^{r/2} \equiv -1 \pmod{n})。

成功概率：设 (n) 有 (k) 个不同的奇素因子。可以证明，随机选取的 (x) 使算法成功的概率至少为 (1 - 1/2^{k-1})。因此，只要 (n) 不是奇素数的幂（这种情况可用经典方法处理），重复数次即可以高概率得到因子。

5.2 量子求阶算法

输入：(x, n)。

步骤：

  1. 选取 (q) 为 2 的幂，使得 (n^2 \le q < 2n^2)。

  2. 将第一个寄存器置于均匀叠加态：

$$
\frac{1}{\sqrt{q}} \sum_{a=0}^{q-1} |a\rangle |0\rangle.
$$

这只需将第一个寄存器的每个比特置于 (\frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)) 态。

  1. 在第二个寄存器中计算 (x^a \bmod n)：

$$
\frac{1}{\sqrt{q}} \sum_{a=0}^{q-1} |a\rangle |x^a \bmod n\rangle.
$$

  1. 对第一个寄存器施加量子傅里叶变换 (A_q)：

$$
\frac{1}{q} \sum_{a=0}^{q-1} \sum_{c=0}^{q-1} \exp!\left(\frac{2\pi i a c}{q}\right) |c\rangle |x^a \bmod n\rangle.
$$

  1. 测量两个寄存器，得到某个 (c) 和 (x^k \bmod n)（其中 (0 \le k < r)）。

5.3 测量概率分析

测得状态 (|c, x^k \bmod n\rangle) 的概率为

$$
P(c, k) = \frac{1}{q} \left| \sum_{\substack{a:, x^a \equiv x^k}} \exp!\left(\frac{2\pi i a c}{q}\right) \right|^2.
$$

因为 (x) 的阶为 (r)，条件 (x^a \equiv x^k \pmod{n}) 等价于 (a \equiv k \pmod{r})。令 (a = br + k)，概率化为

$$
P(c, k) = \frac{1}{q} \left| \sum_{b=0}^{\lfloor(q-k-1)/r\rfloor} \exp!\left(\frac{2\pi i (br+k)c}{q}\right) \right|^2.
$$

提取因子 (\exp(2\pi i k c / q))（模为 1，不影响概率），并记 ({rc}_q) 为 (rc \bmod q) 落在 ((-q/2, q/2]) 中的余数，概率近似为

$$
P(c, k) \approx \frac{1}{q} \left| \sum_{b=0}^{\lfloor q/r \rfloor} \exp!\left(\frac{2\pi i b {rc}_q}{q}\right) \right|^2.
$$

5.4 关键不等式

当 (|{rc}_q| \le r/2) 时，求和中各振幅的相位接近，叠加后得到较大的概率。将求和近似为积分，可证明此时概率的下界为

$$
P(c, k) \ge \frac{4}{\pi^2 r^2} > \frac{1}{3r^2}.
$$

条件 (|{rc}_q| \le r/2) 等价于存在整数 (d) 使得

$$
\left| \frac{c}{q} - \frac{d}{r} \right| \le \frac{1}{2q}.
$$

由于 (q > n^2)，分母不超过 (n) 的分数 (d/r) 若满足上述不等式则是唯一的。

5.5 连分数展开提取阶

测得 (c) 后，已知 (c/q)。利用连分数展开（continued fraction expansion）可以在多项式时间内找到 (c/q) 的最佳有理逼近 (d/r)（分母小于 (n)）。

若 (\gcd(d, r) = 1)，则直接得到 (r)。更一般地，(d/r) 约分后给出的分母是 (r) 的因子。

5.6 成功概率与重复次数

能使算法成功的状态 ((c, x^k)) 的总数至少为 (r \cdot \varphi(r))（其中 (\varphi) 为 Euler 函数），每个的概率至少为 (1/(3r^2))，因此成功概率至少为

$$
\frac{\varphi(r)}{3r}.
$$

由数论定理 (\varphi(r)/r > \delta / \log\log r)，只需重复 (O(\log\log r)) 次实验即可以高概率确定 (r)。

实际优化：

  • 不仅尝试 (c)，还尝试 (c \pm 1, c \pm 2, \ldots)
  • 若 (d/r) 约分后为 (d'/r')，尝试 (r') 的小倍数 (2r', 3r', \ldots) 作为候选阶
  • 若有两个候选阶 (r_1, r_2)，测试 (\operatorname{lcm}(r_1, r_2))

这些优化可将期望实验次数降至 (O(1))。

5.7 具体示例

设 (n = 33)，选 (x = 5)。则 (5) 模 (33) 的阶为 (r = 10)（因为 (5^{10} \equiv 1 \pmod{33})）。

取 (q = 256)（满足 (33^2 = 1089 \le 256) 不成立，此处为演示取较小的 (q)）。测量后，(c) 高概率出现在 (q/r = 256/10 = 25.6) 的整数倍附近，即 (c \approx 0, 26, 51, 77, 102, 128, 154, 179, 205, 230)。

例如测得 (c = 51)，则 (c/q = 51/256 \approx 0.199)。连分数展开给出最佳逼近 (1/5 = 0.2)，于是 (d = 1, r = 5) 或 (r = 10)。验证 (5^5 = 3125 \equiv 3125 - 94 \times 33 = 3125 - 3102 = 23 \not\equiv 1)，而 (5^{10} \equiv 1 \pmod{33})，故 (r = 10)。

进而计算 (\gcd(5^5 - 1, 33) = \gcd(3124, 33) = \gcd(3124, 33))。由于 (3124 = 94 \times 33 + 22)，(\gcd(33, 22) = 11)。因此得到 (33 = 3 \times 11) 的一个非平凡因子 (11)。

6. 算法复杂度总结

步骤 经典复杂度 量子复杂度
模幂运算 (x^a \bmod n) (O(l^3)) (O(l^3)) 时间，(O(l)) 空间
量子傅里叶变换 — (O(l^2)) 门
连分数展开 — (O(l^2))（经典后处理）
总计 指数时间 (O(l^2 \log l \log\log l)) 量子步

其中 (l = O(\log n)) 为输入位数。量子算法的总步数关于 (\log n) 是多项式的，而最好的经典算法（数域筛法）需要 (\exp(O((\log n)^{1/3} (\log\log n)^{2/3}))) 时间——这是指数时间。

7. 离散对数问题

7.1 问题描述

设 (p) 为素数，(g) 为模 (p) 的原根。给定 (x)，求 (r) 使得

$$
g^r \equiv x \pmod{p}.
$$

最快的经典算法（数域筛法的推广）需要 (\exp(O((\log p)^{1/3} (\log\log p)^{2/3}))) 时间。

7.2 量子算法

使用三个量子寄存器。选取 (q) 为接近 (p) 的 2 的幂（(p < q < 2p)）。

  1. 将前两个寄存器置于均匀叠加态，计算第三个寄存器：

$$
\frac{1}{p-1} \sum_{a=0}^{p-2} \sum_{b=0}^{p-2} |a\rangle |b\rangle |g^a x^{-b} \bmod p\rangle.
$$

  1. 对前两个寄存器分别施加量子傅里叶变换 (A_q)，将 (|a\rangle \to |c\rangle)，(|b\rangle \to |d\rangle)：

$$
\frac{1}{(p-1)q} \sum_{a,b=0}^{p-2} \sum_{c,d=0}^{q-1} \exp!\left(\frac{2\pi i(ac+bd)}{q}\right) |c\rangle |d\rangle |g^a x^{-b} \bmod p\rangle.
$$

  1. 测量后，通过分析 ((c, d)) 的值，利用与分解算法类似的连分数方法，可以在多项式时间内恢复出离散对数 (r)。

7.3 推广

该算法不仅适用于 (\mathbb{Z}_p) 上的离散对数，还可推广到：

  • (\mathbb{Z}_{p^\alpha}) 上的离散对数（只要乘法群是循环群）
  • 一般 Abel 群上的离散对数（Boneh–Lipton 1995 的推广）

8. 物理实现与开放问题

8.1 退相干与精度

构建量子计算机的最大障碍是退相干（decoherence）——量子叠加态与环境的相互作用导致信息丢失。Bennett 等人（1994）证明，量子门只需精度 (O(1/t)) 即可完成 (t) 步计算。

8.2 量子纠错

经典计算机可通过冗余和纠错码降低错误率。量子比特的不可克隆定理禁止直接复制量子态，但更复杂的纠错方法已被提出。Shor（1995）给出了减少量子存储器退相干的方案，Bennett 等人（1995）证明了量子比特可以通过噪声信道进行忠实传输。

8.3 更广泛的应用

大数分解与离散对数本身并非广泛有用的问题——它们的重要性主要来自公钥密码学。量子计算机若能解决的唯一困难问题就是这两类，那它将成为一种专用工具。

更有意义的开放问题是：量子计算机能否高效求解 NP 完全问题？目前有一些微弱的迹象表明答案可能是否定的（Bennett et al. 1994），但这一问题尚未定论。另一个有前景的方向是格上短向量问题（short vector in a lattice），它可能适合量子算法求解。

参考文献

[1] Shor P W. Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer[J]. SIAM Journal on Computing, 1997, 26(5): 1484–1509.

[2] Shor P W. Algorithms for quantum computation: Discrete logarithms and factoring[C]// Proceedings of the 35th Annual Symposium on Foundations of Computer Science. IEEE Computer Society Press, 1994: 124–134.

[3] Coppersmith D. An approximate Fourier transform useful in quantum factoring[R]. IBM Research Report RC 19642, 1994.

[4] Simon D. On the power of quantum computation[C]// Proceedings of the 35th Annual Symposium on Foundations of Computer Science. IEEE Computer Society Press, 1994: 116–123.

[5] Nielsen M A, Chuang I L. Quantum Computation and Quantum Information[M]. Cambridge University Press, 2000.