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

はじめに

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

前回の記事では，量子コンピュータ（シミュレータ）で素因数分解しました。
（※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 - 公式サイト
• 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の素因数分解を追記