Shor算法，一种用于分解大整数为素数乘积的量子算法

Shor算法是一种量子算法，用于分解大整数为素数的乘积。这是一个量子计算课程辅导员发布的英文笔记的翻译版本。

• A Probably More Detailed Explanation of the Shor's Algorithm

PKC 公钥密码学

在公钥密码学中，Alice发布一个公钥，Bob使用该公钥来加密消息。然后，Alice使用她的私钥来解密消息。这依赖于密码学的非对称性：使用公钥加密消息很容易，但在没有私钥的情况下解密消息却很困难。这又依赖于单向函数的存在。

单向函数

• 单向函数是指那些正向计算容易，但反向计算困难的函数。例如，RSA公钥加密系统使用因数分解作为单向函数：将两个素数 $p$ 和 $q$ 相乘得到一个合数 $N$ 很容易，但反过来，将一个合数 $N$ 分解为因数 $p$ 和 $q$ 却在数学上非常困难。
• 事实上，已知的最好的经典因数分解算法——数域筛法（number field sieve）需要 $\exp (\mathrm{\Theta }(n^{1/3}{\log }^{2/3}n))$ 次操作，其中 $n=\lceil {\log }_{2}N\rceil$，即 $N$ 的比特数。
• 对于密码学应用来说，单向函数的数学定义至关重要，即它们在平均情况下必须难以反转，而不仅仅是在最坏情况下难以反转。

量子力学因式分解

1994年，Peter Shor发明了一种量子算法，可以在多项式时间内分解整数。这仍然是量子计算最重要和最令人印象深刻的潜在应用之一。Shor的算法建立在之前的查询复杂性算法之上，并基于两个关键见解：

• 量子傅里叶变换算法可以用来解决阶（和周期）查找的数学问题。
• 因数分解问题可以归约为阶/周期查找问题。

周期查找

让我们来看一个周期函数：

$$f(x)=a^x\mathrm{mod}N$$

其中 $a$ 和 $N$ 是正整数，且 $a$ 小于 $N$，并且它们没有公因数。周期或阶（$r$）是满足以下条件的最小（非零）整数：

$$a^r\mathrm{mod}N=1$$

Shor的解决方案是对酉算子使用量子相位估计：

$$U|y⟩\equiv |ay\mathrm{mod}N⟩$$

为了理解这如何有用，让我们看看 $U$ 的特征态可能是什么样子。如果我们从状态 $|1⟩$ 开始，我们可以看到每次应用 $U$ 都会将寄存器的状态乘以 $a$（模 $N$），并且在应用 $r$ 次后，我们将再次到达状态 $|1⟩$。例如，当 $a=3$ 和 $N=35$ 时：

$$\begin{array}{rl}U|1⟩& =|3⟩\\ U^2|1⟩& =|9⟩\\ U^3|1⟩& =|27⟩\\ & ⋮\\ U^{r-1}|1⟩& =|12⟩\\ U^r|1⟩& =|1⟩\end{array}$$

因此，这个循环中的状态叠加态（$|u_0⟩$）将是 $U$ 的特征态：

$$\begin{array}{rl}|u_0⟩& =\frac{1}{\sqrt{r}}\sum _{k=0}^{r-1}|a^k\mathrm{mod}N⟩\\ |u_0⟩& =\frac{1}{\sqrt{12}}(|1⟩+|3⟩+|9⟩+\cdots +|4⟩+|12⟩)\\ U|u_0⟩& =\frac{1}{\sqrt{12}}(U|1⟩+U|3⟩+U|9⟩+\cdots +U|4⟩+U|12⟩)\\ & =\frac{1}{\sqrt{12}}(|3⟩+|9⟩+|27⟩+\cdots +|12⟩+|1⟩)\\ & =|u_0⟩\end{array}$$

这个特征态的特征值为1，这并不是很有趣。一个更有趣的特征态可能是每个计算基态的相位都不同的情况。具体来说，让我们看看在第k个状态的相位与k成正比的情况：

$$\begin{array}{rl}|u_1⟩& =\frac{1}{\sqrt{r}}\sum _{k=0}^{r-1}{e}^{-\frac{2\pi ik}{r}}|a^k\mathrm{mod}N⟩\\ U|u_1⟩& =e^{\frac{2\pi i}{r}}|u_1⟩\end{array}$$

当$a=3$, $N=35$时，我们可以看到这个特征态的形式：

$$\begin{array}{rl}|u_1⟩& =\frac{1}{\sqrt{12}}(|1⟩+e^{-\frac{2\pi i}{12}}|3⟩+e^{-\frac{4\pi i}{12}}|9⟩+\cdots +e^{-\frac{20\pi i}{12}}|4⟩+e^{-\frac{22\pi i}{12}}|12⟩)\\ U|u_1⟩& =\frac{1}{\sqrt{12}}(|3⟩+e^{-\frac{2\pi i}{12}}|9⟩+e^{-\frac{4\pi i}{12}}|27⟩+\cdots +e^{-\frac{20\pi i}{12}}|12⟩+e^{-\frac{22\pi i}{12}}|1⟩)\\ U|u_1⟩& =e^{\frac{2\pi i}{12}}\cdot \frac{1}{\sqrt{12}}(e^{-\frac{2\pi i}{12}}|3⟩+e^{-\frac{4\pi i}{12}}|9⟩+e^{-\frac{6\pi i}{12}}|27⟩+\cdots +e^{-\frac{22\pi i}{12}}|12⟩+e^{-\frac{24\pi i}{12}}|1⟩)\\ U|u_1⟩& =e^{\frac{2\pi i}{12}}|u_1⟩\end{array}$$

（我们可以看到相位的分母中出现了 $r=12$。）

这是一个特别有趣的特征值，因为它包含了 $r$。实际上，必须包含 $r$ 以确保 $r$ 个计算基态之间的相位差是相等的。这并不是唯一具有这种行为的特征态；为了进一步推广这一点，我们可以将一个整数 $s$ 乘以这个相位差，这将在我们的特征值中体现出来：

$$\begin{array}{rl}|u_s⟩& =\frac{1}{\sqrt{r}}\sum _{k=0}^{r-1}{e}^{-\frac{2\pi isk}{r}}|a^k\mathrm{mod}N⟩\\ U|u_s⟩& =e^{\frac{2\pi is}{r}}|u_s⟩\\ |u_s⟩& =\frac{1}{\sqrt{12}}(|1⟩+e^{-\frac{2\pi is}{12}}|3⟩+e^{-\frac{4\pi is}{12}}|9⟩+\cdots +e^{-\frac{20\pi is}{12}}|4⟩+e^{-\frac{22\pi is}{12}}|12⟩)\\ U|u_s⟩& =\frac{1}{\sqrt{12}}(|3⟩+e^{-\frac{2\pi is}{12}}|9⟩+e^{-\frac{4\pi is}{12}}|27⟩+\cdots +e^{-\frac{20\pi is}{12}}|12⟩+e^{-\frac{22\pi is}{12}}|1⟩)\\ U|u_s⟩& =e^{\frac{2\pi is}{12}}\cdot \frac{1}{\sqrt{12}}(e^{-\frac{2\pi is}{12}}|3⟩+e^{-\frac{4\pi is}{12}}|9⟩+e^{-\frac{6\pi is}{12}}|27⟩+\cdots +e^{-\frac{22\pi is}{12}}|12⟩+e^{-\frac{24\pi is}{12}}|1⟩)\\ U|u_s⟩& =e^{\frac{2\pi is}{12}}|u_s⟩\end{array}$$

我们现在对于每个整数 $s$（其中 $0\le s\le r-1$）都有一个唯一的特征态。非常方便的是，如果我们将所有这些特征态相加，不同的相位会抵消掉所有的计算基态，除了 $|1⟩$：

$$\frac{1}{\sqrt{r}}\sum _{s=0}^{r-1}|u_s⟩=|1⟩$$

为此，我们将看一个较小的例子，其中 $a=7$ 且 $N=15$。在这种情况下，$r=4$：

$$\begin{array}{rl}|u_0⟩& =\frac{1}{2}(|1⟩+|7⟩+|4⟩+|13⟩)\\ |u_1⟩& =\frac{1}{2}(|1⟩+e^{-\frac{2\pi i}{4}}|7⟩+e^{-\frac{4\pi i}{4}}|4⟩+e^{-\frac{6\pi i}{4}}|13⟩)\\ |u_2⟩& =\frac{1}{2}(|1⟩+e^{-\frac{4\pi i}{4}}|7⟩+e^{-\frac{8\pi i}{4}}|4⟩+e^{-\frac{12\pi i}{4}}|13⟩)\\ |u_3⟩& =\frac{1}{2}(|1⟩+e^{-\frac{6\pi i}{4}}|7⟩+e^{-\frac{12\pi i}{4}}|4⟩+e^{-\frac{18\pi i}{4}}|13⟩)\end{array}$$$$\frac{1}{2}(|u_0⟩+|u_1⟩+|u_2⟩+|u_3⟩)=|1⟩$$

由于计算基态 $|1⟩$ 是这些特征态的叠加，这意味着如果我们使用状态 $|1⟩$ 对 $U$ 进行量子相位估计（QPE），我们将测量一个相位：

$$\varphi =\frac{s}{r}$$

其中 $s$ 是 $0$ 到 $r-1$ 之间的一个随机整数。我们最终使用算法在 $\varphi$ 上找到 $r$。

从周期查找中分解

• 并不是所有的因数分解问题都很难；我们可以立即发现一个偶数，并知道它的一个因数是2。实际上，选择难以分解的数字有特定的标准，但基本思路是选择两个大素数的乘积。
• 一个通用的因数分解算法首先会检查是否有因数分解的捷径（这个数是偶数吗？这个数是 $N=a^b$ 的形式吗？），然后在最坏的情况下使用 Shor 的周期查找算法。由于我们旨在关注算法的量子部分，我们将直接跳到 $N$ 是两个素数的乘积的情况。

示例：分解 $N=15$

1. 第一步是选择一个介于 $1$ 和 $N-1$ 之间的随机数 $a$。我们可以选择 $11$，它不是 $N$ 的非平凡因数（如果是，则将 $N$ 除以这个 $a$ 并重新开始）。
2. 接下来，我们对 $a=11$ 和 $N=15$ 进行 Shor 的阶数查找算法。记住，我们测量的相位将是 $\frac{s}{r}$，其中$$a^r\mathrm{mod}N=1$$并且 $s$ 是 $0$ 到 $r-1$ 之间的一个随机整数。如果 $r$ 不是偶数，我们不能继续，必须尝试不同的 $a$ 值。
3. 现在我们有了 $r$，我们可能可以用它来找到 $N$ 的一个因数。因为：$$a^r\mathrm{mod}N=1$$那么：$$(a^r-1)\mathrm{mod}N=0$$这意味着 $N$ 必须整除 $a^r-1$。如果 $r$ 也是偶数，那么我们可以写成：$$a^r-1=(a^{r/2}-1)(a^{r/2}+1)$$那么 $N$ 和 $a^{r/2}-1$ 或 $a^{r/2}+1$ 的最大公约数很可能是 $N$ 的一个适当因数。如果在最后一步中，我们通过连分数算法找到 $r=2$。那么我们有：$$a^{r/2}-1={11}^{2/2}-1=10\text{和}{a}^{r/2}+1={11}^{2/2}+1=12$$10 和 15 的（非平凡）最大公约数是 5，12 和 15 的（非平凡）最大公约数是 3。然后我们可以验证 $3\times 5=15$，表明我们已经成功地将 15 分解为 5 和 3。

What is Order Finding?

一些数论基础

为了解释阶数查找问题以及如何使用相位估计来解决它，解释一些数论中的基本概念并引入一些方便的符号将非常有帮助。

首先，对于任何给定的正整数 $N$，我们定义一个集合：

$$Z_N=\{0,1,\dots ,N-1\}$$

例如，$Z_1=\{0\}$，$Z_2=\{0,1\}$，$Z_3=\{0,1,2\}$，等等。

这些是数字的集合，但我们可以将它们看作不仅仅是集合。特别地，我们可以考虑在 $Z_N$ 上的算术运算，例如加法和乘法——如果我们同意总是取模 $N$ 的结果，那么当我们执行这些运算时，我们将始终保持在这个集合内。（加法和乘法这两个特定的运算，都是取模 $N$ 的，将 $Z_N$ 变成一个环，这是代数中一种基本重要的对象。）

例如，3 和 5 是 $Z_7$ 的元素，如果我们将它们相乘，我们得到 $3\cdot 5=15$，除以 7 余 1。

有时我们将其表示如下：

$$3\cdot 5\equiv 1(\mathrm{mod}7)$$

但我们也可以简单地写成 $3\cdot 5=1$，前提是已经明确我们在 $Z_7$ 中工作，以使我们的符号尽可能简单和清晰。

例如，以下是 $Z_6$ 的加法和乘法表：

加法表：

$$\begin{array}{ccccccc}+& 0& 1& 2& 3& 4& 5\\ 0& 0& 1& 2& 3& 4& 5\\ 1& 1& 2& 3& 4& 5& 0\\ 2& 2& 3& 4& 5& 0& 1\\ 3& 3& 4& 5& 0& 1& 2\\ 4& 4& 5& 0& 1& 2& 3\\ 5& 5& 0& 1& 2& 3& 4\end{array}$$

乘法表：

$$\begin{array}{ccccccc}\cdot & 0& 1& 2& 3& 4& 5\\ 0& 0& 0& 0& 0& 0& 0\\ 1& 0& 1& 2& 3& 4& 5\\ 2& 0& 2& 4& 0& 2& 4\\ 3& 0& 3& 0& 3& 0& 3\\ 4& 0& 4& 2& 0& 4& 2\\ 5& 0& 5& 4& 3& 2& 1\end{array}$$

在 $Z_N$ 的 $N$ 个元素中，满足 $gcd(a,N)=1$ 的元素 $a\in Z_N$ 是特殊的。通常包含这些元素的集合用星号表示，如下所示：

$$Z_N^{\ast }=\{a\in Z_N:gcd(a,N)=1\}$$

如果我们将注意力集中在乘法运算上，集合 $Z_N^{\ast }$ 形成一个群——特别是一个阿贝尔群——这是代数中另一种重要的对象。关于这些集合的一个基本事实（实际上关于有限群的一般事实）是，如果我们选择任何元素 $a\in Z_N^{\ast }$ 并反复将 $a$ 自乘，我们最终总会得到数字 1。

举个例子，设 $N=6$。我们有 $5\in Z_6^{\ast }$，因为 $gcd(5,6)=1$，并且如果我们将 5 自乘，我们得到 1（如上表所示）。

$$5^2=1(\text{在 }{Z}_6\text{ 中})$$

再举一个例子，设 $N=21$。如果我们遍历从 0 到 20 的数字，这些是与 21 的最大公约数等于 1 的数字：

$$Z_{21}^{\ast }=\{1,2,4,5,8,10,11,13,16,17,19,20\}$$

对于这些元素中的每一个，都可以将该数字提升到一个正整数幂以得到 1。

以下是这些元素满足条件的最小幂次：

$$\begin{array}{rl}{1}^1& =1\\ 2^6& =1\\ 4^3& =1\\ 5^6& =1\\ 8^2& =1\\ {10}^6& =1\\ {11}^6& =1\\ {13}^2& =1\\ {16}^3& =1\\ {17}^6& =1\\ {19}^6& =1\\ {20}^2& =1\end{array}$$

自然地，我们在所有这些方程中都在 $Z_{21}$ 中工作，我们没有写出来——我们认为这是隐含的，以避免使事情变得混乱。在整个课程中，我们将继续这样做。

你可以检查上面的每一个方程，以及这些是使方程成立的最小正整数幂次，使用以下代码单元（根据需要更改数字）。

N = 21
a = 17
max_power = 12
print("k \t a^k \n  ")
for k in range(1, max_power+1):
    print(f"{k} \t {a**k % N}")

注意，在我们回到 1 之后，循环重复——这是有道理的，因为将 1 乘以 $a$ 会将我们带回 $a$，这是我们开始的地方。

虽然这对于课程的目的并不是必需的，但我们也可以检查当 $gcd(a,N)\ne 1$ 时，我们永远不会回到 1——所以我们依赖于 $a\in Z_N^{\ast }$ 这一事实来使其工作。

问题陈述

现在我们可以陈述阶数查找问题。

输入：正整数 $N$ 和 $a$，满足 $gcd(N,a)=1$

输出：最小的正整数 $r$，使得 $a^r\equiv 1(\mathrm{mod}N)$

或者，用我们刚才介绍的符号来说，我们给定 $a\in Z_N^{\ast }$，并且我们正在寻找最小的正整数 $r$，使得 $a^r=1$。这个数 $r$ 被称为 $a$ 模 $N$ 的阶数。

与 $Z_N^{\ast }$ 中的元素相乘

为了将阶数查找问题与相位估计联系起来，让我们考虑一个系统上的操作，其经典状态对应于 $Z_N$，我们通过一个固定的元素 $a\in Z_N^{\ast }$ 进行乘法。

$$M_a|x⟩=|ax⟩(\text{对于每个 }x\in Z_N)$$

需要明确的是，我们在 $Z_N$ 中进行乘法，所以在方程右侧的态矢量中隐含地取模 $N$。例如，如果 $N=15$ 且 $a=2$，我们有：

$$\begin{array}{rl}{M}_2|0⟩& =|0⟩\\ M_2|1⟩& =|2⟩\\ M_2|2⟩& =|4⟩\\ M_2|3⟩& =|6⟩\\ M_2|4⟩& =|8⟩\\ M_2|5⟩& =|10⟩\\ M_2|6⟩& =|12⟩\\ M_2|7⟩& =|14⟩\\ M_2|8⟩& =|1⟩\\ M_2|9⟩& =|3⟩\\ M_2|10⟩& =|5⟩\\ M_2|11⟩& =|7⟩\\ M_2|12⟩& =|9⟩\\ M_2|13⟩& =|11⟩\\ M_2|14⟩& =|13⟩\end{array}$$

只要 $gcd(a,N)=1$，这就是一个酉操作。它将标准基态 $\{|0⟩,\dots ,|N-1⟩\}$ 的元素重新排列，因此作为矩阵，它是一个置换矩阵。从其定义可以看出，这个操作是确定性的，一个简单的方法来证明它是可逆的就是考虑 $a$ 模 $N$ 的阶数 $r$，并认识到 $M_a$ 的逆是 $M_a^{r-1}$。

$$M_a^{r-1}{M}_a=M_a^r=M_a^r=M_1=I$$

还有另一种不需要任何 $r$ 知识的方法来考虑逆——毕竟，这是我们试图计算的。对于每个元素 $a\in Z_N^{\ast }$，总有一个唯一的元素 $b\in Z_N^{\ast }$ 满足 $ab=1$。我们将这个元素 $b$ 表示为 $a^{-1}$，并且可以高效地计算出来。（欧几里得最大公约数算法的扩展在 $\lg (N)$ 的平方成本下完成。）因此，

$$M_{a^{-1}}{M}_a=M_{a^{-1}a}=M_1=I$$

因此，操作 $M_a$ 是确定性的和可逆的。这意味着它由一个置换矩阵描述，因此是酉的。

乘法运算的特征向量和特征值

现在让我们考虑操作 $M_a$ 的特征向量和特征值，假设 $a\in Z_N^{\ast }$。正如刚才所论证的，这一假设告诉我们 $M_a$ 是酉的。 $M_a$ 有 $N$ 个特征值，可能包括重复多次的相同特征值，通常在选择相应的特征向量时有一定的自由度——但我们不需要担心所有的可能性。让我们从简单的开始，先确定 $M_a$ 的一个特征向量。

$$|{\psi }_0⟩=\frac{1}{\sqrt{r}}(|1⟩+|a⟩+\cdots +|a^{r-1}⟩)$$

这里的 $r$ 是 $a$ 模 $N$ 的阶数——在本文的其余部分也是如此。与此特征向量相关的特征值是 1，因为当我们乘以 $a$ 时它不会改变。

$$M_a|{\psi }_0⟩=\frac{1}{\sqrt{r}}(|a⟩+\cdots +|a^{r-1}⟩+|a^r⟩)=\frac{1}{\sqrt{r}}(|a⟩+\cdots +|a^{r-1}⟩+|1⟩)=|{\psi }_0⟩$$

这是因为 $a^r=1$，所以每个标准基态 $|a^k⟩$ 被移到 $|a^{k+1}⟩$，对于 $k\le r-1$，并且 $|a^{r-1}⟩$ 被移回到 $|1⟩$。非正式地说，就像我们在慢慢搅拌 $|{\psi }_0⟩$，但它已经完全搅拌好了，所以没有任何变化。

这是 $M_a$ 的另一个特征向量。在阶数查找和相位估计的背景下，这个特征向量更有趣。

$$|{\psi }_1⟩=\frac{1}{\sqrt{r}}(|1⟩+{\omega }_r^{-1}|a⟩+\cdots +{\omega }_r^{-(r-1)}|a^{r-1}⟩)$$

或者，我们可以用求和来写这个向量，如下所示：

$$|{\psi }_1⟩=\frac{1}{\sqrt{r}}\sum _{k=0}^{r-1}{\omega }_r^{-k}|a^k⟩$$

这里我们看到复数 ${\omega }_r=e^{2\pi i/r}$ 自然地出现了，这是由于模 $N$ 乘法的底层结构。这次相应的特征值是 ${\omega }_r$。为了看到这一点，我们可以首先这样计算：

$$M_a|{\psi }_1⟩=\sum _{k=0}^{r-1}{\omega }_r^{-k}{M}_a|a^k⟩=\sum _{k=0}^{r-1}{\omega }_r^{-k}|a^{k+1}⟩=\sum _{k=1}^r{\omega }_r^{-(k-1)}|a^k⟩={\omega }_r\sum _{k=1}^r{\omega }_r^{-k}|a^k⟩$$

然后，因为 ${\omega }_r^{-r}=1={\omega }_r^0$ 并且 $|a^r⟩=|1⟩=|a^0⟩$，我们看到

$$\sum _{k=1}^r{\omega }_r^{-k}|a^k⟩=\sum _{k=0}^{r-1}{\omega }_r^{-k}|a^k⟩=|{\psi }_1⟩$$

所以 $M_a|{\psi }_1⟩={\omega }_r|{\psi }_1⟩$。

使用相同的推理，我们可以确定 $M_a$ 的其他特征向量/特征值对。实际上，对于任何选择的 $j\in \{0,\dots ,r-1\}$，我们有

$$|{\psi }_j⟩=\frac{1}{\sqrt{r}}\sum _{k=0}^{r-1}{\omega }_r^{-jk}|a^k⟩$$

是 $M_a$ 的一个特征向量，其相应的特征值是 ${\omega }_r^j$。

$$M_a|{\psi }_j⟩={\omega }_r^j|{\psi }_j⟩$$

还有其他特征向量，例如 $|0⟩$，其特征值为 1，但我们只关心我们刚刚确定的特征向量 $|{\psi }_0⟩,\dots ,|{\psi }_{r-1}⟩$。

应用相位估计

为了对给定的 $a\in Z_N^{\ast }$ 解决阶数查找问题，我们可以将相位估计过程应用于操作 $M_a$。为此，我们不仅需要用量子电路高效地实现 $M_a$，还需要实现 $M_a^2$、$M_a^4$、$M_a^8$ 等，直到获得足够精确的估计。这里我们将解释如何做到这一点，并稍后确定需要多少精度。

让我们从操作 $M_a$ 本身开始。自然地，因为我们使用量子电路模型，我们将使用二进制表示法来编码 $0$ 到 $N-1$ 之间的数字。我们需要编码的最大数字是 $N-1$，所以我们需要的位数是：

$$n=\lfloor {\log }_2(N-1)\rfloor +1$$

例如，如果 $N=21$，我们有 $n=5$。以下是 $Z_{20}$ 的元素作为长度为 5 的二进制字符串的编码：

$$\begin{array}{rl}0& \mapsto 00000\\ 1& \mapsto 00001\\ ⋮\\ 20& \mapsto 10100\end{array}$$

现在，以下是 $M_a$ 作为一个 $n$ 量子比特操作的精确定义：

$$M_a|x⟩=\{\begin{array}{ll}|ax\mathrm{mod}N⟩& \text{如果 }0\le x<N\\ |x⟩& \text{如果 }N\le x<2^n\end{array}$$

关键是，虽然我们只关心 $M_a$ 如何作用于 $|0⟩,\dots ,|N-1⟩$，但我们确实需要指定它如何作用于剩余的 $2^n-N$ 个标准基态——我们需要以仍然给我们一个酉操作的方式来做到这一点。定义 $M_a$ 使其对剩余的标准基态不做任何操作可以实现这一点。

使用上一课中讨论的整数乘法和除法算法，以及它们的可逆、无垃圾实现方法，我们可以构建一个量子电路来执行 $M_a$，对于任何选择的 $a\in Z_N^{\ast }$，其成本为 $O(n^2)$。以下是实现这一点的一种方法。

我们构建一个电路来执行操作：

$$|x⟩|y⟩\mapsto |x⟩|y\oplus f_a(x)⟩$$

其中

$$f_a(x)=\{\begin{array}{ll}ax\mathrm{mod}N& 0\le x<N\\ x& N\le x<2^n\end{array}$$

使用在 Qiskit 教程 中描述的方法。这给我们一个大小为 $O(n^2)$ 的电路。

1. 我们使用 $n$ 个交换门逐个交换量子比特来交换两个 $n$ 量子比特系统。
2. 类似于第一步，我们可以构建一个电路来执行操作：$$|x⟩|y⟩\mapsto |x⟩|y\oplus f_{a^{-1}}(x)⟩$$其中 $a^{-1}$ 是 $a$ 在 $Z_N^{\ast }$ 中的逆。

通过初始化底部的 $n$ 个量子比特并组合这三个步骤，我们得到以下变换：

$$\begin{array}{rl}|x⟩|0^n⟩& \stackrel{\text{step 1}}{\to }|x⟩|f_a(x)⟩\\ & \stackrel{\text{step 2}}{\to }|f_a(x)⟩|x⟩\\ & \stackrel{\text{step 3}}{\to }|f_a(x)⟩|x\oplus f_{a^{-1}}(f_a(x))⟩=|f_a(x)⟩|0^n⟩\end{array}$$

我们得到的电路的总成本是 $O(n^2)$。

为了执行 $M_a^2$、$M_a^4$、$M_a^8$ 等操作，我们可以使用完全相同的方法，只需将 $a$ 替换为 $a^2$、$a^4$、$a^8$ 等作为 $Z_N^{\ast }$ 的元素。也就是说，对于我们选择的任何幂 $k$，我们可以通过计算 $b=a^k\in Z_N^{\ast }$，然后使用 $M_b$ 的电路来创建 $M_a^k$ 的电路，而不是迭代 $M_a$ 的电路 $k$ 次。

计算 $a^k\in Z_N$ 的幂是上一课中提到的模幂运算问题。这个计算可以通过使用模幂运算法（在计算数论中通常称为幂算法）在经典计算机上完成。这次我们不需要用量子电路可逆地实现这个算法，只需要在经典计算机上完成它。

我们很幸运这是可能的。我们实际上将迭代 $M_a$ 大量次的问题（这可能是我们在相位估计中选择的数 $m$ 的指数级）卸载到一个高效的经典计算中。就我们运行的量子电路而言，$M_a$ 迭代 $k$ 次的成本只是 $M_b$ 的成本，对于 $b=a^k$——因此成本是 $O(n^2)$。

对于相位估计问题中任意选择的量子电路，这通常是不可能的，导致相位估计的成本是我们在相位估计中使用的数 $m$ 的指数级。通过计算数论的力量，在当前情况下的成本是 $m$ 的线性级。

方便的特征向量/特征值对

为了理解如何使用相位估计来解决阶数查找问题，让我们从假设我们在操作 $M_a$ 上运行相位估计过程并使用特征向量 $|{\psi }_1⟩$ 开始。事实证明，获得这个特征向量并不容易，所以这不会是故事的结尾——但从这里开始是有帮助的。 $M_a$ 对应于特征向量 $|{\psi }_1⟩$ 的特征值是

$${\omega }_r=e^{2\pi i\frac{1}{r}}$$

也就是说，${\omega }_r=e^{2\pi i\theta }$，其中 $\theta =\frac{1}{r}$。因此，如果我们在 $M_a$ 上运行相位估计过程并使用特征向量 $|{\psi }_1⟩$，我们将得到 $\frac{1}{r}$ 的近似值。通过计算倒数，我们将能够知道 $r$——前提是我们的近似值足够好。

更准确地说，当我们使用 $m$ 个控制量子比特运行相位估计过程时，我们得到一个数字 $y\in \{0,\dots ,2^m-1\}$，我们将 $y/2^m$ 作为 $\theta$ 的猜测，在这种情况下是 $\frac{1}{r}$。为了从这个近似值中找出 $r$，自然的做法是计算我们的近似值的倒数并四舍五入到最接近的整数。

$$\lfloor \frac{2^m}{y}+\frac{1}{2}\rfloor$$

例如，假设 $r=6$，我们使用特征向量 $|{\psi }_1⟩$ 在 $M_a$ 上执行相位估计，使用 $m=5$ 个控制比特。$\frac{1}{r}=\frac{1}{6}$ 的最佳 5 位近似值是 $\frac{5}{32}$，我们有很大的机会（在这种情况下约为 68%）从相位估计中得到结果 $y=5$。我们有

$$\frac{2^m}{y}=\frac{32}{5}=6.4$$

四舍五入到最接近的整数得到 6，这是正确答案。

另一方面，如果我们使用的精度不够，我们可能得不到正确答案。例如，如果我们在相位估计中使用 $m=4$ 个控制量子比特，我们可能得到 $\frac{1}{r}=\frac{1}{6}$ 的最佳 4 位近似值，即 $\frac{3}{16}$。取倒数得到

$$\frac{2^m}{y}=\frac{16}{3}=5.333\dots$$

四舍五入到最接近的整数得到错误答案 5。

我们需要多少精度才能得到正确答案？ 我们知道阶数 $r$ 是一个整数，直观地说，我们需要足够的精度来区分 $\frac{1}{r}$ 与附近的可能值，包括 $\frac{1}{r+1}$ 和 $\frac{1}{r-1}$。我们需要关注的最接近 $\frac{1}{r}$ 的数字是 $\frac{1}{r+1}$，这两个数字之间的距离是

$$\frac{1}{r}-\frac{1}{r+1}=\frac{1}{r(r+1)}$$

因此，如果我们想确保不将 $\frac{1}{r}$ 误认为 $\frac{1}{r+1}$，只需使用足够的精度来保证最佳近似 $y/2^m$ 到 $\frac{1}{r}$ 比到 $\frac{1}{r+1}$ 更接近。如果我们使用足够的精度，使得

$$|\frac{y}{2^m}-\frac{1}{r}|<\frac{1}{2r(r+1)}$$

这样误差小于 $\frac{1}{r}$ 和 $\frac{1}{r+1}$ 之间距离的一半，那么 $y/2^m$ 将比任何其他可能值（包括 $\frac{1}{r+1}$ 和 $\frac{1}{r-1}$）更接近 $\frac{1}{r}$。

我们可以如下双重检查这一点。假设

$$\frac{y}{2^m}=\frac{1}{r}+ϵ$$

对于满足

$$|ϵ|<\frac{1}{2r(r+1)}$$

的 $ϵ$。当我们取倒数时，我们得到

$$\frac{2^m}{y}=\frac{1}{\frac{1}{r}+ϵ}=\frac{r}{1+ϵr}\approx r-ϵr^2$$

通过在分子中取最大值并在分母中取最小值，我们可以如下界定我们与 $r$ 的距离。

$$\frac{ϵr^2}{1+ϵr}\le \frac{r^2}{2r(r+1)}\cdot \frac{1}{1-\frac{r}{2r(r+1)}}=\frac{r}{2r+1}<\frac{1}{2}$$

我们距离 $r$ 小于 $\frac{1}{2}$，因此如预期的那样，我们在四舍五入时会得到 $r$。

不幸的是，因为我们还不知道 $r$ 是多少，我们不能用它来告诉我们需要多少精度。我们可以做的是利用 $r$ 必须小于 $N$ 这一事实来确保我们使用足够的精度。特别地，如果我们使用足够的精度来保证最佳近似 $y/2^m$ 到 $\frac{1}{r}$ 满足

$$|\frac{y}{2^m}-\frac{1}{r}|\le \frac{1}{2N^2}$$

那么当我们取倒数时，我们将有足够的精度来正确确定 $r$。

取 $m=2\log (N)+1$ 可以确保我们有很大的机会使用前面描述的方法获得具有这种精度的估计。（如果我们对成功概率的下限为 40% 感到满意，那么取 $m=2\log (N)$ 就足够了。）

其他特征向量/特征值对

正如我们刚才所看到的，如果我们有 $M_a$ 的特征向量 $|{\psi }_1⟩$，我们将能够通过相位估计来了解 $r$，只要我们使用足够的控制量子比特来获得足够的精度来做到这一点。不幸的是，获得特征向量 $|{\psi }_1⟩$ 并不容易，所以我们需要找出如何继续。

假设我们像上面一样进行，只是用特征向量 $|{\psi }_k⟩$ 代替 $|{\psi }_1⟩$，对于我们选择考虑的任何 $k\in \{0,\dots ,r-1\}$。我们从相位估计过程中得到的结果将是一个近似值

$$\frac{y}{2^m}\approx \frac{k}{r}$$

在假设我们不知道 $k$ 或 $r$ 的情况下，这可能会或可能不会允许我们识别 $r$。例如，如果 $k=0$，我们将得到一个接近 0 的近似值 $y/2^m$，这不幸地告诉我们什么都没有。然而，这是一个不寻常的情况；对于其他 $k$ 值，我们至少能够了解一些关于 $r$ 的信息。

我们可以使用一种称为连分数算法的算法将我们的近似值 $y/2^m$ 转换为附近的分数——如果近似值足够好，包括 $k/r$。我们不会在这里解释连分数算法。相反，这里是关于该算法的一个已知事实的陈述。

给定一个整数 $N\ge 2$ 和一个实数 $\alpha \in (0,1)$，最多有一个整数对 $u,v\in \{0,\dots ,N-1\}$ 满足 $v\ne 0$ 且 $gcd(u,v)=1$，并且满足 $|\alpha -\frac{u}{v}|<\frac{1}{2N^2}$。给定 $\alpha$ 和 $N$，连分数算法找到 $u$ 和 $v$（或者报告它们不存在）。

该算法可以实现为大小为 $O((\log (N){)}^3)$ 的布尔电路。

如果我们有一个非常接近的近似值 $\frac{y}{2^m}$ 到 $\frac{r}{k}$，并且我们为 $N$ 和 $\alpha =\frac{y}{2^m}$ 运行连分数算法，我们将得到 $u$ 和 $v$，如事实中所述。仔细阅读该事实使我们得出结论：

$$\frac{u}{v}=\frac{k}{r}$$

因此，我们不一定能知道 $k$ 和 $r$，我们只知道最简形式的 $\frac{k}{r}$。

例如，正如我们已经注意到的，我们不会从 $k=0$ 中学到任何东西。但这是唯一一个发生这种情况的 $k$ 值。当 $k$ 非零时，它可能与 $r$ 有公因数——但我们从连分数算法中得到的数 $v$ 必须能整除 $r$。

这并不明显，但这是一个已知的事实，如果我们有能力学习 $u$ 和 $v$，对于 $u/v=k/r$，其中 $k\in \{0,\dots ,r-1\}$ 是均匀随机选择的，那么我们很可能在几次采样后恢复 $r$。特别地，如果我们对 $r$ 的猜测是我们观察到的所有 $v$ 值的最小公倍数，那么我们很可能是正确的。一些 $k$ 值不好，因为它们与 $r$ 共享公因数，当我们学习 $u$ 和 $v$ 时，这些公因数对我们是隐藏的。但随机选择的 $k$ 不太可能长时间隐藏 $r$ 的因数，并且我们猜测 $r$ 不正确的概率随着样本数量的增加而指数下降。

不使用特征向量继续

到目前为止，我们还没有解决如何获得 $M_a$ 的特征向量 $|{\psi }_k⟩$ 来运行相位估计过程的问题。事实证明，我们不需要创建它们。我们将做的是在状态 $|1⟩$ 上运行相位估计过程，这里的 $|1⟩$ 是数字 1 的 n 位二进制编码，代替 $M_a$ 的特征向量 $|\psi ⟩$。

到目前为止，我们已经讨论了在特定特征向量上运行相位估计过程，但没有什么能阻止我们在不是 $M_a$ 特征向量的输入状态上运行该过程，这就是我们在这里用状态 $|1⟩$ 所做的。（除非 $a=1$，否则 $a\in Z_N^{\ast }$ 的 $M_a$ 没有 $|1⟩$ 作为特征向量，这是我们将避免的 $a$ 的选择。）

以下方程有助于解释为什么我们选择状态 $|1⟩$ 代替特征向量：

$$|1⟩=\frac{1}{\sqrt{r}}\sum _{k=0}^{r-1}|{\psi }_k⟩$$

通过将右侧与标准基态的内积并使用前面提到的公式可以验证这一点。

更详细地说，假设我们用状态 $|1⟩$ 代替特征向量 $|{\psi }_k⟩$ 运行相位估计过程。在执行量子傅里叶变换后，这使我们得到状态：

$$\frac{1}{\sqrt{r}}\sum _{k=0}^{r-1}|{\psi }_k⟩|{\gamma }_k⟩$$

其中

$$|{\gamma }_k⟩=\frac{1}{2^m}\sum _{y=0}^{2^m-1}\sum _{x=0}^{2^m-1}{e}^{2\pi ix(k/r-y/2^m)}|y⟩$$

表示在执行量子傅里叶变换的逆变换后顶部 $m$ 个量子比特的状态。因此，当测量顶部 $m$ 个量子比特时，我们得到一个近似值 $y/2^m$ 到 $k/r$，其中 $k\in \{0,\dots ,r-1\}$ 是均匀随机选择的。

正如我们已经讨论过的，这使我们在几次独立运行后能够以很高的置信度了解 $r$，这是我们的目标。

总成本

实现每个受控幺正操作 $M_a^k$ 的成本是 $O(n^2)$。有 $m$ 个受控幺正操作，因此受控幺正操作的总成本是 $O(n^3)$。此外，我们有 $m$ 个 Hadamard 门（其成本为 $O(n)$），量子傅里叶变换的成本为 $O(n^2)$。因此，受控幺正操作的成本主导了整个过程的成本——因此总成本为 $O(n^3)$。

除了量子电路本身，还有一些需要执行的经典计算。这包括在 $Z_N$ 中计算幂 $a^k$，其中 $k=2,4,8,\dots ,2^{m-1}$，这些是创建受控幺正门所需的，以及将 $\theta$ 的近似值转换为分数的连分数算法。在这两种情况下，这些计算都可以通过成本为 $O(n^3)$ 的布尔电路来执行。

通常，所有这些界限都可以使用渐近快速算法来改进；这些界限假设我们使用的是基本算术运算的标准算法。

通过阶数查找进行因数分解

最后我们需要讨论的是如何通过解决阶数查找问题来帮助我们进行因数分解。这部分完全是经典的——它与量子计算没有特别的关系。

基本思路是，我们想要分解数字 $N$，并且我们可以递归地进行分解。具体来说，我们可以专注于分裂 $N$ 的任务，这意味着找到任意两个整数 $b,c\ge 2$ 使得 $N=bc$。如果 $N$ 是素数，这是不可能的，但我们可以首先使用素性测试算法高效地测试 $N$ 是否为素数，如果 $N$ 不是素数，我们将尝试分裂它。一旦我们分裂了 $N$，我们可以简单地对 $b$ 和 $c$ 进行递归，直到我们所有的因数都是素数，并且我们获得 $N$ 的素数分解。

分裂偶数很容易：我们只需输出 2 和 $N/2$。 分裂完全幂次也很容易，这意味着数字的形式为 $N=s^j$，其中整数 $s,j\ge 2$，只需近似根 $N^{1/2},N^{1/3},N^{1/4}$ 等，并检查附近的整数作为 $s$ 的嫌疑人。我们不需要超过 $\log (N)$ 步，因为此时根会下降到 2 以下，不会揭示其他候选者。

我们可以做这两件事是很好的，因为阶数查找对偶数或素数幂没有帮助，其中数字 $s$ 恰好是素数。 如果 $N$ 是奇数且不是素数幂，阶数查找允许我们分裂 $N$。

输入一个奇数、合数且不是素数幂的整数 $N$。迭代以下步骤：

1. 随机选择 $a\in \{2,\dots ,N-1\}$。
2. 计算 $d=gcd(a,N)$。
3. 如果 $d>1$，则输出 $b=d$ 和 $c=N/d$ 并停止。否则继续下一步，知道 $a\in Z_N^{\ast }$。
4. 设 $r$ 为 $a$ 模 $N$ 的阶数。（这里我们需要阶数查找。）
5. 如果 $r$ 是偶数：
   1. 计算 $x=a^{r/2}-1(\mathrm{mod}N)$
   2. 计算 $d=gcd(x,N)$。
   3. 如果 $d>1$，则输出 $b=d$ 和 $c=N/d$ 并停止。
6. 如果到达此点，迭代未能找到 $N$ 的因数。

该算法的迭代可能未能找到 $N$ 的因数。具体来说，这在两种情况下发生：

• $a$ 模 $N$ 的阶数是奇数。
• $a$ 模 $N$ 的阶数是偶数且 $gcd(a^{r/2}-1,N)=1$。

使用基本数论可以证明，对于随机选择的 $a$，至少有 1/2 的概率这两种情况都不会发生。事实上，对于 $N$ 的不同素因数的数量 $m$，概率最多为 $2^{-(m-1)}$。这就是为什么假设 $N$ 不是素数幂很重要。假设 $N$ 是奇数也是必要的，这就是为什么（简单的）$N$ 是偶数的情况必须单独处理。

因此，如果我们重复该过程 $t$ 次，每次随机选择 $a$，我们将以至少 $1-2^{-t}$ 的概率成功分裂 $N$。

我们不会详细分析这个过程，但这里是基本思路。如果我们选择的 $a$ 使得 $a$ 模 $N$ 的阶数 $r$ 是偶数，那么考虑数字 $a^{r/2}-1(\mathrm{mod}N)$ 和 $a^{r/2}+1(\mathrm{mod}N)$ 是有意义的。

使用公式 $Z^2-1=(Z+1)(Z-1)$，我们得出结论：

$$(a^{r/2}-1)(a^{r/2}+1)=a^r-1$$

根据阶数的定义，我们知道 $a^r\equiv 1(\mathrm{mod}N)$，这意味着 $N$ 整除 $a^r-1$。因此，$N$ 整除数字 $(a^{r/2}-1)(a^{r/2}+1)$，这意味着 $N$ 的每个素因数必须要么整除 $a^{r/2}-1$，要么整除 $a^{r/2}+1$（或两者都整除）。

对于随机选择的 $a$，我们很可能会有 $N$ 的素因数同时整除这两个项，这使我们能够通过计算 $gcd$ 来分裂 $N$。