利用量子特性來質因數分解。

有錯再麻煩指正，感謝。

目標問題

RSA 是一種常見的加密演算法，使用的是非對稱加密演算法，其中會利用到兩個極大質數相乘，因此 RSA 的可靠度建立於對極大整數因數分解的難度。

Example ref

• Generate keys:
  • Pick two prime numbers $P, Q$, e.g., $P = 5, Q = 11$
  • $N = P\times Q = 55$
  • $Z = (P-1)\times(Q-1) = 40$
  • Pick $e$ such that $\gcd(e, Z) = 1$, e.g., $e = 13$
  • Pick $d$ such that $e\times d\equiv1(\text{mod }Z)$, e.g., $d = 37$
  • $\Rightarrow$ Public key: $(e, N)$; Private key: $(d, N)$
• Usage:
  • Encrypt: $E = M^e(\text{mod }N)\Rightarrow E = M^{13}(\text{mod }55)$
  • Decrypt: $M = E^d(\text{mod }N)\Rightarrow M = E^{37}(\text{mod }55)$

Several lemmas with examples

這邊使用了比較直觀會用到的 Lemma，且使用範例當作說明，不使用證明。

• Example 1: $N = 15$
  $4^2\equiv 1(\text{mod }15), \gcd(15, 3) = 3\mid 15, \gcd(15, 5) = 5 \mid 15$

• Example 2: $N = 21$
  $8^2\equiv 1(\text{mod }21), \gcd(21, 7) = 7\mid 21, \gcd(21, 9) = 3 \mid 21$

• Example 1: $N = 15$
  $2^4\equiv 1(\text{mod }N), \gcd(15, 3) = 3\mid 15, \gcd(15, 5) = 5 \mid 15$

• Example 2: $N = 21$
  $2^6\equiv 1(\text{mod }N), \gcd(21, 7) = 7\mid 21, \gcd(21, 9) = 3 \mid 21$

• Example: $N = 15, g = 2$
  $2^1\equiv 2^5\equiv 2^9\equiv 2(\text{mod }15); 2^2\equiv 2^6\equiv 2^{10}\equiv 4(\text{mod }15)$
  $\Rightarrow p = 4\Rightarrow 2^4\equiv 1(\text{mod }15)$
• Example: $N = 15, g = 7$
  $7^1\equiv 7^5\equiv 7^9\equiv 7(\text{mod }15); 7^2\equiv 7^6\equiv 7^{10}\equiv 4(\text{mod }15)$
  $\Rightarrow p = 4\Rightarrow 2^4\equiv 1(\text{mod }15)$

Solving by Quantum

加速效果 (對 $N$ 作質因數分解):

• 古典最佳: $O(e^{1.9(\log N)^{1/3}(\log\log N)^{2/3}})$ is sub-exponential time
• Shor’s algorithm: $O((\log N)^2(\log\log N)(\log\log\log N))$ is polynomial time

Convert to $U$

$f(x) = g^x (\text{mod } N)$, where $g, N\in\mathbb{Z}, g < N$, and $\gcd(g, N) = 1$
$g^r \equiv 1 (\text{mod } N)$
$U\lvert y\rangle\equiv \lvert g y (\text{mod } N)\rangle$
$U^k\lvert y\rangle\equiv \lvert g^ky (\text{mod } N)\rangle, 0\le k$

• Example: $N = 15, g = 7$
  • $U^0\lvert 1\rangle = \lvert 1\rangle$
  • $U^1\lvert 1\rangle = \lvert 7\rangle$
  • $U^2\lvert 1\rangle = \lvert 4\rangle$
  • $U^3\lvert 1\rangle = \lvert 13\rangle$
  • $U^4\lvert 1\rangle = U^r\lvert 1\rangle = \lvert 1\rangle$
  • $U^5\lvert 1\rangle = \lvert 7\rangle$

Eigenvector

$$\lvert \mu_s\rangle = \frac{1}{\sqrt{r}}\sum_{k = 0}^{r-1}
e^{-\frac{2\pi i s k}{r}}\lvert g^k (\text{mod } N)\rangle, 0\le s < r
$$

$$U\lvert \mu_s\rangle = e^{\frac{2\pi i s}{r}}\lvert \mu_s\rangle$$

則 eigenvalue 可推得為 $e^{\frac{2\pi i s}{r}}$，再經過 IQFT 將相位資訊提出，可得 $\varphi = \frac{s}{r}$。

• Example: $N = 15, g = 7$
  • $\lvert \mu_0\rangle = \frac{1}{2}(\lvert 1\rangle+\qquad\lvert 7\rangle+\qquad
    \lvert 4\rangle+\qquad
    \lvert 13\rangle)$
  • $\lvert \mu_1\rangle = \frac{1}{2}(\lvert 1\rangle+e^{-\frac{2\pi i}{4}}\lvert 7\rangle+e^{-\frac{4\pi i}{4}}
    \lvert 4\rangle+e^{-\frac{6\pi i}{4}}
    \lvert 13\rangle)$
  • $\lvert \mu_2\rangle = \frac{1}{2}(\lvert 1\rangle+e^{-\frac{4\pi i}{4}}
    \lvert 7\rangle+e^{-\frac{8\pi i}{4}}
    \lvert 4\rangle+e^{-\frac{12\pi i}{4}}
    \lvert 13\rangle)$
  • $\lvert \mu_3\rangle = \frac{1}{2}(\lvert 1\rangle+e^{-\frac{6\pi i}{4}}
    \lvert 7\rangle+e^{-\frac{12\pi i}{4}}
    \lvert 4\rangle+e^{-\frac{18\pi i}{4}}
    \lvert 13\rangle)$

Eigenvectors superposition

$$\frac{1}{\sqrt{r}}\sum_{s = 0}^{r-1}\lvert \mu_s\rangle = \lvert 1\rangle
$$

• Example: $N = 15, g = 7$

$$\frac{1}{\sqrt{r}}\sum_{s = 0}^{r-1}\lvert \mu_s\rangle =\lvert \mu_0\rangle +\lvert \mu_1\rangle+\lvert \mu_2\rangle+\lvert \mu_3\rangle
=\lvert 1\rangle$$

Use QPE to get eigenvalue

Eigenvector 輸入使用疊加態 (參考上一篇文章)。輸出時就有這些可能: $\lvert 2^t\frac{0}{r}\rangle,\lvert 2^t\frac{1}{r}\rangle,…, \lvert 2^t\frac{r-1}{r}\rangle$ 共 $r$ 種，$2^t$ 可以看作將小數變為整數的 scaling factor 表達時可以將其轉換成小數的表達方法即可，因此就可以透過統計的方式找出所求 $r$。另外可以看到 $\frac{0}{r} = 0$ 無法得知 $r$ 的資訊，因此有 $\frac{1}{r}$ 的機率會猜到 $0$ 進而無法推算出 $r$。

以上就是利用量子來解決質因數分解的方法。再來進入實作部分。

Qiskit implementation

因為程式已經大到我用 IBM Composer 無法簡單實作了，所以這邊改使用 IBM Qiskit 這個量子 python 工具以及從官方的範例 Shor’s Algorithm - IBM / src 下去修改。並且由於 Qubit 數量已經超過官方提供的 5 位，因此先使用模擬器來進行運算。

差異性: 官方使用的 counting bit 數量是 8-bit (原始版 Shor’s algorithm 需要兩個 register，一個, counting bit, 是 $2n$ bits，另一個, eigenvector, 是 $n$ bits，$n = \log N$。因為 $15$ 是 4-bit 整數，所以 counting bit 使用 8 位)。不過理論上 2-bit 就可以重現結果，但本篇為了說明大部分使用 4-bit。另外一些變數跟函式命名方式有針對本文所使用的符號進行修改。官方在 counting bit 使用 Little Endian 但 $U$ 的運算使用 Big Endian，本篇統一使用 Little Endian (高記憶體位置放 MSB)。

電路架構:

1. U circuit: 實作 $g^x (\text{mod } N)$
2. IQFT
3. QPE: 包含對 t-bit counting bits 做疊加以及 1. 2.

U circuit

這個 circuit 是我個人認為最複雜的區塊，因為它似乎沒有一個直觀的建立方法，且隨著 $g$ 或 $N$ 不一樣就需要重新建立，因此不太像一般古典運算的函式一樣可以修改輸入的變數即可。U circuit 要做的事就是產生一個 $g^x (\text{mod } N)$ 其中 $x$ 隨著 $t$ 增加而增加。為何這邊的 $U$ 需要使用 4-bit，原因是 $\text{mod } 15$ ($0, 1, …, 14$) 的所有數字需要用 4-bit 表示。

另外可以看到這裡的 gmod15_circ 是比較一般化 (general) 的形式，可以同時支援不同的 $g$ 作為輸入，不過仍是需要重新產生不同的電路。有趣的是根據設計可以發現在 $g = 4$ 時可以用到最少的邏輯閘 (只需要 2 個 CSWAP¹ 也就是 6 個 CCX²)。

這裡以 $t = 2$ (counting bit: $c[0], c[1]$), $g = 7$ 為例 (因為 counting bit 兩位就足夠)，真值表可寫成下方表格。$[i]$ 表示二進位第 $i$ 位。請搭配下方圖片查看中間結果。

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

def gmod15_circ(g, power):
    """Controlled multiplication by g mod 15"""
    if g not in [2, 4, 7, 8, 11, 13]:
        raise ValueError("'g' must be 2,4,7,8,11 or 13")
    U = QuantumCircuit(4)
    for _ in range(power):
        if g in [2,13]:
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
        if g in [7,8]:
            U.swap(2,3)
            U.swap(1,2)
            U.swap(0,1)
        if g in [4, 11]:
            U.swap(1,3)
            U.swap(0,2)
        if g in [7, 11, 13]:
            for q in range(4):
                U.x(q)
    U = U.reverse_bits() # reversed because IBM uses Big endian for U
    U = U.to_gate()
    U.name = f"{g}^{{2^{log2(power):.0f}}} mod 15"
    c_U = U.control()
    return c_U

IQFT

這部份我們之前討論過 (請參考前幾篇文章)，所以很容易可以實現。有兩種實作方式，IBM 這裡採用 SWAP 在前的版本。

38
39
40
41
42
43
44
45
46
47
48

def iqft_circ(n):
    """n-qubit QFTdagger (IQFT) the first n qubits in circ"""
    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 = "IQFT"
    return qc

QPE

將上述電路合併在一起。

51
52
53
54
55
56
57
58
59
60
61
62

def qpe_circ(g):
    qc = QuantumCircuit(4 + CBIT_NUM, CBIT_NUM)
    for q in range(CBIT_NUM):
        qc.h(q)  # Initialize counting qubits in state |+>
    qc.x(0 + CBIT_NUM)  # And auxiliary register in state |1>
    for q in range(CBIT_NUM):  # Do controlled-U operations
        qc.append(
            gmod15_circ(g, 2**q), [q] + [i + CBIT_NUM for i in range(4)]
        )
    qc.append(iqft_circ(CBIT_NUM), range(CBIT_NUM))  # Do inverse-QFT
    qc.measure(range(CBIT_NUM), range(CBIT_NUM))
    return qc

轉譯電路 Transpile

103
104
105
106
107
108
109
110
111
112
113
114

# Assign the backend
backend = Aer.get_backend("aer_simulator")

# Build circuit
qc = qpe_circ(g)

# Transpile circuit
t_qc = transpile(qc, backend)
t_qc.metadata = {"backend": backend}

# Circuit visualization
print(qc.draw(fold=-1, output="text"))

執行電路

要執行一個量子電路需要先進行 assemble (將資訊打包並產生 task ID) 接著再使用對應的 backend 執行。執行過一次 QPE (如下方程式碼) 算一次 Attempt。可以看到每個 Attempt 只用到一個 shot。若是用 1024 shot 統計可以發現 counting bit 只會出現四種結果 (‘0000’, ‘0100’, ‘1000’, ‘1100’)。

63
64
65
66
67
68
69
70
71
72

def run_qpe(t_qc):
    # Simulate Results
    qobj = assemble(t_qc, shots=1)
    # Setting memory=True below allows us to see a list of each sequential reading
    result = t_qc.metadata["backend"].run(qobj, memory=True).result()
    readings = result.get_memory()
    print("Register Reading: " + readings[0])
    phase = int(readings[0], 2) / (2**CBIT_NUM)
    print("Corresponding Phase: %f" % phase)
    return phase

運氣好只要一次 Attempt 就可以得到質因數結果，反之，要重複執行。因此會有一個迴圈去重複嘗試 run_qpe (第 81 行) 直到得到結果。

75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96

def solve(t_qc, N, g):
    factor_found = False
    attempt = 0
    while not factor_found:
        attempt += 1
        print("\nAttempt %i:" % attempt)
        phase = run_qpe(t_qc)  # Phase = s/r
        frac = Fraction(phase).limit_denominator(
            N
        )  # Denominator should (hopefully!) tell us r
        r = frac.denominator
        print("Result: r = %i" % r)
        if phase != 0:
            # Guesses for factors are gcd(x^{r/2} ±1 , 15)
            guesses = [gcd(g ** (r // 2) - 1, N), gcd(g ** (r // 2) + 1, N)]
            print("Guessed Factors: %i and %i" % (guesses[0], guesses[1]))
            for guess in guesses:
                if (
                    guess not in [1, N] and (N % guess) == 0
                ):  # Check to see if guess is a factor
                    print("*** Non-trivial factor found: %i ***" % 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
108
109
110
111
112
113
114
115
116
117
118
119
120

import numpy as np
from qiskit import QuantumCircuit, Aer, transpile, assemble
from qiskit.converters import circuit_to_dag
from math import gcd, log2
from fractions import Fraction

print("Imports Successful")
CBIT_NUM = 4  # The number of counting bits


def gmod15_circ(g, power):
    """Controlled multiplication by g mod 15"""
    if g not in [2, 4, 7, 8, 11, 13]:
        raise ValueError("'g' must be 2,4,7,8,11 or 13")
    U = QuantumCircuit(4)
    for _ in range(power):
        if g in [2, 13]:
            U.swap(0, 1)
            U.swap(1, 2)
            U.swap(2, 3)
        if g in [7, 8]:
            U.swap(2, 3)
            U.swap(1, 2)
            U.swap(0, 1)
        if g in [4, 11]:
            U.swap(1, 3)
            U.swap(0, 2)
        if g in [7, 11, 13]:
            for q in range(4):
                U.x(q)
    U = U.reverse_bits()  # reversed because IBM uses Big endian for U
    U = U.to_gate()
    U.name = f"{g}^{{2^{log2(power):.0f}}} mod 15"
    c_U = U.control()
    return c_U


def iqft_circ(n):
    """n-qubit QFTdagger (IQFT) the first n qubits in circ"""
    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 = "IQFT"
    return qc


def qpe_circ(g):
    qc = QuantumCircuit(4 + CBIT_NUM, CBIT_NUM)
    for q in range(CBIT_NUM):
        qc.h(q)  # Initialize counting qubits in state |+>
    qc.x(0 + CBIT_NUM)  # And auxiliary register in state |1>
    for q in range(CBIT_NUM):  # Do controlled-U operations
        qc.append(gmod15_circ(g, 2**q), [q] + [i + CBIT_NUM for i in range(4)])
    qc.append(iqft_circ(CBIT_NUM), range(CBIT_NUM))  # Do inverse-QFT
    qc.measure(range(CBIT_NUM), range(CBIT_NUM))
    return qc


def run_qpe(t_qc):
    # Simulate Results
    qobj = assemble(t_qc, shots=1)
    # Setting memory=True below allows us to see a list of each sequential reading
    result = t_qc.metadata["backend"].run(qobj, memory=True).result()
    readings = result.get_memory()
    print("Register Reading: " + readings[0])
    phase = int(readings[0], 2) / (2**CBIT_NUM)
    print("Corresponding Phase: %f" % phase)
    return phase


def solve(t_qc, N, g):
    factor_found = False
    attempt = 0
    while not factor_found:
        attempt += 1
        print("\nAttempt %i:" % attempt)
        phase = run_qpe(t_qc)  # Phase = s/r
        frac = Fraction(phase).limit_denominator(
            N
        )  # Denominator should (hopefully!) tell us r
        r = frac.denominator
        print("Result: r = %i" % r)
        if phase != 0:
            # Guesses for factors are gcd(x^{r/2} ±1 , 15)
            guesses = [gcd(g ** (r // 2) - 1, N), gcd(g ** (r // 2) + 1, N)]
            print("Guessed Factors: %i and %i" % (guesses[0], guesses[1]))
            for guess in guesses:
                if (
                    guess not in [1, N] and (N % guess) == 0
                ):  # Check to see if guess is a factor
                    print("*** Non-trivial factor found: %i ***" % guess)
                    factor_found = True


def main():
    N = 15
    g = 7

    # Assign the backend
    backend = Aer.get_backend("aer_simulator")

    # Build circuit
    qc = qpe_circ(g)

    # Transpile circuit
    t_qc = transpile(qc, backend)
    t_qc.metadata = {"backend": backend}

    # Circuit visualization
    print(qc.draw(fold=-1, output="text"))

    solve(t_qc, N, g)


if __name__ == "__main__":
    main()

$ cat /proc/version
Linux version 5.10.16.3-microsoft-standard-WSL2 (oe-user@oe-host) (x86_64-msft-linux-gcc (GCC) 9.3.0, GNU ld (GNU Binutils) 2.34.0.20200220) #1 SMP Fri Apr 2 22:23:49 UTC 2021

$ uname -a
Linux DESKTOP-M72RH5C 5.10.16.3-microsoft-standard-WSL2 #1 SMP Fri Apr 2 22:23:49 UTC 2021 x86_64 x86_64 x86_64 GNU/Linux

$ python --version
Python 3.8.10

$ pip list
Package              Version
-------------------- -----------
certifi              2022.5.18.1
cffi                 1.15.0
charset-normalizer   2.0.12
click                8.1.3
cryptography         37.0.2
Deprecated           1.2.13
dill                 0.3.5.1
idna                 3.3
mpmath               1.2.1
networkx             2.8.3
ntlm-auth            1.5.0
numpy                1.22.4
pbr                  5.9.0
pip                  20.0.2
pkg-resources        0.0.0
ply                  3.11
psutil               5.9.1
pycparser            2.21
python-constraint    1.4.0
python-dateutil      2.8.2
python-tsp           0.2.1
qiskit               0.36.2
qiskit-aer           0.10.4
qiskit-ibmq-provider 0.19.1
qiskit-ignis         0.7.1
qiskit-terra         0.20.2
requests             2.27.1
requests-ntlm        1.1.0
retworkx             0.11.0
scipy                1.8.1
setuptools           44.0.0
six                  1.16.0
stevedore            3.5.0
symengine            0.9.2
sympy                1.10.1
tabulate             0.8.9
tsplib95             0.7.1
tweedledum           1.1.1
urllib3              1.26.9
websocket-client     1.3.2
websockets           10.3
wheel                0.34.2
wrapt                1.14.1

執行結果

• 運氣較好:

Attempt 1:
Register Reading: 0100
Corresponding Phase: 0.250000
Result: r = 4
Guessed Factors: 3 and 5
*** Non-trivial factor found: 3 ***
*** Non-trivial factor found: 5 ***

• 運氣較差:

Attempt 1:
Register Reading: 0000
Corresponding Phase: 0.000000
Result: r = 1

Attempt 2:
Register Reading: 0000
Corresponding Phase: 0.000000
Result: r = 1

Attempt 3:
Register Reading: 0000
Corresponding Phase: 0.000000
Result: r = 1

Attempt 4:
Register Reading: 1000
Corresponding Phase: 0.500000
Result: r = 2
Guessed Factors: 3 and 1
*** Non-trivial factor found: 3 ***

結語

有些部分還是不能完全理解。尤其是 $U$ 的理論及實作，這裡留下一些我個人的疑問等待後續找到答案:

• $\mu$ 的部分，輸出似乎不是 $\mu$ 本身，而是四種狀態 (1, 4, 7, 13) 隨機出現。
• 無法直觀看出 QPE 中 phase kickback 的影響。

上述兩點已解決，請參考 Shor's Algorithm Part2。

References

• RSA (cryptosystem) - wiki
• How Quantum Computers Break Encryption | Shor’s Algorithm Explained - minutephysics youtube
• Shor’s algorithm - wiki
• Shor’s Algorithm - IBM
• QC — Period finding in Shor’s Algorithm - Jonathan Hui
• The seed_simulator doesn’t work for list of circuits #585
• Large Integer Factorization: Number Field Sieve vs Shor - Vladimir Silva