Shor

(Updated: )

我大概一辈子也不会忘记这个算法了

直接把报告贴上来罢

Shor 算法

Shor算法简介

RSA 难以破解的原因是大素数因子分解的困难性：经典素数分解算法（如试除法、费马分解法等）在处理大整数时需要指数级时间复杂度。

Shor算法在量子计算机上仅需多项式时间完成因数分解，其时间复杂度为 ，因此可以极大加速破解速度，使得 RSA 失效

量子傅里叶变换

DFT

经典的傅里叶变换形如：

其中

QFT

QFT

在量子算法中，量子傅里叶变换实现以下转换:

即：

相对应的酉矩阵为

具体为

整理 QFT 得出量子电路易实现的形式

将 展开，并将 代入

将 展开

这一步将 展开

电路实现

Hadamard 门可以实现 的功能

除此之外，一个 门为 , 是可控 门

其中

我们对两个量子比特 施加 门，有如下效果：

所以最终的电路为如下：

alt text

运行电路的具体过程如下：

先施加一个 门：

施加 门后得到:

施加所有的 门后得到:

由于

代入得到

类似地运行完剩余电路，对后 位进行变换，得到：

注意到此时得到的状态与上面我们推导出的状态在比特位上相反，所以要交换所有的比特来得到最终状态

模指

量子加法器

alt text

上面的 Conditional Phase Shift 模块就是之前的 门

那么我们可以如此构造一个电路，注意我们要先做 QFT, 这里的 是 QFT 后的第 位

这个电路可以记为:

alt text

所以，利用 QFT 和 IQFT， 我们可以得到一个量子加法器：

alt text

量子模加器

使用上面的量子加法器，可以将 QFT 变换后的状态 转化为

这样我们就可以来构造模加器的电路

具体地，我们先将值 加上一位 , 并整体做 QFT, 得到状态 ；之后施加一个加法器 ，得到状态为

注意这个 由于我们在做 QFT 前量子位多给了一位，所以不会溢出

之后施加一个反门 ，得到状态

如果 , 那么我们得到了想要的状态；

如果 , 那么我们需要将多减去的 再加回来

这时我们就可以先施加 IQFT 得到 ，再将最高位作为控制位输出到最下面的辅助位中，再做 QFT 恢复原本的状态

如果最高位为 , 说明 , 那么我们就施加一个可控加法器 得到状态

这里可控加法器的构造与使用 构造 的方法类似

此时还剩下最下面的辅助位没有还原，如果不还原辅助位，那么我们的门操作是不可逆的

所以我们按照前一部分类似的方法，将辅助位还原即可

最终得到了模加器的电路

量子模加器的电路图如下：

alt text

量子可控模乘器

有了模加器，我们就可以使用 个模加器构造如下电路：

alt text

这里 是控制位

如果 , 那么所有模加器都收到使能，可以正常运作；同时如果 , 那么第 个加法器两位都使能，就会将状态转换为加上 后的结果

由于 , 那么每一位 接在模加器 的控制位上，如果 , 就将结果加上 并对 取模

最后得到的结果即为

模指模块

有了模乘器，我们用两个模乘器和一个交换器即可完成最后的模指模块：

alt text

整个电路运行过程为：

这样在原本的位上得到了模指后的结果，并且将辅助位还原为

对于 模块，只需要重复 模块 次即可

对于 ，可以二分构造 模块，然后再合并，可以从合并 次变成合并 次，从而降低构造的复杂度

量子相位估计 QPE

设有一个酉算符 及其本征态 ，满足：

其中 是一个未知的实数，称为本征相位（eigenphase）

我们的目的是求出这个

首先准备初始态：

对于初态作用 门得到

依次对工作寄存器中的第 个比特（从低位到高位）控制地施加 到数据寄存器上，使系统变为：

由于
，代入后得到：

在工作寄存器上执行逆量子傅里叶变换 IQFT，能够将相位信息 转化为整数表示。

具体地，将 作用 IQFT, 由于 QFT 与 IQFT 表示成：

那么在前 位作用 IQFT 后有:

注意到求和项

当 , 此时 , 则

当 , 此时不同的指数项干涉，相位相消，结果为

因此最终状态趋近于 , 我们就可以测量得到 , 再根据这个 来估计

Shor 算法

约化素数因子分解

素数因子分解问题可以约化为 “求函数 的周期” 问题

具体地，假设一个合数 , 求

现随机取一个整数 , 且

定义阶数 是使得 满足的最小正整数

如果 是奇数，就换一个 ;

如果 是偶数，那么 可以变为

于是可以设

则解为：

如果这组解是一组平庸解，即

那么就换一个 , 否则就求出了

于是问题转化为求解 的周期

求解 的周期

有了上面的 QFT 模块和模指模块，我们就可以来构造求 的周期的电路了

具体地，对一个状态 施加模指电路后得到：

注意到 可以视为一个量子门，即对应着一个矩阵

我们想要知道这个矩阵的特征值

举一个 的例子：

可以看到 的阶数 , 而状态也在 次变换后回到

所以构造一个状态

则

所以 是矩阵的一个特征向量，对应的一个特征值为

进一步，我们构造 , 可以得到

这些特征值/特征向量十分特殊，因为他们包含了我们要求的周期：

这时我们将所有 相加，可以使得除了 的所有状态的相位抵消, 即

具体地，对于 :

由于 是这些特征向量的叠加态，所以，如果我们将 作为 做 QPE (量子相位估计)， 那么我们会得到一个相位估计值：

其中 是一个随机值

这样，我们对 做连分数即可估计出 和

如果 , 那么我们就得出了最终结果 ;

否则，结果为 , 不是我们想要的结果

但是，根据数论的一些定理，结果正确的概率为 , 所以期望在 步得出正确的

至此，Shor 算法的主要流程便介绍完毕了

电路实现

为了保证精确度，我们可以使用 个量子比特来做 QFT 和 QPE， 个量子比特来做模指

具体地，先作用 门，得到初态

然后将 个模指模块的控制位接到 个量子位，经过这一步原先状态变为

这里

由于 , 并且

所以原状态可以展开为：

那么在前 位作用 IQFT 后有:

注意到求和项

当 , 此时 , 则

当 , 此时不同的指数项干涉，相位相消，结果为

所以最终表达式为

这说明最终的状态趋近于 个 叠加，我们只需要测量就可以以 的等概率得到任意一个 , 进而可以得到

最终的线路图如下：

alt text

实验操作

实验方法

由于大数分解对应的 Shor 算法中的模指模块较难实现，并且真正的大数在 IBM 的量子计算机上分解比较困难；

相对来说，对于 的小样例模指模块较容易实现，并且真机上运行效果会更好，所以本次课题选择在真机上提交 的小样例来运行，同时在经典计算机上模拟 的小样例观察模拟结果

实验步骤

仿真

构建模指模块

模指模块 需要两个参数:

对于 , 我们选择的 需要与 互质，所以要先判断两者是否互质

由于 , 这与 位二进制数的最大值相同，所以对 取模实质上就是二进制数自然溢出的过程，所以只需要构造如下代码即可：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def U_a(a, N, power):
    if a == 1 or gcd(a, N) != 1:
        raise ValueError("gcd(a, N) is not 1")
    U = QuantumCircuit(QUBIT)
    for _iteration in range(power):
        if a in [2, 13]:
            U.swap(2, 3)
            U.swap(1, 2)
            U.swap(0, 1)
        if a in [7, 8]:
            U.swap(0, 1)
            U.swap(1, 2)
            U.swap(2, 3)
        if a in [4, 11]:
            U.swap(1, 3)
            U.swap(0, 2)
        if a in [7, 11, 13]:
            for q in range(4):
                U.x(q)
    U = U.to_gate()
    U.name = f"{a}^{power} mod N"
    c_U = U.control()
    return c_U

逆量子傅里叶变换 IQFT

对于 位的 IQFT, 先将所有位施加交换门 ，然后对于第 位，施加 个 门，最后加上 门，即可完成 QFT 电路的构造

1
2
3
4
5
6
7
8
9
10
def QFT(n):
    qc = QuantumCircuit(n)
    for qubit in range(n // 2):
        qc.swap(qubit, n - qubit - 1)
    for j in range(n):
        for m in range(j):
            qc.cp(-np.pi / float(2 ** (j - m)), m, j)
        qc.h(j)
    qc.name = "QFT†"
    return qc

量子相位估计 QPE

定义两个常数

创建一个 位的量子电路，并给前 位作用 门；

之后为后 位作用模指门，每个们控制位连接到前 位上；

再对前 位做 IQFT, 最后测量即可得到相位 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
QUBIT = 4
N_COUNT = 8

def get_phase(job):
    counts = job.result().get_counts()
    count_list = list(counts.items())
    count_list.sort(key=lambda x: x[1], reverse=True)
    print(count_list)
    print("Register Reading: " + count_list[0][0])
    phase = int(count_list[0][0], 2) / (2 ** N_COUNT) # phi = s / r
    print(f"Corresponding Phase: {phase}")
    return phase


def simulate(qc):
    aer_sim = Aer.Aer.get_backend('aer_simulator')
    job = aer_sim.run(transpile(qc, aer_sim), memory=True)
    return job
 

def QPE(a, N):
    qc = QuantumCircuit(QUBIT + N_COUNT, N_COUNT)
    for q in range(N_COUNT):
        qc.h(q)  # Initialize counting qubits in state |+>
    qc.x(N_COUNT)  # And auxiliary register in state |1>
    for q in range(N_COUNT):  # Do controlled-U operations
        qc.append(U_a(a, N, 2 ** q),
                  [q] + [i + N_COUNT for i in range(QUBIT)])
    qc.append(QFT(N_COUNT), range(N_COUNT))  # Do inverse-QFT
    qc.measure(range(N_COUNT), range(N_COUNT))
    # Simulate Results
    job = simulate(qc)
    return get_phase(job)

连分数计算周期

在得到相位 后，就可以对 做连分数得到

随后求出两个因子

如果两个因子是平凡的 , 就重新选取 ; 否则就得到了答案

代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def Shor(N):
    a = 1
    FACTOR_FOUND = False
    ATTEMPT = 0
    while not FACTOR_FOUND:
        a = a % (N - 1) + 1
        if a == 1 or gcd(a, N) != 1:
            continue
        ATTEMPT += 1
        print(f"\nATTEMPT {ATTEMPT}: a = {a}")
        phase = QPE(a, N)  # Phase = s/r
        frac = Fraction(phase).limit_denominator(N)
        r = frac.denominator
        print(f"Result: r = {r}")


        if phase != 0:
            guesses = [gcd(a ** (r // 2) - 1, N), gcd(a ** (r // 2) + 1, N)]
            print(f"Guessed Factors: {guesses[0]} and {guesses[1]}")
            for guess in guesses:
                if guess not in [1, N] and (N % guess) == 0:
                    print(f"*** Non-trivial factor found: {guess} ***")
                    FACTOR_FOUND = True

完整代码

最后的完整代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
import numpy as np
from qiskit import QuantumCircuit, transpile
import qiskit_aer as Aer
from math import gcd
from fractions import Fraction

QUBIT = 4
N_COUNT = 8

def U_a(a, N, power):
    if a == 1 or gcd(a, N) != 1:
        raise ValueError("gcd(a, N) is not 1")
    U = QuantumCircuit(QUBIT)
    for _iteration in range(power):
        if a in [2, 13]:
            U.swap(2, 3)
            U.swap(1, 2)
            U.swap(0, 1)
        if a in [7, 8]:
            U.swap(0, 1)
            U.swap(1, 2)
            U.swap(2, 3)
        if a in [4, 11]:
            U.swap(1, 3)
            U.swap(0, 2)
        if a in [7, 11, 13]:
            for q in range(4):
                U.x(q)
    U = U.to_gate()
    U.name = f"{a}^{power} mod N"


    c_U = U.control()
    return c_U


def QFT(n):
    qc = QuantumCircuit(n)
    for qubit in range(n // 2):
        qc.swap(qubit, n - qubit - 1)
    for j in range(n):
        for m in range(j):
            qc.cp(-np.pi / float(2 ** (j - m)), m, j)
        qc.h(j)
    qc.name = "QFT†"
    return qc


def get_phase(job):
    counts = job.result().get_counts()
    count_list = list(counts.items())
    count_list.sort(key=lambda x: x[1], reverse=True)
    print(count_list)
    print("Register Reading: " + count_list[0][0])
    phase = int(count_list[0][0], 2) / (2 ** N_COUNT)
    print(f"Corresponding Phase: {phase}")
    return phase


def simulate(qc):
    aer_sim = Aer.Aer.get_backend('aer_simulator')
    job = aer_sim.run(transpile(qc, aer_sim), memory=True)
    return job


def QPE(a, N):
    qc = QuantumCircuit(QUBIT + N_COUNT, N_COUNT)
    for q in range(N_COUNT):
        qc.h(q)  # Initialize counting qubits in state |+>
    qc.x(N_COUNT)  # And auxiliary register in state |1>
    for q in range(N_COUNT):  # Do controlled-U operations
        qc.append(U_a(a, N, 2 ** q),
                  [q] + [i + N_COUNT for i in range(QUBIT)])
    qc.append(QFT(N_COUNT), range(N_COUNT))  # Do inverse-QFT
    qc.measure(range(N_COUNT), range(N_COUNT))
    # Simulate Results
    job = simulate(qc)
    return get_phase(job)


def Shor(N):
    a = 1
    FACTOR_FOUND = False
    ATTEMPT = 0
    while not FACTOR_FOUND:
        a = a % (N - 1) + 1
        if a == 1 or gcd(a, N) != 1:
            continue
        ATTEMPT += 1
        print(f"\nATTEMPT {ATTEMPT}: a = {a}")
        phase = QPE(a, N)  # Phase = s/r
        frac = Fraction(phase).limit_denominator(N)
        r = frac.denominator
        print(f"Result: r = {r}")
        if phase != 0:
            # Guesses for factors are gcd(x^{r/2} ±1 , 15)
            guesses = [gcd(a ** (r // 2) - 1, N), gcd(a ** (r // 2) + 1, N)]
            print(f"Guessed Factors: {guesses[0]} and {guesses[1]}")
            for guess in guesses:
                if guess not in [1, N] and (N % guess) == 0:
                    # Guess is a factor!
                    print(f"*** Non-trivial factor found: {guess} ***")
                    FACTOR_FOUND = True


if __name__ == "__main__":
    Shor(15)

提交到真机

在 python 中绑定 IBM 用户与 api-token

注册 IBM 用户后获取 api-token, 即可提交任务

绑定用户并运行电路的代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
api_token = "..." # 在 IBM 官网上获取到 api-token

def account():
    QiskitRuntimeService.save_account(channel="ibm_quantum", 
        token=api_token, overwrite=True, set_as_default=True)
    service = QiskitRuntimeService(channel="ibm_quantum")
    backend = service.least_busy(operational=True, simulator=False, min_num_qubits=2)
    print(backend.name)
    return backend


def run(qc, backend):
    job = backend.run(transpile(qc, backend))
    print("任务已提交，任务 ID:", job.job_id())
    while True:
        job_status = job.status()
        # 如果任务完成，获取结果
        if not job.running():
            break
    return job

完整代码

最终的完整代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
import numpy as np
import time
from qiskit import QuantumCircuit, transpile
from qiskit_ibm_runtime import QiskitRuntimeService
import qiskit_aer as Aer
from math import gcd
from fractions import Fraction

QUBIT = 4
N_COUNT = 8
api_token = "546d02ef0ab61d9bcb1f6ae7e0c138059d122b1319057b6e1eb7fcea029cb
    04ccd072f9792890991fee88f2f0c66942f25809413fb95597a54e84ede989e7923"


def QFT(n):
    qc = QuantumCircuit(n)
    for qubit in range(n // 2):
        qc.swap(qubit, n - qubit - 1)
    for j in range(n):
        for m in range(j):
            qc.cp(-np.pi / float(2 ** (j - m)), m, j)
        qc.h(j)
    qc.name = "QFT†"
    return qc   


def U_a(a, N, power):
    if a == 1 or gcd(a, N) != 1:
        raise ValueError("gcd(a, N) is not 1")
    U = QuantumCircuit(QUBIT)
    for _iteration in range(power):
        if a in [2, 13]:
            U.swap(2, 3)
            U.swap(1, 2)
            U.swap(0, 1)
        if a in [7, 8]:
            U.swap(0, 1)
            U.swap(1, 2)
            U.swap(2, 3)
        if a in [4, 11]:
            U.swap(1, 3)
            U.swap(0, 2)
        if a in [7, 11, 13]:
            for q in range(4):
                U.x(q)
    U = U.to_gate()
    U.name = f"{a}^{power} mod N"
    c_U = U.control()
    return c_U

def get_phase(job):
    counts = job.result().get_counts()
    count_list = list(counts.items())
    count_list.sort(key=lambda x: x[1], reverse=True)
    print(count_list)
    print("Register Reading: " + count_list[0][0])
    phase = int(count_list[0][0], 2) / (2 ** N_COUNT)
    print(f"Corresponding Phase: {phase}")
    return phase



def simulate(qc):
    aer_sim = Aer.Aer.get_backend('aer_simulator')
    job = aer_sim.run(transpile(qc, aer_sim), memory=True)
    return job


def run(qc, backend):
    job = backend.run(transpile(qc, backend))
    print("任务已提交，任务 ID:", job.job_id())
    while True:
        job_status = job.status()
        # 如果任务完成，获取结果
        if not job.running():
            break
    return job

def QPE(a, N, backend):

    qc = QuantumCircuit(QUBIT + N_COUNT, N_COUNT)
    for q in range(N_COUNT):
        qc.h(q)  # Initialize counting qubits in state |+>
    qc.x(N_COUNT)  # And auxiliary register in state |1>
    for q in range(N_COUNT):  # Do controlled-U operations
        qc.append(U_a(a, N, 2 ** q),
                  [q] + [i + N_COUNT for i in range(QUBIT)])
    qc.append(QFT(N_COUNT), range(N_COUNT))  # Do inverse-QFT
    qc.measure(range(N_COUNT), range(N_COUNT))
    # Simulate Results
    # job = simulate(qc)
    job = run(qc, backend)
    return get_phase(job)


def Shor(N, backend):
    a = 1
    FACTOR_FOUND = False
    ATTEMPT = 0
    while not FACTOR_FOUND:
        a = a % (N - 1) + 1
        if a == 1 or gcd(a, N) != 1:
            continue
        ATTEMPT += 1
        print(f"\nATTEMPT {ATTEMPT}: a = {a}")
        phase = QPE(a, N, backend)  # Phase = s/r
        frac = Fraction(phase).limit_denominator(N)
        r = frac.denominator
        print(f"Result: r = {r}")
        if phase != 0:
            # Guesses for factors are gcd(x^{r/2} ±1 , 15)
            guesses = [gcd(a ** (r // 2) - 1, N), gcd(a ** (r // 2) + 1, N)]
            print(f"Guessed Factors: {guesses[0]} and {guesses[1]}")
            for guess in guesses:
                if guess not in [1, N] and (N % guess) == 0:
                    # Guess is a factor!
                    print(f"*** Non-trivial factor found: {guess} ***")
                    FACTOR_FOUND = True


def account():
    QiskitRuntimeService.save_account(channel="ibm_quantum", 
        token=api_token, overwrite=True, set_as_default=True)
    service = QiskitRuntimeService(channel="ibm_quantum")
    backend = service.least_busy(operational=True, simulator=False, min_num_qubits=2)
    print(backend.name)
    return backend


if __name__ == "__main__":
    backend = account()
    Shor(15, backend)

实验结果

仿真运行结果

仿真的运行结果如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
ATTEMPT 1: a = 2
[('00000000', 268), ('11000000', 257), ('01000000', 257), ('10000000', 242)]
Register Reading: 00000000
Corresponding Phase: 0.0
Result: r = 1

ATTEMPT 2: a = 4
[('00000000', 522), ('10000000', 502)]
Register Reading: 00000000
Corresponding Phase: 0.0
Result: r = 1

ATTEMPT 3: a = 7
[('01000000', 273), ('11000000', 269), ('10000000', 249), ('00000000', 233)]
Register Reading: 01000000
Corresponding Phase: 0.25
Result: r = 4
Guessed Factors: 3 and 5
*** Non-trivial factor found: 3 ***
*** Non-trivial factor found: 5 ***

从结果可以看到，第一次选择的模指的 , 即底数

此时测得相位 , 无法得到 ，因此更换

第二次 依然得不到结果，更换

此时测得相位 , 得到

由于 , 所以已经得到了两个素数因子，算法结束，结果符合预期

另一次测量结果如下：

1
2
3
4
5
6
7
8
ATTEMPT 1: a = 2
[('01000000', 263), ('11000000', 262), ('10000000', 250), ('00000000', 249)]
Register Reading: 01000000
Corresponding Phase: 0.25
Result: r = 4
Guessed Factors: 3 and 5
*** Non-trivial factor found: 3 ***
*** Non-trivial factor found: 5 ***

可以看到，这次测量中 算出了正确周期 , 直接得到了因子 和

真机运行结果

提交到真机后运行结果如下：

alt text

alt text

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
ibm_brisbane

ATTEMPT 1: a = 2
任务已提交，任务 ID: cwyqh74tdtng0087r2v0
[('00010100', 36), ('00000110', 36), ('00100110', 34), ('00011100', 34), 
('00011110', 33), ('01000100', 32), ('00001010', 31), ('00000111', 30), 
('01111101', 29), ('00110001', 29), ('00001100', 29), ('00000010', 28), 
('00010110', 28), ('01011000', 28), ('11010100', 28), ('00011101', 27), 
('00000011', 27), ('00111001', 27), ('00111101', 27), ('00111110', 27), 
('01011100', 26), ('00100001', 26), ('00001000', 25), ('01000110', 25), 
('00010101', 25), ('01001110', 24), ('11111101', 24), ('11111010', 24), 
('00010000', 24), ('00010010', 24), ('00001110', 24), ('01001000', 24), 
('00101111', 24), ('11110111', 24), ('00010001', 24), ('00101001', 24), 
('01011110', 24), ('11101010', 23), ('00100011', 23), ('01010000', 23), 
('01001010', 23), ('11001010', 23), ('00101010', 23), ('00011000', 22), 
('00101011', 22), ('11111111', 22), ('11110110', 22), ('00001111', 22), 
('00110101', 22), ('00111000', 22), ('11000100', 22), ('11111011', 22), 
('11011101', 22), ('00100101', 22), ('00011111', 21), ('01101011', 21), 
('01001111', 21), ('01000111', 21), ('01001101', 21), ('00011001', 21), 
('00000000', 21), ('01110101', 21), ('11011000', 20), ('00001011', 20), 
('11101101', 20), ('01001100', 20), ('01010110', 20), ('00110011', 20), 
('11110101', 20), ('01000011', 20), ('00111111', 20), ('01011101', 20), 
('00101100', 19), ('00100111', 19), ('00101101', 19), ('00000100', 19), 
('00100000', 19), ('01111110', 19), ('11010110', 19), ('01100010', 19), 
('00010111', 19), ('00110010', 19), ('10010000', 18), ('11011100', 18), 
('11110011', 18), ('00000001', 18), ('00101110', 18), ('00100100', 18), 
('00011011', 18), ('01100111', 18), ('00110000', 18), ('11011010', 18), 
('10110011', 17), ('01010101', 17), ('11111110', 17), ('01000010', 17), 
('10110101', 17), ('00101000', 17), ('01110001', 17), ('10011110', 17), 
('11101011', 16), ('11101111', 16), ('11110100', 16), ('01011001', 16), 
('10010100', 16), ('01111100', 16), ('11001110', 16), ('11010001', 16), 
('11110001', 16), ('11010000', 16), ('00110110', 16), ('00110111', 16), 
('00011010', 16), ('10101011', 16), ('11101110', 16), ('11000010', 16), 
('01111001', 16), ('11100110', 15), ('00111100', 15), ('01011010', 15), 
('01100001', 15), ('11111100', 15), ('10111100', 15), ('01010010', 15), 
('01010011', 15), ('11100100', 15), ('00110100', 15), ('00111011', 15), 
('01000000', 15), ('00111010', 15), ('11111001', 15), ('01110110', 15), 
('11001100', 15), ('01101110', 15), ('11110000', 15), ('01111000', 14), 
('01101111', 14), ('11000110', 14), ('00010011', 14), ('10010111', 14), 
('11010010', 14), ('10100011', 14), ('10011100', 14), ('10010110', 13), 
('11010101', 13), ('01100011', 13), ('10000111', 13), ('11011001', 13), 
('11100111', 13), ('10111101', 13), ('01111111', 13), ('11011110', 13), 
('10110001', 13), ('01110011', 13), ('01000001', 13), ('11101001', 13), 
('10011001', 13), ('10000110', 13), ('10000000', 13), ('11110010', 13), 
('11011111', 13), ('01010100', 13), ('10111001', 13), ('11001101', 12),
('10101101', 12), ('11001111', 12), ('00000101', 12), ('01011111', 12), 
('01110111', 12), ('11100010', 12), ('01101100', 12), ('10111010', 12), 
('11101000', 12), ('11101100', 11), ('11000000', 11), ('10111110', 11), 
('10111011', 11), ('10001100', 11), ('10111111', 11), ('11010111', 11), 
('01001011', 11), ('10101110', 11), ('10100000', 11), ('10001000', 11), 
('10001011', 11), ('01110100', 11), ('10011010', 11), ('10001010', 11), 
('01000101', 11), ('01100000', 10), ('11011011', 10), ('10000101', 10), 
('01101101', 10), ('00001101', 10), ('01110000', 10), ('01010001', 10), 
('00001001', 10), ('01100110', 10), ('11100000', 10), ('10110100', 10), 
('10000011', 10), ('10000010', 10), ('01001001', 9), ('10110010', 9), 
('10010010', 9), ('10011101', 9), ('01011011', 9), ('11000111', 9), 
('01100100', 9), ('10100100', 9), ('10110111', 9), ('00100010', 9), 
('11000011', 9), ('11100001', 9), ('11111000', 9), ('10001111', 8), 
('01010111', 8), ('11001000', 8), ('10100010', 8), ('01110010', 8), 
('01101010', 8), ('11010011', 8), ('10101010', 8), ('10011000', 8), 
('01111010', 8), ('10010101', 8), ('10110110', 8), ('11001001', 8), 
('11100011', 8), ('10100001', 8), ('10101111', 8), ('10011011', 8), 
('10111000', 8), ('11100101', 7), ('10011111', 7), ('01100101', 7), 
('10001101', 7), ('10010001', 7), ('10100111', 7), ('01111011', 7), 
('01101001', 7), ('10101000', 7), ('10110000', 7), ('10001110', 6), 
('10001001', 6), ('10101001', 6), ('10100110', 6), ('10010011', 6), 
('10000100', 5), ('10000001', 5), ('10101100', 4), ('11001011', 4), 
('01101000', 4), ('11000101', 4), ('11000001', 3), ('10100101', 1)]
Register Reading: 00010100
Corresponding Phase: 0.078125
Result: r = 13
Guessed Factors: 3 and 5
*** Non-trivial factor found: 3 ***
*** Non-trivial factor found: 5 ***

可以看到，这次测量得到了错误的周期 , 但是最后也得到了正确答案

这说明真机上误差较大，容易计算出错误的周期

另外一次测量结果如下:


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
ibm_kyiv

ATTEMPT 1: a = 2
任务已提交，任务 ID: cxk8f3mfdnwg008sqdh0
[('11000000', 35), ('10110110', 33), ('11000011', 31), ('10111110', 31), 
('10100000', 30), ('11101011', 30), ('11000101', 29), ('10101001', 28), 
('11100111', 28), ('10101000', 27), ('10111000', 27), ('10100001', 27), 
('10110001', 27), ('10110100', 27), ('11001110', 27), ('10100010', 27), 
('10000100', 26), ('11100011', 26), ('11100101', 26), ('10101101', 26), 
('11001001', 26), ('11010100', 26), ('10100111', 26), ('10100100', 26), 
('11101001', 25), ('11000010', 25), ('10100110', 25), ('10010101', 25), 
('11110100', 24), ('11010011', 24), ('11001000', 24), ('00001001', 24), 
('11110010', 24), ('10000111', 24), ('10101010', 24), ('10000001', 24), 
('11100001', 24), ('10111100', 23), ('10110101', 23), ('11010101', 23), 
('11100110', 23), ('10110010', 22), ('00001010', 22), ('11101101', 22), 
('10011000', 22), ('10101100', 22), ('10000110', 21), ('11100010', 21), 
('11000001', 21), ('00000010', 21), ('10101110', 21), ('01110110', 21), 
('10101111', 21), ('11001010', 21), ('11011011', 21), ('10100011', 21), 
('11010001', 21), ('10110011', 21), ('10101011', 21), ('10111011', 20), 
('10111001', 20), ('01100001', 20), ('01111000', 20), ('00001000', 20), 
('00100010', 20), ('11000100', 20), ('10111111', 20), ('10010100', 19), 
('01010100', 19), ('10111010', 19), ('01101011', 19), ('11101111', 19), 
('10010011', 19), ('00100001', 19), ('11110101', 19), ('11001011', 18), 
('10011100', 18), ('11011001', 18), ('10001000', 18), ('10011101', 18), 
('11000110', 18), ('10000000', 18), ('11011100', 18), ('10111101', 18), 
('11000111', 18), ('10001010', 18), ('10011111', 17), ('10011010', 17), 
('11010010', 17), ('00101111', 17), ('11011010', 17), ('11110111', 17), 
('01100100', 17), ('11001111', 17), ('01100110', 17), ('11101000', 17), 
('01001110', 17), ('10001110', 17), ('10001111', 17), ('11100000', 17), 
('00100111', 17), ('01110100', 17), ('01001111', 17), ('10000011', 17), 
('00100110', 17), ('10000101', 17), ('10010010', 16), ('10011011', 16), 
('11010000', 16), ('10010000', 16), ('10011110', 16), ('01100101', 16), 
('11111111', 16), ('00100100', 16), ('01111110', 16), ('10110111', 16), 
('11011111', 16), ('01111101', 16), ('11011101', 16), ('11101110', 16), 
('11011110', 16), ('00000000', 16), ('01111001', 15), ('01110000', 15), 
('10001011', 15), ('00100000', 15), ('00101110', 15), ('11001100', 15), 
('10001100', 15), ('00100011', 15), ('00101100', 15), ('10110000', 15), 
('01100111', 15), ('01101101', 15), ('11011000', 15), ('10010110', 15), 
('00101011', 14), ('00000111', 14), ('11110011', 14), ('01110010', 14), 
('11111000', 14), ('10001001', 14), ('00101000', 14), ('00000110', 14), 
('00000011', 14), ('11111010', 14), ('01011001', 14), ('10010111', 14), 
('10011001', 14), ('01101111', 14), ('01001101', 14), ('00010101', 14), 
('11110000', 14), ('01111011', 14), ('01111010', 13), ('00011111', 13), 
('10100101', 13), ('01101110', 13), ('01000101', 13), ('11111100', 13), 
('00001100', 13), ('00000101', 13), ('01000011', 13), ('11110001', 13), 
('00100101', 13), ('01001010', 13), ('01101100', 13), ('01000110', 13), 
('01011111', 13), ('01001100', 13), ('00001101', 12), ('11100100', 12), 
('10010001', 12), ('11101100', 12), ('11001101', 12), ('01100000', 12), 
('01101001', 12), ('11101010', 12), ('01100010', 12), ('00101001', 12), 
('11111101', 12), ('01110111', 12), ('00000100', 12), ('11110110', 12), 
('01110001', 12), ('00010000', 12), ('01000111', 12), ('10000010', 12), 
('01010101', 11), ('01110101', 11), ('10001101', 11), ('11111011', 11), 
('01011101', 11), ('01000010', 11), ('01010001', 11), ('11010110', 11), 
('01111100', 11), ('00001111', 11), ('01101000', 11), ('11010111', 11), 
('00010001', 11), ('01001000', 11), ('01010111', 10), ('11111001', 10), 
('00000001', 10), ('01000000', 10), ('00110111', 10), ('11111110', 10), 
('00010010', 10), ('00111110', 10), ('01001001', 10), ('00001110', 10), 
('00011100', 10), ('00111011', 9), ('01011100', 9), ('01101010', 9), 
('00111100', 9), ('01011110', 9), ('01100011', 9), ('01111111', 9), 
('01000100', 9), ('00011101', 9), ('00111101', 8), ('01010110', 8), 
('00110001', 8), ('00101010', 8), ('00110101', 8), ('00101101', 8), 
('01001011', 8), ('00110100', 7), ('01010011', 7), ('01010000', 7), 
('01110011', 7), ('00011110', 7), ('00011000', 7), ('00010110', 7), 
('00010111', 7), ('01011010', 7), ('00010100', 6), ('00110110', 6), 
('00011001', 6), ('00111111', 6), ('00110010', 5), ('00011010', 5), 
('00111000', 5), ('00110011', 5), ('01011000', 5), ('01000001', 5), 
('00110000', 5), ('00001011', 4), ('00111010', 4), ('01011011', 4), 
('01010010', 4), ('00011011', 2), ('00010011', 2), ('00111001', 2)]
Register Reading: 11000000
Corresponding Phase: 0.75
Result: r = 4
Guessed Factors: 3 and 5
*** Non-trivial factor found: 3 ***
*** Non-trivial factor found: 5 ***

可以看到, 这次运行结果在 11000000 这个状态处得到了最大频率

于是得到了正确的周期 , 并求出了两个因子

未来展望

  • 安全性

  • 有关 问题 (Bounded-error Quantum Polynomial-time)

    质因数分解属于 , 从经典的角度看, Shor 算法在多项式时间解决了 问题

    那么研究 或许对于 有帮助?

参考文献

[1]: P. Shor (1997), Polynomial-time algorithms for prime factorization
and discrete logarithms on a quantum computer, SIAM J. Comp., 26,
pp. 1484-1509. Also on quant-ph/9508027

[2]: St´ephane Beauregard, Circuit for Shor’s algorithm using 2n+3 qubits. Also on quant-ph/0205095

[3]: Unathi Skosana and MarkTame, Demonstration of Shor’s factoring algorithm for N = 21 on IBM quantum processors, Scientific Reports volume 11, Article number: 16599 (2021)

[4]: Qiskit textbook: quantum-fourier-transform. Also on Github

[5]: Qiskit textbook: shor. Also on Github