Edited at

量子コンピュータ(シミュレータ)で素因数分解する【QISKit編】

More than 1 year has passed since last update.


はじめに

(注)私は量子情報の専門家でも,情報の専門家でもありません。

不適切な表現等ありましたらご指摘いただけると幸いです。

前回の記事では,量子コンピュータ(シミュレータ)で素因数分解しました。

(※2018/01/24付の前回記事中,逆量子フーリエ変換の実装で,不適切な箇所があり訂正しました。申し訳ありません。)

ここまでいくつかの計算例を記事にしてきましたが,Composerをいじるのに少々飽きてきました。

そこで今回は,QISKitおよびPython3.6を用いて,改めて実装してみようと思います。

QiitaでもQISKit関連の記事は数えるほどしかございませんので,今後増えていけばいいなあと個人的には感じています。

$$

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

$$


QISKit


概要

QISKitは量子情報ソフトウェア開発キットです。なんと公式ページに日本語版があります。

QISKitでの実装は@hiroyuki827さんのこの記事や,@shiibassさんのこの記事等があります。

他にもQISKitはこの記事この記事等でも少し言及されています。

IBM Qと同様にIBMが提供しているということで,今回はこれを採用しました。


インストール他

pipを使用してインストールします。

$ pip install qiskit

また,個人的にPython-fireを用いて開発中なこともあり,それもインストールします。

Python-fireに概要ついてはこちらの記事が端的かつ分かりやすかったです。

$ pip install fire


QISKitでの実装


使用するレジスタ

前回の記事と同様に,q[6],a[7],c[7]の3つのレジスタを考えます。

qp = QuantumProgram()

#--- localで使用するのでAPIの設定割愛
backend = 'local_qasm_simulator'

q = qp.create_quantum_register("q", 6)
a = qp.create_quantum_register("a", 7)
c = qp.create_classical_register("c", 7)
qcirc = qp.create_circuit("qcirc", [q, a], [c])

また,qには各量子ビットにアダマール変換ゲートを適用しておきます。

def input_state(circ, q, n):

for i in range(n):
circ.h(q[i])


逆量子フーリエ変換

$n$量子ビットの入力$\ket{\mathbf{Q}}=\ket{q_{n-1} q_{n-2}\cdots q_1 q_0}$($q_i$は2進数表記)に対する逆量子フーリエ変換は次のように実装されます。なお,各変換後は量子ビットの並び順が反転しますので,SWAPゲートを利用して,変換前の順序に合わせておきます。

def rev_qreg(circ, q, n):

for i in range(int(n/2)):
circ.swap(q[i], q[(n-1)-i])

def iqft(circ, q, n):
# 上位ビットから処理
for j in range(n-1,-1,-1):
circ.h(q[j])
for k in reversed(range(j)):
circ.cu1((-1)*math.pi/float(2**(j-k)), q[j], q[k])
rev_qreg(circ, q, n)

ここで,引数(circ,q,n)は(量子回路のオブジェクト,量子レジスタのオブジェクト,量子ビット数)です。


量子オラクル

前回の記事を参考に,コンパイル版の回路を組みます。

#--- 量子オラクルの設定

qcirc.cx(q[0], a[2])
qcirc.cx(q[0], a[5])
qcirc.x(a[0])


測定とヒストグラム表示

あとはこれらを組み合わせてZ基底で測定し,データをまとめるのみです。

次のようになります。

import fire

import math # for pi
from qiskit import QuantumCircuit, QuantumProgram
# Import basic plotting tools
from qiskit.tools.visualization import plot_histogram

def input_state(circ, q, n):
for i in range(n):
circ.h(q[i])

def rev_qreg(circ, q, n):
for i in range(int(n/2)):
circ.swap(q[i], q[(n-1)-i])

def iqft(circ, q, n):
# 上位ビットから処理
for j in range(n-1,-1,-1):
circ.h(q[j])
for k in reversed(range(j)):
circ.cu1((-1)*math.pi/float(2**(j-k)), q[j], q[k])
rev_qreg(circ, q, n)

def shor_57():
#--- 量子プログラムの開始
qp = QuantumProgram()
# localで使用するのでAPIの設定割愛
backend = 'local_qasm_simulator'

q = qp.create_quantum_register("q", 6)
a = qp.create_quantum_register("a", 7)
c = qp.create_classical_register("c", 7)
qcirc = qp.create_circuit("qcirc", [q, a], [c])

#--- 初期化
input_state(qcirc, q, 6)
#--- 量子オラクルの設定
qcirc.cx(q[0], a[2])
qcirc.cx(q[0], a[5])
qcirc.x(a[0])
#--- 逆量子フーリエ変換
iqft(qcirc, q, 6)

#--- 測定
for i in range(6):
qcirc.measure(q[i], c[i])

# Start Simulation
simulate = qp.execute(["qcirc"], backend=backend, shots=128, timeout=600)
print(simulate.get_counts("qcirc"))
# ヒストグラム表示
plot_histogram(simulate.get_counts("qcirc"))

if __name__ == '__main__':
fire.Fire()

同ディレクトリからこのファイル(hoge.py)を次のように実行するとヒストグラムが表示されます。

$ python .\hoge.py shor_57

何度か計算しましたが,約50%の確率で観測値$s=32$が得られます。ローカルにて計算しているので時間が結構かかります。shots,timeoutを調整しています。


素因数分解の結果を出力

上記に追加して,計算結果から位数を計算する部分を下に追記します。

def reduction(p, q):

common = gcd(p, q)
return (p // common, q // common)

# 2進数表記 -> 10進数表記

integer = []
for i in simulate.get_counts("qcirc"):
integer.append(int(i,2))
print(integer)

# 最大公約数の計算
for o in integer:
if o:
# 約分
p, q = reduction(o, 2**6)
# 有理数で表現
a = Fraction(p, q)
# 分母がNの場合 次へ
if a.denominator == 64: continue
print(o, a)
# 分母を取得しそれが偶数ならば最大公約数を求める
if 1 & a.denominator:
pass
else:
print(int(gcd(37**(a.denominator/2)-1,57)), int(gcd(37**(a.denominator/2)+1,57)))
else:
# 0のとき解なし
pass

最終的に,完成版は次のようになります。

import fire

import math # for pi
from fractions import gcd, Fraction # for shor
from qiskit import QuantumCircuit, QuantumProgram
# Import basic plotting tools
from qiskit.tools.visualization import plot_histogram

def input_state(circ, q, n):
for i in range(n):
circ.h(q[i])

def rev_qreg(circ, q, n):
for i in range(int(n/2)):
circ.swap(q[i], q[(n-1)-i])

def iqft(circ, q, n):
# 上位ビットから処理
for j in range(n-1,-1,-1):
circ.h(q[j])
for k in reversed(range(j)):
circ.cu1((-1)*math.pi/float(2**(j-k)), q[j], q[k])
rev_qreg(circ, q, n)

def reduction(p, q):
common = gcd(p, q)
return (p // common, q // common)

def shor_57():
#--- 量子プログラムの開始
qp = QuantumProgram()
# localで使用するのでAPIの設定割愛
backend = 'local_qasm_simulator'

q = qp.create_quantum_register("q", 6)
a = qp.create_quantum_register("a", 7)
c = qp.create_classical_register("c", 7)
qcirc = qp.create_circuit("qcirc", [q, a], [c])

#--- 初期化
input_state(qcirc, q, 6)
#--- 量子オラクルの設定
qcirc.cx(q[0], a[2])
qcirc.cx(q[0], a[5])
qcirc.x(a[0])
#--- 逆量子フーリエ変換
iqft(qcirc, q, 6)

#--- 測定
for i in range(6):
qcirc.measure(q[i], c[i])

# Start Simulation
simulate = qp.execute(["qcirc"], backend=backend, shots=128, timeout=600)
print(simulate.get_counts("qcirc"))
# ヒストグラム表示,今回は表示させません
#plot_histogram(simulate.get_counts("qcirc"))

# 2進数表記 -> 10進数表記
integer = []
for i in simulate.get_counts("qcirc"):
integer.append(int(i,2))
print(integer)

# 最大公約数の計算
for o in integer:
if o:
# 約分
p, q = reduction(o, 2**6)
# 有理数で表現
a = Fraction(p, q)
# 分母がNの場合 次へ
if a.denominator == 64: continue
print(o, a)
# 分母を取得しそれが偶数ならば最大公約数を求める
if 1 & a.denominator:
pass
else:
print(int(gcd(37**(a.denominator/2)-1,57)), int(gcd(37**(a.denominator/2)+1,57)))
else:
# 0のとき解なし
pass

if __name__ == '__main__':
fire.Fire()

最後にもう一度計算させてみました。出力結果は下記の通りです。求めたい解(3,19)が約50%の確率で得られました。

$ python .\hoge.py shor_57

{'0100000': 61, '0000000': 67}
[32, 0]
32 1/2
3 19


おまけ:1679の素因数分解

ある程度大きな数字だとどれくらいできるだろう,ということで1679を素因数分解してみました。コードは下記の通りです。関数としては,$f(x)=392^x \,\,\,\text{mod}\,\,\,1679$を選択しました。

import fire

import math # for pi
from fractions import gcd, Fraction # for shor
from qiskit import QuantumCircuit, QuantumProgram
# Import basic plotting tools
from qiskit.tools.visualization import plot_histogram

def input_state(circ, q, n):
for i in range(n):
circ.h(q[i])

def rev_qreg(circ, q, n):
for i in range(int(n/2)):
circ.swap(q[i], q[(n-1)-i])

def iqft(circ, q, n):
# 上位ビットから処理
for j in range(n-1,-1,-1):
circ.h(q[j])
for k in reversed(range(j)):
circ.cu1((-1)*math.pi/float(2**(j-k)), q[j], q[k])
rev_qreg(circ, q, n)

def reduction(p, q):
common = gcd(p, q)
return (p // common, q // common)

def shor_1679():
#--- 量子プログラムの開始
qp = QuantumProgram()
# localで使用するのでAPIの設定割愛
backend = 'local_qasm_simulator'

q = qp.create_quantum_register("q", 6)
a = qp.create_quantum_register("a", 11)
c = qp.create_classical_register("c", 11)
qcirc = qp.create_circuit("qcirc", [q, a], [c])

#--- 初期化
input_state(qcirc, q, 6)
#--- 量子オラクルの設定
qcirc.x(a[0])
qcirc.cx(q[0], a[0])
qcirc.cx(q[0], a[3])
qcirc.cx(q[0], a[7])
qcirc.cx(q[0], a[8])
qcirc.cx(q[1], a[0])
qcirc.cx(q[1], a[1])
qcirc.cx(q[1], a[3])
qcirc.cx(q[1], a[5])
qcirc.cx(q[1], a[6])
qcirc.cx(q[1], a[8])
qcirc.cx(q[1], a[9])
qcirc.ccx(q[0], q[1], a[0])
qcirc.ccx(q[0], q[1], a[1])
qcirc.ccx(q[0], q[1], a[2])
qcirc.ccx(q[0], q[1], a[8])
qcirc.ccx(q[0], q[1], a[9])

#--- 逆量子フーリエ変換
iqft(qcirc, q, 6)

#--- 測定
for i in range(6):
qcirc.measure(q[i], c[i])

# Start Simulation
simulate = qp.execute(["qcirc"], backend=backend, shots=32, timeout=1800)
print(simulate.get_counts("qcirc"))
# ヒストグラム表示,今回は表示させません
plot_histogram(simulate.get_counts("qcirc"))

# 2進数表記 -> 10進数表記
integer = []
for i in simulate.get_counts("qcirc"):
integer.append(int(i,2))
print(integer)

# 最大公約数の計算
for o in integer:
if o:
# 約分
p, q = reduction(o, 2**6)
# 有理数で表現
a = Fraction(p, q)
# 分母がNの場合 次へ
if a.denominator == 2**6: continue
print(o, a)
# 分母を取得しそれが偶数ならば最大公約数を求める
if 1 & a.denominator:
pass
else:
print(int(gcd(392**(a.denominator/2)-1,1679)), int(gcd(392**(a.denominator/2)+1,1679)))
else:
# 0のとき解なし
pass

if __name__ == '__main__':
fire.Fire()

ローカルで計算すると時間がかなりかかりますが,下記のような出力になります。

(23, 73)が得られていることが確認できますね。

$ python .\hoge2.py shor_1679

{'00000000000': 7, '00000110000': 10, '00000100000': 7, '00000010000': 8}
[0, 48, 32, 16]
48 3/4
23 73
32 1/2
23 1
16 1/4
23 73


おわりに

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


  • QISKitを用いて57を素因数分解

  • QISKitを用いて1679を素因数分解

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


参考文献


  • 量子アルゴリズム

  • クラウド量子計算入門

  • クラウド量子計算 量子アセンブラ入門


来歴

2018.01.25 新規作成

2018.01.26 計算例に1679の素因数分解を追記