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

11 min read読了の目安(約10000字

現在,千葉大学大学院で符号暗号を研究しているM1です.
月1ぐらいのペースで,研究関連?の分野を勉強したことをまとめていこうと思います.
理論と実装どっちもやりますが,基本的には実装メインで書いたQiitaの記事の理論部分の補強をメインで書いていきます.たまに自分の研究分野の理論を簡単に実装を詳細に紹介するスタイルをとると思います.
*そんなこんなでガバとかあったら指摘お願いします・・・

今回は,Shorの素因数分解アルゴリズムについてです.
下記のQiitaの記事の理論部分を詳しくしたって感じです.

https://qiita.com/0917Laplace/items/982f171d023ef2fe6320

(アルゴリズムの正当性とか上で書かんくても良かったな…)

追記
原論文をより深く理解したい方はこちらを

https://github.com/awlaplace/Quantum/blob/master/Shor_factorization/Shor.pdf

概要

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 整数 ab の最大公約数
\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 のうち最小なもの
xN が互いに素であれば,x\mathrm{mod} \ N での位数が存在することは次のように分かります
(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 以下の実数全体の集合(abは実数)
( a, b )a より大きくb より小さい実数全体の集合,などを表すこともありますが,今回は使いませんし,文脈で分かるかと

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

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

i). ランダムに x \in \mathbb{Z} / (N - 1) \mathbb{Z} \wedge x \neq 0 を選ぶ
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). に戻る

i) は x = 0 だと議論が無意味なので省きます.
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} + 1N の倍数となりますが,このことは 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 で,721 の素因数ゆえ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でも1時間ぐらい).

iii)の方針

実装はpythonのqiskitというモジュールを使います
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 を作用させるような量子ゲートを作る

理論

まず, 2N^2 \leq q < 4N^2 の範囲で2べきの形で表せる q をとります.
q = 2^{\lceil 2 \mathrm{log}_2 N \rceil + 1} なので,このような q は一意的に存在します.
今回なら, N = 35 なので, 2,450 = 2N^2 \leq q \leq 4N^2 = 4,900 の範囲で2べきのものを考えると, q = 4,096 となります.

このとき, a \in \mathbb{Z}/ q \mathbb{Z} を取って,

\displaystyle \frac{1}{q^{1/2}} \sum_{a = 0}^{q - 1} | a \rangle

なる状態を考えます.
a \in [0, q - 1] で各状態の係数は1なので, a を動かしたときの重ね合わせを考えると,振幅に更に q^{1/2} の重みがかかります
次に, a \in [0, q - 1] で別の状態

\displaystyle \frac{1}{q^{1/2}} \sum_{a = 0}^{q - 1} | x^a \ \mathrm{mod} \ N \rangle

を考えます.
そして,上記の状態の合成を考えると,

\displaystyle \frac{1}{q^{1/2}} \sum_{a = 0}^{q - 1} | a \rangle \otimes | x^a \ \mathrm{mod} \ N \rangle

となります.

実装

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変換を作用させる

理論

状態 | a \rangle に対して量子Fourier変換を作用させると,

\displaystyle \frac{1}{q} \sum_{c = 0}^{q - 1} \left( \sum_{a = 0}^{q - 1} \mathrm{exp} \left( 2 \pi \sqrt{-1} \frac{ac}{q} \right) \right) | c \rangle \otimes | x^a \ \mathrm{mod} \ N \rangle

になります.

実装

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.上記で得られた値を測定することで位数を得る

理論

状態 (|c \rangle, | x^k \ \mathrm{mod} \ N \rangle) を観測する確率を考えます.詳細は省略しますが,この確率は少なくとも \displaystyle \frac{1}{3 r^2} 以上です.
そして, r を得る確率は(これも詳細は省略しますが)少なくとも \displaystyle \frac{\phi(r)}{3 r} 以上になります.
*上2つの詳細は論文に書かれています(行間もある上に細々とした議論ですが)

実は

\displaystyle \underset{r \rightarrow \infty}{\mathrm{lim}} \frac{\phi(r) \mathrm{log} \ \mathrm{log} r}{r} = e^{- \gamma}

であることが知られています.ここで, \gammaEuler定数です.
よって, O(\mathrm{log} \ \mathrm{log} r) の計算量で求めることができます.

実装

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 ***

まとめ

符号暗号(っていうか耐量子暗号)をやるなら,ある程度は敵(量子計算機)のことも知っておかないとな〜と思い,勉強してまとめてみました.ちなみに,今やっている研究にSimonのアルゴリズムなどが関係してくるっぽいので,そちらも勉強したらまとめると思います(どこに?).
今回は実装が先にあったので,そこに合わない原論文の理論的な部分(と議論が面倒な部分)は省略してしまいましたが,興味のある方はぜひ論文を読んでみてください.

参考文献

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