$$
\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}}
$$
はじめに
前回に引き続き、算術演算シリーズです。今回は「剰余加算」です。アルゴリズムを説明した後、自作の量子計算シミュレータqlazyで、動作の確認をします。
参考にさせていただいた論文・記事は以下です。
- V. Vedral,A. Barenco,A. Ekert; "Quantum Networks for Elementary Arithmetic Operations" (arXiv:quant-ph/9511018)
- 量子コンピュータ（シミュレータ）でモジュロ加算器を作る(@converghub)
- 量子コンピュータ（シミュレータ）でモジュロ加算器を作る【QISKit編】(@converghub)
剰余加算の実現方法
参考論文に全体の回路図が出ているので、まずそれを掲載します。
論文には、
\ket{a,b} \rightarrow \ket{a,a+b \space mod \space N}
を計算する回路と説明されています(ただし、$0 \leqq a,b < N$です)。これで本当に剰余加算が実現できるのでしょうか。図を参照しながら1ステップずつ見ていきます。
まずはじめに、$a,b,N$用のレジスタおよび補助的なレジスタを用意します。図の左側に「$a,b,N,0$」と記号が書かれています。各々何ビットのレジスタにするかは、注意深く決める必要があります。適当に決めてしまうと、入力の値によって途中の計算でオーバーフローしておかしな結果を出してしまう可能性があるからです(言うまでもないことですが)。$a,b$は加算する値で、その結果は$b$のレジスタ(以後、レジスタ<b>と書くことにします)に入りますので、$b$のビット数は$a$よりも1つ多くしておく必要があります。また、今回の前提は、$0 \leqq a,b < N$だったので、$N$のビット数は計算したいビット数よりも1つ多くしておく必要があります。さらに、後で説明するようにレジスタ\とレジスタ\は途中スワップさせますので、同じビット数にしておいた方が良いです。ということを、すべて合わせて考慮すると、計算したいビット数$n$に対して、レジスタ\には$n+1$ビット、レジスタ\には$n+2$ビット、レジスタ\には$n+1$ビットを割り当てておけば良いことになります。それから図には明示的に記載されていないのですが、ADDERを実現するためには、[前回の記事](https://qiita.com/SamN/items/814efecbd9b6aa72d174)で説明したように、桁上げ情報を格納しておくためのレジスタ\も必要です。レジスタ\と同じ数の量子ビットを用意しておいて最上位ビットは\と共有しますので、さらに追加で$n+1$ビット必要になります論文には、このあたり逐一説明されていないのですが、そういうことになるのかなと思っています。間違っていたらご指摘ください)。この考え方に基づくと、結局、$n$ビット同士の剰余加算を実行するためには、$(n+1)+(n+2)+(n+1)+(n+1)+1=4n+6$ビット必要になります。
さて、それでは実際の量子計算について見ていきます。図中の番号(1)から(10)に至るまで、状態がどう変化するかを以下の表に示しました。ここで、上の行から順番にレジスタ<a>,<b>,<N>,<t>の状態を表しています(<t>は図中一番下の補助量子ビットを表しています)。また、この表は状態の変化を表しているので、ケット記号を使って、例えば$\ket{a}$等と書くべきなのですが、面倒なので単に$a$とのみ書くことにしました。
| 初期状態 | (1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) |
|---|---|---|---|---|---|---|---|---|---|---|
| $a$ | $a$ | $N$ | $N$ | [↑]$0$ [↓]$N$ |
[↑]$0$ [↓]$N$ |
$N$ | $a$ | $a$ | $a$ | $a$ |
| $b$ | $a+b$ | $a+b$ | $a+b-N$ | [↑]$a+b-N$ [↓]$a+b-N$ |
[↑]$a+b-N$ [↓]$a+b$ |
[↑]$a+b-N$ [↓]$a+b$ |
[↑]$a+b-N$ [↓]$a+b$ |
[↑]$b-N$ [↓]$b$ |
[↑]$b-N$ [↓]$b$ |
[↑]$a+b-N$ [↓]$a+b$ |
| $N$ | $N$ | $a$ | $a$ | $a$ | $a$ | $a$ | $N$ | $N$ | $N$ | $N$ |
| $t=0$ | $0$ | $0$ | $0$ | [↑]$1$ [↓]$0$ |
[↑]$1$ [↓]$0$ |
[↑]$1$ [↓]$0$ |
[↑]$1$ [↓]$0$ |
[↑]$1$ [↓]$0$ |
$0$ | $0$ |
まず、(1)は、単なるADDERなので、レジスタ<b>に$a+b$が入るだけです。
(2)は、レジスタ<a>とレジスタ<N>のスワップです。$a$と$N$が入れ替わります。
(3)は、ADDERの逆(引き算)です。レジスタ<b>の値$a+b$からレジスタ<a>の値$N$を引いたもの、すなわち$a+b-N$が、レジスタ<b>に入ります。ここで、$a+b$と$N$の大小関係によって、レジスタ<b>の最上位ビットが1になるか0になるかが変わります。$a+b-N$が負の数になる場合、アンダーフロー状態となり、レジスタ<b>の最上位ビットは1になります。そうでない場合は0です。以後、どっちのパターンかによって、場合分けして考えていきます。
(4)は、2つのことを行っています。まず、レジスタ<b>に対する操作についてです。レジスタ<b>の最上位ビットが0だった場合(つまり、$a+b \geqq N$)、レジスタ<t>が1に変化し、そうでない場合(つまり、$a+b < N$)、レジスタは0のままです。表の中で、[↑]と記載したのは$a+b \geqq N$の場合、[↓]と記載したのは、$a+b < N$の場合です。次に、レジスタ\に対する操作についてです。ここで、レジスタ\に矢印がささったような制御ユニタリ系のちょっと見かけない記号が出てきました。論文によると、制御ビット(レジスタ)が1の場合、レジスタ\を$\ket{0}$にし、そうでない場合そのまま何もしない演算子とのことです(ここでは、「アロー演算子」と呼ぶことにします)。結局、$a+b \geqq N$の場合([↑])、レジスタ\は$0$、レジスタ\は$a+b-N$(アンダーフローなし)、$a+b < N$の場合([↓])、レジスタ\は$N$、レジスタ\は$a+b-N$(アンダーフローあり)となります。
(5)は、単なるADDERです。[↑]/[↓]の場合分けによって、レジスタ<b>が異なる変化をします。[↑]の場合、$0$を足すので、$a+b-N$のまま、[↓]の場合、$N$を足すので、アンダーフロー状態が解消されて、$a+b$となります。
(6)は、レジスタ<a>に対するアロー演算子です。レジスタ<t>の制御ビットがどっちであったとしても、レジスタ<a>には$N$が入ることになります。
(7)は、単なるスワップです。レジスタ<a>とレジスタ<N>の値が入れ替わります。
(8)は、逆ADDER(引き算)です。ちょっとトリッキーなのですが、これの意図はレジスタ<t>を$0$に戻すことです。レジスタ<b>から$a$を引くことで、[↑]$b-N$、[↓]$b$となります。すると、[↑]の場合、レジスタ<b>はアンダーフロー状態になり、最上位ビットに1が立ちます。[↓]の場合、最上位ビットは0のままです。この状態にしておいて、次のステップに行きます。
(9)は、レジスタ<b>の最上位ビットによって、レジスタ<t>の状態が変わる回路です。前段でのレジスタ<b>の最上位ビットが1の場合のみ、Xゲートが働くので、レジスタ<t>の値は0に戻ることになります。
(10)は、ADDERです。(8)(9)でトリッキーなことをやったので、レジスタ<b>をもとに戻します。[↑]の場合、$a+b-N$、[↓]の場合、$a+b$となります。これで、$a+b \space mod \space N$が実現できたことになります。
シミュレータで動作確認
さて、それではシミュレータで、この剰余加算の動作を確認してみます。今回想定したのは、3ビット同士の剰余加算です。具体的には、b=5,N=9を固定にして、aを取りうる値すべて(0,1,2,...,7)の重ね合わせとして、一気に剰余加算をやってみます。全体のPythonコードは以下の通りです。
【2021.9.5追記】qlazy最新版でのソースコードはここに置いてあります。
from qlazypy import QState
def swap(self,id_0,id_1):
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 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 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 arrow(self,N,q,id):
for i in range(len(id)):
if (N>>i)%2 == 1:
self.cx(q,id[i])
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 encode(self,decimal,id):
for i in range(len(id)):
if (decimal>>i)%2 == 1:
self.x(id[i])
return self
def create_register(digits):
num = 0
id_a = [i for i in range(digits+1)]
num += len(id_a)
id_b = [i+num for i in range(digits+2)]
num += len(id_b)
id_c = [i+num for i in range(digits+2)]
id_c[digits+1] = id_b[digits+1] # share the qubit id's
num += (len(id_c)-1)
id_N = [i+num for i in range(digits+1)]
num += len(id_N)
id_t = [i+num for i in range(1)]
num += len(id_t)
id_r = id_b[:-1] # to store the result
return (num,id_a,id_b,id_c,id_N,id_t,id_r)
def superposition(self,id):
for i in range(len(id)-1):
self.h(id[i])
return self
def result(self,id_a,id_b):
# measurement
id_ab = id_a + id_b
iid_ab = id_ab[::-1]
freq = self.m(id=iid_ab).frq
# set results
a_list = []
r_list = []
for i in range(len(freq)):
if freq[i] > 0:
a_list.append(i%(2**len(id_a)))
r_list.append(i>>len(id_a))
return (a_list,r_list)
if __name__ == '__main__':
# add methods
QState.swap = swap
QState.sum = sum
QState.i_sum = i_sum
QState.carry = carry
QState.i_carry = i_carry
QState.arrow = arrow
QState.plain_adder = plain_adder
QState.i_plain_adder = i_plain_adder
QState.modular_adder = modular_adder
QState.encode = encode
QState.superposition = superposition
QState.result = result
# create registers
digits = 3
num,id_a,id_b,id_c,id_N,id_t,id_r = create_register(digits)
# set iput numbers
b = 5
N = 9
# initialize quantum state
qs = QState(num)
qs.superposition(id_a) # set superposition of |0>,|1>,..,|7> for |a>
qs.encode(b,id_b)
qs.encode(N,id_N)
# execute modular adder
qs.modular_adder(N,id_a,id_b,id_c,id_N,id_t)
a_list,r_list = qs.result(id_a,id_r)
for i in range(len(a_list)):
print("{0:}+{1:} mod {2:} -> {3:}".format(a_list[i],b,N,r_list[i]))
qs.free()
何をやっているか順に説明します。
def swap(self,id_0,id_1):
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
は、2つのレジスタのスワップです。swapゲートを対象となる量子ビット数分繰り返しています。
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
で、ADDERの逆演算を定義しています。関数plain_adderの処理を逆に実施しているだけです。
def arrow(self,N,q,id):
for i in range(len(id)):
if (N>>i)%2 == 1:
self.cx(q,id[i])
return self
で、アロー演算子を定義しています。実は、この実装をどうすればよいか、ちょっと悩みました。@converghubさんの記事では、Nの値を格納しておくための補助量子ビットを追加で用意した上で、Toffoliゲートを複数使って、実現されていました。が、論文には、複数のCNOTゲートで実現できると書いてあるので、補助量子ビットなしで何とかなるんだろと思い、考えてみました。結局、Nはこの量子計算を行うにあたり、はじめから与えられている値なので、Nのビット配列のどこに1が立っているのかを見て、CNOTゲートのターゲットとなる量子ビット番号を決めて、それを複数適切に並べれば良いということに思い至りました。それを上のように実装してみました。
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
で、剰余加算を定義しています。最初に示した回路図をそのまま実現しています。
その他の関数sum、i_sum、carry、i_carry、plain_adder、encode、superposition、resultについては、前回の記事で説明したので、省略します。関数create_registerは、前回と実装が違っていますが、今回の計算に合わせてレジスタを用意するものです。こちらも説明省略します。
というわけで、プログラム本体の説明です。
# create registers
digits = 3
num,id_a,id_b,id_c,id_N,id_t,id_r = create_register(digits)
レジスタ<a>,<b>,<N>,<t>、ADDERで必要になる桁上げ用補助レジスタ<c>、剰余加算のための補助レジスタ<t>、および結果を取得するレジスタ<r>(=測定するレジスタ=レジスタ<b>)を決め、全体で必要になる量子ビット数も合わせて出力し、各変数に格納しておきます。
# set iput numbers
b = 5
N = 9
# initialize quantum state
qs = QState(num)
qs.superposition(id_a) # set superposition of |0>,|1>,..,|7> for |a>
qs.encode(b,id_b)
qs.encode(N,id_N)
今回、固定入力するb,Nの値を設定し、aについては取りうる値すべての重ね合わせ状態とします。
# execute modular adder
qs.modular_adder(N,id_a,id_b,id_c,id_N,id_t)
剰余加算を実行します。
a_list,r_list = qs.result(id_a,id_r)
for i in range(len(a_list)):
print("{0:}+{1:} mod {2:} -> {3:}".format(a_list[i],b,N,r_list[i]))
測定によって得られた結果から、レジスタ<a>の状態とレジスタ<r>(=レジスタ<b>)の状態をセットにして表示します。
実行結果は、以下の通りです。
4+5 mod 9 -> 0
5+5 mod 9 -> 1
6+5 mod 9 -> 2
7+5 mod 9 -> 3
0+5 mod 9 -> 5
1+5 mod 9 -> 6
2+5 mod 9 -> 7
3+5 mod 9 -> 8
というわけで、剰余加算が正しくできることが確認できました。
おわりに
前回「加算」、今回「剰余加算」が確認できましたので、次回は「制御剰余乗算」(論文には"Controlled-multiplier modulo N"と書いてありましたが、日本語これでいい？)に挑戦です。いままでの回路を部品にして、組み合わせていくようなので、だんだん量子ビット数的に厳しくなっていく予感がありますが、できる範囲で頑張ってみたいと思います。
以上