量易简

秀尔算法

背景介绍

整数分解（Integer Factorization）指将一个正整数 $N$ 分解为若干更小（大于 $1$ 的）正整数的乘积，比如 $45 = 5 \times 9$ 。若这些更小的整数均为素数，这种分解又称为素因数分解（Prime Factorization），比如 $45 = 3 \times 3 \times 5 = 3^{2} \times 5$ 。由

定理 1（算术基本定理，Fundamental Theorem of Arithmetic）. 任何一个大于 $1$ 的正整数 $N$ 可以唯一分解成有限个素数的乘积 $N = p_{1}^{m_{1}} p_{2}^{m_{2}} \dots p_{s}^{m_{s}}$ ，这里 $p_{1} < p_{2} < \dots < p_{s}$ 均为素数，指数 $m_{i}$ 是正整数。

正整数的素因数分解唯一，继而正整数的素因数分解问题方才有意义。整数分解问题历史悠久，上述算数基本定理可考的证明最早可以追溯到古希腊数学家欧几里得的著作《几何原本》。

现代密码学中最重要的公开密钥密码体系（Public Key Cryptosystem）之一 RSA 便是基于素因数分解问题提出的，该公钥体系的安全性正是基于整数分解问题的困难性。目前已知的渐近运行时间最好的经典算法是普通数域筛选法（General Number Field Sieve $^{[1]}$ ），其运行时间是

\begin{matrix} exp ({(\frac{64 log N}{9})}^{\frac{1}{3}} (log log N)^{\frac{2}{3}}) . \\ (1.1) \end{matrix}

秀尔算法（Shor's Algorithm），是秀尔在 1994 年提出的量子算法 $^{[2]}$ ，可以在多项式时间

\begin{matrix} O ((log N)^{2} (log log N) (log log log N)) \\ (1.2) \end{matrix}

给出 $N$ 的一个整数分解。进而表明整数分解问题属于复杂度类 BQP，是量子计算机可以高效求解的一个问题。由于 RSA 密码体系的安全性依赖于整数分解问题的困难性，秀尔算法的提出使得这些密码体系在量子计算机面前不再安全。为了让更多的读者更好地了解量子算法特别是秀尔算法的魅力，我们将先介绍整体的算法流程与框架，再给出基于量易伏平台的具体用例，将更详细的技术细节后置于附录，以供想要深入研究的同学进一步学习。

章节目录

我们在第 2 段介绍秀尔算法的经典部分和一些需要引入的简单知识，并在算法 11 中给出完整的秀尔算法流程，大量的证明和细节留到第 6 段供感兴趣的读者深入阅读；
我们在第 3 段介绍秀尔算法的量子部分——量子求阶算法，部分证明也留到第 6 段再补充；
我们在第 4 段中使用量子编程平台量易伏给出两个例子——分解 $15$ 和 $63$ ，供读者自己体验一下量子编程的乐趣；
我们在第 6 段中补充定理的证明及算法的一些细节，不同部分需要的知识储备不尽相同，这里将其分为 5 个部分，供读者自行抉择：
- 6.1 子段补充量子求阶算法中 $U_{a, N}$ 门的本征值本征向量相关内容，需要一点点线性代数的知识背景；
- 6.2 子段补充量子求阶算法中连分数展开这一步骤的正确性的证明，不需要知识背景但证明略繁琐需要读者的耐心；
- 6.3 子段补充量子求阶算法成功概率相关内容，需要一点点初等数论的知识背景；
- 6.4 子段补充秀尔算法搜索 $a$ 具有“良性”的概率，更紧的下界需要一些抽象代数的知识背景；
- 6.5 子段补充 $U_{a, N}$ 门的分解，这是物理实现必要的步骤，这里给出了参考文献供读者延伸阅读。

算法总览

阅读本节后续内容需要一些基本的数学背景知识如下：

记整数集合为 $Z$ ，正整数集合为 $N^{*}$ ；
给定整数 $N, a$ ，若 $N / a$ 也是整数，称 $a$ 整除（Divide） $N$ ， $a$ 是 $N$ 的 因数（Factor），记为 $a | N$ ；否则称 $a$ 不整除 $N$ ，记为 $a ∤ N$ ； $N$ 的平凡因数指 $\pm 1$ 和 $\pm N$ ；
素数（Prime） 指没有非平凡因数的正整数，正 合数（Composite） 指既不是素数也不是 $1$ 的正整数；
偶数必然有素因数 $2$ ，负数的因数与其相反数无异；
两个整数的 最大公因数（Greatest Common Divisor） 定义为同时整除这两个整数的最大的正整数，记为 $g c d (\cdot, \cdot)$ ，可以使用欧几里得辗转相除法高效计算；称最大公因数为 $1$ 的两个整数互素；
给定整数 $a, b, N$ ，若 $N | (a - b)$ ，称 $a$ 与 $b$ 模 $N$ 同余，记为 $a \equiv b (mod N)$ ；
给定互素的整数 $a$ 和 $N$ ，若 $r$ 是使得 $a^{r} \equiv 1 (mod N)$ 的最小正整数，称 $r$ 为 $a$ 模 $N$ 的阶，记为 $r = o r d (a, N)$ 。

我们重申整数分解问题和素因数分解问题，并引入素数判定问题如下：

问题 2（整数分解）. 给定一个正整数 $N > 1$ ，若 $N$ 是合数，求 $N$ 的一个不平凡因数。

问题 3（素因数分解）. 给定一个正整数 $N > 1$ ，求 $N$ 的素因数分解。

问题 4（素数判定）. 给定一个正整数 $N > 1$ ，判断 $N$ 是否是素数。

关于素数判定问题，目前已知的最好的经典算法是确定性的多项式时间算法—— AKS 素数测试算法 $^{[3, 4]}$ 。

因为素数判定问题可以有效解决，所以这里关于整数分解和素因数分解的问题描述，便可改为更一般地表达，去掉输入为合数这一个条件。并且我们可以声明：素因数分解问题可以多项式归约到整数分解问题，也就是说如果存在求解整数分解问题的多项式时间算法，那么存在求解素因数分解问题的多项式时间算法：比如我们要求 $N$ 的素因数分解，通过整数分解算法找到了一个非平凡因子 $a$ ，那么 $N = a \times (\frac{N}{a})$ 就是 $N$ 的一个整数分解。我们分别对 $a$ 和 $\frac{N}{a}$ 递归地做整数分解算法，直到每个因数都是素数（通过素数判定算法检验），整个过程就是一个素因数分解算法。这个素因数分解算法的复杂度为： $($ 素数判定的复杂度 $+$ 整数分解的复杂度 $) \times$ 迭代次数，这里迭代次数 $< N$ 的各素因数重数和的两倍 $< 2 {log}_{2} N$ ，这里 $N$ 关于素因子 $p$ 的重数指使得 $p^{m} | N$ 的最大整数 $m$ 。故而上述归约是多项式归约。

故我们只需要考虑整数分解问题。

给定输入正整数 $N$ 后我们可以调用 AKS 素数测试算法，判断 $N$ 是否是素数，如果它是素数，返回 $\emptyset$ 表示 $N$ 没有非平凡因数。如果它不是素数，那么它一定是合数。有两类合数很容易地判断出来，同时也能找到其非平凡因数：偶数和整数幂。偶数的判定是简单的，我们用计算 $N / 2$ ，若它是整数， $N$ 便有非平凡因数 $2$ ；否则 $N$ 就不是偶数，即奇数。

整数幂是指形如 $a^{m}$ 的一类合数，其中 $a$ 和 $m$ 是大于 $1$ 的整数，例如 $3^{2}$ ， $4^{2}$ 和 $6^{6}$ 等等。我们可以用以下算法判断正合数 $N$ 是不是整数幂：

算法 5（整数幂判定）.

输入： 正合数 $N > 1$ 。

输出： [因数 $a$ , 对应的幂次 $m$ ]，若 $N = a^{m}$ 是整数幂；否则输出 False。

for m in range(2, int(math.log(N,2)) + 1):
    a = int(N ** (1 / m))
    if N == a ** m
        return [a,m]
return False

复制

这里 int(N ** (1 / m)) 可以用数值方法计算，确保多项式时间内可以完成整数幂的判定。特别地，当已知 $N$ 是奇数时，int(math.log(N,2)) 可以改为 int(math.log(N,3))。如果 $N$ 是整数幂，算法 5 返回的因数 $a$ 便是 $N$ 的一个不平凡因数。

经过素数判定、偶数判定和整数幂判定之后，我们便可假设正整数 $N$ 是奇合数且不是整数幂，也就是说 $N$ 是至少包含两个不同的奇素因数的奇合数（否则 $N$ 是奇数且只含有一种素因数可以推出 $N$ 是某个素数的幂，与 $N$ 不是整数幂矛盾）。至此才进入狭义的整数分解问题：

问题 6（狭义的整数分解）. 给定一个不是整数幂的正奇合数 $N$ ，求 $N$ 的一个不平凡因数。

对于这一问题有着非常多的经典求解办法，秀尔基于其中一种思路给出了量子版本的算法：基于量子求阶算法的整数分解算法。

我们再来回顾一下阶的概念：

定义 7（模 $N$ 的阶）. 给定整数 $a$ 和 $N$ ，若存在 $r$ 是使得 $a^{r} \equiv 1 (mod N)$ 的最小正整数，称 $r$ 为 $a$ 模 $N$ 的阶，记为 $r = o r d (a, N)$ ；否则称 $a$ 模 $N$ 的阶不存在。在不引发歧义的时候简称为 $a$ 的阶。

如果我们知道整数 $a$ 和其模 $N$ 的阶 $r$ ，便有 $N | (a^{r} - 1)$ 。若追加假设 $r$ 是偶数，那么由平方差公式有 $N | (a^{r / 2} - 1) (a^{r / 2} + 1)$ ，由阶定义中的最小性知 $N ∤ (a^{r / 2} - 1)$ 。若再追加假设 $N ∤ (a^{r / 2} + 1)$ 即 $a^{r / 2} \equiv̸ - 1 (mod N)$ ，有 $N ∤ (a^{r / 2} \pm 1)$ 但 $N ∣ (a^{r / 2} - 1) (a^{r / 2} + 1)$ 。故而 $g c d (N, a^{r / 2} \pm 1)$ 均为 $N$ 的不平凡因数。更详细准确的证明可以参照下述定理 10。

定义 8（“良性”）. 给定互素的整数 $a$ 和 $N$ ，设 $r = o r d (a, N)$ 是 $a$ 模 $N$ 的阶，若 $r$ 是偶数且 $a^{r / 2} \equiv̸ - 1 (mod N)$ ，我们称 $a$ 模 $N$ 是“良性”的；否则称为“非良性”的。

这里的“良性”是为了简化表述给出的称谓，无通用正式的术语。我们依旧在不引发歧义时略去“模 $N$ ”以简化表述。

我们已经知道“良性”的 $a$ 及其阶 $r$ 可以诱导出 $N$ 的非平凡因子，但如何找到“良性”的 $a$ 和如何计算 $a$ 的阶都是尚未解决的问题。但是如果我们有对任意的 $a$ 和 $N$ 可以计算阶 $o r d (a, N)$ 的算法，便可随机抽选整数 $a$ ，然后计算其阶 $r$ ，再判断 $a$ 是不是“良性”的，若其是“良性”的，便可返回 $g c d (N, a^{r / 2} \pm 1)$ 为 $N$ 的不平凡因数；否则重新抽选 $a$ ，并重复上述流程。特别地，如果我们抽选到的 $a$ 与 $N$ 不互素，虽然 $a$ 的阶不存在，但 $o r d (a, N)$ 已经是 $N$ 的一个非平凡因数了，此时也没必要进入求阶步骤了。

例 9. 考虑 $N = 315$ ， $a = 2$ ，可以验证 $o r d (2, 315) = 12$ 及 $2^{6} \equiv 64 \equiv̸ - 1 (mod 315)$ ，进而 $2$ 是“良性”的， $g c d (2^{6} - 1, 315) = 63$ 和 $g c d (2^{6} + 1, 315) = 5$ 均为 $315$ 的非平凡因数。特别地，恰有 $315 = 63 \times 5$ 。

在例 9 中我们指出若 $a$ 模 $N$ 是“良性”的，其便能诱导 $N$ 的一个不平凡的分解，这里更准确地表述成定理如下：

定理 10 给定正整数 $a, N$ ，设 $r = ord (a, N)$ 。若 $N$ 是奇数且 $a$ 模 $N$ 是“良性”的，则有 $N$ 的一个不平凡的分解： $N = gcd (a^{r / 2} - 1, N) \cdot gcd (a^{r / 2} + 1, N)$ 。

证明. 设 $b = a^{r / 2}$ ，易验证 $ord (b, N) = 2$ 及 $N | (b^{2} - 1) = (b - 1) (b + 1)$ 。设 $gcd (b - 1, N) = x_{1}$ ， $gcd (b + 1, N) = x_{2}$ ， $b - 1 = x_{1} y_{1}$ 及 $b + 1 = x_{2} y_{2}$ 。考虑到 $N$ 为奇数，进而有

\begin{matrix} gcd (x_{1}, x_{2}) = gcd (b + 1, b - 1, N) = gcd (b + 1, b - 1, 2, N) = 1, \\ (2.1) \end{matrix}

故而可设 $N = x_{1} x_{2} z$ ，这里 $gcd (\cdot, \dots, \cdot)$ 表示一组数的最大公因数。由

\begin{matrix} \begin{matrix} \begin{matrix} N | (b - 1) (b + 1) ⟺ x_{1} x_{2} z | x_{1} y_{1} x_{2} y_{2} ⟺ z | y_{1} y_{2} gcd (b - 1, N) = gcd (x_{1} y_{1}, x_{1} x_{2} z) = x_{1} ⟺ gcd (y_{1}, x_{2} z) = 1 ⟹ gcd (y_{1}, z) = 1 gcd (b + 1, N) = gcd (x_{2} y_{2}, x_{1} x_{2} z) = x_{2} ⟺ gcd (y_{2}, x_{1} z) = 1 ⟹ gcd (y_{2}, z) = 1 \end{matrix} ⎫ ⎪ ⎬ ⎪ ⎭ ⟹ z = 1, \end{matrix} \\ (2.2) \end{matrix}

此时 $N = x_{1} x_{2}$ 。由 $b \equiv̸ - 1 (mod N)$ 有 $x_{2} \neq N$ ；由 $ord (b, N) \neq 1$ 有 $x_{1} \neq N$ 。故而 $N = x_{1} x_{2}$ 是 $N$ 的一个不平凡的分解， $x_{1}$ 和 $x_{2}$ 都是 $N$ 的非平凡因数。 $□$

秀尔算法是对整数分解问题的一个经典求解算法的量子改进：

两者均将整数分解问题转化为求阶问题（Order-Finding）；
秀尔在求阶问题上获得了突破，提出了量子求阶算法（Quantum Order-Finding Algorithm）；
将量子求阶算法替换经典算法中求阶的步骤，便可得到量子版本的整数分解算法，即秀尔算法。

关于求阶这一步骤我们将在下一段中展开，这也是秀尔算法中唯一量子的部分。作为本段的最后，我们总结一下秀尔算法的流程如下：

算法 11（秀尔算法）.

输入： 正整数 $N$ 。

输出： $N$ 的一个非平凡的正因数，或 $\emptyset$ 表示 $N$ 没有非平凡的正因数。

若 $N = 1$ ，返回 $\emptyset$ ；
使用 AKS 素数测试算法判定 $N$ 是否为素数，若是返回 $\emptyset$ ；
若 $N$ 是偶数，返回 $2$ ；
用数值计算的方法依次判断 $\sqrt{N}, \sqrt[3]{N}, \dots, \sqrt[⌊ {log}_{3} (N) ⌋]{N}$ 是否为整数，若为整数则返回这个整数；
随机选取整数 $a$ 使得 $2 \leq a \leq N - 2$ ；
使用欧几里得辗转相除法计算最大公因数 $g c d (a, N)$ ；
若 $g c d (a, N) \neq 1$ ，返回 $g c d (a, N)$ ；
否则，使用（量子）求阶算法计算 $r = o r d (a, N)$ ；
若 $r$ 为奇数或者 $a^{r / 2} \equiv - 1 (mod N)$ ，转移到步骤 5；
否则，返回 $g c d (a^{r / 2} + 1, N)$ （或 $g c d (a^{r / 2} - 1, N)$ ）。

关于秀尔算法复杂度的分析需要 6.3 子段分析算法 11 步骤 8 的成功概率，以及 6.4 子段分析算法 11 步骤 9 发生循环的概率作为支撑。

量子求阶算法

因为算法旨在求解 $a$ 模 $N$ 的阶 $o r d (a, N)$ ，我们称其为量子求阶算法。考虑到

\begin{matrix} f (x) = a^{x} mod N \\ (3.1) \end{matrix}

为周期函数，其（最小正）周期恰是 $o r d (a, N)$ ，量子求阶算法可作为 量子周期求解算法（Period-Finding）的一个特例，这里 $x mod y := x - y ⌊ x / y ⌋ \in {0, \dots, y - 1}$ 是 $x$ 除以 $y$ 的余数。从本质上讲，量子相位估计算法也是量子周期求解算法的一个特例，但是从周期求解算法的方式来理解量子求阶算法和量子相位估计算法相对要困难一些，感兴趣的同学可以自行尝试了解 $^{[5]}$ 。

秀尔的量子求阶算法基于量子相位估计算法，将阶这一信息编码到某个量子门（或量子子电路）的本征相位上去，借由相位估计算法估计出该相位，进而可以恢复出阶的值。

将阶编码到本征相位

给定 $a \in {2, \dots, N - 1}$ 与 $N$ 互素，现要求解 $a$ 模 $N$ 的阶 $r = o r d (a, N)$ 。使用 $L = ⌈ {log}_{2} (N + 1) ⌉$ 个量子比特来编码，设有量子门

\begin{matrix} U_{a, N} = N - 1 \sum x = 0 | a x mod N ⟩ ⟨ x | + 2^{L} - 1 \sum x = N | x ⟩ ⟨ x | . \\ (3.2) \end{matrix}

习题 1. 验证 $U_{a, N}$ 为酉矩阵。

我们可以计算（证明参见附录 6.1 子段。）

引理 12. $U_{a, N}$ 有本征值形如

\begin{matrix} λ_{s} = exp (\frac{2 π i s}{r}), \\ (3.3) \end{matrix}

且有对应的本征态

\begin{matrix} | u_{s} ⟩ = \frac{1}{\sqrt{r}} r - 1 \sum k = 0 exp (\frac{- 2 π i s k}{r}) | a^{k} mod N ⟩ . \\ (3.4) \end{matrix}

为了让读者们更清晰的看到 $U_{a, N}$ 及其性质，我们给出如下例子：

例 13. 令 $N = 7$ ， $a = 2$ ，我们有 $r = o r d (a, N) = 3$ ， $L = ⌈ {log}_{2} (N + 1) ⌉ = 3$ ，

\begin{matrix} U_{2, 7} = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ \begin{matrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . \\ (3.5) \end{matrix}

$U_{2, 7}$ 部分本征值及本征态如下：

\begin{matrix} λ_{0} = 1, & | u_{0} ⟩ = \frac{1}{\sqrt{3}} (| 1 ⟩ & + & | 2 ⟩ + & | 4 ⟩); λ_{1} = exp (\frac{2 π i}{3}), & | u_{1} ⟩ = \frac{1}{\sqrt{3}} (| 1 ⟩ & + & exp (\frac{- 2 π i}{3}) & | 2 ⟩ + & exp (\frac{- 4 π i}{3}) & | 4 ⟩); λ_{2} = exp (\frac{4 π i}{3}), & | u_{2} ⟩ = \frac{1}{\sqrt{3}} (| 1 ⟩ & + & exp (\frac{- 4 π i}{3}) & | 2 ⟩ + & exp (\frac{- 8 π i}{3}) & | 4 ⟩) . \\ (3.6) \end{matrix}

当我们对算符 $U_{a, N}$ 及其本征态 $| u_{s} ⟩$ 执行量子相位估计算法时便可得到相位参数 $s / r$ 的一个估计，借助连分数展开算法便有概率恢复出其分母 $r$ 。

考虑到在 $r$ 未知的前提下，我们无法制备量子态 $| u_{s} ⟩$ ，但我们有如下恒等式（证明参见附录 6.1 子段）：

引理 14.

\begin{matrix} \frac{1}{\sqrt{r}} r - 1 \sum s = 0 | u_{s} ⟩ = | 1 ⟩ . \\ (3.7) \end{matrix}

也就是说，计算基中的 $| 1 ⟩$ 可以看做本征态 ${| u_{s} ⟩}$ 的一个同概率幅叠加。

当我们对算符 $U_{a, N}$ 及量子态 $| 1 ⟩$ 执行量子相位估计算法时，可以得到相位参数 ${s / r | s = 0, 1, \dots, r - 1}$ 的估计值的一个叠加态，记为混态 $ρ = {1 / r, | 2^{t} ˜ s / r ⟩}$ （相位估计算法输出的是两个寄存器，当我们只关心辅助寄存器时，它的状态就是一个混态）。关于混态的系综表示可以参考 1.4 节，这里可以简单理解为对该叠加态的辅助寄存器进行测量时如同先以概率 $1 / r$ 抽选到量子态 $| 2^{t} ˜ s / r ⟩$ ，再对 $| 2^{t} ˜ s / r ⟩$ 进行测量。这里相位估计算法的精度取 $2^{- 2 L - 1}$ ，如果在辅助寄存器使用 $t = 2 L + 1 + ⌈ {log}_{2} (1 + \frac{2}{ϵ π^{2}}) ⌉$ 个量子比特，相位估计算法的成功概率不少于 $1 - ϵ$ （符号和理论的正确性参照 4.3 节）。

连分数展开解码相位为阶

当我们对辅助寄存器 $ρ$ 用计算基进行测量后，设测量结果为 $| m ⟩$ ，我们有至少 $1 - ϵ$ 的把握认为存在 $s$ 使得 $2^{- t} m$ 是 $s / r$ 的一个 $2^{- 2 L - 1}$ -估计。对 $2^{- t} m$ 做连分数展开，并依次计算其截断 $p / q$ 的分母 $q$ ，找到满足条件 $q < N$ 的最大分母 $q$ 。

这里举例演示一下连分数展开算法，更多细节参考篇末的附录 6.2 子段：

例 15. 考虑 $N = 63$ 时，对 $2^{- t} m = 2729 / 8192$ 执行连分数展开：

\begin{matrix} \frac{2729}{8192} = & \frac{1}{\frac{8192}{2729}} = \frac{1}{3 + \frac{5}{2729}} = \frac{1}{3 + \frac{1}{545 + \frac{4}{5}}} = & 0 + \frac{1}{3 + \frac{1}{545 + \frac{1}{1 + \frac{1}{4}}}} =: [0, 3, 545, 1, 4], \\ (3.8) \end{matrix}

可以看出连分数展开的过程就是欧几里得辗转相除法，每一次迭代将分数表示为带分数的形式，然后将带分数的分数部分的分子化为 $1$ ，将新的分母作为下一次迭代的输入，循环至带分数为整数止。分别称式 $(3.8)$ 中的 $[0] = 0$ ， $[0, 3] = 0 + \frac{1}{3} = \frac{1}{3}$ ， $[0, 3, 545] = 0 + \frac{1}{3 + \frac{1}{545}} = \frac{545}{1636}$ 和 $[0, 3, 545, 1] = 0 + \frac{1}{3 + \frac{1}{545 + \frac{1}{1}}} = \frac{546}{1639}$ 为连分数 $[0, 3, 545, 1, 4]$ 的第 $0$ ，第 $1$ ，第 $2$ 和第 $3$ 截断。考虑其中满足条件分母 $< 63$ 的最大分母为 $3$ 。

注 16. 连分数展开截断总能恢复出 $s / r$ 的分母，但 $s / r$ 的分母不一定是 $r$ ：当 $g c d (s, r) \neq 1$ 时， $s / r$ 的分母是 $r / g c d (s, r)$ 。 $q = r$ 的概率即 $g c d (s, r) = 1$ 的概率 $φ (r) / r$ ，这里 $q$ 是连分数展开各截断的分母中满足 $< N$ 条件的最大值， $φ (r)$ 是欧拉函数（不超过 $r$ 且与 $r$ 互素的正整数数量，参见附录 6.1 子段定义 28）。

这里可以恢复出 $r$ 的概率 $φ (r) / r$ 没有良好的下界，也就是说存在极端的情况 $φ (r) / r$ 趋近于 $0$ ，为了提高恢复出 $r$ 的概率，目前最好的办法是执行两次求阶算法并计算分母的最小公倍数：

定理 17. 通过执行两次量子求阶算法，可以得到两个分母 $q$ 和 $q^{'}$ ，计算 $q$ 和 $q^{'}$ 的最小公倍数 $~ r := q q^{'} / gcd (q, q^{'})$ ，则 $~ r = r$ 的概率 $\geq 6 (1 - ϵ)^{2} / π^{2} \geq 60.7 % \cdot (1 - ϵ)^{2}$ ，这里 $ϵ$ 是相位估计算法的失败概率。

本定理的证明参阅附录 6.3 子段。

总结

本段的最后，我们再总结一下量子求阶算法，并将在下段给出两个具体的例子。

算法 18（量子求阶算法）.

输入： 整数 $a \in {2, \dots, N - 1}$ 与 $N$ 互素，相位估计算法的失败概率 $ϵ$ 。

输出： $o r d (a, N)$ 。

令 $L = ⌈ {log}_{2} (N + 1) ⌉$ ；
执行两次：使用 $t = 2 L + 1 + ⌈ {log}_{2} (1 + \frac{2}{ϵ π^{2}}) ⌉$ 个量子比特作为辅助寄存器对量子门 $U_{a, N}$ 和量子态 $| 1 ⟩$ 执行量子相位估计算法并对辅助寄存器进行测量。分别得到 $| m ⟩$ 和 $| m^{'} ⟩$ ；
分别对 $2^{- t} m$ 和 $2^{- t} m^{'}$ 做连分数展开，分别计算它们各个截断的分母中满足 $< N$ 条件的最大整数 $q_{j}$ 和 $q_{k}^{'}$ ；
计算 $r_{j k} = q_{j} q_{k}^{'} / g c d (q_{j}, q_{k}^{'})$ ，若 $a^{r_{j k}} \equiv 1 (mod N)$ 返回 $o r d (a, N) = r_{j k}$ ，否则转移到步骤 2 重新计算。

算法 18 步骤 4 的成功率不小于 $60.7 % \cdot (1 - ϵ)^{2}$ 。特别地，取 $t = 2 L + 1 + 3$ 时， $ϵ = 2.895 %$ ，算法 18 步骤 4 的成功概率不小于 $57.3 %$ ，此时执行至多 $4$ 次步骤 4 可以使计算出 $o r d (a, N)$ 的概率不少于 $96.6 %$ 。关于量子门 $U_{a, N}$ 的制备可以参阅附录 6.5 子段；更多的细节以及证明也请参考篇末的附录，其需要更多的数学背景知识。

简例与量易伏平台演示

在本小节中，我们使用量易伏平台来演示秀尔算法。我们首先使用 QComposer 分解一个简单的例子—— $15 = 3 \times 5$ ，并给出对应电路图。然后考虑一个更复杂的例子—— $63 = 3^{2} \times 7$ ，并给出基于混合语言编程的代码。

分解 $15$

依次做素数判定、偶数判定和整数幂判定，知 $15$ 满足狭义的整数分解的输入条件。随机选取 $a \in [2, 13]$ 得到 $a = 2$ （我们这里先随机拿到了 $9$ 和 $5$ ，为了更好地演示算法重新随机得到了这个值）。验证知 $g c d (2, 15) = 1$ ，进入量子求阶算法求 $o r d (2, 15)$ 。（ $g c d (9, 15) = 3$ ， $g c d (5, 15) = 5$ 均可经典地得到 $15$ 的分解。）

我们需要 $L = ⌈ {log}_{2} (15 + 1) ⌉ = 4$ 个量子比特来编码算符 $U_{2, 15}$ （本子段中也简记为 $U$ ）：

\begin{matrix} U_{2, 15} = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ \begin{matrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . \\ (4.1) \end{matrix}

我们可以用如下电路实现 $C (U_{2, 15})$ ，第 $0$ 比特控制作用在第 $1, 2, 3, 4$ 比特上的一个 $U_{2, 15}$ 门的一种电路分解如下图 1 所示：

在量易伏 QComposer 上实现 \operatorname{controlled^0-}U_{2,15}^{1,2,3,4} — 图 1：在量易伏 QComposer 上实现 $c o n t r o l l e d^{0} - U_{2, 15}^{1, 2, 3, 4}$

对应 OPENQASM 代码如下：

OPENQASM 2.0;
include "qelib1.inc";
qreg q[6];
creg c[6];
ccx q[0], q[4], q[1];
ccx q[0], q[4], q[2];
ccx q[0], q[4], q[3];
ccx q[0], q[3], q[4];
ccx q[0], q[4], q[3];
ccx q[0], q[2], q[3];
ccx q[0], q[3], q[2];
ccx q[0], q[1], q[2];
ccx q[0], q[2], q[1];

复制

相位估计算法中还需要用到其它 $C (U_{2, 15}^{2^{j}})$ 门:

\begin{matrix} C (U_{2, 15}^{2}) & = {CSWAP}_{0, 1, 3} \cdot {CSWAP}_{0, 2, 4}; C (U_{2, 15}^{2^{j}}) & = I, for j \geq 3, \\ (4.2) (4.3) \end{matrix}

这里 ${CSWAP}_{j, k, l}$ 为第 $j$ 比特为控制比特、以第 $k$ 和 $l$ 比特为受控比特的 $CSWAP$ 门。使用 $2 L + 1 = 9$ 个量子比特估计量子态 $| 0001 ⟩$ 的相位，QComposer 中的电路如下图 2 所示（这里略去了可以忽略的 $I$ 门）：

在量易伏 QComposer 上求解 \textrm{Ord}(2,15)" — 图 2：在量易伏 QComposer 上求解 $Ord (2, 15)$

其对应的 OPENQASM 代码如下：

OPENQASM 2.0;
include "qelib1.inc";
qreg q[13]; // q[0],...,q[8] 为辅助寄存器，q[9],...,q[12] 为工作寄存器
creg c[13];
x q[12]; // 制备初始状态|0>|1>，下面开始执行相位估计算法
h q[0];
h q[1];
h q[2];
h q[3];
h q[4];
h q[5];
h q[6];
h q[7];
h q[8]; // 相位估计的辅助步，对辅助寄存器取平均叠加完毕
ccx q[8], q[12], q[9];
ccx q[8], q[12], q[10];
ccx q[8], q[12], q[11];
ccx q[8], q[11], q[12];
ccx q[8], q[12], q[11];
ccx q[8], q[10], q[11];
ccx q[8], q[11], q[10];
ccx q[8], q[9], q[10];
ccx q[8], q[10], q[9]; // 上一注释至此实现了 C(U) 门，q[8] 控制工作寄存器
cswap q[7], q[9], q[11];
cswap q[7], q[10], q[12]; //上一注释至此实现了 C(U^2) 门，q[7] 控制工作寄存器
// 略去为 I 的其他 C(U^(2^j)) 门，准备开始对辅助寄存器做逆傅里叶变换
swap q[0], q[8];
swap q[1], q[7];
swap q[2], q[6];
swap q[3], q[5];
h q[8];
cu(0, 0, -pi/2) q[8], q[7];
h q[7];
cu(0, 0, -pi/4) q[8], q[6];
cu(0, 0, -pi/2) q[7], q[6];
h q[6];
cu(0, 0, -pi/8) q[8], q[5];
cu(0, 0, -pi/4) q[7], q[5];
cu(0, 0, -pi/2) q[6], q[5];
h q[5];
cu(0, 0, -pi/16) q[8], q[4];
cu(0, 0, -pi/8) q[7], q[4];
cu(0, 0, -pi/4) q[6], q[4];
cu(0, 0, -pi/2) q[5], q[4];
h q[4];
cu(0, 0, -pi/32) q[8], q[3];
cu(0, 0, -pi/16) q[7], q[3];
cu(0, 0, -pi/8) q[6], q[3];
cu(0, 0, -pi/4) q[5], q[3];
cu(0, 0, -pi/2) q[4], q[3];
h q[3];
cu(0, 0, -pi/64) q[8], q[2];
cu(0, 0, -pi/32) q[7], q[2];
cu(0, 0, -pi/16) q[6], q[2];
cu(0, 0, -pi/8) q[5], q[2];
cu(0, 0, -pi/4) q[4], q[2];
cu(0, 0, -pi/2) q[3], q[2];
h q[2];
cu(0, 0, -pi/128) q[8], q[1];
cu(0, 0, -pi/64) q[7], q[1];
cu(0, 0, -pi/32) q[6], q[1];
cu(0, 0, -pi/16) q[5], q[1];
cu(0, 0, -pi/8) q[4], q[1];
cu(0, 0, -pi/4) q[3], q[1];
cu(0, 0, -pi/2) q[2], q[1];
h q[1];
cu(0, 0, -pi/256) q[8], q[0];
cu(0, 0, -pi/128) q[7], q[0];
cu(0, 0, -pi/64) q[6], q[0];
cu(0, 0, -pi/32) q[5], q[0];
cu(0, 0, -pi/16) q[4], q[0];
cu(0, 0, -pi/8) q[3], q[0];
cu(0, 0, -pi/4) q[2], q[0];
cu(0, 0, -pi/2) q[1], q[0];
h q[0]; // 逆傅里叶变换结束，相位估计算法进入测量阶段
measure q[8] -> c[8];
measure q[7] -> c[7];
measure q[6] -> c[6];
measure q[5] -> c[5];
measure q[4] -> c[4];
measure q[3] -> c[3];
measure q[2] -> c[2];
measure q[1] -> c[1];
measure q[0] -> c[0];

复制

某次运行结果如下：

{
    "000000011": 237,
    "000000010": 257,
    "000000000": 269,
    "000000001": 261
}

复制

测量结果大约为等频率的 $| 000000000 ⟩$ ， $| 010000000 ⟩$ ， $| 100000000 ⟩$ 和 $| 110000000 ⟩$ ，分别对应真分数 $0$ ， $1 / 4$ ， $1 / 2$ 和 $3 / 4$ 。这里仅 $3 / 4$ 可以做连分数展开：

\begin{matrix} \frac{3}{4} = \frac{1}{1 + \frac{1}{3}}, \\ (4.4) \end{matrix}

其第 $1$ 截断为 $1$ 。

故而每次相位估计算法有约 $1 / 2$ 的概率可以得到正确结果 $o r d (2, 15) = 4$ ，两次实验有约 $3 / 4$ 的概率可以得到正确结果。由 $2^{4 / 2} = 4 \equiv̸ - 1 (mod 15)$ 知 $15 | (4 - 1) (4 + 1)$ 诱导了 $15$ 的一个不平凡的分解：

\begin{matrix} 15 = gcd (15, 4 - 1) \times gcd (15, 4 + 1) = 3 \times 5. \\ (4.5) \end{matrix}

至此完成了秀尔算法分解 $15$ 的演示。

后话：代码演示中“略去为 $I$ 的其他 $C (U^{2^{j}})$ 门”有循环论证之嫌，正确的演示应当将这些门逐个分解，而不是利用理论证明的方法计算得出 $C (U^{2^{j}}) = I$ ，在下一子段中，我们在演示 $63$ 的分解时避开了这一问题。

分解 $63$

在本子段中，我们基于混合语言编程给出秀尔算法的演示代码，推荐读者在本地使用 QCompute SDK 进行模拟实验。

选取 $a = 2$ ，验证知 $g c d (2, 63) = 1$ ，进入量子求阶算法求 $o r d (2, 63)$ 。我们需要 $L = ⌈ {log}_{2} (63 + 1) ⌉ = 6$ 个量子比特来编码算符 $U_{2, 63}$ （本子段中也简记为 $U$ ），使用 $2 L + 1 = 13$ 个辅助量子比特来估计相位。考虑到 $j > 0$ 时，

\begin{matrix} C (U^{2^{j}}) = C (U^{2^{j m o d 2}}), \\ (4.6) \end{matrix}

这里 $j mod 2$ 表示 $j$ 除以 $2$ 的余数 $\in {0, 1}$ 。故而我们只需要分解门 $C (U)$ ， $C (U^{2})$ 和 $C (U^{4})$ （细节详见下述代码）。量易伏中的代码如下：

import numpy as np
from QCompute import *


# 初始化量子环境，选择使用本地自研模拟器
# 注：量易伏也支持云端模拟器和真机后端，访问 https://quantum-hub.baidu.com/ 获取更多信息。
env = QEnv()
env.backend(BackendName.LocalBaiduSim2)

pi = np.pi
L = 6  # 编码 U 门需要的量子比特数、工作寄存器的量子比特数
N = 3 * L + 1  # 算法需要使用的总量子比特数
t = 2 * L + 1  # 相位估计算法所需的辅助比特数、辅助寄存器的量子比特数

q = [env.Q[i] for i in range(0, N)]  # 量子比特的列表

X(q[N - 1])  # 在工作寄存器制备初始状态 |1>，下面开始执行相位估计算法

# 相位估计的第一步，对辅助寄存器取平均叠加完毕，下面开始将相位估计算法的转存步骤
for i in range(0, t):
    H(q[i])

# 实现门 C(U)，辅助寄存器的最后一个量子比特控制工作寄存器
CCX(q[2 * L], q[t + 5], q[t + 0])
CCX(q[2 * L], q[t + 5], q[t + 1])
CCX(q[2 * L], q[t + 5], q[t + 2])
CCX(q[2 * L], q[t + 5], q[t + 4])
CCX(q[2 * L], q[t + 5], q[t + 3])
CCX(q[2 * L], q[t + 4], q[t + 5])
CCX(q[2 * L], q[t + 5], q[t + 4])
CCX(q[2 * L], q[t + 3], q[t + 4])
CCX(q[2 * L], q[t + 4], q[t + 3])
CCX(q[2 * L], q[t + 2], q[t + 3])
CCX(q[2 * L], q[t + 3], q[t + 2])
CCX(q[2 * L], q[t + 1], q[t + 2])
CCX(q[2 * L], q[t + 2], q[t + 1])
CCX(q[2 * L], q[t + 0], q[t + 1])
CCX(q[2 * L], q[t + 1], q[t + 0])

# 依次作用其他 C(U^(2^j)) 门
s = 2 * L - 1
while s >= 0:
    if s % 2 == 1:
        CCX(q[s], q[t + 5], q[t + 1])
        CCX(q[s], q[t + 5], q[t + 3])
        CCX(q[s], q[t + 4], q[t + 0])
        CCX(q[s], q[t + 4], q[t + 2])
        CCX(q[s], q[t + 3], q[t + 5])
        CCX(q[s], q[t + 5], q[t + 3])
        CCX(q[s], q[t + 2], q[t + 4])
        CCX(q[s], q[t + 4], q[t + 2])
        CCX(q[s], q[t + 1], q[t + 3])
        CCX(q[s], q[t + 3], q[t + 1])
        CCX(q[s], q[t + 0], q[t + 2])
        CCX(q[s], q[t + 2], q[t + 0])  # 本条件下的门组合后即门 C(U^2)
    else:
        CCX(q[s], q[t + 5], q[t + 1])
        CCX(q[s], q[t + 5], q[t + 3])
        CCX(q[s], q[t + 4], q[t + 0])
        CCX(q[s], q[t + 4], q[t + 2])
        CCX(q[s], q[t + 3], q[t + 5])
        CCX(q[s], q[t + 3], q[t + 1])
        CCX(q[s], q[t + 5], q[t + 3])
        CCX(q[s], q[t + 2], q[t + 4])
        CCX(q[s], q[t + 2], q[t + 0])
        CCX(q[s], q[t + 4], q[t + 2])
        CCX(q[s], q[t + 1], q[t + 5])
        CCX(q[s], q[t + 5], q[t + 1])
        CCX(q[s], q[t + 0], q[t + 4])
        CCX(q[s], q[t + 4], q[t + 0])  # 本条件下的门组合后即门 C(U^4)

# 逆傅里叶变换中的 swap 步骤
for i in range(0, t // 2):
    SWAP(q[i], q[t - i - 1])

# 逆傅里叶变换中的控制旋转步骤
for i in range(0, t - 1):
    H(q[t - i - 1])
    for j in range(0, i + 1):
        CRZ(-pi / 2 ** (i - j + 1))(q[t - j - 1], q[t - i - 2])
        RZ(-pi / 2 ** (i - j + 2))(q[t - j - 1])
H(q[0])

# 完成逆傅里叶变换，完成相位估计算法，等待测量
MeasureZ(q[:L], range(0, L))  # 只测量辅助寄存器即前 L 个比特
env.commit(8192, downloadResult=False)

复制

执行算法 $8192$ 次，忽略总频次 $5 %$ 的一些小概率结果，其余结果经数据处理后如下，其中 "4096": 1419 表示 $| 4096 ⟩$ 出现了 $1419$ 次：

{
    "4096": 1419,
    "0": 1287,
    "5461": 955,
    "2731": 955,
    "1365": 940,
    "6827": 896,
    "6826": 252,
    "1366": 250,
    "2730": 246,
    "5462": 219,
    "2732": 56,
    "5460": 53,
    "6828": 52,
    "1364": 49,
    "2729": 47,
    "1367": 41,
    "5463": 36,
    "6825": 34,
}

复制

分别对这些结果对应的分数做连分数展开并截断，得：

{
    {{0, 1/2}: 1419},
    {{0}: 1287},
    {{0, 1, 1/2, 2/3, 5461/8192}: 955},
    {{0, 1/2, 1/3, 2731/8192}: 955},
    {{0, 1/6, 682/4093, 1365/8192}: 940}, 
    {{0, 1, 5/6, 3411/4093, 6827/8192}: 896}, 
    {{0, 1, 4/5, 5/6, 1704/2045, 3413/4096}: 252}, 
    {{0, 1/5, 1/6, 341/2045, 683/4096}: 250}, 
    {{0, 1/3, 1365/4096}: 246}, 
    {{0, 1, 2/3, 2731/4096}: 219}, 
    {{0, 1/2, 1/3, 683/2048}: 56}, 
    {{0, 1, 1/2, 2/3, 1365/2048}: 53}, 
    {{0, 1, 5/6, 851/1021, 1707/2048}: 52},
    {{0, 1/6, 170/1021, 341/2048}: 49}, 
    {{0, 1/3, 545/1636, 546/1639, 2729/8192}: 47}, 
    {{0, 1/5, 1/6, 136/815, 137/821, 410/2457, 1367/8192}: 41},
    {{0, 1, 2/3, 1091/1636, 1093/1639, 5463/8192}: 36},
    {{0, 1, 4/5, 5/6, 679/815, 684/821, 2047/2457, 6825/8192}: 34}
}

复制

考虑各行第一个分量中分母不超过 $N = 63$ 且分母最大的分数的分母（这里仅有 $1, 2, 3, 6$ ），其中可以得到的分母为 $2$ ， $3$ 和 $6$ 的频次分别为 $1419$ ， $2567$ 和 $2514$ 。考虑到 $o r d (2, 63) = 6$ ，故而每次相位估计算法有约 $2514 / 8192 \approx 30.7 %$ 的概率可以得到正确结果 $o r d (2, 63) = 6$ ，两次实验有约

\begin{matrix} 1 - {(1 - \frac{2514}{8192})}^{2} + 2 \times \frac{1419}{8192} \times \frac{2567}{8192} \approx 62.8 % \\ (4.7) \end{matrix}

的概率可以得到正确结果 $o r d (2, 63) = 6$ ，这里最后一个加项表示两次测量中一次结果为 $2$ 而另一次为 $3$ 。由 $2^{6 / 2} = 8 \equiv̸ - 1 (mod 63)$ 知 $63 | (8 - 1) (8 + 1)$ 诱导了 $63$ 的一个不平凡的分解：

\begin{matrix} 63 = gcd (63, 8 - 1) \times gcd (63, 8 + 1) = 7 \times 9. \\ (4.8) \end{matrix}

参考资料

[1] Pomerance, Carl. "A tale of two sieves." Notices Amer. Math. Soc., 1996.

[2] Shor, Peter W. "Algorithms for quantum computation: discrete logarithms and factoring." Proceedings 35th Annual Symposium on Foundations of Computer Science, IEEE, 1994.

[3] Agrawal, Manindra, Neeraj Kayal, and Nitin Saxena. "PRIMES is in P." Annals of Mathematics (2004): 781-793.

[4] Lenstra Jr, H. W., and Carl Pomerance. "Primality testing with Gaussian periods." Archived version 20110412, Dartmouth College, US, 2011.

[5] Nielsen, Michael A., and Isaac L. Chuang. "Quantum Computation and Quantum Information: 10th Anniversary Edition." Cambridge University Press, 2010.

[6] Vedral, Vlatko, Adriano Barenco, and Artur Ekert. "Quantum networks for elementary arithmetic operations." Physical Review A 54.1(1996): 147.

[7] Beauregard, Stephane. "Circuit for Shor's algorithm using 2n+3 qubits." arXiv preprint , 2002.

附录

附录中会补充一些定理的证明，以及更多的细节，建议打开两个窗口平行阅读。附录中有些公式会有两个编号，此时第一个编号是该公式上一次在本节中出现时的编号。

本段目录重申如下：

6.1 子段补充量子求阶算法中 $U_{a, N}$ 门的本征值本征向量相关内容，需要一点点线性代数的知识背景
6.2 子段补充量子求阶算法中连分数展开这一步骤的正确性的证明，不需要知识背景，但证明略繁琐需要读者的耐心；
6.2 子段补充量子求阶算法成功概率相关内容，需要一点点初等数论的知识背景；
6.4 子段补充秀尔算法搜索 $a$ 具有“良性”的概率，更紧的下界需要一些抽象代数的知识背景；
6.5 子段补充 $U_{a, N}$ 门的分解，这是物理实现必要的步骤，这里给出了参考文献供读者延伸阅读，另外读者也可以参考 QCompute SDK 中的秀尔算法代码样例和教程。

$U_{a, N}$ 门的本征

本子段旨在补全引理 12 和引理 14 的证明，需要的知识仅有线性代数，不熟练 Dirac 符号的读者可以再回顾下 1.2 节。这里先重申一下 $U_{a, N}$ 门的定义：

\begin{matrix} U_{a, N} = N - 1 \sum x = 0 | a x mod N ⟩ ⟨ x | + 2^{L} - 1 \sum x = N | x ⟩ ⟨ x | . \\ (3.2; 6.1) \end{matrix}

引理 12. $U_{a, N}$ 有本征值形如

\begin{matrix} λ_{s} = exp (\frac{2 π i s}{r}), \\ (3.3; 6.2) \end{matrix}

且有对应本征态

\begin{matrix} | u_{s} ⟩ = \frac{1}{\sqrt{r}} r - 1 \sum k = 0 exp (\frac{- 2 π i s k}{r}) | a^{k} mod N ⟩ . \\ (3.4; 6.3) \end{matrix}

证明. 可以验证

\begin{matrix} U_{a, N}^{r} = N - 1 \sum x = 0 | a^{r} x mod N ⟩ ⟨ x | + 2^{L} - 1 \sum x = N | x ⟩ ⟨ x | = I_{2^{L}}, \\ (6.4) \end{matrix}

进而有 $U_{a, N}$ 的任意本征值 $λ$ 满足 $λ^{r} = 1$ ，故而存在 $s \in {0, 1, \dots, r - 1}$ ，使得 $λ = exp (\frac{2 π i s}{r})$ 。设 $| u_{s} ⟩$ 如式 $(6.3)$ ，考虑到

\begin{matrix} U_{a, N} | u_{s} ⟩ = & ⎛ ⎝ N - 1 \sum x = 0 | a x mod N ⟩ ⟨ x | + 2^{L} - 1 \sum x = N | x ⟩ ⟨ x | ⎞ ⎠ \cdot \frac{1}{\sqrt{r}} r - 1 \sum k = 0 exp (\frac{- 2 π i s k}{r}) | a^{k} mod N ⟩ = & \frac{1}{\sqrt{r}} N - 1 \sum x = 0 r - 1 \sum k = 0 exp (\frac{- 2 π i s k}{r}) | a x mod N ⟩ ⟨ x | a^{k} mod N ⟩ + 0 = & \frac{1}{\sqrt{r}} r - 1 \sum k = 0 exp (\frac{- 2 π i s k}{r}) | a^{k + 1} mod N ⟩ = & \frac{1}{\sqrt{r}} r \sum k^{'} = 1 exp (\frac{- 2 π i s (k^{'} - 1)}{r}) | a^{k^{'}} mod N ⟩ = & exp (\frac{2 π i s}{r}) \cdot \frac{1}{\sqrt{r}} r - 1 \sum k^{'} = 0 exp (\frac{- 2 π i s k^{'}}{r}) | a^{k^{'}} mod N ⟩ = exp (\frac{2 π i s}{r}) \cdot | u_{s} ⟩, \\ (6.5) (6.6) \end{matrix}

有 $| u_{s} ⟩$ 是本征值 $λ_{s}$ 的一个本征向量，这里 $(6.5)$ 处使用了变量代换 $k^{'} = k + 1$ 。

再来验算 $| u_{s} ⟩$ 是单位向量，进而 $| u_{s} ⟩$ 是 $U_{a, N}$ 的一个本征态：

\begin{matrix} \begin{matrix} ⟨ u_{s} | u_{s} ⟩ = & \frac{1}{\sqrt{r}} r - 1 \sum k^{'} = 0 exp (\frac{2 π i s k^{'}}{r}) ⟨ a^{k^{'}} mod N | \cdot \frac{1}{\sqrt{r}} r - 1 \sum k = 0 exp (\frac{- 2 π i s k}{r}) | a^{k} mod N ⟩ = & \frac{1}{r} r - 1 \sum k^{'} = 0 r - 1 \sum k = 0 exp (\frac{2 π i s (k^{'} - k)}{r}) ⟨ a^{k^{'}} mod N | a^{k} mod N ⟩ = & \frac{1}{r} r - 1 \sum k = 0 exp (\frac{2 π i s (k - k)}{r}) = 1, \end{matrix} \\ (6.7) \end{matrix}

这里 $k, k^{'} \in {0, \dots, r - 1}$ 蕴含 $a^{k^{'}} \equiv a^{k} (mod N)$ 当且仅当 $k^{'} = k$ 。 $□$

引理 14.

\begin{matrix} \frac{1}{\sqrt{r}} r - 1 \sum s = 0 | u_{s} ⟩ = | 1 ⟩ . \\ (3.7; 6.8) \end{matrix}

证明.

\begin{matrix} \frac{1}{\sqrt{r}} r - 1 \sum s = 0 | u_{s} ⟩ = & \frac{1}{\sqrt{r}} r - 1 \sum s = 0 \frac{1}{\sqrt{r}} r - 1 \sum k = 0 exp (\frac{- 2 π i s k}{r}) | a^{k} mod N ⟩ = & \frac{1}{r} r - 1 \sum k = 0 (r - 1 \sum s = 0 exp (\frac{- 2 π i s k}{r})) | a^{k} mod N ⟩ = & \frac{1}{r} r - 1 \sum k = 0 r δ_{0, k} | a^{k} mod N ⟩ = | 1 ⟩ . \\ (6.9) \end{matrix}

其中 $\sum_{s = 0}^{r - 1} exp (\frac{- 2 π i s k}{r}) = δ_{0, k}$ 为等比数列求和， $δ_{0, k}$ 为克罗内克 $δ$ 函数，定义为 $δ_{0, k} := ⟨ 0 | k ⟩ = {\begin{matrix} 0, & 若 k \neq 0; 1, & 若 k = 0. \end{matrix}$ $□$

连分数展开

本子段讲述算法 18 的步骤 3 连分数展开相关内容的正确性，所必要的连分数相关知识封闭，全部呈现如下。这里没有引入太多高深的数学知识与工具，只要读者有足够的耐心，定能从中所得。对证明过程并不感兴趣的读者可以只阅读各个定理的内容，直至定理 28，其支撑起了量子求阶算法连分数展开这一步骤的正确性。另外，更本质的量子周期求解算法中也需要用到连分数展开这一步骤。

定义 19（连分数，Continued Fraction）. 给定非负整数 $a_{0}$ 和 $M$ 个正整数 $a_{1}, \dots, a_{M}$ ，设正有理数

\begin{matrix} a = a_{0} + \frac{1}{a_{1} + \frac{1}{a_{2} + \frac{⋱}{\dots + \frac{1}{a_{M}}}}} . \\ (6.10) \end{matrix}

称式 $(6.10)$ 为正有理数 $a$ 的连分数表示，简记为 $a = [a_{0}, \dots, a_{M}]$ 。特别地， $a$ 确定且 $a_{M} > 1$ 时， $a$ 连分数表示唯一。

鉴于正有理数的连分数表示唯一，我们有如下简单的算法计算其连分数表示，称为连分数展开算法。

算法 20（连分数展开）.

输入： 正整数对 $[p, q]$ 表示正有理数 $a = p / q$ 。

输出： 正有理数 $a$ 的连分数表示 $[a_{0}, \dots, a_{M}]$ ，其中 $a_{M} > 1$ 。

def Euclidean_Algorithm(p, q, L=[]):
    """
    :param p: 分子，整数
    :param q: 分母，正整数
    :param L: 连分数展开数组的储存器，整数数组，可以理解为使用了数据结构队列
    :return: 将 p/q 连分数展开，各项构成的整数数组
    """
    # L 为连分数表示结果的寄存器，首次被调用时初始化为 []
    Quotient = p // q # 计算商
    Remainder = p - q * Quotient # 计算余数
    L.append(Quotient) # 将商添加到连分数表示队列尾，即 a_n
    if Remainder != 0: # 在余数不为 0 时
        # 递归调用，注意上一步的分母是这一步的分子，上一步的余数是这一步的分母
        return Euclidean_Algorithm(q, Remainder, L) 
    else:
        return L # 返回值 L=[a_0,...,a_M] 即输出

复制

在流程上，上述算法就是欧几里得辗转相除法，特别的是欧几里得辗转相除法关心最后一次余数为 $0$ 时的除数是多少，而连分数展开算法关心每一步的商是多少。

定义 21（截断）. 给定连分数 $[a_{0}, \dots, a_{M}]$ 及 $n = 0, \dots, M - 1$ ，称连分数 $[a_{0}, \dots, a_{n}]$ 为连分数 $[a_{0}, \dots, a_{M}]$ 的第 $n$ 截断（ $n$ -th convergent）。

这里因为此处 convergent 的中文译名不常见，直译第 $n$ “收敛”又会让人感到费解，笔者特选择“截断”一词作为意译。如例 15 所示对 $2729 / 8192 = [0, 3, 545, 1, 4]$ 执行连分数展开有：

\begin{matrix} \begin{matrix} \frac{2729}{8192} = & \frac{1}{\frac{8192}{2729}} = \frac{1}{3 + \frac{5}{2729}} = \frac{1}{3 + \frac{1}{545 + \frac{4}{5}}} = & 0 + \frac{1}{3 + \frac{1}{545 + \frac{1}{1 + \frac{1}{4}}}} =: [0, 3, 545, 1, 4], \end{matrix} \\ (3.8; 6.11) \end{matrix}

可以看出连分数展开的过程就是欧几里得辗转相除法，每一次迭代将分数表示为带分数的形式，然后将带分数的分数部分的分子化为 $1$ ，将新的分母作为下一次迭代的输入，循环至带分数为整数止。而将有理数表示为带分数就是“相除”，分子化为 $1$ 就是 “辗转”。式 $(6.11)$ 中的 $[0] = 0$ ， $[0, 3] = 0 + \frac{1}{3} = \frac{1}{3}$ ， $[0, 3, 545] = 0 + \frac{1}{3 + \frac{1}{545}} = \frac{545}{1636}$ 和 $[0, 3, 545, 1] = 0 + \frac{1}{3 + \frac{1}{545 + \frac{1}{1}}} = \frac{546}{1639}$ 即为连分数 $[0, 3, 545, 1, 4]$ 的第 $0$ ，第 $1$ ，第 $2$ 和第 $3$ 截断。

引入截断是为了更好解释连分数的一些性质，这里我们首先给出各截断与连分数表示之间的关系，如下：

定理 22. 给定连分数 $[a_{0}, \dots, a_{M}]$ ，令 $p_{0} = a_{0}, q_{0} = 1, p_{1} = 1 + a_{0} a_{1}, q_{1} = a_{1}$ ，对 $n = 2, \dots, M$ 设

\begin{matrix} \begin{matrix} p_{n} = & a_{n} p_{n - 1} + p_{n - 2}; q_{n} = & a_{n} q_{n - 1} + q_{n - 2} . \end{matrix} \\ (6.12) \end{matrix}

则对 $n = 0, 1, \dots, M$ ，我们有连分数

\begin{matrix} [a_{0}, \dots, a_{n}] = p_{n} / q_{n} . \\ (6.13) \end{matrix}

特别地， ${p_{n}}, {q_{n}}$ 都是严格单调递增的正整数列。

证明. 可以直接验证定理 22 在 $n = 1, 2$ 时成立。对 $k \geq 3$ 归纳，假设定理 22 在 $n = 1, \dots, k - 1$ 时成立。当 $n = k$ 时，我们需要验证定理 22 仍成立。考虑连分数表示 $[a_{0}, \dots, a_{n - 2}, a_{n - 1} + 1 / a_{n}]$ ，这里尽管 $a_{n - 1} + 1 / a_{n}$ 不是整数，但上述的连分数的定义仍可推广至此，我们可以验证定理 22 在 $n = k - 1$ 的情形对该新连分数仍成立。设 $[a_{0}, \dots, a_{n - 2}, a_{n - 1} + 1 / a_{n}]$ 在定理 22 中诱导出的序列 ${p_{n}}, {q_{n}}$ 为 ${{~ p}_{n}}, {{~ q}_{n}}$ ，那么我们有 ${~ p}_{n - 2} = p_{n - 2}, {~ q}_{n - 2} = q_{n - 2}, {~ p}_{n - 3} = p_{n - 3}$ 和 ${~ q}_{n - 3} = q_{n - 3}$ 。然后由式 $(6.12)$ 有：

\begin{matrix} \begin{matrix} \frac{{~ p}_{n - 1}}{{~ q}_{n - 1}} = & \frac{(a_{n - 1} + 1 / a_{n}) {~ p}_{n - 2} + {~ p}_{n - 3}}{(a_{n - 1} + 1 / a_{n}) {~ q}_{n - 2} + {~ q}_{n - 3}} = \frac{(a_{n - 1} + 1 / a_{n}) p_{n - 2} + p_{n - 3}}{(a_{n - 1} + 1 / a_{n}) q_{n - 2} + q_{n - 3}} = & \frac{(a_{n - 1} p_{n - 2} + p_{n - 3}) + p_{n - 2} / a_{n}}{(a_{n - 1} q_{n - 2} + q_{n - 3}) + q_{n - 2} / a_{n}} = \frac{p_{n - 1} + p_{n - 2} / a_{n}}{q_{n - 1} + q_{n - 2} / a_{n}} = \frac{p_{n}}{q_{n}} . \end{matrix} \\ (6.14) \end{matrix}

由 $\forall n \in {1, \dots, M}, a_{n}, p_{n}, q_{n} \geq 1$ 及式 $(6.12)$ 有 ${p_{n}}, {q_{n}}$ 都是严格单调递增的正整数列。 $□$

除了单调之外，截断还有下述额外性质，其蕴含了 $p_{n} / q_{n}$ 是既约分数，或者说 $p_{n}$ 与 $q_{n}$ 互素：

推论 23（裴蜀问题）. 符号定义同定理 22，对整数 $n = 1, \dots, M$ 有：

\begin{matrix} q_{n} p_{n - 1} - p_{n} q_{n - 1} = {(- 1)}^{n} . \\ (6.15) \end{matrix}

特别地 $q_{n} p_{n - 1} - p_{n} q_{n - 1} = {(- 1)}^{n}$ 还蕴含了 $p_{n}$ 与 $q_{n}$ 互素。

证明. 由定义验证可知 $q_{1} p_{0} - p_{1} q_{0} = - 1$ ，即 $n = 1$ 时式 $(6.15)$ 成立。当 $n \geq 2$ 时，将式 $(6.12)$ 代入式 $(6.15)$ 左边有：

\begin{matrix} q_{n} p_{n - 1} - p_{n} q_{n - 1} = (a_{n} q_{n - 1} + q_{n - 2}) p_{n - 1} - (a_{n} p_{n - 1} + p_{n - 2}) q_{n - 1} = - (q_{n - 1} p_{n - 2} - p_{n - 1} q_{n - 2}) . \\ (6.16) \end{matrix}

对 $n$ 归纳可知对任意正整数 $n$ 有 $q_{n} p_{n - 1} - p_{n} q_{n - 1} = {(- 1)}^{n}$ 。

考虑最大公因数 $gcd (p_{n}, q_{n}) | (q_{n} p_{n - 1} - p_{n} q_{n - 1}) = {(- 1)}^{n}$ ，故而 $gcd (p_{n}, q_{n}) = 1$ ，即 $p_{n}$ 与 $q_{n}$ 互素。 $□$

除了整数列 ${p_{n}}$ 和 ${q_{n}}$ 单调以外， $∣ ∣ \frac{p_{M}}{q_{M}} - \frac{p_{m}}{q_{m}} ∣ ∣$ 关于 $m$ 也是单调的，即推论 26，但这个单调性的证明就没那么好证了。笔者不才，下面给出一种比较繁琐的证明方式，需要途经引理 24 和 25 方才得到推论 26。

引理 24. 符号定义同定理 22，对任意的非负整数 $n$ 有

\begin{matrix} \begin{matrix} \frac{p_{n + 1}}{q_{n + 1}} - \frac{p_{n}}{q_{n}} = & \frac{{(- 1)}^{n}}{q_{n} q_{n + 1}} \frac{p_{n + 2}}{q_{n + 2}} - \frac{p_{n}}{q_{n}} = & \frac{{(- 1)}^{n} a_{n + 2}}{q_{n} q_{n + 2}} \end{matrix} \\ (6.17) \end{matrix}

更一般地，对正整数 $m$ 有

\begin{matrix} \frac{p_{n + m}}{q_{n + m}} - \frac{p_{n}}{q_{n}} = {(- 1)}^{n} ⎛ ⎝ ⌊ m / 2 ⌋ - 1 \sum i = 0 \frac{a_{n + 2 i + 2}}{q_{n + 2 i} q_{n + 2 i + 2}} + \frac{δ_{- 1} ({(- 1)}^{m})}{q_{n + m - 1} q_{n + m}} ⎞ ⎠, \\ (6.18) \end{matrix}

其中 $δ_{- 1} ({(- 1)}^{m}) = 1$ 若 $m$ 为奇数， $δ_{- 1} ({(- 1)}^{m}) = 0$ 若 $m$ 为偶数。特别地，

\begin{matrix} sgn (\frac{p_{n + m}}{q_{n + m}} - \frac{p_{n}}{q_{n}}) = {(- 1)}^{n}, \\ (6.19) \end{matrix}

这里 $sgn (x) \in {- 1, 0, 1}$ 表示实数 $x$ 的符号。

证明. 由推论 23 有

\begin{matrix} \frac{p_{n + 1}}{q_{n + 1}} - \frac{p_{n}}{q_{n}} = \frac{p_{n + 1} q_{n} - p_{n} q_{n + 1}}{q_{n} q_{n + 1}} = \frac{{(- 1)}^{n}}{q_{n} q_{n + 1}}, \\ (6.20) \end{matrix}

再由式 $(6.12)$ 有

\begin{matrix} \begin{matrix} \frac{p_{n + 2}}{q_{n + 2}} - \frac{p_{n}}{q_{n}} = & (\frac{p_{n + 2}}{q_{n + 2}} - \frac{p_{n + 1}}{q_{n + 1}}) + (\frac{p_{n + 1}}{q_{n + 1}} - \frac{p_{n}}{q_{n}}) = \frac{{(- 1)}^{n + 1}}{q_{n + 1} q_{n + 2}} + \frac{{(- 1)}^{n}}{q_{n} q_{n + 1}} = & \frac{{(- 1)}^{n} (- q_{n} + q_{n + 2})}{q_{n} q_{n + 1} q_{n + 2}} = \frac{{(- 1)}^{n} (a_{n + 2} q_{n + 1})}{q_{n} q_{n + 1} q_{n + 2}} = \frac{{(- 1)}^{n} a_{n + 2}}{q_{n} q_{n + 2}}, \end{matrix} \\ (6.21) \end{matrix}

完成了式 $(6.17)$ 的证明。更一般地，对正整数 $m$ 有

\begin{matrix} \frac{p_{n + m}}{q_{n + m}} - \frac{p_{n}}{q_{n}} = & ⌊ m / 2 ⌋ - 1 \sum i = 0 (\frac{p_{n + 2 i + 2}}{q_{n + 2 i + 2}} - \frac{p_{n + 2 i}}{q_{n + 2 i}}) + (\frac{p_{n + m}}{q_{n + m}} - \frac{p_{n + 2 ⌊ m / 2 ⌋}}{q_{n + 2 ⌊ m / 2 ⌋}}) = & ⌊ m / 2 ⌋ - 1 \sum i = 0 (\frac{p_{n + 2 i + 2}}{q_{n + 2 i + 2}} - \frac{p_{n + 2 i}}{q_{n + 2 i}}) + δ_{- 1} ({(- 1)}^{m}) (\frac{p_{n + m}}{q_{n + m}} - \frac{p_{n + m - 1}}{q_{n + m - 1}}) = & ⌊ m / 2 ⌋ - 1 \sum i = 0 \frac{{(- 1)}^{n + 2 i} a_{n + 2 i + 2}}{q_{n + 2 i} q_{n + 2 i + 2}} + δ_{- 1} ({(- 1)}^{m}) (\frac{{(- 1)}^{n + m - 1}}{q_{n + m - 1} q_{n + m}}) = & {(- 1)}^{n} ⎛ ⎝ ⌊ m / 2 ⌋ - 1 \sum i = 0 \frac{a_{n + 2 i + 2}}{q_{n + 2 i} q_{n + 2 i + 2}} + \frac{δ_{- 1} ({(- 1)}^{m})}{q_{n + m - 1} q_{n + m}} ⎞ ⎠ . \\ (6.22) (6.23) (6.24) \end{matrix}

这里式 $(6.22)$ 处将第二个加项采用 $δ_{- 1} ({(- 1)}^{m})$ 重新表示，式 $(6.23)$ 处代入了式 $(6.17)$ 。因为式 $(6.24)$ 第二个乘项 $> 0$ ，故而 $\frac{p_{n + m}}{q_{n + m}} - \frac{p_{n}}{q_{n}}$ 的符号与 ${(- 1)}^{n}$ 相同。 $□$

以引理 24 为工具，我们先考虑 ${∣ ∣ \frac{p_{M}}{q_{M}} - \frac{p_{m}}{q_{m}} ∣ ∣}$ 作为 $m$ 的数列相邻两项的大小关系：

引理 25. 符号定义同引理 23，给定两个正整数 $s > m$ ，有

\begin{matrix} ∣ ∣ ∣ \frac{p_{s}}{q_{s}} - \frac{p_{m}}{q_{m}} ∣ ∣ ∣ < ∣ ∣ ∣ \frac{p_{s}}{q_{s}} - \frac{p_{m - 1}}{q_{m - 1}} ∣ ∣ ∣ . \\ (6.25) \end{matrix}

证明. 考虑式 $(6.25)$ 左边与右边差的符号，带入式 $(6.19)$ 有：

\begin{matrix} sgn (∣ ∣ ∣ \frac{p_{s}}{q_{s}} - \frac{p_{m}}{q_{m}} ∣ ∣ ∣ - ∣ ∣ ∣ \frac{p_{s}}{q_{s}} - \frac{p_{m - 1}}{q_{m - 1}} ∣ ∣ ∣) = & sgn ({(- 1)}^{m} (\frac{p_{s}}{q_{s}} - \frac{p_{m}}{q_{m}}) - {(- 1)}^{m - 1} (\frac{p_{s}}{q_{s}} - \frac{p_{m - 1}}{q_{m - 1}})) = & {(- 1)}^{m} sgn ((\frac{p_{s}}{q_{s}} - \frac{p_{m}}{q_{m}}) + (\frac{p_{s}}{q_{s}} - \frac{p_{m - 1}}{q_{m - 1}})) = & {(- 1)}^{m} sgn (2 (\frac{p_{s}}{q_{s}} - \frac{p_{m + 1}}{q_{m + 1}}) + (\frac{p_{m + 1}}{q_{m + 1}} - \frac{p_{m}}{q_{m}}) + (\frac{p_{m + 1}}{q_{m + 1}} - \frac{p_{m - 1}}{q_{m - 1}})) = & {(- 1)}^{m} sgn (2 (\frac{p_{s}}{q_{s}} - \frac{p_{m + 1}}{q_{m + 1}}) + (\frac{{(- 1)}^{m}}{q_{m} q_{m + 1}}) + (\frac{{(- 1)}^{m + 1} a_{m + 1}}{q_{m - 1} q_{m + 1}})) = & {(- 1)}^{m} sgn (2 (\frac{p_{s}}{q_{s}} - \frac{p_{m + 1}}{q_{m + 1}}) + {(- 1)}^{m} \frac{q_{m - 1} - a_{m + 1} q_{m}}{q_{m - 1} q_{m} q_{m + 1}}) = & {(- 1)}^{m} sgn (2 {(- 1)}^{m + 1} ∣ ∣ ∣ \frac{p_{s}}{q_{s}} - \frac{p_{m + 1}}{q_{m + 1}} ∣ ∣ ∣ + {(- 1)}^{m + 1} ∣ ∣ ∣ \frac{q_{m - 1} - a_{m + 1} q_{m}}{q_{m - 1} q_{m} q_{m + 1}} ∣ ∣ ∣) = & - sgn (2 ∣ ∣ ∣ \frac{p_{s}}{q_{s}} - \frac{p_{m + 1}}{q_{m + 1}} ∣ ∣ ∣ + ∣ ∣ ∣ \frac{q_{m - 1} - a_{m + 1} q_{m}}{q_{m - 1} q_{m} q_{m + 1}} ∣ ∣ ∣) = & - 1, \\ (6.26) (6.27) (6.28) (6.29) \end{matrix}

其中式 $(6.26)$ 处带入了式 $(6.12)$ ； $(6.27)$ 处 $sgn$ 函数内第一个加项再次带入了式 $(6.19)$ ，第二个加项中由数列 ${q_{n}}$ 的严格单调递增性（定理 22）可以推出 $q_{m - 1} - a_{m + 1} q_{m} < (1 - a_{m + 1}) q_{m - 1} \leq 0$ ，故而加上绝对值符号后外面的系数发生了变号；式 $(6.28)$ 处提出了因数 $(- 1)^{m + 1}$ ；式 $(6.29)$ 中 $q_{m - 1} - a_{m + 1} q_{m}$ 严格大于 $0$ 诱导该 $sgn$ 函数值为 $1$ （不是 $0$ ）。故式 $(6.25)$ 得证。 $□$

将引理 25 推广到一般情形有：

\begin{matrix} ∣ ∣ ∣ \frac{p_{s}}{q_{s}} - \frac{p_{s - 1}}{q_{s - 1}} ∣ ∣ ∣ < ∣ ∣ ∣ \frac{p_{s}}{q_{s}} - \frac{p_{s - 2}}{q_{s - 2}} ∣ ∣ ∣ < \dots < ∣ ∣ ∣ \frac{p_{s}}{q_{s}} - \frac{p_{0}}{q_{0}} ∣ ∣ ∣ . \\ (6.30) \end{matrix}

推论 26. 符号定义同引理 23，给定三个正整数 $s > m_{2} > m_{1}$ ，有

\begin{matrix} ∣ ∣ ∣ \frac{p_{s}}{q_{s}} - \frac{p_{m_{2}}}{q_{m_{2}}} ∣ ∣ ∣ < ∣ ∣ ∣ \frac{p_{s}}{q_{s}} - \frac{p_{m_{1}}}{q_{m_{1}}} ∣ ∣ ∣, \\ (6.31) \end{matrix}

即连分数 $p_{s} / q_{s}$ 的第 $m$ 截断与该连分数的距离关于 $m$ 严格单调递减。

下面的定理 27 讲的是连分数截断的必然性，若一个有理数 $p / q$ 是另一个有理数 $x$ 满足一定条件的逼近，那么 $p / q$ 必然是 $x$ 的连分数表示的某个截断：

定理 27. 给定有理数 $x$ ，若另一个（既约）有理数 $p / q$ 满足

\begin{matrix} ∣ ∣ ∣ \frac{p}{q} - x ∣ ∣ ∣ \leq \frac{1}{2 q^{2}} . \\ (6.32) \end{matrix}

那么 $p / q$ 是 $x$ 的连分数表示的某个截断，也就是说

\begin{matrix} x := [a_{0}, a_{1}, \dots, a_{M}] ⟹ \exists n \in {0, \dots, M}, s.t. p / q = [a_{0}, \dots, a_{n}], \\ (6.33) \end{matrix}

证明. 令 $p / q = [a_{0}, \dots, a_{n}]$ ，如式 $(6.12)$ 定义正整数列 ${p_{j}}$ 和 ${q_{j}}$ ，设 $δ = 2 q_{n}^{2} x - 2 p_{n} q_{n}$ ，那么我们有

\begin{matrix} x = \frac{p_{n}}{q_{n}} + \frac{δ}{2 q_{n}^{2}} . \\ (6.34) \end{matrix}

再由 $p = p_{n}$ ， $q = q_{n}$ 及式 $(6.32)$ 有 $| δ | \leq 1$ 。设

\begin{matrix} λ = 2 (\frac{q_{n} p_{n - 1} - p_{n} q_{n - 1}}{δ}) - \frac{q_{n - 1}}{q_{n}}, \\ (6.35) \end{matrix}

将式 $(6.15)$ 代入，上式 $(6.35)$ 化简为

\begin{matrix} λ = \frac{2 (- 1)^{n}}{δ} - \frac{q_{n - 1}}{q_{n}} . \\ (6.36) \end{matrix}

经验算可知

\begin{matrix} x = \frac{λ p_{n} + p_{n - 1}}{λ q_{n} + q_{n - 1}}, \\ (6.37) \end{matrix}

和（由定理 22）

\begin{matrix} x = [a_{0}, a_{1}, \dots, a_{n}, λ] . \\ (6.38) \end{matrix}

（a）当 $sgn (δ) = (- 1)^{n}$ 时，由式 $(6.38)$ 有

\begin{matrix} λ = \frac{2}{| δ |} - \frac{q_{n - 1}}{q_{n}} > 2 - 1 = 1. \\ (6.39) \end{matrix}

此时可记 $λ = [a_{n + 1}, \dots, a_{M}]$ ，那么我们有 $x = [a_{0}, \dots, a_{M}]$ 。

（b）当 $sgn (δ) = - (- 1)^{n}$ 时，我们不能使用相同的办法得到如式 $(6.39)$ 的结果。先做如下修改：

\begin{matrix} p / q = [a_{0}, a_{1}, \dots, a_{n - 1}, a_{n} - 1, 1], \\ (6.40) \end{matrix}

这里由算法 20 的流程可知 $a_{n} \neq 1$ ，否则 $p / q = [a_{0}, a_{1}, \dots, a_{n - 2}, a_{n - 1} + 1]$ 。然后以新的连分数表示计算数列 ${p_{j}}$ ， ${q_{j}}$ ，并改记为 ${{~ p}_{j}}$ ， ${{~ q}_{j}}$ 。类似式 $(6.36)$ 和式 $(6.39)$ ，我们有

\begin{matrix} λ = \frac{2 (- 1)^{n + 1}}{δ} - \frac{{~ q}_{n}}{{~ q}_{n + 1}} = \frac{2}{| δ |} - \frac{{~ q}_{n}}{{~ q}_{n + 1}} > 2 - 1 = 1. \\ (6.41) \end{matrix}

此时可记 $λ = [a_{n + 2}, \dots, a_{M}]$ ，由式 $(6.39)$ ，我们有 $x = [a_{0}, \dots, a_{n - 1}, a_{n} - 1, 1, a_{n + 2}, \dots, a_{M}]$ 。

综（a）（b）所述， $p / q$ 是 $x$ 的连分数表示的某个截断。 $□$

有了推论 26 和定理 27 还不足以证明连分数展开步骤的正确性，所缺失的最后一块拼图是有理数的间隙有多大：

定理 28. 当相位估计算法返回了 $s / r$ 的一个 $2^{- 2 L - 1}$ -估计 $2^{- t} m$ 时， $2^{- t} m$ 的连分数表示有唯一截断 $p / q$ 满足 $q < N$ 且 $| p / q - 2^{- t} m | \leq 2^{- 2 L - 1}$ 。此时 $p / q$ 恰是 $2^{- t} m$ 的连分数表示的截断 $p_{n} / q_{n}$ 中满足 $q_{n} < N$ 且 $q_{n}$ 最大的那个截断。

证明. 由 $2^{- t} m$ 是 $s / r$ 的 $2^{- 2 L - 1}$ -估计及 $r < N < 2^{L}$ ，有

\begin{matrix} | s / r - 2^{- t} m | \leq 2^{- 2 L - 1} \leq \frac{1}{2 r^{2}}, \\ (6.42) \end{matrix}

由定理 27 有 $s / r$ 是 $2^{- t} m$ 的连分数表示的某个截断。考虑分母不超过 $N$ 的有理数全体构成的集合 $Q_{(N)}$ ，其中的任意两个元素的距离不少于

\begin{matrix} \frac{1}{N - 1} - \frac{1}{N} = \frac{1}{N^{2} - N} > 2^{- 2 L} \geq 2 | s / r - 2^{- t} m | . \\ (6.43) \end{matrix}

故而集合 $Q_{(N)}$ 中至多一个元素与 $2^{- t} m$ 距离 $\leq | s / r - 2^{- t} m |$ ，也就是说 $s / r$ 是集合 $Q_{(N)}$ 中与 $2^{- t} m$ 距离最近（且唯一）的元素。由推论 26 知， $2^{- t} m$ 的各截断与 $2^{- t} m$ 的距离随截断序数增加而减少，故 $s / r$ 恰是满足 $q_{n} < N$ 且 $q_{n}$ 最大的那个截断 $p_{n} / q_{n}$ 。

量子求阶算法的成功概率

本子段旨在给出量子求阶算法的成功概率，即定理 17 的证明，需要引入一些初等数论相关的知识：主要扩展模 $N$ 的阶这一概念。

首先我们引入欧拉函数（Euler's Totient Function）和费马-欧拉定理：

定义 29（欧拉函数）. 给定正整数 $N$ ，记 ${1, 2, \dots, N}$ 中与 $N$ 互素的元素数量为欧拉函数 $φ (N)$ 。设 $N = \prod_{i = 1}^{s} p_{i}^{m_{i}}$ 为 $N$ 的素因数分解，其中 $p_{i}$ 为互不相同的素数， $m_{i}$ 为正整数，有欧拉函数

\begin{matrix} φ (N) = N s \prod i = 1 (1 - \frac{1}{p_{i}}) = s \prod i = 1 (p_{i}^{m_{i}} - p_{i}^{m_{i} - 1}) . \\ (6.44) \end{matrix}

特别地， $φ (1) = 1$ ，对素数 $p$ 有 $φ (p) = p - 1$ 。

定理 30（费马-欧拉定理）. 任给定互素的正整数 $a$ 和 $N$ ，有

\begin{matrix} a^{φ (N)} \equiv 1 (mod N) . \\ (6.45) \end{matrix}

由费马-欧拉定理可知 $a$ 与 $N$ 互素时 $a$ 模 $N$ 的阶 $ord (a, N)$ 存在。另外阶 $ord (a, N)$ 与欧拉函数 $φ (N)$ 之间还有整除关系，更一般的形式即下述习题 2，为了证明习题 2 我们还要引入两个数的最大公约数与这两个数之间的关系，即裴蜀恒等式（Bézout's Identity）。

定理 31（裴蜀恒等式）. 给定整数 $a, b$ ，存在整数 $x, y$ ，使得

\begin{matrix} a x + b y = gcd (a, b) . \\ (6.46) \end{matrix}

由裴蜀恒等式便可推导出：

习题 2. 任给定正整数 $a, b, N$ ，若正整数 $b$ 满足 $a^{b} \equiv 1 (mod N)$ ，则有 $a$ 模 $N$ 的阶 $ord (a, N)$ 存在，且 $ord (a, N) | b$ 。特别地， $ord (a, N) | φ (N)$ ， $ord (a, N) | gcd (b, φ (N))$ 。

最后我们再来证明量子求阶算法的成功概率，即定理 17。

定理 17. 通过执行两次量子求阶算法，可以得到两个分母 $q$ 和 $q^{'}$ ，计算 $q$ 和 $q^{'}$ 的最小公倍数 $~ r := q q^{'} / gcd (q, q^{'})$ ，则 $~ r = r$ 的概率 $\geq \frac{6 (1 - ϵ)^{2}}{π^{2}} \geq 60.7 % \cdot (1 - ϵ)^{2}$ ，这里 $ϵ$ 是相位估计算法的失败概率。

证明. 先假设两次相位估计都达到了精度 $2^{- 2 L - 1}$ 。

由注 16，可设 $q = r / gcd (s, r)$ ， $q = r / gcd (s^{'}, r)$ ，有

\begin{matrix} \begin{matrix} ~ r = & \frac{q q^{'}}{gcd (q, q^{'})} = \frac{r^{2}}{gcd (s, r) gcd (s^{'}, r) gcd (r / gcd (s, r), r / gcd (s^{'}, r))} = & \frac{r^{2}}{gcd (r \cdot gcd (s^{'}, r), r \cdot gcd (s, r)))} = \frac{r^{2}}{r \cdot gcd (gcd (s^{'}, r), gcd (s, r)))} = \frac{r}{gcd (s, s^{'}, r)}, \end{matrix} \\ (6.47) \end{matrix}

这里用到了最大公因数的一条性质 $a \cdot gcd (b, c) = gcd (a \cdot b, a \cdot c)$ 。所以 $~ r = r$ 当且仅当 $gcd (s, s^{'}, r) = 1$ 。由引理 14， $s$ 和 $s^{'}$ 独立地服从 ${0, 1, \dots, r - 1}$ 上的均匀分布，那么 $gcd (s, s^{'}) = 1$ 的概率为

\begin{matrix} \prod p prime < r (1 - 1 / p^{2}) \geq 1 - \sum p prime < r 1 / p^{2} \geq 1 - \infty \sum j = 2 1 / j^{2} = 1 - (π^{2} / 6 - 1) \geq 35.5 %, \\ (6.48) \end{matrix}

这里 $p prime < r$ 表示 $p$ 遍历所有 $< r$ 的素数，使用了级数 $\sum_{j = 1}^{\infty} j^{- 2} = π^{2} / 6$ 这一结果（最常见的证明需要傅里叶分析）。

通过更细致地分析，借由黎曼 $ζ$ 函数的欧拉乘积公式（Euler Product Formula for the Riemann zeta Function）：

\begin{matrix} \infty \sum j = 1 \frac{1}{j^{s}} = \prod p prime \frac{1}{1 - p^{- s}}, \\ (6.49) \end{matrix}

可以得到

\begin{matrix} \prod p prime < r (1 - 1 / p^{2}) \geq \prod p prime (1 - 1 / p^{2}) = {(\infty \sum j = 1 \frac{1}{j^{2}})}^{- 1} = \frac{6}{π^{2}} \geq 60.7 %, \\ (6.50) \end{matrix}

这里 $p prime$ 表示遍历全部素数 $p$ ，计算累积。

考虑到量子相位估计算法的成功概率为 $1 - ϵ$ 及 $gcd (s, s^{'}) = 1$ 蕴含 $gcd (s, s^{'}, r) = 1$ ，综上所述，本引理得证。 $□$

由注 16 有：

推论 32. 若只执行一次量子求阶算法，得到的分母 $q = r$ 的概率不少于 $(1 - ϵ) \cdot φ (r) / r$ ，这里 $ϵ$ 是相位估计算法的失败概率， $φ$ 为定义 29 中的欧拉函数。

“良性”的比例

本子段需要一些抽象代数（或称近世代数）的知识，即环论的中国剩余定理（或称孙子定理），用于估算秀尔算法中搜索的 $a$ 具有“良性”的概率。阅读本节需要知道有限群的阶，群中元素的阶，商环等概念。

整数环的商环的中国剩余定理可描述如下，为了合理地简化问题，我们这里假设 $N$ 是奇数，即没有素因子 $2$ ：

定理 33（中国剩余定理）. 设 $N = \prod_{i = 1}^{s} p_{i}^{m_{i}} > 1$ 为正奇数 $N$ 的素因数分解，其中 $p_{i}$ 为互不相同的奇素数， $m_{i}$ 为正整数，那么有环同构：

\begin{matrix} Z / N Z ≅ s \prod i = 1 Z / p_{i}^{m_{i}} Z . \\ (6.51) \end{matrix}

其诱导了环的乘法群的群同构：

\begin{matrix} (Z / N Z)^{*} ≅ s \prod i = 1 (Z / p_{i}^{m_{i}} Z)^{*} ≅ s \prod i = 1 (C_{p_{i}^{m_{i} - 1}} \times C_{p_{i} - 1}) . \\ (6.52) \end{matrix}

这里 $C_{n}$ 表示 $n$ 阶循环群。记 $p_{i} - 1 = 2^{b_{i}} \cdot q_{i}$ 使得 $q_{i}$ 为奇数，且 $b_{i}$ 为正整数，有：

\begin{matrix} (Z / N Z)^{*} ≅ s \prod i = 1 (C_{p_{i}^{m_{i} - 1}} \times C_{q_{j}}) \times s \prod i = 1 C_{2^{b_{i}}} =: A_{Q} \times t \prod l = 1 C_{2^{l}}^{n_{l}} =: A_{Q} \times A_{C}, \\ (6.53) \end{matrix}

这里 $t$ 是各 $b_{i}$ 中的最大值, $n_{l}$ 是与 $l$ 相等的 $b_{i}$ 的数量。特别地，这里 $A_{Q}$ 为一个奇数阶交换群， $A_{C}$ 为一个 $2$ 的幂次阶交换群（ $| A_{C} | = \prod_{l = 1}^{t} (2^{l})^{n_{l}} = 2^{\sum_{l = 1}^{t} l n_{l}}$ ，这里用 $| \cdot |$ 表述群的元素数量，称群的阶）。

考虑 $a = (\to b, \to c) \in A_{Q} \times A_{C}$ 具有偶数阶 $2^{r} \cdot q$ ，这里 $\to b \in A_{Q}$ 和 $\to c \in A_{C}$ 分别称为 $a$ 在 $A_{Q}$ 和 $A_{C}$ 中的投影， $r$ 为正整数， $q$ 为正奇数，那么

\begin{matrix} a^{2^{r - 1} q} = 2^{r - 1} q (\to b, \to c) = (2^{r - 1} q \to b, 2^{r - 1} q \to c), \\ (6.54) \end{matrix}

这里向量表示的环中元素改用加法表示交换群上的乘法，改用数乘表示群上的整数幂运算。考虑到 $A_{Q}$ 为一个奇数阶交换群， $A_{C}$ 为一个 $2$ 的幂次阶交换群，我们有（由群中元素的阶整除群的阶，类似习题 2）

\begin{matrix} a^{2^{r} q} = (2^{r} q \to b, 2^{r} q \to c) = (0, 0) ⟹ q \to b = 0 \in A_{Q} \land 2^{r} \to c = 0 \in A_{C} . \\ (6.55) \end{matrix}

将上式 $(6.55)$ 带入式 $(6.54)$ 可以得到 $a^{2^{r - 1} q} = (0, 2^{r - 1} q \to c)$ 。所以， $a^{2^{r - 1} q} \neq - 1 \in (Z / N Z)^{*} ⟺ (0, 2^{r - 1} q \to c) \neq - 1 \in (Z / N Z)^{*}$ ，与 $a$ 在 $A_{Q}$ 中的投影 $\to b$ 无关。故而当我们考虑元素 $a$ 是否具有“良性”时，可以不失一般性地假设 $A_{Q} = {0}$ 为平凡群，只考虑另一部分：交换群

\begin{matrix} A_{C} := t \prod l = 1 C_{2^{l}}^{n_{l}} = {x = (x_{11}, x_{12}, \dots, x_{t n_{t}}) | x_{i j} \in C_{2^{l}} = {0, 1, \dots, 2^{l} - 1}} . \\ (6.56) \end{matrix}

仍记 $- 1 \in (Z / N Z)^{*}$ 在 $A_{C}$ 上的投影为 $- 1$ ，有：

\begin{matrix} - 1 = (n_{1}      1, \dots, 1, n_{2}      2, \dots, 2, 2^{2}, \dots, 2^{t - 1}) \in A_{C} . \\ (6.57) \end{matrix}

假设

\begin{matrix} \exists x = (x_{11}, x_{12}, \dots, x_{t n_{t}}) \in A_{C}, \exists l \geq 1 \in Z, \\ (6.58) \end{matrix}

使得

\begin{matrix} 2^{l - 1} x = (2^{l - 1} x_{11}, 2^{l - 1} x_{12}, \dots, 2^{l - 1} x_{t n_{t}}) = - 1, \\ (6.59) \end{matrix}

即

\begin{matrix} \forall i \in {1, \dots, t}, \forall j \in {1, \dots, n_{i}}, 2^{l - 1} x_{i j} = 2^{i - 1} . \\ (6.60) \end{matrix}

由 $x_{i j} \in C_{2^{i}}$ 有

\begin{matrix} \forall i \leq l - 1, 2^{l - 1} x_{i j} = 0 \in C_{2^{i}} . \\ (6.61) \end{matrix}

上式 $(6.61)$ 与式 $(6.60)$ 矛盾，这迫使

\begin{matrix} n_{1} = n_{2} = \dots n_{l - 1} = 0, \\ (6.62) \end{matrix}

亦有式 $(6.60)$ 可化简为

\begin{matrix} \forall i \in {l, \dots, t}, \forall j \in {1, \dots, n_{i}}, 2^{l - 1} x_{i j} = 2^{i - 1} \in C_{2^{i}} . \\ (6.63) \end{matrix}

此时相当于 $x_{i j} \in {0, 1, \dots, 2^{i} - 1}$ 且 $2^{l - 1} x_{i j} \equiv 2^{i - 1} (mod 2^{i})$ ，解之得

\begin{matrix} \exists k_{i j} \in {1, \dots, 2^{l - 1}}, x_{i j} = 2^{i - l} (2 k_{i j} - 1) . \\ (6.64) \end{matrix}

故而满足式 $(6.63)$ 的 $x = (x_{l 1}, x_{l 2}, \dots, x_{t n_{t}})$ 的数量为

\begin{matrix} ∣ ∣ ∣ ∣ t \prod i = l n_{i} \prod j = 1 {x_{i j} = 2^{i - l} (2 k_{i j} - 1) | k_{i j} = 1, \dots, 2^{l - 1}} ∣ ∣ ∣ ∣ = t \prod i = l n_{i} \prod j = 1 2^{l - 1} = 2^{(l - 1) \sum_{i = l}^{t} n_{i}} =: 2^{(l - 1) s}, \\ (6.65) \end{matrix}

其中 $s := \sum_{i = l}^{t} n_{i}$ 为 $N$ 所含不同素因数的数量。

设 $u = min {i | n_{i} \neq 0}$ ，由式 $(6.62)$ 我们有式 $(6.59)$ 蕴含 $l \leq u$ 。故由式 $(6.65)$ 有

\begin{matrix} \begin{matrix} ∣ ∣ {x \in A_{C} | \exists l \in N^{*}, 2^{l - 1} x = - 1 \in A_{C}} ∣ ∣ = & u \sum l = 1 ∣ ∣ {x \in A_{C} | 2^{l - 1} x = - 1 \in A_{C}} ∣ ∣ = & u \sum l = 1 2^{(l - 1) s} + t \sum l = u + 1 0 = \frac{2^{u s} - 1}{2^{s} - 1} . \end{matrix} \\ (6.66) \end{matrix}

那么 $a$ 具备“良性”的概率

\begin{matrix} \begin{matrix} p (2 | ord (a, N) \land a^{ord (a, N) / 2} \equiv̸ - 1 (mod N)) = & \frac{∣ ∣ {x \in A_{C} | x \neq 0 \land \forall l \in N^{*}, 2^{l - 1} x \neq - 1 \in A_{C}} ∣ ∣}{| A_{C} |} = & \frac{| A_{C} | - | {0} | - ∣ ∣ {x \in A_{C} | \exists l \in N^{*}, 2^{l - 1} x = - 1 \in A_{C}} ∣ ∣}{| A_{C} |} = & 1 - \frac{1 + \frac{2^{u s} - 1}{2^{s} - 1}}{| A_{C} |} = 1 - \frac{2^{s} - 1 + 2^{u s} - 1}{(2^{s} - 1) 2^{\sum_{l = u}^{t} l n_{l}}} \geq & 1 - \frac{2^{u s} + 2^{s} - 2}{(2^{s} - 1) 2^{u s}} \geq 1 - \frac{2^{s} + 2^{s} - 2}{(2^{s} - 1) 2^{s}} = 1 - \frac{2}{2^{s}} . \end{matrix} \\ (6.67) \end{matrix}

整理结论得到下述定理，这个结果已经过度放缩了，式 $(6.67)$ 中蕴含了更紧的下界。

定理 34. 随机抽取 $a \in {2, 3, \dots, N - 2}$ ，其具有“良性”即满足

\begin{matrix} 2 | ord (a, N) \land a^{ord (a, N) / 2} \equiv̸ - 1 (mod N) \\ (6.68) \end{matrix}

的概率 $p \geq 1 - \frac{2}{2^{s}}$ ，其中 $s$ 是 $N$ 互不相同的素因数的数量。

$C (U_{a, N})$ 门的实现

在量易伏实现中，我们使用代数手段将 $C (U_{a, N})$ 分解为 $C^{2} (X)$ 门的组合，并借由经典计算简化了 $C (U_{a, N}^{2^{j}})$ 门。这在逻辑上其实是说不通的，无论是门的分解的过程，还是省略计算的过程都有循环论证之嫌：涉及到或者蕴含到已知 $N$ 的分解这一事实。若想引入一个常规的泛用的 $C (U_{a, N})$ 门的分解可以参照文献 $[6, 7]$ 等。另外读者也可以参考 QCompute SDK 中的秀尔算法代码样例和教程。

除了这种为了控制辅助量子比特数量而设计的复杂的方法之外，我们可以借助经典的可逆计算的途径来实现 $C (U_{a, N})$ 。

从（经典）可逆计算开始

什么是经典的可逆计算？大家应该都知道经典计算机是基于 2 进制数及逻辑运算实现的，逻辑运算一般用三种基本的门“与”门 $\land$ 、“或”门 $\lor$ 和“非”门 $\neg$ 来实现，除去“非”门外，“与”门和“或”门都是两个输入一个输出：

\begin{matrix} \forall a, b \in {0, 1}, a \land b = a b, a \lor b = \neg (\neg a \land \neg b) = 1 - (1 - a) (1 - b), \neg a = 1 - a . \\ (6.69) \end{matrix}

这种无法从输出完全恢复出输入信息的计算就被称为不可逆计算，反之称为 可逆计算。可以看出对可逆计算而言，输入与输出的比特数量应当相等，比如上述的“非”门就是可逆计算。

可以预料到当我们保留输入信息，并在额外的寄存器上实现计算时，实现的计算必然是可逆计算，比如：

\begin{matrix} (a, b, 0) \to (a, b, f (a, b)), \\ (6.70) \end{matrix}

这里 $f (a, b)$ 是任意一个 ${0, 1}^{2} \to {0, 1}$ 的函数。

最后我们来讲，通过引入辅助量子比特，量子计算中的 $X$ 门、 $C (X)$ 门和 Toffoli 门 $C^{2} (X)$ ，便可模拟所有的可逆计算，其中：

\begin{matrix} \begin{matrix} X | a ⟩ = & | \neg a ⟩, C^{2} (X) | a ⟩ | b ⟩ | 0 ⟩ = & | a ⟩ | b ⟩ | a \land b ⟩, C (X) | a ⟩ | 0 ⟩ = & | a ⟩ | a ⟩ =: FANOUT | a ⟩ \end{matrix} \\ (6.71) \end{matrix}

这里“扇出” $FANOUT$ 是将一个经典比特复制为两份，这在量子计算对计算基能有效实现。另外由“或”门的分解，我们也可以得到其量子版本：

\begin{matrix} (X \otimes X \otimes X) \cdot C^{2} (X) \cdot (X \otimes X \otimes I) | a ⟩ | b ⟩ | 0 ⟩ = | a ⟩ | b ⟩ | a \lor b ⟩ . \\ (6.72) \end{matrix}

回到 $U_{a, N}$ 门

\begin{matrix} U_{a, N} = N - 1 \sum x = 0 | a x mod N ⟩ ⟨ x | + 2^{L} - 1 \sum x = N | x ⟩ ⟨ x | . \\ (3.2; 6.73) \end{matrix}

给定一个输入 $x \in {0, 1, \dots, 2^{L} - 1}$ ，我们总能经典地计算判断 $x \leq N$ ，并在 $x \leq N$ 时计算 $a x mod N$ ，这个判断和计算的过程终归是用逻辑运算实现的，即最终被分解为“与”门、“非”门和“扇出”门的组合，进而可被量子电路等效实现。最终我们获得了一个由若干 $X$ 门、 $C (X)$ 门和 Toffoli 门组成的量子电路 $C$ ，其实现了

\begin{matrix} C | x ⟩ | 0 ⟩^{\otimes t} = (U_{a, N} | x ⟩) \otimes | 0 ⟩^{\otimes t}, \\ (6.74) \end{matrix}

其中 $t$ 是整个电路需要引入的辅助量子比特的数目。根据可逆计算的性质，上述电路中所需量子门的数目和 $t$ 都不超过经典计算中逻辑运算次数的两倍，这保证了 $C$ 可以被有效实现。

最后更新时间：2021年12月20日

背景介绍
算法总览
量子求阶算法
简例与量易伏平台演示
参考资料
附录