Shor算法

Shor算法是一种量子计算算法，由美国数学家彼得·肖尔（Peter Shor）在1994年提出。这个算法能够高效地分解大整数，是量子计算领域的重要突破之一。在经典计算机上，大整数的分解被认为是一个非常困难的问题，特别是当涉及到用作现代加密技术基础的大质数时。Shor算法的提出表明，一旦量子计算机足够强大，它能够破解包括RSA加密在内的多种基于大整数分解难题的加密系统。

在本章中，我们将重点关注Shor算法的量子部分，它实际上解决了寻找周期（period finding）的问题。由于分解问题可以在多项式时间内转化为寻找周期问题，因此一个有效的寻找周期算法也可以用来有效地分解整数。由于寻找周期本身就是一个值得研究的问题，我们将首先解决这个问题，然后再讨论如何用它来进行分解。

1. 问题：寻找周期

让我们看一下这个周期函数：

f(x) = a^x \bmod{N}

模运算（缩写为”mod”）意味着求一个数除以另一个数的余数。例如:

17 \bmod 5 = 2

因为 17 \div 5 = 3 余 2 （即 17 = (3\times 5) + 2 ）。在Python中，模运算用 \% 符号表示。

这种行为在模运算中使用，其中数字在达到某个值（模数）后会”回绕”。使用模运算，我们可以写：

17 = 2 \pmod 5

注意，这里 \pmod 5 适用于整个等式（因为它在括号中），而上面的等式中它只适用于等式的左边。

其中 a 和 N 是正整数， a 小于 N ，并且它们没有公共因子。周期或阶（ r ）是最小的（非零）整数，使得：

a^r \bmod N = 1

2. 解决方案

Shor的解决方案是在幺正算符上使用量子相位估计：

U|y\rangle \equiv |ay \bmod N \rangle

为了看出这怎么有帮助，让我们来计算U的一个本征态可能是什么样子。如果我们从态 |1\rangle 开始，我们可以看到U的每次成功应用都会将我们寄存器的状态乘以 a \pmod N ，经过 r 次应用后，我们将再次到达态 |1\rangle 。例如， a=3 和 N=35 时：

\begin{aligned}
U|1\rangle &= |3\rangle&\\
U^2|1\rangle &= |9\rangle \\
U^3|1\rangle &= |27\rangle \\
& \vdots \\
U^{(r-1)}|1\rangle &= |12\rangle \\
U^r|1\rangle &= |1\rangle
\end{aligned}

所以这个循环中态的叠加（ |u_0\rangle ）将是 U 的一个本征态：

|u_0\rangle = \tfrac{1}{\sqrt{r}}\sum_{k=0}^{r-1}{|a^k \bmod N\rangle}

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

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

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

现在对于每个整数值 s ，其中 0 \leq s \leq r-1 ，我们都有一个独特的本征态。非常方便的是，如果我们把所有这些本征态加起来，不同的相位会抵消掉除 |1\rangle 之外的所有计算基态：

\tfrac{1}{\sqrt{r}}\sum_{s=0}^{r-1} |u_s\rangle = |1\rangle

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

由于计算基态 |1\rangle 是这些本征态的叠加，这意味着如果我们在 U 上使用态 |1\rangle 做QPE，我们将测量一个相位：

\phi = \frac{s}{r}

其中 s 是 0 到 r-1 之间的随机整数。我们最终使用连分数算法对 \phi 进行处理以找到 r 。

接下来，我们将使用 DeepQuantum 来演示Shor算法。在这个演示中，我们将提供 U 的线路而不作解释，但在第4节中，我们将讨论如何有效地构建 U^{2^j} 的电路。

3. DeepQuantum 实现

在这个例子中，我们将解决 a=7 和 N=15 时的周期查找问题。我们提供了 U 的线路，其中：

U|y\rangle = |ay\bmod 15\rangle

在这里不做具体解释。为了创建 U^x，我们只需重复该线路 x 次即可。在下一节中，我们将讨论一种创建这些线路的通用高效方法。函数 c_amod15 返回 a 的受控 U 门，重复 power 次。

from typing import Any, List, Optional, Union
import torch 
import deepquantum as dq
import numpy as np

class c_amod15(dq.Ansatz):
    def __init__(
        self,
        nqubit: int,
        a: int,
        power: int,
        minmax: Optional[List[int]] = None,
        controls: Union[int, List[int], None]=None,
        reverse: bool = False,
        init_state: Any = 'zeros',
        den_mat: bool = False,
        mps: bool = False,
        chi: Optional[int] = None,
        show_barrier: bool = False
    ) -> None:
        super().__init__(nqubit=nqubit,wires=None, minmax=minmax, ancilla=None, controls=controls,
                         init_state=init_state, name='c_amod15', den_mat=den_mat,
                         mps=mps, chi=chi)


        if a not in [2,4,7,8,11,13]:
            raise ValueError("'a' must be 2,4,7,8,11 or 13")

        for _ in range(power):
            if a in [2, 13]:
                self.swap([self.wires[2],self.wires[3]], controls)
                self.swap([self.wires[1],self.wires[2]], controls)
                self.swap([self.wires[0],self.wires[1]], controls)
            if a in [7, 8]:
                self.swap([self.wires[0],self.wires[1]], controls)
                self.swap([self.wires[1],self.wires[2]], controls)
                self.swap([self.wires[2],self.wires[3]], controls)
            if a in [4, 11]:
                self.swap([self.wires[1],self.wires[3]], controls)
                self.swap([self.wires[0],self.wires[2]], controls)
            if a in [7, 11, 13]:
                for q in range(4):
                    self.x(self.wires[q], controls)

cir=c_amod15(nqubit=6,a=2,power=1,minmax=[2, 5],controls=0)
cir.draw() 

c:\Users\HP\.conda\envs\dq\lib\site-packages\qiskit\visualization\circuit\matplotlib.py:266: FutureWarning: The default matplotlib drawer scheme will be changed to "iqp" in a following release. To silence this warning, specify the current default explicitly as style="clifford", or the new default as style="iqp".
  self._style, def_font_ratio = load_style(self._style)

我们将使用 8 个计数量子位:

n_count = 8  
a = 2

我们还导入了 QFT 线路：

class QuantumFourierTransform(dq.Ansatz):
    def __init__(
        self,
        nqubit: int,
        minmax: Optional[List[int]] = None,
        reverse: bool = False,
        init_state: Any = 'zeros',
        den_mat: bool = False,
        mps: bool = False,
        chi: Optional[int] = None,
        show_barrier: bool = False
    ) -> None:
        super().__init__(nqubit=nqubit, wires=None, minmax=minmax, ancilla=None, controls=None,
                         init_state=init_state, name='QuantumFourierTransform', den_mat=den_mat,
                         mps=mps, chi=chi)
        self.reverse = reverse
        for i in self.wires:
            self.qft_block(i)
            if show_barrier:
                self.barrier(self.wires)
        if not reverse:
            for i in range(len(self.wires) // 2):
                self.swap([self.wires[i], self.wires[-1 - i]])

    def qft_block(self, n):
        self.h(n)
        k = 2
        for i in range(n, self.minmax[1]):
            self.cp(i + 1, n, torch.pi / 2 ** (k - 1))
            k += 1

有了这些构建模块，我们就可以轻松地构建 Shor 算法的线路：

cir = dq.QubitCircuit(12)

for q in range(n_count):
    cir.h(q)

cir.x(n_count)

for q in range(n_count):
    cir.add(c_amod15(nqubit=12,a=a, power=2**q, minmax=[8, 11],controls=q))

cir.add(QuantumFourierTransform(nqubit=12, minmax=[0, 7],reverse=True).inverse())

cir.measure(wires=list(range(n_count)))

让我们看看测量得到的结果：

cir()
res=cir.measure(wires=list(range(n_count)))
print(res)


import matplotlib.pyplot as plt

# 准备x轴和y轴的数据
keys = list(res.keys())
values = list(res.values())
counts = res
# 创建条形图
plt.bar(keys, values)

# 设置图表的标题和轴标签
plt.xlabel('Keys')
plt.ylabel('Values')

# 设置x轴的刻度标签，以便它们更容易阅读
plt.xticks(rotation=90)

# 显示图表
plt.show()

{'10000000': 246, '11000000': 239, '01000000': 275, '00000000': 264}

由于我们有 8 个量子位，这些结果对应的测量相位为：

import pandas as pd

rows, measured_phases = [], []
for output in counts:
    decimal = int(output, 2)  # 将二进制字符串转换为十进制
    phase = decimal / (2**n_count)  # 计算对应的本征值
    measured_phases.append(phase)
    # 将这些值添加到我们表格的行中：
    rows.append([f"{output}(bin) = {decimal:>3}(dec)", 
                 f"{decimal}/{2**n_count} = {phase:.2f}"])
# 在表格中打印行
headers=["输出", "相位"]
df = pd.DataFrame(rows, columns=headers)
print(df)

                         输出              相位
0  10000000(bin) =    128(dec)        128/256 = 0.50
1  11000000(bin) =    192(dec)        192/256 = 0.75
2  01000000(bin) =     64(dec)        64/256 = 0.25
3  00000000(bin) =      0(dec)        0/256 = 0.00

我们现在可以使用连分数算法来尝试找到 s 和 r 。Python 内置了此功能：我们可以使用 fractions 模块将浮点数转换为 Fraction 对象，我们可以使用 .limit_denominator() 方法获得最接近我们浮点数的分数。

from fractions import Fraction

# 获取最接近0.666的分数，分母小于15
approx_fraction = Fraction(0.666).limit_denominator(15)
print(approx_fraction)

2/3

r 必须小于 N，因此我们将最大分母设置为 15：

rows = []
for phase in measured_phases:
    frac = Fraction(phase).limit_denominator(15)
    rows.append([phase, f"{frac.numerator}/{frac.denominator}", frac.denominator])
# 以表格形式打印
headers=["相位", "分数", "猜测的r值"]
df = pd.DataFrame(rows, columns=headers)
print(df)

     相位   分数  猜测的r值
0    0.50   1/2      2
1    0.75   3/4      4
2    0.25   1/4      4
3    0.00   0/1      1

我们可以看到，两个测量特征值为我们提供了正确的结果： r=4，而且我们可以看到 Shor 算法有失败的可能性。这些不好的结果是因为 s = 0，或者 s 和 r 不是互质的，我们得到了 r 的一个因子而不是 r 本身。解决这个问题最简单的方法就是重复实验，直到我们得到一个令人满意的 r 结果。

4. 模幂运算

您可能已经注意到，通过重复 U 来创建 U^{2^j} 门的方法随 j 呈指数增长，并且不会产生多项式时间算法。我们希望有一种方法来创建算子：

U^{2^j}|y\rangle = |a^{2^j}y \bmod N \rangle

它随 j 多项式增长。幸运的是，高效计算：

a^{2^j} \bmod N

是可能的。经典计算机可以使用称为重复平方的算法来计算指数。在我们的例子中，由于我们只处理形式为 2^j 的指数，因此重复平方算法变得非常简单：

def a2jmodN(a, j, N):
    """ a^{2^j} (mod N) """
    for i in range(j):
        a = np.mod(a**2, N)
    return a
a2jmodN(7, 2049, 53)

47

如果可以在 Python 中实现高效算法，那么我们就可以在量子计算机上使用相同的算法。不幸的是，尽管模幂运算电路关于 j 呈多项式增长，但它们并不简单，而且是 Shor 算法的瓶颈。

5. 从周期查找到因式分解

并非所有因式分解问题都很困难；我们可以立即发现偶数，并知道它的一个因子是 2。事实上，有特定的标准来选择难以因式分解的数字，但基本思想是选择两个大质数的乘积。

通用的因式分解算法首先会检查是否有因式分解整数的捷径（数是偶数吗？数的形式是否为 N = a^b?），然后在最坏情况下使用 Shor 的周期查找算法。由于我们的目标是关注算法的量子部分，因此我们将直接跳到 N 是两个质数乘积的情况。

示例：分解 15

为了在少量量子位上看到因式分解的示例，我们将分解 15，我们都知道它是两个不太大的质数 3 和 5 的乘积。

N = 15

第一步是在 1 和 N-1 之间选择一个随机数 a ：

a=2

接下来我们快速检查它是否已经是 N 的非平凡因子：

from math import gcd 
gcd(a, N)

1

接下来，我们对 a = 2 和 N = 15 执行 Shor 的阶查找算法。请记住，我们测量的相位将是 s/r，其中:

a^r \bmod N = 1

并且 s 是 0 到 r-1 之间的随机整数。

print(df)

     相位   分数  猜测的r值
0    0.50   1/2      2
1    0.75   3/4      4
2    0.25   1/4      4
3    0.00   0/1      1

我们可以很容易地找到 r 的猜测值为4：

现在我们有了 r，我们可能可以用它来找到 N 的一个因子。因为：

a^r \bmod N = 1

那么：

(a^r – 1) \bmod N = 0

这意味着 N 必须整除 a^r-1。而且如果 r 也是偶数，那么我们可以写：

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

（如果 r 不是偶数，我们就无法进一步了，必须使用不同的 a 值重试）。那么 N 和 a^{r/2}-1 或 a^{r/2}+1 的最大公约数很可能是 N 的一个真因子 [2]：

r=4
guesses = [gcd(a**(r//2)-1, N), gcd(a**(r//2)+1, N)]
print(guesses)

[3, 5]

6. 参考文献

[1] Beauregard S. Circuit for Shor’s algorithm using 2n+ 3 qubits[J]. arXiv preprint quant-ph/0205095, 2002.

[2] Nielsen M A, Chuang I L. Cambridge Series on Information and the Natural Sciences[J]. Quantum Computation and Quantum Information, 2000

文档链接

案例代码