[量子计算]Shor算法

秀尔算法（Shor's algorithm）是一种非常著名的量子算法，甚至可以说是量子计算的一块金字招牌了。Shor's algorithm可以在多项式时间内完成大整数质因数分解。所以Shor's algorithm从诞生之时，就和以RSA算法为根基的加密技术形成了不可调和的矛盾。2001年，IBM的研究小组使用Shor's algorithm完成了15=3×515 = 3×515=3×5这个整数分解运算 。时至今日，快20年过去了，量子计算机没有摧毁RSA。未来这一天来临的时候，希望我们已经准备好了。下图是，经典大数分解和秀尔算法的复杂度对比。

Shor算法的整体流程

完整的Shor算法是需要经典计算机和量子计算机协作完成的。其中量子计算机实现一个周期查找的函数，经典计算机负责整个算法流程的控制，以及调用量子算法。

首先、还是表达一下Shor's algorithm的问题描述：给定一个合成数N，找到整数p在1和N之间且不包含1和N，并且N整除于p。

1. 选择任意数字a<Na < Na<N
2. 计算gcd(a,N)gcd(a, N)gcd(a,N)。可以使用辗转相除法来计算。
3. 若gcd(a,N)≠1gcd(a, N) \not = 1gcd(a,N)​=1(不等于，这里不等号渲染有问题)，则我们有了一个NNN非平凡的因数，因此这部分结束了。否则，利用下面的周期查找子程序（Period-finding subroutine，使用量子计算机完成）来找出下面这个函数的周期rrr：

f(x)=axmodNf(x)=a^{x} {mod} Nf(x)=axmodN

换句话说，找出aaa在ZN{Z} _{N}ZN​里面的阶 rrr，或者最小的正整数rrr令 f(x+r)=f(x)f ( x + r ) = f ( x )f(x+r)=f(x)。

1. 若rrr是奇数，回到第一步。
2. 若ar/2≡−1(modN)a^{ r /2} ≡ -1 (mod N)ar/2≡−1(modN)，回到第一步。否则，计算gcd(ar/2+1,N)gcd(a^{r/2}+1, N)gcd(ar/2+1,N)与gcd(ar/2−1,N)gcd(a^{r/2}-1, N)gcd(ar/2−1,N)，并验证他们是不是分别是N非平凡的因数。成功则分解完成。

第三个步骤中，周期查找函数是使用量子计算机完成的。这里暂且不分析这一部分，而关注于整体流程。第一个步骤就是遍历小于N的整数；第二步其实就是碰碰运气，如果a和N的最大公约数不是1，那就是运气好，算法直接结束了；如果最大公约数是1，则调用周期查找子程序计算f(x)=axmodNf(x)=a^{x} {mod} Nf(x)=axmodN这个函数的阶数rrr，modmodmod NNN就是以NNN为模的取模运算。

因为f(x+r)=f(x)f ( x + r ) = f ( x )f(x+r)=f(x)，而且f(0)=1f ( 0 ) = 1f(0)=1 modmodmod NNN，因此可得：

f(0)=f(0+r)=armodN=1modNf ( 0 ) = f (0+ r )=a^{r} mod N =1mod Nf(0)=f(0+r)=armodN=1modN

ar−1=0modNa^{r}-1=0 modNar−1=0modN

如果rrr为偶数，而且ar/2a^{ r /2}ar/2不等于-1，则可得：

(ar/2+1)(ar/2−1)=0modN(a^{r/2}+1)(a^{r/2}-1)=0 modN(ar/2+1)(ar/2−1)=0modN

观察上面的等式可以很容易发现，(ar/2+1)(ar/2−1)(a^{r/2}+1)(a^{r/2}-1)(ar/2+1)(ar/2−1)和NNN存在不是1的最大公约数，也即(ar/2+1)(a^{r/2}+1)(ar/2+1)或者(ar/2−1)(a^{r/2}-1)(ar/2−1)与NNN非1的最大公约数。通过辗转相除法即可计算出来，然后验证这两个最大公约数是不是NNN的因子即可。如果是的话，算法就成功了。

整个算法的充分性，其实是比较明显的，因为最终都会验证结果是不是NNN的因子。算法的必要性，即是不是一定能找到质因子，我其实还没理解清楚，这里先略过了。接下来就是Shor算法的量子部分，即前文说的Period-finding subroutine。

周期查找算法（Period-finding algorithm）

周期查找算法的作用就是求f(x)=axmodNf(x)=a^{x} {mod} Nf(x)=axmodN这个函数的周期。查找周期的前提是周期存在，这个函数的周期存在吗？欧拉对这个问题做过证明，有如下结论：gcd(y,N)=1gcd ( y, N ) = 1gcd(y,N)=1 ，按 EulerEulerEuler 定理,周期 rrr 一定存在。

周期查找算法其实本质上是特殊情况的相位估计算法变种（如果对相位估计算法，可以看看我前面相位估计的文章），整体的算法流程如下：

1. ∣0⟩∣0⟩\vert 0 \rangle\vert 0 \rangle∣0⟩∣0⟩; 初始化状态
2. ⟶12t∑x=02t−1∣j⟩∣0⟩\longrightarrow\dfrac{1}{\sqrt{2^t}} \displaystyle\sum_{x=0}^{2^t-1} {\vert j \rangle} {\vert 0 \rangle}⟶2t​1​x=0∑2t−1​∣j⟩∣0⟩; 使用Hadamard门生成叠加态
3. ⟶12t∑j=02t−1∣x⟩∣f(x)⟩≈1r2t∑l=0r−1∑x=02t−1e2πilx/r∣x⟩∣f^(l)⟩\longrightarrow\dfrac{1}{\sqrt{2^t}} \displaystyle\sum_{j=0}^{2^t-1} {\vert x \rangle} {\vert f(x) \rangle} \approx \dfrac{1}{\sqrt{r2^t}} \displaystyle \sum_{l=0}^{r-1} \sum_{x=0}^{2^t-1} {e^{2\pi i lx / r}} {\vert x \rangle} {\vert \hat f(l) \rangle}⟶2t​1​j=0∑2t−1​∣x⟩∣f(x)⟩≈r2t​1​l=0∑r−1​x=0∑2t−1​e2πilx/r∣x⟩∣f^​(l)⟩;
   这一步主要是通过"控制UUU变换"，其中U∣x⟩∣y⟩=∣x⟩∣y⨁f(x)⟩U{\vert x \rangle}{\vert y \rangle} = {\vert x \rangle}{\vert y \bigoplus f(x) \rangle}U∣x⟩∣y⟩=∣x⟩∣y⨁f(x)⟩。
4. ⟶1r∑l=0r−1∣l/r~⟩∣f^(l)⟩\longrightarrow \dfrac{1}{\sqrt{r}} \sum_{l=0}^{r-1} {\vert \widetilde{l / r} \rangle} {\vert \hat f(l) \rangle}⟶r​1​∑l=0r−1​∣l/r​⟩∣f^​(l)⟩;
   这一步主要是应用傅里叶反变换。
5. ⟶l/r~\longrightarrow { \widetilde{l / r} }⟶l/r​；
   这一步就是测量上部分寄存器的值。
6. ⟶r\longrightarrow {r}⟶r；
   通过连分数算法计算r的值。

其中，第三步使用了控制U门，而且引入新的表示态：

∣f^(l)⟩=1r∑x=0r−1e−2πilx/r∣x⟩{\vert \hat f(l) \rangle} = \dfrac{1}{\sqrt{r}} \displaystyle \sum_{x=0}^{r-1} {e^{-2\pi i lx / r}} {\vert x \rangle}∣f^​(l)⟩=r​1​x=0∑r−1​e−2πilx/r∣x⟩

其实就是f(x)f(x)f(x)傅里叶变换之后的态，由于2t2^t2t不一定是周期rrr的整倍数，所以步骤三有一个约等号。此外，由于随机性，lll是有很多可能取值的。因此有人可能会质疑如果lll是rrr的整数倍，则连分数算出的rrr是有问题的。但这个其实不用担心，和rrr互质的会更多，所以从概率来讲，正确的rrr会很高概率地被发现。

本来是打算每个算法都做一个线路实现的，可是Shor算法很难搞，这里就只能略过了。这篇水完，量子计算方面应该会暂时告一段落了，毕竟还有太多事情要做。而且写的这些东西，我也不满意，毕竟是现学现卖了，准备之后有时间修改一下。

【1】.《量子计算与量子信息》Michael A. Nielsen
【2】.《量子力学》列文.朗道
【3】.《量子力学原理》P.A.M. Dirac
【4】. ETHZ讲义