Factoring（ショアの素因数分解） 〜Quantum Zooやっていく〜

July 17, 2021, 9:16 a.m. edited Feb. 2, 2023, 11:53 a.m.

#Quantum Zoo  #量子情報  #量子力学 

$$ \def\bra#1{\mathinner{\left\langle{#1}\right|}} \def\ket#1{\mathinner{\left|{#1}\right\rangle}} \def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}} $$

なんだかんだ量子¹アルゴリズム、 Grover と Shor しか知らないので、 Quantum Zoo をやっていく突貫企画。 多分すぐ更新は途絶える

ということで初回はお馴染み Factoring （素因数分解）

ここに訳を、と思ったが、なんか昔の翻訳記事は翻訳に許可とったりしてて、私はそんな手間をとりたくないので、英語読んで。

2023-01-31 追記：以下の項目に（非推奨）が多くあるのは、位数発見アルゴリズムが Shor の素因数分解アルゴリズムのベースとなっていることを知らない頃に書いた項目であるため。現在は（非推奨）のついていない項目だけ読めば良いはず。

背景

\(n\) ビットの整数を素因数分解する最速の古典アルゴリズムは \(2^{\tilde{O}(n^{1/3})}\) であったが、 Peter Shor が開発した量子アルゴリズムにより \(\tilde{O}(n^3)\) となった²。すごい

問題

ある整数 \(N\) が与えられたとき、それを割り切ることのできる整数 \(p\in[1,N-1]\) を見つける。

アルゴリズム

位数発見アルゴリズム以外は古典アルゴリズムのみとなっている。

1. そもそも \(N\) が偶数なら 2 を返して終わり
2. ランダムに整数 \(a\in[2,N-1]\) を選ぶ。そして、最大公約数 \(\gcd(a, N)\) を計算⁴する。これでもしも 1 ではない最大公約数が見つかれば、その値を返して終わり
3. \(a^r\equiv 1\pmod N\) を満たす位数 \(r\) を探す位数発見アルゴリズムを実行する
4. 得られた \(r\) がもしも奇数であれば、再び 2. からやり直し。また、 \(a^{r/2}\equiv -1\pmod N\) を満たしてしまっていても、 2. からやり直し。問題なければ、 \(\gcd(a^{r/2}\pm 1, N)\) が答えなので、これを返す

これで動く理由は「なぜ」を参照のこと。

Qiskit による実装

位数発見アルゴリズムの order_finding.pyを利用しているので非常にシンプル。コードの全体は GitHub にある。

import random
import math
from typing import Optional

import numpy as np
import matplotlib.pyplot as plt
from qiskit import QuantumCircuit, Aer, transpile, assemble
from qiskit.visualization import plot_histogram
from sympy import Rational
from sympy.ntheory.continued_fraction import continued_fraction, continued_fraction_convergents

from qft import qft
from elementary import ax_modM
from order_finding import order_finding


def shor(N: int, show_hist: bool = True):
    """Shor's factoring algorithm: given $N\in\mathbb{Z}$, it finds a prime factor of $N$.
    Args:
        N (int): $N$
    Returns:
        A prime factor of $N$
    """

    if N % 2 == 0:
        return 2

    for _ in range(N):
        a = random.randint(2, N - 1)
        gcd = math.gcd(a, N)
        if gcd != 1:
            return gcd

        r = order_finding(x=a, N=N, show_hist=show_hist)

        if r % 2 == 1:
            continue
        if (x ** (r // 2) + 1) % N == 0:
            continue

        factor_candidate = math.gcd(x ** (r // 2) - 1, N)
        if N % factor_candidate == 0:
            return factor_candidate
        factor_candidate = math.gcd(x ** (r // 2) + 1, N)
        if N % factor_candidate == 0:
            return factor_candidate

    raise Exception('A prime factor is Not found!')


if __name__ == '__main__':
    print(shor(N=9))

実際に \(N=9\) で動かすと、 3 と返してくれる。以上。

（位数発見アルゴリズムがベースであることを知らない頃に書いた量よりもずっと少なく書けてびっくりしてる）

アルゴリズム (非推奨)

2023-01-31 追記：ここから先は今見るとちょっと微妙なので↑までの方法が良いと思う。ただ、「なぜ」と「計算量」は使えるはず

[P. Shor, 1996] のアルゴリズムに従いたかったが、提案当初、という感じでなんか難しかったので、 ここ に載ってるものにしたがう³（ただ、このサイト、一部で \(a^r\) が \(ar\) になっちゃってるね・・）。

古典パート（非推奨）

ランダムに整数 \(a\in[2,N-1]\) を選ぶ。そして、最大公約数 \(\gcd(a, N)\) を計算⁴する。これでもしも 1 ではない最大公約数が見つかれば古典だけで終わり。そうでない場合は、以下の量子パートに進む。

量子パート（非推奨）

$$f(x):=a^x \pmod N$$

と定義する⁵。このとき、

$$f(x+r)=f(x)$$

を満たす \(r\) を探す。はじめに

$$\frac{1}{\sqrt{N}}\sum_{x=0}^{N-1}\ket{x}\ket{0}$$

という状態を用意する。それから \(f\) を計算する。

$$\frac{1}{\sqrt{N}}\sum_{x=0}^{N-1}\ket{x}\ket{f(x)}$$

ここに量子フーリエ変換⁶ \(U_{QFT}\) を適用する。

$$\frac{1}{N}\sum_{x=0}^{N-1}\sum_{y=0}^{N-1}e^{2\pi ixy/N}\ket{y}\ket{f(x)}$$

ここで、 2 つのレジスタ（\(\ket{y}\) と \(\ket{f(x)}\)）を測定して \(y\) および \(f(x_0)\) (\(x_0\in[0,N-1]\)) が得られる確率は

$$\frac{1}{N^2}\left|\sum_{x\ \mathrm{where}\ f(x)=f(x_0)}e^{2\pi ixy/N}\right|^2=\frac{1}{N^2}\left|\sum_b e^{2\pi i(x_0+br)y/N}\right|^2$$

となる。つまり、 \(f\) は \(r\) ごとに同じ値をとる周期関数なので、共通の \(y\) および \(f(x)=f(x_0)\) となる項が複数存在し、それらは \(r\) を何回掛けて \(x_0\) に足すかという形で整数 \(b\) を用いて右辺のように表される。このとき、 \(ry/N\) が整数であれば

$$e^{2\pi i(x_0+br)y/N}=e^{2\pi i(x_0+(b+1)r)y/N}$$

となることから確率は大きくなる（つまり、打ち消し合わないで済むということ）ので、大体それを満たす \(y\) が観測される。さらに、この観測された \(y\) を用いて \(y/N\) を計算できるが、これは整数 \(m\) を用いると

$$\frac{ry}{N}=m$$ $$\frac{y}{N}=\frac{m}{r}$$

となるため、つまり、 \(y/N\) を既約分数にすると、その分母に \(r\) が出てくることがあるということ。すごい！！・・とは言っても、そんな偶然に頼りたくないので、今回 (\(y=y_1\)) は \(m=m_1\) 、もう一回 (\(y=y_2\)) やったら \(m=m_2\) が得られたとする。すると、

$$y_1=\frac{Nm_1}{r},\ y_2=\frac{Nm_2}{r}$$

より、 \(m_1\), \(m_2\) が互いに素であれば⁷ \(y_1\) と \(y_2\) の最大公約数が \(\frac{N}{r}\) となるので \(r\) が求まる。

再びの古典パート（非推奨）

こうして得られた \(r\) がもしも奇数だったら、もう一回最初の古典パートからやり直し。また、 \(a^{r/2}\equiv -1\pmod N\) を満たしてしまっていても、もう一回最初の古典パートからやり直し。

大丈夫であれば、

$$\gcd(a^{r/2}\pm 1, N)$$

が答えとなる。

なぜ

詳しいことは上に挙げた wiki 見て。

まず、

$$a^r\equiv 1\pmod N$$

となる \(r\) が存在する（proof: \(a^x\pmod N\) により構成される集合の大きさは有限なので \(a^\mu\equiv a^\nu\pmod N\) を満たす \(\mu,\nu\in\mathbb{N}\), \(\mu>\nu\) が存在する。 \(\gcd(a,N)=1\) より、両辺を \(a\) で割ることができるので、 \(\nu\) 回割ると \(a^{\mu-\nu}\equiv 1\pmod N\) が得られ、 \(r=\mu-\nu\) は存在する ■）。すると、 \(a^r-1=(a^{r/2}+1)(a^{r/2}-1)\) は \(N\) で割り切れるということになる（ここで割る 2 できる必要があるため \(r\) は奇数であってはならない。それと、割り切れるだけではまだ \(\gcd(a^{r/2}\pm 1, N)\neq 1\) とは言い切れないのだが、上に挙げた wiki ではこれを背理法により示している）。

そして、この \(r\) は量子パートで用いている \(r\) と同じなの、と疑問も出るが、これは

$$f(x+r)=f(x)$$ $$a^{x+r}\pmod N=a^x\pmod N$$

つまり、

$$a^{x+r}\equiv a^x\pmod N$$

両辺を \(a^x\) で割って

$$a^r\equiv 1\pmod N$$

より、示される。

計算量

QFT の計算量が \(O((\log N)^2)=O(n^2)\) であり、かつ \(m_1\), \(m_2\) が互いに素になる確率も定数である⁷ことから、クエリ計算量としては \(O(n^2)\) でできそうである。しかしながら、時間計算量としては \(f(x)\) は \(O(1)\) ではできず、後述する実装に用いる論文 p.5 では \(n^3\) 必要とあるので、諸々のゲート操作に必要な log 計算量を無視した \(\tilde{O}\) 表記を用いると \(\tilde{O}(n^3)\) ということになるのだろう。

ただし、 Wikipedia を見ると少し違う計算量に思える。。また、余談だが、 Shor 本人が計算量を答えているページがある。

例（非推奨）

\(N=8\) とする。答えは \(p=2\) だけれど、これを得られるか確認。

まずは古典でランダムに \(a=3\) が得られたとする。すると \(\gcd(3, 8)=1\) なので量子パートへ突入。

\(f(x)=3^x\pmod 8\) について \(f(x+r)=f(x)\) を満たす \(r\) を探す。はじめに

$$ \begin{align} &\frac{1}{\sqrt{8}}\sum_{x=0}^{8-1}\ket{x}\ket{0}\\ =&\frac{1}{\sqrt{8}}(\ket{0}+\ket{1}+\ket{2}+\ket{3}+\ket{4}+\ket{5}+\ket{6}+\ket{7})\ket{0} \end{align} $$

という状態を用意し、 \(f\) を計算。

$$\frac{1}{\sqrt{8}}(\ket{0}\ket{1}+\ket{1}\ket{3}+\ket{2}\ket{1}+\ket{3}\ket{3}+\ket{4}\ket{1}+\ket{5}\ket{3}+\ket{6}\ket{1}+\ket{7}\ket{3})$$

量子フーリエ変換。

$$ \begin{align} &\frac{1}{8}\sum_{x=0}^{8-1}\sum_{y=0}^{8-1}e^{2\pi ixy/8}\ket{y}\ket{f(x)}\\ =&\frac{1}{8}\{(\ket{0}+\ket{1}+\ket{2}+\ket{3}+\ket{4}+\ket{5}+\ket{6}+\ket{7})\ket{1}\\ &+(\ket{0}+e^{2\pi i/8}\ket{1}+e^{4\pi i/8}\ket{2}+e^{6\pi i/8}\ket{3}+e^{8\pi i/8}\ket{4}+e^{10\pi i/8}\ket{5}+e^{12\pi i/8}\ket{6}+e^{14\pi i/8}\ket{7})\ket{3}\\ &+(\ket{0}+e^{4\pi i/8}\ket{1}+e^{8\pi i/8}\ket{2}+e^{12\pi i/8}\ket{3}+e^{16\pi i/8}\ket{4}+e^{20\pi i/8}\ket{5}+e^{24\pi i/8}\ket{6}+e^{28\pi i/8}\ket{7})\ket{1}\\ &+(\ket{0}+e^{6\pi i/8}\ket{1}+e^{12\pi i/8}\ket{2}+e^{18\pi i/8}\ket{3}+e^{24\pi i/8}\ket{4}+e^{30\pi i/8}\ket{5}+e^{36\pi i/8}\ket{6}+e^{42\pi i/8}\ket{7})\ket{3}\\ &+(\ket{0}+e^{8\pi i/8}\ket{1}+e^{16\pi i/8}\ket{2}+e^{24\pi i/8}\ket{3}+e^{32\pi i/8}\ket{4}+e^{40\pi i/8}\ket{5}+e^{48\pi i/8}\ket{6}+e^{56\pi i/8}\ket{7})\ket{1}\\ &+(\ket{0}+e^{10\pi i/8}\ket{1}+e^{20\pi i/8}\ket{2}+e^{30\pi i/8}\ket{3}+e^{40\pi i/8}\ket{4}+e^{50\pi i/8}\ket{5}+e^{60\pi i/8}\ket{6}+e^{70\pi i/8}\ket{7})\ket{3}\\ &+(\ket{0}+e^{12\pi i/8}\ket{1}+e^{24\pi i/8}\ket{2}+e^{36\pi i/8}\ket{3}+e^{48\pi i/8}\ket{4}+e^{60\pi i/8}\ket{5}+e^{72\pi i/8}\ket{6}+e^{84\pi i/8}\ket{7})\ket{1}\\ &+(\ket{0}+e^{14\pi i/8}\ket{1}+e^{28\pi i/8}\ket{2}+e^{42\pi i/8}\ket{3}+e^{56\pi i/8}\ket{4}+e^{70\pi i/8}\ket{5}+e^{84\pi i/8}\ket{6}+e^{98\pi i/8}\ket{7})\ket{3}\}\\ =&\frac{1}{8}\{(4\ket{0}+(1+e^{\pi i/2}+e^{\pi i}+e^{3\pi i/2})\ket{1}+(1+e^{\pi i}+e^{2\pi i}+e^{3\pi i})\ket{2}+(1+e^{3\pi i/2}+e^{3\pi i}+e^{9\pi i/2})\ket{3}+(1+e^{2\pi i}+e^{4\pi i}+e^{6\pi i})\ket{4}+(1+e^{5\pi i/2}+e^{5\pi i}+e^{15\pi i/2})\ket{5}+(1+e^{3\pi i}+e^{6\pi i}+e^{9\pi i})\ket{6}+(1+e^{7\pi i/2}+e^{7\pi i}+e^{21\pi i/2})\ket{7})\ket{1}\\ &+(4\ket{0}+(e^{\pi i/4}+e^{3\pi i/4}+e^{5\pi i/4}+e^{7\pi i/4})\ket{1}+(e^{\pi i/2}+e^{3\pi i/2}+e^{5\pi i/2}+e^{7\pi i/2})\ket{2}+(e^{3\pi i/4}+e^{9\pi i/4}+e^{15\pi i/4}+e^{21\pi i/4})\ket{3}+(e^{\pi i}+e^{3\pi i}+e^{5\pi i}+e^{7\pi i})\ket{4}+(e^{5\pi i/4}+e^{15\pi i/4}+e^{25\pi i/4}+e^{35\pi i/4})\ket{5}+(e^{3\pi i/2}+e^{9\pi i/2}+e^{15\pi i/2}+e^{21\pi i/2})\ket{6}+(e^{7\pi i/4}+e^{21\pi i/4}+e^{35\pi i/4}+e^{49\pi i/4})\ket{7})\ket{3}\}\\ =&\frac{1}{8}\{(4\ket{0}+4\ket{4})\ket{1}\\ &+(4\ket{0}-4\ket{4})\ket{3}\} \end{align} $$

2 つのレジスタを測定し、 \(4\), \(1\) が得られたとすると、

$$\frac{y}{N}=\frac{4}{8}=\frac{1}{2}$$

より、 \(r=2\) を得る（本当は上述したように何回かやってちゃんと確かめる。いや、確かめずとも先へ進んでもそう問題ないのかも）。この \(r\) は奇数ではなく、かつ \(a^{r/2}=3^{2/2}=3\) なので \(a^{r/2}\equiv -1\pmod N\) も満たしていないので問題ない。したがって、

$$\gcd(a^{r/2}\pm 1, N)=\gcd(3\pm 1, 8)=4,2$$

より、確かに素因数分解できた⁸。

Qiskit による実装（非推奨）

ショアのアルゴリズムを Qiskit により実装したものを以下に示す。実装の全体は GitHub に置いているので、そこからたどってほしい（ただし、リンク先は後述する内容を反映したコードになっている）。

import random
import math
from typing import Optional

import numpy as np
import matplotlib.pyplot as plt
from qiskit import QuantumCircuit, Aer, transpile, assemble
from qiskit.visualization import plot_histogram

from qft import qft
from elementary import ax_modN


def shor(N: int, a: Optional[int] = None):
    if a is None:
        random.randint(2, N - 1)
    gcd = math.gcd(a, N)
    if gcd != 1:
        print(f'Answer is found in only classical computation: {gcd}')
        return

    N_len = int(np.ceil(np.log2(N)))
    qc = QuantumCircuit(10 * N_len - 2, 2 * N_len)

    qc.h(range(N_len))

    qc.append(ax_modN(a=a, N=N), range(10 * N_len - 2))

    qc.append(qft(n=N_len), range(N_len))

    qc.measure(range(2 * N_len), range(2 * N_len))

    backend = Aer.get_backend('aer_simulator_matrix_product_state')#('aer_simulator')
    qc = transpile(qc, backend)
    job = backend.run(qc, shots=1000)
    hist = job.result().get_counts()
    plot_histogram(hist)
    plt.show()


if __name__ == '__main__':
    shor(N=8, a=3)

重要なのは色を塗った部分で、まずアダマールにより重ね合わせ状態を作り、それからオラクル \(f(x)=a^x\pmod N\) を適用し、そして量子フーリエ変換をする。このうち最も難しいのは実はオラクル計算であり、ここではシンプルに ax_modN() を呼び出すだけにしているが、実際は elementary.py にあるようにとても複雑な量子回路となっている。これを実装するために論文 Quantum Networks for Elementary Arithmetic Operations およびその日本語解説記事を大いに参考にさせていただいた。

また、 Aer の backend として 'aer_simulator_matrix_product_state' を用いているが、これは今回のケースではこちらのほうが高速であったためである（参考）。

\(N=8\), \(a=3\) の場合の実行結果を示す（試行回数 1000 回により得られた分布）。

若干表記がわかりづらいが、例えば左端の 001000 は \(\ket{0}\ket{1}\) を表している。ゆえに、この結果は上述した例において手計算した結果と合致している。

（長かった。。。）

↓

↓

これで終わり、ではなかった。

もう 1 つの例（連分数展開による近似の活用）（非推奨）

\(f(x)=3^x\pmod 7\) の例を考えてみる。 \(7\) は素因数分解できないので、今回は単に周期発見問題として \(r=6\) を求めるだけ（\(f(x+6)=f(x)\)）。そして \(7\) より大きい最小の \(2\) のべき乗として \(N=2^3=8\) を用いる（ここに \(N\) という文字を使うので以降は \(\bmod\) の方に \(M\) の文字を使う：\(M=7\)）。そうしてヒストグラムを描くと

さっきより棒が多い（なお、第 2 レジスタは省略している）。ピークは \(\ket{y}=\ket{0},\ket{4}\) であるが、念のためすべての \(y\in[0,7]\) について \(\frac{y}{N}\) を計算すると

$$\frac{1}{8},\ \frac{1}{4},\ \frac{3}{8},\ \frac{1}{2},\ \frac{5}{8},\ \frac{3}{4},\ \frac{7}{8}$$

となり、さっきみたいに分母に \(r=6\) なんてない。というか、そもそも最大公約数で出てきてほしい値が \(\frac{N}{r}=\frac{4}{3}\) で整数じゃない。

つまり、今回の例のようにうまく求まらないものが存在する。ここで、もう少しヒストグラムを詳しく見てみる。実はピークの 2 本（\(y=0,4\)）より小さいものの \(1,3,5,7\) が \(2,6\) よりも大きいのには理由がある。上に出てきた \(y=\frac{Nm}{r}\) に \(m=0,1,2,\cdots\) を代入すると、

\(m\)	\(y=\frac{Nm}{r}\)
0	0
1	1.333
2	2.666
3	4
4	5.333
5	6.666

となる。つまり、 1.333 は 2 よりも 1 に近いので \(y=1\) の方が \(y=2\) よりも大きくなっている。 3 以降も同様である。したがって、実際には小数である \(y\) を分子に置いた \(\frac{y}{N}\) で \(N\) を大きくすれば \(\frac{m}{r}\) を近似して求められそうである⁹。

そこで連分数展開を用いる。

まず、 \(N\gg M\) として適当に \(N=2^6=64>49=M^2\gg M\) をとる。これでヒストグラムを描くと、

となる。このうち最も大きいピーク 2 つは \(y=0,32\) で 0 は論外として 32 を見ると \(\frac{32}{64}=\frac{1}{2}\) で連分数展開どころではない。ゆえにその次点で大きいピークのうち一番左（目盛りと少しずれててわかりづらい）の \(001011=11\) で連分数展開をおこなってみる。

$$\frac{11}{64}=0+\frac{1}{5+\frac{1}{1+\frac{1}{4+\frac{1}{2}}}}$$

したがって、この近似を順番に列挙していくと、

$$0+\frac{1}{5}=\frac{1}{5}$$
$$0+\frac{1}{5+\frac{1}{1}}=\frac{1}{6}$$

あっ、見つかった。したがって、 \(\frac{m}{r}=\frac{1}{6}\) より \(r=6\) が得られた。このあたりの話はこのスライド（p.12 から）がわかりやすい。

ここまでの内容を反映した Qiskit による実装コードは次のようになる。

import random
import math
from typing import Optional

import numpy as np
import matplotlib.pyplot as plt
from qiskit import QuantumCircuit, Aer, transpile, assemble
from qiskit.visualization import plot_histogram
from sympy import Rational
from sympy.ntheory.continued_fraction import continued_fraction, continued_fraction_convergents

from qft import qft
from elementary import ax_modM


def shor(M: int, a: Optional[int] = None, N_len: Optional[int] = None, show_hist: bool = True, use_only_period: bool = False):
    if a is None:
        random.randint(2, M - 1)
    gcd = math.gcd(a, M)
    if gcd != 1:
        print(f'Answer is found in only classical computation: {gcd}')
        return gcd

    if N_len is None:
        N_len = int(np.ceil(np.log2(M ** 2)))
    N = 2 ** N_len
    qc = QuantumCircuit(10 * N_len - 2, N_len)

    qc.h(range(N_len))

    qc.append(ax_modM(a=a, M=M, N_len=N_len), range(10 * N_len - 2))

    qc.append(qft(n=N_len), range(N_len))

    qc.measure(range(N_len), range(N_len))

    backend = Aer.get_backend('aer_simulator_matrix_product_state')#('aer_simulator')
    qc = transpile(qc, backend)
    job = backend.run(qc, shots=10000)
    hist = job.result().get_counts()

    if show_hist:
        plot_histogram(hist)
        plt.show()

    print(sorted(hist.items(), key=lambda x: x[1], reverse=True))
    y_list = []
    for measured_key, _ in sorted(hist.items(), key=lambda x: x[1], reverse=True):
        y = int(measured_key[-N_len:], 2)
        if y == 0:
            continue

        for fraction in continued_fraction_convergents(continued_fraction(Rational(y, N))):
            maybe_r = fraction.denominator
            if maybe_r != 1 and maybe_r < M and a ** maybe_r % M == 1:
                if use_only_period:
                    return maybe_r

                if maybe_r % 2 == 1 or (a ** (maybe_r // 2) + 1) % M == 0:
                    continue
                gcd = math.gcd(a ** (maybe_r // 2) + 1, M)
                if gcd != 1:
                    return gcd
                gcd = math.gcd(a ** (maybe_r // 2) - 1, M)
                if gcd != 1:
                    return gcd

    print('gcd was not found?!')

if __name__ == '__main__':
    print(shor(M=8, a=3))
    print(shor(M=7, a=3, use_only_period=True))

shor(M=8, a=3) とすると \(f(x)=3^x\pmod 8\) での周期 \(r\) を見つけそれを用いて 8 の素因数 4 を返す。一方、 shor(M=7, a=3, use_only_period=True) とすると \(f(x)=3^x\pmod 7\) での周期 \(r=6\) を見つけそれを返す。双方とも用いる \(N\) は \(N=64=8^2>7^2\) となる。

以上！！