量子アルゴリズム① Shorの素因数分解アルゴリズム

現在，千葉大学大学院で符号暗号を研究しているM1です．
隔週ぐらいのペースで，研究関連？の分野を勉強したことをまとめていこうと思います．
＊そんなこんなでガバとかあったら指摘お願いします・・・

今回は，Shorの素因数分解アルゴリズムについてです．

概要

Shor, Peter W. "Algorithms for quantum computation: discrete logarithms and factoring." Proceedings 35th annual symposium on foundations of computer science. Ieee, 1994.

こちらは，現状の公開鍵暗号系でよく用いられるRSA暗号やElgamal暗号に基づくNP困難問題（素因数分解・離散対数問題）は，量子コンピュータによって多項式時間で解読されることを解説した論文です．

素因数分解のみに絞って，深入りできるところはして，そうでないところは実装して理解を深めます．

追記
大雑把に理論を知りたい方は　こちら
原論文を詳細に理解したい方は　こちら

記号など

＊は知っている方向けの補足です
$ \bullet $ $(a, b)$ $ \cdots $ 整数 $a$ と $b$ の最大公約数
$ \bullet $ $ p $ $ \cdots $ 奇素数
$ \bullet $ $ \mathbb{Z} $ $ \cdots $ 整数全体の集合
$ \bullet $ $ \mathbb{Z}/N \mathbb{Z} $ $ \cdots $ 0以上 $N - 1$ 以下の整数全体の集合で，和と積は $\mathrm{mod} \ N$ 下での通常の和と積に定義されているもの
＊$N$ を法とする剰余類環
$ \bullet $ $ 位数 $ $ \cdots $ $ x $ を $ \mathbb{Z}/N \mathbb{Z} $ の元としたとき， $ x^r \equiv 1 \ \mathrm{mod} \ N $ なる $ \mathbb{Z}/N \mathbb{Z} $ の元 $r$ が存在すれば，そのような $r$ のうち最小なもの
＊$(x, N) = 1$ （つまり，互いに素）であれば，Eulerの定理から { $ 1 \leq r \leq r - 1 $ $ | $ $ x^r \equiv 1 \ \mathrm{mod} \ N $ } なる集合は空でないことが分かり，この集合は有限（特に整列集合）であることから最小値が存在します（その値を$r$とすればいい）．
$ \bullet $ $ \lfloor x \rfloor $ $ \cdots $ $x$ 以下の最大の整数（$x$は実数）
$ \bullet $ $ [ a, b ] $ $ \cdots $ $ a $ 以上 $b$ 以下の実数全体の集合（$a$と$b$は実数）
＊$ ( a, b ) $ で$a$ より大きく$b$ より小さい実数全体の集合，などを表すこともありますが，今回は使いませんし，文脈で分かるかと

Shorの素因数分解アルゴリズム

Input : $N$($2$より大きい自然数)
Output : $N$のある素因数

i). ランダムに $ x \in \mathbb{Z} / (N - 1) \mathbb{Z} $ を選ぶ
ii). $ (x, N) \geq 2 $ かどうか then $ (x, N) $ を出力して終了 else iii). へ
iii). $ \mathrm{mod} \ N $ における $x$ の位数 $r$ を求める
iv). $r$ が偶数かつ $ x^{r/2} \not \equiv N - 1 \ \mathrm{mod} \ N $ かどうか
then $ (x^{r/2} - 1, N) $ を出力して終了
else i). に戻る

iii) が量子アルゴリズムを用いる箇所になります，それ以外の例えば $ (x, N) $ を求める部分はEuclidの互除法（古典）で求められます．

ちなみにですが，このアルゴリズムの正当性は次のように示されます．
$ (x^{r/2} - 1, N) $ が $N$ の素因数になっているかどうか，を判定すればいいですが，これは，$ 2 \leq (x^{r/2} - 1, N) $ を示せばOKです．
まず，$(x^{r/2} - 1, N)$ が $ N $ 以下の値として存在することはEuclidの互除法の構成から保証されます．
$ (x^{r/2} - 1, N) \not \equiv 0 \ \mathrm{mod} \ N $ について，そうでないとすると，$r$ の最小性に矛盾します．
$ (x^{r/2} - 1, N) \not \equiv 1 \ \mathrm{mod} \ N $ について，そうでないとすると，$ x^r - 1 = (x^{r/2} - 1)(x^{r/2} + 1) \equiv \ 0 \ \mathrm{mod} \ N $ ゆえ， $ x^{r/2} + 1 $ は $N$ の倍数となりますが，このことは $ x^{r/2} \not \equiv N - 1 \ \mathrm{mod} \ N $ という仮定に矛盾します．
以上よりOKです．
停止性は，特に任意の$x$でiii) が有限ステップか（ある$x$で無限ループにならないか）どうか気になるところですが，（一般的な構成をしていないので，今回は証明は省略しますが）これもOKです．

例

ここでは $ N = 35 $ とします．
(1) $ x = 8 $
$ 8^4 \equiv 1 \ \mathrm{mod} \ 35 $ ゆえ，$8$ の $\mathrm{mod} \ 35$ の位数は $4$ となる．特に偶数．
また，$ (8^2 - 1, 21) = (63, 21) = 7 $ で，$7$ は $21$ の素因数ゆえOK．
(2) $ x = 11 $
$ {11}^3 \equiv 1 \ \mathrm{mod} \ 35 $ ゆえ，$11$ の $\mathrm{mod} \ 35$ の位数は $3$ となるが，これは奇数ゆえ不適．
(3) $ x = 19 $
$ {19}^6 \equiv 1 \ \mathrm{mod} \ 35 $ ゆえ，$19$ の $\mathrm{mod} \ 35$ の位数は $6$ となる．特に偶数．
しかし，$ ({19}^3 + 1, 35) = (6860, 35) = 35 $ ゆえ，不適．

※厳密には，例えば，(1)で「$8$ の $\mathrm{mod} \ 35$ の位数は $4$ となる．」と言うためには，「$ 8^1 \not \equiv 1 \ \mathrm{mod} \ 35$ かつ $ 8^2 \not \equiv 1 \ \mathrm{mod} \ 35$ かつ $ 8^3 \not \equiv 1 \ \mathrm{mod} \ 35$」を示す必要がありますが，ここでは省略します（実際に計算すればいいので）．

$x$ の選びようによっては，上手くいったりいかなかったりって感じです．
実は，Nielsen-Chunang で，$N$の相異なる素因数の数を$m$個とするとき $　\displaystyle 1 - \frac{1}{2^m} $ 以上の確率で「上手くいく」$x$ が取れることが知られています．
＊RSA暗号では，大きい素数$2$つの合成数を考えることが多いので，その場合では $ \frac{3}{4} $ 以上の確率で「上手くいく」ことになります

モジュールを用いた実装

実装がめんどくさい方は，実は次のコード（ここでは $N = 21$）でi).〜iv).全てを実行できます．

Shor.py

from qiskit import BasicAer
from qiskit.aqua import QuantumInstance
from qiskit.aqua.algorithms import Shor

N = 21
shor = Shor(N)
backend = BasicAer.get_backend('qasm_simulator')
quantum_instance = QuantumInstance(backend, shots=1024)
result = shor.run(quantum_instance)
print(f"The list of factors of {N} as computed by the Shor's algorithm is {result['factors'][0]}.")

出力

The list of factors of 21 as computed by the Shor's algorithm is [3, 7].

実験とかするときに用いられることを想定しているようです．→参考
ただ，$ 4 (\lfloor \mathrm{log}_2 N \rfloor + 1) + 2 $ 量子ビットを用いるので，実行時間はかなり遅いです．
$ N = 15, 21 $ ぐらいであれば，普通に $ [2, \lfloor \sqrt{N} \rfloor ] $ の範囲の整数で素因数があるかを探索する方がよっぽど早いですし、手元の環境だと $N = 35$ は返ってくる気配がありませんでした・・・（$N = 21$でも１時間ぐらい）．

iii)の方針

理論的なことは飛ばして，実装メインです．Qiskit Textbook を参考にしています．
＊和訳
また，以下では，「上手くいく」 $ (N, x) = (35, 8)$ を例とします．
目標は以下のUnitary演算子で量子位相推定を行い， $r$ を算出することです．

方針

1. $ U | y \rangle \equiv | 8y \ \mathrm{mod} \ 35 \rangle $ として，$ y = 1 $ なる初期状態に，$ U $ を繰り返し作用させるような量子ゲートを作る
2. 固有状態を得るために，量子Fourier変換を作用させる
3. 上記で得られた値を測定することで位数を得る

って感じです．それでは順にコードを見ていきます．

初めに以下をインポートしています．

from qiskit import QuantumCircuit, Aer, transpile, assemble
# from numpy.random import randint
from math import gcd
import math
import numpy as np
from fractions import Fraction
print("Imports Successful")

出力

Imports Successful

1.$ U | y \rangle \equiv | 8y \ \mathrm{mod} \ 35 \rangle $ として，$ y = 1 $ なる初期状態に，$ U $ を作用させるような量子ゲートを作る

def c_xmod35(x: int, N: int, power: int) -> list:
    """Controlled multiplication by a mod 35"""
    gate_number = int(math.log2(N))
    U = QuantumCircuit(gate_number)        
    for iteration in range(power):
        U.swap(0, 1)
        U.swap(1, 2)
        U.swap(2, 3)
        U.cx(0, 1)
        U.cx(2, 0)
        U.cx(2, 3)
        U.cx(2, 4)
        U.cx(1, 0)
        U.cx(1, 4)
    U = U.to_gate()
    U.name = "%i^%i mod 35" % (x, power)
    c_U = U.control()
    return c_U

2.固有状態を得るために，量子Fourier変換を作用させる

def qft_dagger(n: int) -> list:
    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

3.上記で得られた値を測定することで位数を得る

def qpe_xmod35(x: int, N: int) -> int:
    n_count = 3
    gate_number = int(math.log2(N))
    qc = QuantumCircuit(gate_number+n_count, n_count)
    for q in range(n_count):
        qc.h(q)     
    qc.x(3+n_count) 
    for q in range(n_count): 
        qc.append(c_xmod35(x, N, 2**q), 
                 [q] + [i+n_count for i in range(gate_number)])
    qc.append(qft_dagger(n_count), range(n_count)) 
    qc.measure(range(n_count), range(n_count))
    qasm_sim = Aer.get_backend('qasm_simulator')
    t_qc = transpile(qc, qasm_sim)
    obj = assemble(t_qc, shots=1)
    result = qasm_sim.run(assemble(t_qc), memory=True).result()
    readings = result.get_memory()
    print("Register Reading: " + readings[0])
    phase = int(readings[0],2)/(2**n_count)
    print("Corresponding Phase: %f" % phase)
    return phase

今までのを全てまとめると次のようになります．

N = 35
factor_found = False
attempt = 0
while not factor_found:
    x = 8 # 本当は x = randint(2, N - 1)
    if gcd(x, N) == 1:
        attempt += 1
        print("\nAttempt %i:" % attempt)
        phase = qpe_xmod35(x, N) # Phase = s/r
        frac = Fraction(phase).limit_denominator(N) 
        r = frac.denominator
        if phase != 0 :
            guesses = [gcd(x**(r//2)-1, N), gcd(x**(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: 
                    print("*** Non-trivial factor found: %i ***" % guess)
                    factor_found = True

出力（↓がそのまま出力されるとは限りません）

Attempt 1:
Register Reading: 011
Corresponding Phase: 0.375000
Guessed Factors: 35 and 1

Attempt 2:
Register Reading: 110
Corresponding Phase: 0.750000
Guessed Factors: 7 and 5
*** Non-trivial factor found: 7 ***
*** Non-trivial factor found: 5 ***

まとめ

符号暗号（っていうか耐量子暗号）をやるなら，ある程度は敵（量子計算機）のことも知っておかないとな〜と思い，勉強してまとめてみました．
素因数分解なのに，位数を経由する，という発想が面白かったです．周期を求めるのに強い，という量子計算機の性質が上手く活かされているアルゴリズムだと感じました．
また，別の量子アルゴリズムも勉強しようと思います．

参考文献

$ \bullet $ 縫田光司. 耐量子計算機暗号. 森北出版 株式会社, 2020.