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

July 17, 2021, 9:16 a.m. edited Nov. 27, 2021, 9:17 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}} $$

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

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

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

背景

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

問題

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

アルゴリズム

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

古典パート

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

量子パート

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

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

$$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)}$$

ここに量子フーリエ変換6 \(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\) が互いに素であれば7 \(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\) が互いに素になる確率も定数である7ことから、クエリ計算量としては \(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$$

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

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}\) を近似して求められそうである9

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

まず、 \(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\) となる。

以上!!


  1. Mac の Google 日本語入力、いくら量子アルゴリズムと入力してもすぐに忘れて漁師アルゴリズムになるのなんとかならんの?(手動登録なんてしたくないよ。がんばって) 

  2. \(\tilde{O}\) ってなんぞ、となったが、 log を無視するビッグオーのことらしい 

  3. なんか文献によって載ってる Shor のアルゴリズムの内容、違うんだよね。。 

  4. ユークリッドの互除法により \(O(\log N)\) で計算できる 

  5. \(\bmod\) ってなんだっけとなりがちだが、要は余りのことで、 \(a\equiv b\pmod c\) だと \(a-b\) は \(c\) で割り切れるということを表す(合同式)。ここからわかるように合同式では負の値も使える 

  6. \(U_{QFT}\ket{x}=\frac{1}{\sqrt{N}}\sum_{y=0}^{N-1}e^{2\pi ixy/N}\ket{y}\) 

  7. 互いに素にならない確率は \(1-\pi^2/6\) 未満(量子情報科学入門 p.76 に素数の自乗の逆数和を用いた証明が載っている) 

  8. なお、この \(a\) で問題ないかは特に確かめずに LaTeX を書き進めていたので、もし満たさなかったらどうしようとビクビクしていた() 

  9. この辺のアイデアのもう少し詳しいところ