量子コンピュータを使用する素因数分解アルゴリズム

量子コンピュータを使うと素因数分解が高速にできる。

が、実際はどうやって素因数分解をするんでしょう?
そして具体的に何が高速化するんでしょう?これを知るのが本稿の目的です。

  1. 素因数分解アルゴリズム(古典的)
  2. 実際のアルゴリズム(量子的)
    1. 初期状態の準備
    2. レジスタ1の状態に量子フーリエ変換を施す
    3. レジスタ2に、レジスタ1を元にした関数を作用させる
    4. レジスタ2を観測する
    5. レジスタ1を量子フーリエ変換する
    6. レジスタ1の観測
    7. 連分数展開によって周期\(r\)を見つける
  3. 最大公約数を求めるアルゴリズム(古典)

素因数分解アルゴリズム


量子的に高速に出来るのであれば、そのアイデア自体は古典的に考えられるはずです。

素因数分解する古典的アルゴリズムは以下のように記述されます。
\(m\)を素因数分解したい数とし、\(x\)をmと互いに素な整数(\(\text{gcd}(x,m)=1\)を満たす数)と表記します。

ここで、\(\text{gcd}(a,b)\)は、\(a\)と\(b\)の最大公約数(greatest common divisor)を出力する関数です。

  1. \(a(=0,1,\cdots )\)を引数として、関数\(f_m(a)=x^a~ \text{mod}~ m\)を計算
  2. \(f_m(a)\)の周期\(r\)を見つける。実は\(f_m(a)\)は周期\(r\)の周期関数(周期\(r\)は実際に計算してみるまで分からない。数式で表せば、\(f_n(a)=f_n(a+r)\))。
  3. 周期\(r\)が判明したら、\(p=\text{gcd}(x^{r/2}+1,m),~q=\text{gcd}(x^{r/2}-1,m)\)を計算
  4. \(p,q\)は\(m\)の素因数となっている

という風に素因数分解を行うことが出来ます。

具体的に\(m=21\)を素因数分解してみましょう。\(m\)と互いに素な整数として、\(x=11\)を選びます。
実際に横軸に\(a\)をとり、関数\(f_m(a)\)を計算すると以下のようになります。

この計算結果を見ますと、周期\(r\)が6だと分かります。

周期が\(6\)が分かったので、\(p=\text{gcd}(x^{r/2}+1,m),~q=\text{gcd}(x^{r/2}-1,m)\)を計算します。すると、
\(
\begin{align}
p&=\text{gcd}(x^{r/2}+1,m)=\text{gcd}(11^{6/2}+1,21)=3 \\
q&=\text{gcd}(x^{r/2}-1,m)=\text{gcd}(11^{6/2}-1,21)=7
\end{align}
\)

と分かります。\(21\)は\(3\times 7\)と書けるので、素因数分解を行うことが出来ました。

さて、量子的に素因数分解を行う場合、上のアルゴリズムはShorのアルゴリズムと呼ばれます。
上の計算は、

  • 最大公約数を求めるアルゴリズム(Euclidの互除法, 計算量\(O(\text{log} n)\))
  • 周期\(r\)を見つけるアルゴリズム(
    古典的:\(O(\exp[(\text{log}n)^{1/3}(\text{ln}\text{ln}n)^{2/3}])\),
    量子的:\(O((\text{log}n)^{2}(\text{log}\text{log}n)(\text{log}\text{log}\text{log}n))\)

です。最も時間が掛かる部分は、周期\(r\)を見つける部分で、これを見つけるために量子コンピュータを用いるのです。Shorのアルゴリズムは、この部分だけ量子コンピュータを使うのです。

量子的に高速に計算する、という部分は、
\(f_m(a)\)の周期\(r\)を見つける
という部分です。

実際のアルゴリズム


ここで紹介するモデルは、量子ビット(キュービット)を出さないで説明するものです。

キュービットの考えはShorのアルゴリズムの本質ではなく、実現できるか?の部分で重要です。
すなわち、ある状態を指定するために
・1つの系で、非常に高い量子番号に属する状態まで制御する(\(1[つの系]\times N[量子状態]\))
ではなく、
・複数の系で、基底状態と1つの励起状態を制御する(\(Q[つの系]\times 2[量子状態]\))
の方が現実で実現しやすいという考えです。

1. 初期状態の準備



2つの系(2つの記憶素子、レジスタ)を用意し、初期状態として各々の基底状態を表します。
この時点で、2つの系は独立に存在するもので、系同士に相互作用や、もつれ等は存在しません。

それぞれの系を添え字\(1,2\)で表すと、系全体の状態\(|\psi\rangle\)は(直積の記号を省略して)
\(
|\psi\rangle=|0\rangle_1|0\rangle_2
\)

と書き表せます。
ここで、\(|a\rangle\)は量子数\(a\)に属する固有状態を表しています。

図示すると、上の図のようになり、左はレジスタ1, 右はレジスタ2の量子状態を表していて、それぞれ、基底状態であることを表しています。

2. レジスタ1の状態に量子フーリエ変換を施す


レジスタ1に量子フーリエ変換を施すことで、レジスタ1の全状態の確率振幅を等確率にします。
量子フーリエ変換は量子状態に対して施す変換で、

\(
\begin{align}
\mathcal{F}|a\rangle&=\frac{1}{\sqrt{N}}\sum_{b=0}^{N-1}e^{-2\pi i a b/N}|b\rangle \\
\mathcal{F}^{-1}|a\rangle&=\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}e^{2\pi i a b/N}|a\rangle
\end{align}
\)

で定義します。基底状態\(|0\rangle\)に対して量子フーリエ変換を施すと、考えたい全励起状態を等確率にすることが出来ます。
\(
\begin{align}
\mathcal{F}|0\rangle&=\frac{1}{\sqrt{N}}\sum_{b=0}^{N-1}e^{-2\pi i 0 b/N}|b\rangle \\
&=\frac{1}{\sqrt{N}}\sum_{b=0}^{N-1}|b\rangle
\end{align}
\)

量子フーリエ変換を作用させると、上図のようになります。

3. レジスタ2に、レジスタ1を元にした関数を作用させる


レジスタ2に、レジスタ1を参照して計算した結果をレジスタ2に出力する、という演算を行います。
この操作を行うことで、レジスタ1と2をもつれさせます(もつれ、という言葉が正しいかは分かりません。ですが、言いたいことはレジスタ1と2を結びつける、という操作を行います)。

この演算を行う事が出来る装置があるかはわかりませんが、あると仮定します。この ”あるかどうか分からない演算を行う装置” はオラクル(神託装置)といいます。

このオラクルを用いて、レジスタ1の量子数\(a\)の値を引数として、レジスタ2の量子状態を\(|f_m(a)\rangle\)にしていきます。

この作用を\(\hat{U}_{1\to2}\)という演算子で表すと、レジスタ1の結果をレジスタ2に格納するので、

\(
\begin{align}
\hat{U}_{1\to2} \hat{\mathcal{F}}_1|0\rangle_1|0\rangle_2
&=\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}\hat{U}_{1\to2}|a\rangle_1|0\rangle_2 \\
&=\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}|a\rangle_1|x^a~ \text{mod}~ m\rangle_2
\end{align}
\)

4. レジスタ2を観測する


レジスタ2を観測し、レジスタ2の量子数を知ります。
この観測は以下の2つの影響を系に与えます。

・レジスタ2の量子状態を確定させる
・レジスタ1が取り得る量子状態に制限を加える

1回の”観測”の操作はただ1点だけを与えますが、数式上扱いづらいです。
なので、数式では何度も観測を行った時に得られる期待値を数式で表すと分かりやすくなります。
レジスタ2の状態を量子数\(k\)の状態に見出す確率を知るために、左から\(\langle k|_2\)を作用させます。
すると、以下のように式変形をすることが出来ます。

\(
\begin{align}
\langle k|\hat{U}_{1\to 2} \hat{\mathcal{F}}_1|0\rangle_1|0\rangle_2
&=\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}|a\rangle_1|x^a ~\text{mod}~ m\rangle_2 \\
&\approx\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}|a\rangle_1\langle k|x^a ~\text{mod}~ m\rangle_2 \\
&=\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}|a\rangle_1\left(\frac{1}{\sqrt{J}}\sum_{j=0}^J \delta_{a,a_j^{(k)}}\right)
\end{align}
\)
ここで、\(a_j^{(k)}\)は、
\(k=x^{a_{j}^{(k)}}~ \text{mod} ~m,~~(j=0,1,\cdots,J-1)\)を満たす\(J\)個の\(a\)です。

この観測により、レジスタ1の状態を特定の状態を\(|a_j^{(k)}\rangle,~~(j=0,1,\cdots, J-1)\)のみ存在させることになります。レジスタ1のとれる状態は周期的で周期\(r\)をもち、
\(
\displaystyle a_j^{(k)}=a_0^{(k)}+j\cdot r,~~(j=0,1,\cdots, J-1)
\)

という形で表すことが出来ます。

レジスタ1の状態を図示するとこのようになっているはずです。

係数をまとめると分かりやすくなります。
単なる式変形ですが、状態に適当な係数が掛かった状態の重ね合わせで表現されるので、なじみ深い形…だと私は思います。
\(
\begin{align}
\langle k|\hat{U}_{1\to 2} \hat{\mathcal{F}}_1|0\rangle_1|0\rangle_2
\approx=\sum_{a=0}^{N-1}\left(\frac{1}{\sqrt{NJ}}\frac{1}{\sqrt{J}}\sum_{j=0}^{J-1} \delta_{a,a_j^{(k)}}\right)|a\rangle_1
\end{align}
\)

5. レジスタ1を量子フーリエ変換する



レジスタ1の状態に対して量子フーリエ変換を行うことで、レジスタ1の持つ周期を知ることが出来ます。

量子フーリエ変換を行う前は量子の属する状態そのものに周期\(r\)の情報は含まれていませんが、量子フーリエ変換を行った後は量子数そのものに周期\(r\)に関係する値が含まれるようになります。

実際に観測を行うと、

\(
\begin{align}
\hat{\mathcal{F}}^{-1}_1\langle k|\hat{U}_{1\to 2} \hat{\mathcal{F}}_1|0\rangle_1|0\rangle_2
&\approx\sum_{a=0}^{N-1}\left(\frac{1}{\sqrt{NJ}}\frac{1}{\sqrt{J}}\sum_{j=0}^J \delta_{a,a_j^{(k)}}\right)\hat{\mathcal{F}}^{-1}_1|a\rangle_1 \\
&=\sum_{b=0}^{N-1}\left(\frac{1}{N\sqrt{J}}\sum_{j=0}^{J-1} e^{2\pi i a_j^{(k)}b/N}\right)|b\rangle_1 \\
&=\sum_{b=0}^{N-1}\left[\frac{1}{N\sqrt{J}}e^{2\pi i a_0^{(k)}b/N}\sum_{j=0}^{J-1} e^{2\pi i jrb/N}\right]|b\rangle_1
\end{align}
\)

レジスタ1の量子数\(b’\)に見出す確率は、左から\(\langle b’|_1\)を作用させて絶対値の2乗をとると分かります。

まず、左から\(\langle b’|_1\)を作用させると、
\(
\begin{align}
C'(b’)=\langle b’|_1\hat{\mathcal{F}}^{-1}_1\langle k|\hat{U}_{1\to 2} \hat{\mathcal{F}}_1|0\rangle_1|0\rangle_2
&\approx\frac{1}{N\sqrt{J}}e^{2\pi i a_0^{(k)}b’/N}\sum_{j=0}^{J-1} e^{2\pi i jrb’/N}
\end{align}
\)

となります。確率は上記\(C'(b)\)の値の絶対値の2乗で与えられるので(表記を\(b’\to b\)に変更)、
\(
\begin{align}
|C'(b)|^2&=\left|\frac{1}{N\sqrt{J}}e^{2\pi i a_0^{(k)}b/N}\sum_{j=0}^{J-1} e^{2\pi i jrb/N}\right|^2 \\
&=\frac{1}{N^2 J}\left|\sum_{j=0}^{J-1} e^{2\pi i jrb/N}\right|^2
\end{align}
\)
となります。確率密度は
\(
\displaystyle \left|\sum_{j=0}^{J-1} e^{i(2\pi jrb/N)}\right|
\)
に依存します。この形は離散フーリエ変換の場合に似ています。
つまり、この項が大きくなる時の状態がレジスタ1の量子状態として観測されることが期待され、
観測される状態とは、\(j\)が変化しても位相が一致しているとき、すなわち

\(2\pi rb/N=2\pi \times (整数)\)

を満たすような\(r\)の時に存在確率密度が大きくなる、ということになります。
整数をsと表記すると、

\(rb/N= s~\to~ b/N= s/r\)

と表されます。

6. レジスタ1の観測



確率密度分布は\(b=a_0^{(k)}+j\cdot r\)の量子状態を中心に分布しますが、あくまで、確率分布ですので、観測を行った場合は若干ずれた所に見出したりします。

実際にレジスタ1の観測を行った際に得られた量子数を\(A\)とします。すると、関係式

\(\frac{A}{N}\approx \frac{s}{r}\)

が近似的に成立していると考えることが出来ます。
今、左辺は完全に分かっている量ですが、右辺の分子・分母は整数を持つ、ということくらいしかわかりません。これを求めるために連分数展開を用います。

7. 連分数展開によって周期\(r\)を見つける


レジスタ1の観測によって得られた状態の量子数を\(A\)とします。今、
\(\displaystyle
\frac{A}{N}=\frac{s}{r}
\)

ここで、\(A\)は観測された量子数、
\(N\)はレジスタの全量子状態数、
\(s (\lt N)\)は任意の整数、
\(r\)は求めたい周期です。

これから行いたい操作は、\(\frac{A}{N}\)の分子分母を出来るだけ小さい数にする、という作業を行いたいのです。
これは単純に割るだけで済む話ではありません。なぜなら、\(\frac{A}{N}\)は綺麗に割り切れる数ではない可能性があるからです。

例えば、\(k=184, N=243\)という値であった場合、直感的に\(\frac{184}{243}\approx \frac{3}{4}\)だということが分かります。

これを機械的に行うためには連分数展開を利用します。
連分数展開とは、ある実数\(x\)を整数\(a_n\)を用いて以下のように展開することです。
\(\displaystyle
x=a_0+\frac{1}{a_1+\frac{1}{a_2+\frac{1}{\cdots}}}
\)

ここで、\(a_n~(n=0,1\cdots)\)は、
\(
\begin{align}
&a_0=\lfloor x \rfloor,~~b_0=\frac{1}{x-a_0} \\
&a_n=\lfloor b_{n-1} \rfloor,~~b_n=\frac{1}{b_{n-1}-a_n}
\end{align}
\)
と得られます。ここで\(\lfloor x \rfloor\)はガウスの記号で、実数の整数部分を表しています。
\(x\)が有理数であれば、連分数展開は有限の項数で終わります。

連分数展開後、\(x\)の分数による近似は
\(\displaystyle
x\approx\frac{d_n}{r_n}
\)

ここで、
\(
\begin{align}
& d_0=a_0,~d_1=1+a_0a_1,~d_n=a_n d_{n-1}+d_{n-2},\\
& r_0=1,~r_1=a_1,~r_n=a_n r_{n-1}+r_{n-2},
\end{align}
\)

として順々に求める事が出来ます。
また、連分数の打切りに関しては定理

「\(\frac{q}{p}\)を\(\left|\frac{q}{p}-x\right|\lt \frac{1}{2q^2}\)を満たす任意の有理数とする時、\(\frac{q}{p}\)は\(x\)の近似分数である。さらに、その近似分数は\(p,q\)の最大公約数は1である。」

を用いて、条件式\(\left|\frac{q}{p}-x\right|\lt \frac{1}{2q^2}\)から決めます。

連分数展開によって分子・分母を出力するプログラムはこんな感じで実装できます。

最大公約数を求めるアルゴリズム(古典)


最大公約数を求めるアルゴリズムはEuclidのアルゴリズムと呼ばれ、他のアルゴリズムと比べても早く終わります。そのプログラムは以下で実装できます。

function gcd(a,b)
  implicit none
  integer,intent(in)::a,b
  integer::gcd
  !compute the greatest common divisor
  integer::x,r
 
  x=a
  gcd=b
  r=mod(x,gcd)
  do while( r .gt. 0)
     x=gcd
     gcd=r
     r=mod(x,gcd)
  enddo
 
  return
end function gcd

参考文献

Elisa Bäumer, Jan-Grimo Sobez, Stefan Tessarini, Shor’s Algorithm
https://qudev.phys.ethz.ch/content/QSIT15/Shors%20Algorithm.pdf

C. P. ウィリアムズ、S. H. クリアウォータ著、西野哲郎、新井隆、渡邉昇訳「量子コンピューティング」(springer, 2000)

exp(-ikx)/x の無限区間に渡る積分を数値的に計算したい。

問題:以下の積分を数値的に計算せよ。
\( \displaystyle
-\frac{1}{2\pi i}\int_{-\infty}^{\infty} \frac{e^{-ikx}}{x}dx,~~(kは任意の実数)
\)


1/xのフーリエ変換である
\( \displaystyle
\int_{-\infty}^{\infty} \frac{e^{-ikx}}{x}dx,~~(kは任意の実数)
\)

数値的に計算します。

この積分値は一意には決まりません。
複素関数論によると、この積分値は\(x=0\)周りに存在する特異点を上周りに迂回して積分するか、下周りに迂回して積分するか、で答えが異なります。

特異点を上周りに回る場合、答えは
\(
\begin{eqnarray}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{e^{-ikx}}{x+i\varepsilon}dx =
\left\{
\begin{aligned}
-2\pi i~~~(k\gt 0)\\
0~~~(k\lt 0)
\end{aligned}
\right.
\end{eqnarray}
\)
であり、特異点を下周りに回る場合、
\(
\begin{eqnarray}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{e^{-ikx}}{x-i\varepsilon}dx =
\left\{
\begin{aligned}
0~~~(k\gt 0)\\
2\pi i~~~(k \lt 0)
\end{aligned}
\right.
\end{eqnarray}
\)
となります。
私が疑問に思うのは、積分すると本当にステップ状の関数が現れるのか?ということです。にわかには信じられません。

ここでは、上を回った時の計算にのみ注目します。
答えに出てくる係数を減らしたいため、被積分関数
\( \displaystyle
-\frac{1}{2\pi i}\int_{-\infty}^{\infty} \frac{e^{-ikx}}{x}dx
\)

を考えます。

数値計算方法


さて、問題を整理しましょう。今取り組む問題は、

問題:以下の積分を数値的に計算せよ。
\( \displaystyle
-\frac{1}{2\pi i}\int_{-\infty}^{\infty} \frac{e^{-ikx}}{x}dx,~~(kは任意の実数)
\)

です。読み取れる情報は

・始点は、実軸上のマイナス無限大
・終点は、実軸上のプラス無限大
・\(x=0\)に特異点が存在する

という情報です。これらは理由なく勝手に変えてはいけません。
被積分関数が収束するならば、xがどんなに大きくても無限遠まで振動していきますが、フレネル積分のように積分結果は収束することを期待しましょう。

さて、\(x=0\)に特異点があるために、実軸上のマイナス無限からプラス無限まで積分していこうとすると\(x=0\)を通過しようとする時に無限大になるため、オーバーフローします。
ということで、積分経路を特異点周りで迂回します。

始点と終点は決められていますが、積分経路の指定はありませんので何も問題ありません。
積分経路は以下の図のように取ります。

このようにとると、積分経路上に特異点はないですので計算することが出来ます。
始点と終点は十分大きい値、ここでは(±2000,0)という値にします。

横軸をkとして、積分を計算すると以下のようになります。

複素関数論から、特異点を複素平面上方を迂回した場合に得られる解は
\(
\begin{eqnarray}
\lim_{\varepsilon\to +0} -\frac{1}{2\pi i}\int_{-\infty}^{\infty}\frac{e^{-ikx}}{x+i\varepsilon}dx =
\left\{
\begin{aligned}
1~~~(k\gt 0)\\
0~~~(k\lt 0)
\end{aligned}
\right.
\end{eqnarray}
\)
です。数値計算によって得られた解は、\(k\lt 0\)で値0, \(k\gt 0\)で値1であるので、数値的にも正しそうです。
しかし、厳密にそれぞれの値になっているわけではないことに気が付きます。なぜならば、\(k=0\)周りでは値は収束していないように見えます。
この原因は始点と終点としてとった、大きい値(±2000,0)が十分ではない事を意味しています。kが小さいと、被積分関数の\(\exp(-ikx)\)による振動はゆっくりになるため、収束が遅くなっている、と考えられます。

ジョルダンの補題の利用


この振動を無くすためには被積分関数が漸近領域で収束するように積分経路を変えればいいのです。

始点と終点は変わりません。図は\(k\)が負の時の経路です。

経路の取り方は無数にあり、上の経路は一つ前の計算で選んだ経路を出来る限り変更しないように経路をとっているだけにすぎません。

ジョルダンの補題を利用すれば、被積分関数が収束する領域であれば積分の寄与はゼロなので、数値的に積分する場合の視点、終点を複素平面上へ移すことが出来ます。
この考えの下実際に計算すると積分結果はこうなります。

綺麗なエッジが作られました。これならば解析解と数値解が一致した、と言えるでしょう。複素平面上の数値計算は一筋縄ではいかないことが実感できました。

数値計算プログラム


プログラムは以下のものを用いています。

数値積分はquadpackのガウス=ルジャンドル求積法の自動積分コードを複素関数用に変更して用いています。
実際に使用しているquadpackコードは以下のもので、
https://slpr.sakura.ne.jp/qp/supplement_data/quadpack_cqag.f90
これと、下のはmainプログラムです。

program main
  implicit none
  integer::ier,Nt,i,j,sn
  double precision::eps
  complex(kind(0d0))::a,b,s,s0
  complex(kind(0d0))::xa(1:5),xb(1:5)
  complex(kind(0d0)),external::g
     
  eps=1d-10
  ier=1

  Nt=1000
  do i=1,Nt
     t=i*(1d0-(-1d0))/dble(Nt)-1d0
     s=dcmplx(0d0,0d0)
     
     sn=1
     if(t.gt.0d0)sn=-1
     if(t.eq.0)sn=0
     
     xa(1)=dcmplx(-2000d0,sn*2000d0)
     xb(1)=dcmplx(-1d0,0d0)

     xa(2)=dcmplx(-1d0,0d0)
     xb(2)=dcmplx(-1d0,1d0)

     xa(3)=dcmplx(-1d0,1d0)
     xb(3)=dcmplx( 1d0,1d0)

     xa(4)=dcmplx( 1d0,1d0)
     xb(4)=dcmplx( 1d0,0d0)

     xa(5)=dcmplx( 1d0,0d0)
     xb(5)=dcmplx( 2000d0,sn*2000d0)  

     do j=1,5
        s0=dcmplx(0d0,0d0)
        call cqag_sk(g,xa(j),xb(j),eps,s0,3,ier)
        s=s+s0
     enddo
     write(10,*)t,dble(s),dimag(s)
  enddo  
     
  stop
end program main

function g(z)
  implicit none
  complex(kind(0d0))::g
  complex(kind(0d0)),intent(in)::z

  double precision::pi=dacos(-1d0)
 
  g=exp(-dcmplx(0d0,1d0)*z*t)/(z)
  g=g*(-1d0/(2d0*pi*dcmplx(0d0,1d0)))

  return
end function g

オリジナルのquadpackコードは
http://www.netlib.org/quadpack/(パブリックドメイン)です。