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

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

[前回の記事](https://qiita.com/converghub/items/13e4022995b6668aba87)では，量子コンピュータ（シミュレータ）で素因数分解しました。
<b>（※2018/01/24付の前回記事中，逆量子フーリエ変換の実装で，不適切な箇所があり訂正しました。申し訳ありません。）<b>

ここまでいくつかの計算例を記事にしてきましたが，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 - 公式サイト](https://www.qiskit.org/#/ja)
* [QISKit - マニュアル](https://www.qiskit.org/documentation/ja/)

QISKitでの実装は@hiroyuki827さんの[この記事](https://qiita.com/hiroyuki827/items/77b9922ed8acba96df31)や，@shiibassさんの[この記事](https://qiita.com/shiibass/items/315921c68875a5a1ba33)等があります。
他にもQISKitは[この記事](https://qiita.com/tsujishin/items/eaad0beedff048cb7491)や[この記事](https://qiita.com/kyamaz/items/eebd0f7734c8a33de2ad)等でも少し言及されています。

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



## インストール他

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

```
$ pip install qiskit
```

また，個人的にPython-fireを用いて開発中なこともあり，それもインストールします。
Python-fireに概要ついては[こちらの記事](https://qiita.com/koga1020/items/0c51e365612689615c91)が端的かつ分かりやすかったです。

```
$ pip install fire
```

# QISKitでの実装
## 使用するレジスタ
前回の記事と同様に，q[6]，a[7]，c[7]の3つのレジスタを考えます。

```python3
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には各量子ビットにアダマール変換ゲートを適用しておきます。

```python3
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ゲートを利用して，変換前の順序に合わせておきます。

```python3
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）は（量子回路のオブジェクト，量子レジスタのオブジェクト，量子ビット数）です。

## 量子オラクル
前回の記事を参考に，コンパイル版の回路を組みます。

```python3
#--- 量子オラクルの設定
qcirc.cx(q[0], a[2])
qcirc.cx(q[0], a[5])
qcirc.x(a[0])
```

## 測定とヒストグラム表示
あとはこれらを組み合わせてZ基底で測定し，データをまとめるのみです。
次のようになります。

```python3
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
```

<img src="https://qiita-image-store.s3.amazonaws.com/0/104268/2f72de51-b0a2-045f-c85b-85356fbcbded.png" width=50%>

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

## 素因数分解の結果を出力
上記に追加して，計算結果から位数を計算する部分を下に追記します。

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

```python3
# 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
```

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

```python3
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$を選択しました。

```python3
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の素因数分解を追記