量子アルゴリズムの基本：算術演算の確認（べき剰余）

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

## はじめに

[加算](https://qiita.com/SamN/items/814efecbd9b6aa72d174)、[剰余加算](https://qiita.com/SamN/items/39016dd3e0567df81b21)、[制御-剰余乗算](https://qiita.com/SamN/items/fc350338c2222561488b)という具合に見てきましたが、それは「Shorのアルゴリズム」で重要な役割を果たす「べき剰余」を実行するための準備として、という位置づけでした。今回はその「べき剰余」について見ていきます。アルゴリズムを説明した後、自作の量子計算シミュレータ[qlazy](https://github.com/samn33/qlazy)で、動作を確認します。

参考にさせていただいた論文は以下です。

* [V. Vedral,A. Barenco,A. Ekert; "Quantum Networks for Elementary Arithmetic Operations" (arXiv:quant-ph/9511018)](https://arxiv.org/abs/quant-ph/9511018)

## べき剰余の実現方法

参考論文に全体の回路図が出ているので、まずそれを掲載します。


<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391058/1ad92cf4-e701-96ed-cf92-4dc05b0ad021.png" width=150%>

$a^x \space mod \space N$という「べき剰余」を実行する回路です(図には$x^a \space mod \space N$と書いてありますが、$a^x \space mod \space N$の間違いではないかと思われます)。

図を見てわかる通り、[前回の記事](https://qiita.com/SamN/items/fc350338c2222561488b)で確認した「制御-剰余乗算(以下CMMと省略)」とスワップ(以下SWAPと表記)とCMMの逆演算(以下iCMMと表記)をセットにして繰り返し適用していくだけなので、考え方はそれほど難しくないです。

説明のため、レジスタに名前をつけます。図の左側の上にxと記載してあるレジスタを\<x\>、1と記載してあり最終結果が格納されるレジスタを\<xx\>、一番下の0と記載されている レジスタを\<y\>と名付けておきます。aを格納するレジスタが記載されていないですが、この量子計算を実現するために外から与えるパラメータのようなもので、あえて言うなら古典レジスタに格納されていると思えば良いと思います。その他、CMMを実現するために必要になるレジスタもありますが、詳細は[前回の記事](https://qiita.com/SamN/items/fc350338c2222561488b)を参照してください。

では、左から順に、レジスタ\<xx\>とレジスタ\<y\>の状態変化を表にして示します。状態なのでケット記号を使って$\ket{A}$のように記載すべきなのですが、簡単のため$A$とのみ記載します。また、各段階での計算は$mod \space N$前提で実行しているものと思ってください。xの値によって、各段階のCMMやSWAPやiCMMがONになったりOFFになったりするので、場合分けをして表にしてみました。

* $x=0$ ($x_{0}=0,x_{1}=0$)の場合

|レジスタ|初期状態|0番目のCMM|SWAP|0番目のiCMM|1番目のCMM|SWAP   |1番目のiCMM|
|--------|--------|----------|----|-----------|----------|-------|-----------|
|\<xx\>  |$1$     |$1$       |$1$ |$1$        |$1$       |$1$    |$1$        |
|\<y\>   |$0$     |$0$       |$0$ |$0$        |$0$       |$0$    |$0$        |

* $x=1$ ($x_{0}=1,x_{1}=0$)の場合

|レジスタ|初期状態|0番目のCMM|SWAP|0番目のiCMM|1番目のCMM|SWAP   |1番目のiCMM|
|--------|--------|----------|----|-----------|----------|-------|-----------|
|\<xx\>  |$1$     |$a$       |$1$ |$a$        |$a$       |$a$    |$a$        |
|\<y\>   |$0$     |$1$       |$a$ |$0$        |$0$       |$0$    |$0$        |

* $x=2$ ($x_{0}=0,x_{1}=1$)の場合

|レジスタ|初期状態|0番目のCMM|SWAP|0番目のiCMM|1番目のCMM|SWAP   |1番目のiCMM|
|--------|--------|----------|----|-----------|----------|-------|-----------|
|\<xx\>  |$1$     |$1$       |$1$ |$1$        |$1$       |$a^{2}$|$a^{2}$    |
|\<y\>   |$0$     |$0$       |$0$ |$0$        |$a^{2}$   |$1$    |$0$        |

* $x=3$ ($x_{0}=1,x_{1}=1$)の場合

|レジスタ|初期状態|0番目のCMM|SWAP|0番目のiCMM|1番目のCMM|SWAP   |1番目のiCMM|
|--------|--------|----------|----|-----------|----------|-------|-----------|
|\<xx\>  |$1$     |$a$       |$1$ |$a$        |$a$       |$a^{3}$|$a^{3}$    |
|\<y\>   |$0$     |$1$       |$a$ |$0$        |$a^{3}$   |$a$    |$0$        |


というわけで、レジスタ\<xx\>に$a^{x} \space mod \space N$の値が入ることがわかります。ここで1点注意しておきたいのがSWAPの動作です。回路図のSWAPには制御ビットとつながる線が書かれていませんが、CMM-SWAP-iCMMが一体となって制御ビットとつながっていることを意味しています。つまり、SWAPは厳密に言うと、制御SWAP(正式名称は知りません。フレドキンゲートの拡張みたいですがちょっと違うかな?)です。

さて、"CMM-SWAP-iCMM"のセットを2回繰り返したところまで書けたので、あとはこれを延長すれば、任意のべき剰余が計算できるというわけです。

数式で説明すると、以下のような量子計算を$i=0,1,2,...$と繰り返していると思えば良いです。

```math
\begin{align}
& \ket{a^{2^0 x_0 + 2^1 x_1 + \cdots 2^{i-1} x_{i-1}}, 0} \\
& \rightarrow \ket{a^{2^0 x_0 + 2^1 x_1 + \cdots 2^{i-1} x_{i-1}},a^{2^0 x_0 + 2^1 x_1 + \cdots 2^{i} x_{i}}} \\
& \rightarrow \ket{a^{2^0 x_0 + 2^1 x_1 + \cdots 2^{i} x_{i}},a^{2^0 x_0 + 2^1 x_1 + \cdots 2^{i-1} x_{i-1}}} \\
& \rightarrow \ket{a^{2^0 x_0 + 2^1 x_1 + \cdots 2^{i} x_{i}}, 0}
\end{align}
```

ここで、1行目から2行目はCMM、2行目から3行目はSWAP、3行目から4行目はiCMMの演算を表しています。これを$n$回繰り返すと、最大で$x=2^{n}-1$までのべき剰余$a^{x} \space mode \space N$が実行できることがわかります。

## シミュレータで動作確認

それでは、シミュレータで動作確認をしていくのですが、その前に各レジスタに必要となるビット数を見積もってみます。上の図だけ見ると、とてもシンプルに見えるのですが、CMMやiCMMを実現するためには内部で「剰余加算」を実行していて、「剰余加算」を実行するために内部で「加算」を実行していて、各々を実行するために多くの補助レジスタが必要になります。シミュレータでの実行は限定されたメモリ範囲内でやるしかないため、今回はとても小さい数での動作確認になります。具体的には、$a=2,N=3,x=0,1,2$で動作させることを前提にします。

まず、$a=2$については外部パラメータということで量子レジスタとして割り当てることはしません。$x$については2ビットで十分なので、レジスタ\<x\>には2ビットを割り当てます。レジスタ\<xx\>には、$mod \space 3$の結果が格納されるので2ビット割り当てることにします。CMMを実行するために、内部で剰余加算を実行しますが、それへの入力レジスタが別途必要になり、今回レジスタ\<xxx\>と名付けることにします(前回の制御-剰余乗算の記事では\<xx\>と呼んでいました)。\<xxx\>には最大で$a^2=4$が入るので、3ビット割り当てます。\<y\>は、CMMの内部で剰余加算をするもう一方の入力レジスタなので、\<xxx\>よりも1つ多いビット数を用意しておく必要があるため、4ビットとします。加算を実行するためには桁上げ情報格納のためのレジスタ\<c\>がさらに必要になります。それには\<y\>と同じ数のビットを用意しておき、最上位の1ビットを\<y\>と共有するため、追加で3ビット必要になります。$N=3$を格納するためのレジスタも、上の回路図には書かれていないですが、必要です。剰余加算回路の内部で\<xxx\>とスワップするため同じビット数、すなわち3ビットを割り当てます。その他、剰余加算でアンダーフロー記述のための補助レジスタが1ビット必要になります。以上、すべてを合わせると「18量子ビット」となります。ちょっとややこしい説明になってしまいましたが、前回までの、[加算](https://qiita.com/SamN/items/814efecbd9b6aa72d174)、[剰余加算](https://qiita.com/SamN/items/39016dd3e0567df81b21)、[制御-剰余乗算](https://qiita.com/SamN/items/fc350338c2222561488b)の記事をよく見ればわかると思います(と言いつつ、間違っていたらすみません、汗)。

それではシミュレータで、「べき剰余」の動作を確認してみます。先程説明したように、$a=2,N=3,x=0,1,2$で動作させます。x=3の場合は、今回ビット数の関係で諦めましたので、xのすべての値に対して重ね合わせることはせず、$x=0,1,2$でforループを回して結果を逐一表示するようにしました。全体のPythonコードは以下の通りです。

【2021.9.5追記】[qlazy](https://github.com/samn33/qlazy)最新版でのソースコードは[ここ](https://github.com/samn33/qlazy/blob/master/example/py/Arithmetic/modular_exponentiation.py)に置いてあります。

```python
from qlazypy import QState

def swap(self,id_0,id_1):
Arithmetic/ctr_modular_multiplier_1.py
    dim = min(len(id_0),len(id_1))
    for i in range(dim):
        self.cx(id_0[i],id_1[i]).cx(id_1[i],id_0[i]).cx(id_0[i],id_1[i])
    return self
    
def ctr_swap(self,id_ctr,id_0,id_1):

    dim = min(len(id_0),len(id_1))
    for i in range(dim):
        self.ccx(id_ctr[0],id_0[i],id_1[i])
        self.ccx(id_ctr[0],id_1[i],id_0[i])
        self.ccx(id_ctr[0],id_0[i],id_1[i])
    return self
    
def sum(self,q0,q1,q2):

    self.cx(q1,q2).cx(q0,q2)
    return self
    
def i_sum(self,q0,q1,q2):

    self.cx(q0,q2).cx(q1,q2)
    return self
    
def carry(self,q0,q1,q2,q3):

    self.ccx(q1,q2,q3).cx(q1,q2).ccx(q0,q2,q3)
    return self

def i_carry(self,q0,q1,q2,q3):

    self.ccx(q0,q2,q3).cx(q1,q2).ccx(q1,q2,q3)
    return self

def plain_adder(self,id_a,id_b,id_c):

    depth = len(id_a)
    for i in range(depth):
        self.carry(id_c[i],id_a[i],id_b[i],id_c[i+1])
    self.cx(id_a[depth-1],id_b[depth-1])
    self.sum(id_c[depth-1],id_a[depth-1],id_b[depth-1])
    for i in reversed(range(depth-1)):
        self.i_carry(id_c[i],id_a[i],id_b[i],id_c[i+1])
        self.sum(id_c[i],id_a[i],id_b[i])
    return self

def arrow(self,N,q,id):

    for i in range(len(id)):
        if (N>>i)%2 == 1:
            self.cx(q,id[i])
    return self

def arrow2(self,i,a,id_ctr,id_x,id_xx):

    aa = a * 2**i
    for j in range(len(id_xx)):
        if (aa>>j)%2 == 1:
            self.ccx(id_ctr[0],id_x[i],id_xx[j])
    return self
    
def arrow3(self,id_ctr,id_x,id_y):

    for j in range(len(id_x)):
        self.ccx(id_ctr[0],id_x[j],id_y[j])
    return self
    
def i_plain_adder(self,id_a,id_b,id_c):

    depth = len(id_a)
    for i in range(depth-1):
        self.i_sum(id_c[i],id_a[i],id_b[i])
        self.carry(id_c[i],id_a[i],id_b[i],id_c[i+1])
    self.i_sum(id_c[depth-1],id_a[depth-1],id_b[depth-1])
    self.cx(id_a[depth-1],id_b[depth-1])
    for i in reversed(range(depth)):
        self.i_carry(id_c[i],id_a[i],id_b[i],id_c[i+1])
    return self

def modular_adder(self,N,id_a,id_b,id_c,id_N,id_t):

    self.plain_adder(id_a,id_b,id_c)
    self.swap(id_a,id_N)
    self.i_plain_adder(id_a,id_b,id_c)
    self.x(id_b[len(id_b)-1])
    self.cx(id_b[len(id_b)-1],id_t[0])
    self.x(id_b[len(id_b)-1])
    self.arrow(N,id_t[0],id_a)
    self.plain_adder(id_a,id_b,id_c)
    self.arrow(N,id_t[0],id_a)
    self.swap(id_a,id_N)
    self.i_plain_adder(id_a,id_b,id_c)
    self.cx(id_b[len(id_b)-1],id_t[0])
    self.plain_adder(id_a,id_b,id_c)
    return self

def i_modular_adder(self,N,id_a,id_b,id_c,id_N,id_t):

    self.i_plain_adder(id_a,id_b,id_c)
    self.cx(id_b[len(id_b)-1],id_t[0])
    self.plain_adder(id_a,id_b,id_c)
    self.swap(id_a,id_N)
    self.arrow(N,id_t[0],id_a)
    self.i_plain_adder(id_a,id_b,id_c)
    self.arrow(N,id_t[0],id_a)
    self.x(id_b[len(id_b)-1])
    self.cx(id_b[len(id_b)-1],id_t[0])
    self.x(id_b[len(id_b)-1])
    self.plain_adder(id_a,id_b,id_c)
    self.swap(id_a,id_N)
    self.i_plain_adder(id_a,id_b,id_c)
    
    return self

def ctr_modular_multiplier(self,a,N,id_ctr,id_x,id_xx,id_y,id_c,id_N,id_t):

    depth = len(id_x)
    for i in range(depth):
        self.arrow2(i,a,id_ctr,id_x,id_xx)
        self.modular_adder(N,id_xx,id_y,id_c,id_N,id_t)
        self.arrow2(i,a,id_ctr,id_x,id_xx)
    self.x(id_ctr[0])
    self.arrow3(id_ctr,id_x,id_y)
    self.x(id_ctr[0])
    return self
    
def i_ctr_modular_multiplier(self,a,N,id_ctr,id_x,id_xx,id_y,id_c,id_N,id_t):

    depth = len(id_x)
    self.x(id_ctr[0])
    self.arrow3(id_ctr,id_x,id_y)
    self.x(id_ctr[0])
    for i in reversed(range(depth)):
        self.arrow2(i,a,id_ctr,id_x,id_xx)
        self.i_modular_adder(N,id_xx,id_y,id_c,id_N,id_t)
        self.arrow2(i,a,id_ctr,id_x,id_xx)
    return self

def modular_exponentiation(self,a,N,id_x,id_xx,id_xxx,id_y,id_c,id_N,id_t):

    depth = len(id_x)
    id_ctr = [0]
    for i in range(depth):
        aa = a**(2**i)
        id_ctr[0] = id_x[i]
        self.ctr_modular_multiplier(aa,N,id_ctr,id_xx,id_xxx,id_y,id_c,id_N,id_t)
        self.ctr_swap(id_ctr,id_xx,id_y)
        self.i_ctr_modular_multiplier(aa,N,id_ctr,id_xx,id_xxx,id_y,id_c,id_N,id_t)
    return self
    
def encode(self,decimal,id):

    for i in range(len(id)):
        if (decimal>>i)%2 == 1:
            self.x(id[i])
    return self

def decode(self,id):

    iid = id[::-1]
    return self.m(id=iid,shots=1).lst

def create_register():

    num = 0
    id_x = [i for i in range(2)]
    num += len(id_x)
    id_xx = [i+num for i in range(2)]
    num += len(id_xx)
    id_xxx = [i+num for i in range(3)]
    num += len(id_xxx)
    id_y = [i+num for i in range(4)]
    num += len(id_y)
    id_c = [i+num for i in range(4)]
    id_c[len(id_c)-1] = id_y[len(id_y)-1]
    num += (len(id_c)-1)
    id_N = [i+num for i in range(3)]
    num += len(id_N)
    id_t = [i+num for i in range(1)]
    num += len(id_t)
    id_r = id_xx

    return (num,id_x,id_xx,id_xxx,id_y,id_c,id_N,id_t,id_r)

if __name__ == '__main__':

    # add methods
    QState.swap = swap
    QState.ctr_swap = ctr_swap
    QState.sum = sum
    QState.i_sum = i_sum
    QState.carry = carry
    QState.i_carry = i_carry
    QState.arrow = arrow
    QState.arrow2 = arrow2
    QState.arrow3 = arrow3
    QState.plain_adder = plain_adder
    QState.i_plain_adder = i_plain_adder
    QState.modular_adder = modular_adder
    QState.i_modular_adder = i_modular_adder
    QState.ctr_modular_multiplier = ctr_modular_multiplier
    QState.i_ctr_modular_multiplier = i_ctr_modular_multiplier
    QState.modular_exponentiation = modular_exponentiation
    QState.encode = encode
    QState.decode = decode
    
    # set input numbers
    a = 2
    N = 3
    x_list = [0,1,2]

    # create registers
    num,id_x,id_xx,id_xxx,id_y,id_c,id_N,id_t,id_r = create_register()

    for x in x_list:
        
        # initialize quantum state
        qs = QState(num)
        qs.encode(x,id_x)
        qs.encode(N,id_N)
        qs.encode(1,id_xx)
            
        # execute controlled modular multiplier
        qs.modular_exponentiation(a,N,id_x,id_xx,id_xxx,id_y,id_c,id_N,id_t)
        res = qs.decode(id_r)
        print("{0:}^{1:} mod {2:} -> {3:}".format(a,x,N,res))

        qs.free()
```

前回までの記事([加算](https://qiita.com/SamN/items/814efecbd9b6aa72d174)、[剰余加算](https://qiita.com/SamN/items/39016dd3e0567df81b21)、[制御-剰余乗算](https://qiita.com/SamN/items/fc350338c2222561488b))で、大半の関数は説明済なので、今回の差分だけを説明します。

```python
def ctr_swap(self,id_ctr,id_0,id_1):

    dim = min(len(id_0),len(id_1))
    for i in range(dim):
        self.ccx(id_ctr[0],id_0[i],id_1[i])
        self.ccx(id_ctr[0],id_1[i],id_0[i])
        self.ccx(id_ctr[0],id_0[i],id_1[i])
    return self
```

で、制御ビットつきのスワップを定義しています。通常のスワップは制御NOTを3つ組み合わせますが、さらに別の制御ビットによってON/OFF制御するため、制御-制御NOT、つまりToffoliゲート3つに置き換えています。

```python
def modular_exponentiation(self,a,N,id_x,id_xx,id_xxx,id_y,id_c,id_N,id_t):

    depth = len(id_x)
    id_ctr = [0]
    for i in range(depth):
        aa = a**(2**i)
        id_ctr[0] = id_x[i]
        self.ctr_modular_multiplier(aa,N,id_ctr,id_xx,id_xxx,id_y,id_c,id_N,id_t)
        self.ctr_swap(id_ctr,id_xx,id_y)
        self.i_ctr_modular_multiplier(aa,N,id_ctr,id_xx,id_xxx,id_y,id_c,id_N,id_t)
    return self
```

で、べき剰余を定義しています。上で示した回路図をそのまま実行しています。

```python
def create_register():

    num = 0
    id_x = [i for i in range(2)]
    num += len(id_x)
    id_xx = [i+num for i in range(2)]
    num += len(id_xx)
    id_xxx = [i+num for i in range(3)]
    num += len(id_xxx)
    id_y = [i+num for i in range(4)]
    num += len(id_y)
    id_c = [i+num for i in range(4)]
    id_c[len(id_c)-1] = id_y[len(id_y)-1]
    num += (len(id_c)-1)
    id_N = [i+num for i in range(3)]
    num += len(id_N)
    id_t = [i+num for i in range(1)]
    num += len(id_t)
    id_r = id_xx

    return (num,id_x,id_xx,id_xxx,id_y,id_c,id_N,id_t,id_r)
```

で、各レジスタに割り当てる量子ビット番号と全体の量子ビット数を計算してリターンする関数を定義しています。

というわけで、プログラム本体の説明です。

```python
# set input numbers
a = 2
N = 3
x_list = [0,1,2]

# create registers
num,id_x,id_xx,id_xxx,id_y,id_c,id_N,id_t,id_r = create_register()
```

$a,N$の値と$x$の値を入れるリストを設定して、各レジスタおよび必要となる全体量子ビット数を適当な変数に格納します。

```python
for x in x_list:
    
    # initialize quantum state
    qs = QState(num)
    qs.encode(x,id_x)
    qs.encode(N,id_N)
    qs.encode(1,id_xx)
        
    # execute controlled modular multiplier
    qs.modular_exponentiation(a,N,id_x,id_xx,id_xxx,id_y,id_c,id_N,id_t)
    res = qs.decode(id_r)
    print("{0:}^{1:} mod {2:} -> {3:}".format(a,x,N,res))
```

xの値を変えながら、べき剰余の計算を実行し、その結果を逐次表示します。

実行結果は、以下の通りです。

```
2^0 mod 3 -> 1
2^1 mod 3 -> 2
2^2 mod 3 -> 1
```

というわけで、べき剰余が正しく計算できることが確認できました。


## おわりに

結局、とても小さい数でしか実行できませんでしたが、とりあえず確認できたことにします。今回18量子ビットという見積りのもとで実行しましたが、参考論文には$N$のビット数$n$に対して、$7n+1$量子ビットで良いという記載があります。だとすると、本当は15量子ビットで実行可能なはず、ということでしょうか？ちょっとわかりませんが、もしかすると何らかの誤まった解釈から余分な量子ビットを使ってしまっているのかもしれません(が、すみません。今回はここまでとします。みなさま、お使いのシミュレータで実行してみて、量子ビット削減できたよ！ということが、もしあれば、後で教えて下さい)。

「べき剰余」のアルゴリズムについては、今回参考にした論文が発表された1990年代半ば以降、研究の進展がとてもあるようなので、今後余裕があれば見ておこうと思います。

次回は、予定では「Shorのアルゴリズム」なのですが、べき剰余の部分をどうしようか、考え中です。

以上