量子コンピュータ（シミュレータ）で素因数分解する

はじめに

（注）私は量子情報の専門家でも，情報の専門家でもありません。
不適切な表現等ありましたらご指摘いただけると幸いです。

前回の記事では量子コンピュータ（シミュレータ）上でモジュロ加算器を作りました。

制御モジュロ乗算器・モジュロ累乗器も作成しようとしましたが，量子ビット数の上限が20のため，一旦諦めました。

そこで，コンパイル版の量子関数を使用し，Shorのアルゴリズム中，素因数分解を量子コンピュータ（シミュレータ）で実行してみました。

この記事では，コンパイル版とは関数表を用いて量子オラクルに量子ゲートを組み込むことを意味するものとします。周期が容易に推定できる問題に限り有効であるのかな，と個人的に感じています。

以下，備忘録になりますがその経緯をまとめます。できるだけ数式変形の詳細に触れず，オブラートに概要を残せればと考えています。

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

前準備

関数$f(x)=a^x ,,\text{mod},, N$を考えるために，$N$と互いに素な数$a(<N)$を選ばないといけません。今回は$N=57$，$a=37$としました。つまり，考えるべき数式は次のようになります。

y = f(x) = 37^x \,\,\,\text{mod}\,\, 57

関数表の作成

$N=57(<2^7)$より$y$を7量子ビット，$x$を$6$量子ビットとしました。

$x$の量子ビット数は関数表にどれくらい書き下すかによりますが，今回はかなり少ないので結果的には足りました。

関数表は次の通りです。ここで，$x_i$，$y_j$は$x$，$y$の2進数表記としています。

$x_5$	$x_4$	$x_3$	$x_2$	$x_1$	$x_0$	$y_6$	$y_5$	$y_4$	$y_3$	$y_2$	$y_1$	$y_0$
0	0	0	0	0	0	0	0	0	0	0	0	1
0	0	0	0	0	1	0	1	0	0	1	0	1
0	0	0	0	1	0	0	0	0	0	0	0	1
0	0	0	0	1	1	0	1	0	0	1	0	1
0	0	0	1	0	0	0	0	0	0	0	0	1
0	0	0	1	0	1	0	1	0	0	1	0	1
0	0	0	1	1	0	0	0	0	0	0	0	1
0	0	0	1	1	1	0	1	0	0	1	0	1
$\vdots$	$\vdots$	$\vdots$	$\vdots$	$\vdots$	$\vdots$	$\vdots$	$\vdots$	$\vdots$	$\vdots$	$\vdots$	$\vdots$	$\vdots$

上表より， $y$の値は$x_0$に依存しそれより上位のビットは影響を及ぼしていないことがわかりました。

したがって，上表を満たすように量子回路を組みます。それが結果的に量子オラクルにあたります。量子回路は次に示す通り，バリアに挟まれている部分です。
（この時点で周期はわかったようなものですが，量子アルゴリズムでも計算できることを検証するため，そのまま進めます）

画像の左端では，$x$の各量子ビットにアダマール変換ゲートを通して重ね合わせ状態にしています。

これにより，先ほど作った量子回路に「全ての場合の$x$」の入力で量子計算し，yq[0]～yq[6]を測定することで結果を確認することができます。なお，最上位ビットがそれぞれxq[0]，yq[0]として計算しています。

計算結果は次のようになりました。$\ket{0100101_2}=\ket{37}$と$\ket{0000001_2}=\ket{1}$の重ね合わせ状態であることがわかりました。

Shorのアルゴリズムにしたがうと，この後，$\ket{\text{xq[0]xq[1]xq[2]xq[3]xq[4]xq[5]}}$に逆量子フーリエ変換を適用することにより干渉が起き，欲しい解の確率振幅を増やすことになります。

逆量子フーリエ変換ゲート・SWAPゲートの準備

量子フーリエ変換・逆量子フーリエ変換では，入力に対し出力の量子ビットの順序が逆になるという特徴があります。ですので，量子フーリエ変換の直後にSWAPゲートを通して，量子ビットの順番を変更しておくと混乱せずに計算を進められます。

SWAPゲートはOpenQASMの量子アセンブリでは次のようになります。

gate swap swap1,swap2 {
cx swap1,swap2;
cx swap2,swap1;
cx swap1,swap2;
}

続いて，逆量子フーリエ変換ゲートは下記のようになります。
量子フーリエ変換は，以前の記事で少し触りましたが，逆量子フーリエ変換の定義式は下記になります。ここで，$\omega=e^{2\pi i/N'}$，$N'=2^n(n$は量子ビット数)です。時間短縮のため前回記事のゲートの配置の規則性を利用して，6量子ビットの場合の逆量子フーリエ変換のサブルーチンを作成しました。

\ket{j} \xrightarrow{QFT^{-1}} \frac{1}{\sqrt{N'}}\sum_{k=0}^{N'-1}\omega^{-jk}\ket{k}

gate invqft32 q0,q1,q2,q3,q4,q5 {
h q0;
cu1 (-pi/2) q0,q1;
cu1 (-pi/4) q0,q2;
cu1 (-pi/8) q0,q3;
cu1 (-pi/16) q0,q4;
cu1 (-pi/32) q0,q5;
h q1;
cu1 (-pi/2) q1,q2;
cu1 (-pi/4) q1,q3;
cu1 (-pi/8) q1,q4;
cu1 (-pi/16) q1,q5;
h q2;
cu1 (-pi/2) q2,q3;
cu1 (-pi/4) q2,q4;
cu1 (-pi/8) q2,q5;
h q3;
cu1 (-pi/2) q3,q4;
cu1 (-pi/4) q3,q5;
h q4;
cu1 (-pi/2) q4,q5;
h q5;
}

位数発見のゲート

上述した内容を踏まえて，周期発見のゲートのOpenQASM・量子回路は次のようになります。

include "qelib1.inc";
qreg xq[6];
qreg yq[7];
creg c[6];
gate swap swap1,swap2 {
cx swap1,swap2;
cx swap2,swap1;
cx swap1,swap2;
}
gate invqft32 q0,q1,q2,q3,q4,q5 {
h q0;
cu1 (-pi/2) q0,q1;
cu1 (-pi/4) q0,q2;
cu1 (-pi/8) q0,q3;
cu1 (-pi/16) q0,q4;
cu1 (-pi/32) q0,q5;
h q1;
cu1 (-pi/2) q1,q2;
cu1 (-pi/4) q1,q3;
cu1 (-pi/8) q1,q4;
cu1 (-pi/16) q1,q5;
h q2;
cu1 (-pi/2) q2,q3;
cu1 (-pi/4) q2,q4;
cu1 (-pi/8) q2,q5;
h q3;
cu1 (-pi/2) q3,q4;
cu1 (-pi/4) q3,q5;
h q4;
cu1 (-pi/2) q4,q5;
h q5;
}

h xq[0];
h xq[1];
h xq[2];
h xq[3];
h xq[4];
h xq[5];

barrier xq[0],xq[1],xq[2],xq[3],xq[4],xq[5],yq[0],yq[1],yq[2],yq[3],yq[4],yq[5],yq[6];
x yq[6];
cx xq[5],yq[1];
cx xq[5],yq[4];
barrier xq[0],xq[1],xq[2],xq[3],xq[4],xq[5],yq[0],yq[1],yq[2],yq[3],yq[4],yq[5],yq[6];

invqft32 xq[0],xq[1],xq[2],xq[3],xq[4],xq[5];
swap xq[0],xq[5];
swap xq[1],xq[4];
swap xq[2],xq[3];

// get results
measure xq[0] -> c[5];
measure xq[1] -> c[4];
measure xq[2] -> c[3];
measure xq[3] -> c[2];
measure xq[4] -> c[1];
measure xq[5] -> c[0];

計算結果

計算結果は下記のとおりです（Shotsは100）。

約50%の確率で$\ket{100000_2}=\ket{32}$が得られました。

もう一つ，Shotsが8192のときの計算結果は下記のとおりです。

いずれの場合も，$\ket{100000_2}=\ket{32}$の得られる確率が最大で約50%でした。

位数発見

先ほどの計算結果から，観測値$s$として，$100000_2=32$，$000000_2=0$が得られました。これは，整数$k=(0,1,2,...,63)$に対して，位数$r$は$\phi=k/r$の6量子ビットの推定として，$2^6 (k/r)$を得たことになります。これにより，位数を計算します。

s=2^6\left( \frac{k}{r} \right) \to \frac{s}{64}=\frac{k}{r}

• s = 32のとき：
  $\frac{s}{64}=\frac{32}{64}=\frac{1}{2}$より$r=2$
  
  このとき，$r$が偶数なので最大公約数を計算できます。

>>> from math import gcd
>>> print(gcd(37**int(2/2)-1,57),gcd(37**int(2/2)+1,57))
3 19

• s = 0のとき：位数は求まりません

以上より，$s=2$のとき，57(=$3\times 19$)を素因数分解できました。

おわりに

今回は，以下のことを行いました。

• 量子コンピュータ（シミュレータ）を用いて57を素因数分解

ここまで読んでくださってありがとうございました。

参考文献

• 量子アルゴリズム
• クラウド量子計算入門
• クラウド量子計算　量子アセンブラ入門

来歴

2018.01.24　新規作成。
2018.01.25　逆量子フーリエ変換のアセンブリ訂正。訂正に伴い文章，結果修正。