可逆論理での箱玉系シミュレーション

量子計算は可逆計算なので、量子計算機上では可逆セル・オートマトンを実現できるはずです。
そう思って検索しているうちに、箱玉系を量子計算で実装する方法について書かれた[論文](https://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/41295/1/1221_05.pdf)に行き当たったので実装してみました。

## 可逆計算の土台を作る

量子計算シミュレータだと扱えるビット数に制限があったり必要 CPU コストが大きくなりがちなうえ、量子回路のデバッグも大変です。
古典可逆計算だと扱えるビット数が多く CPU コストもかからないうえ、回路のデバッグも print 文を適宜挿入することで簡単にできます。

今回使う道具は制御ノットゲート、トフォリゲート、スワップゲートだけなので量子計算シミュレータ上で動かしたからといって特に古典論と違う結果が出るわけではないということもあり、Python 上で古典可逆計算システムを作ってその上で動かすことにします。

古典可逆計算として実装できれば、原理的には同じ手順で量子計算も可能です。

### レジスタ

箱玉系の実装では系のビット列をひとまとめにして扱うのが便利なので、ビット列をひとまとめにして扱うクラス `Register` をつくります。

`Register` インスタンスの一つのビット列上をビットが移動していく形で箱玉系を実現するのが今回のゴールです。

```python
# 現在使用中のレジスタの数
REGCOUNT = 0

class Register(object):
    u'''ビット列をひとつのレジスタとして扱うクラス。
    '''
    NBITS = 32  # レジスタのビット幅は固定
    def __init__(self, title: str) -> None:
        u'''ビット列を確保する。
        全ビットの値は 0。
        '''
        self.title = title
        self.v : int = 0  # このレジスタが保持するビット列
        self.lock_count = 0  # ロック回数
        self.is_deallocated = False
        global REGCOUNT
        REGCOUNT += 1
    def deallocate(self) -> None:
        u'''当該ビット列を返却する。
        全ビットの値は 0 でなければならない。
        '''
        self.check()
        assert(self.v == 0)
        self.is_deallocated = True
        global REGCOUNT
        REGCOUNT -= 1
    def __str__(self) -> str:
        u'''ビット列の内容を表示する。
        '''
        return '{0}({1})'.format(
            'DEALLOCATED' if self.is_deallocated else
            ''.join(reversed([('0' if (self.v & (1 << i)) == 0 else '1')
                              for i in range(Register.NBITS)])),
            self.title)
    def must_be_zero(self) -> None:
        u'''検証用: 当該ビット列の全ビットの値が 0 であることを確認する。
        '''
        self.check()
        assert(self.v == 0)
    def set_locked(self, b: bool) -> None:
        u'''検証用: 引数に True を与えた場合、このレジスタの値の変更を禁止する。
        '''
        if b:
            self.lock_count += 1
        else:
            self.lock_count -= 1
    def is_locked(self) -> bool:
        u'''検証用: このレジスタの値の変更が禁止されているなら True を返す。
        '''
        return self.lock_count > 0
    def check(self) -> None:
        u'''検証用: 解放済みでないか確認する。
        '''
        if self.is_deallocated:
            raise RuntimeError('Register {0} already DEALOOCATED'.format(
                self.title))
    def check_not_locked(self) -> None:
        u'''検証用: 解放済みでなく、行進禁止でもないことを確認する。
        '''
        self.check()
        if self.is_locked():
            raise RuntimeError('Register {0} is locked'.format(self.title))
```

一時変数として扱うビット列を含めて、`Register` インスタンス作成時に値 0 のビット列をまとめて確保し、使い終わった後は 0 クリアしてから `deallocate` メソッドで解放する、という使い方です。

一時変数として使った `Register` をきちんと解放できているかを確認するために、グローバル変数 `REGCOUNT` を設けています。
可逆計算としては、計算を進めて `REGCOUNT` が無限に増大し続けないことが必要です。

また、デバッグ目的・動作検証目的で以下のメソッドを用意しています:

- `must_be_zero` : ビット列の値を 0 にきちんと戻しているか確認する
- `set_locked` / `is_locked` : ビット列を更新禁止とマークすることで、ゲート操作あるいは回路操作での意図しない変更を検知する
- `check` : ゲート実装において、解放済みのビット列を参照・操作したりしていないか確認する
- `check_not_locked` : ゲート実装において、解放済みでないこと、および更新が禁止されていないことを確認する

`set_locked` での更新禁止状態・許可状態切り替えを簡便に間違いなく行うために、ユーティリティとして with 構文で使う `Lock` クラスも準備しました:

```python
class Lock(object):
    u'''検証用: with 区間でレジスタを更新禁止状態にする。
    呼び出した関数でレジスタが破壊されていないことを確認するために用いる。
    '''
    def __init__(self, r: Register) -> None:
        u'''更新禁止対象のレジスタを引数に与える。
        '''
        self.r = r
    def __enter__(self) -> 'Lock':
        self.r.set_locked(True)
        return self
    def __exit__(self,
                 exc_type: Optional[Type[BaseException]],
                 exc_value: Optional[BaseException],
                 traceback: Optional[TracebackType]) -> None:
        self.r.set_locked(False)
```

### 可逆ゲート

`Register` が表すビット列を操作する可逆ゲートとして、量子計算でもおなじみのノットゲート(`X`)、制御ノットゲート(`CX`)およびトフォリゲート(`CCX`)を用意します。

```python
def X(b: Register, index: int) -> None:
    u'''ビット列 b の index 番目の位置のビットを反転する。
    '''
    b.check_not_locked()
    b.v = b.v ^ (1 << index)
def CX(c: Register, b: Register) -> None:
    u'''ビット列 c を制御ビットとして、ビット列 b を反転する。
    '''
    c.check()
    b.check_not_locked()
    b.v = c.v ^ b.v
def CCX(c1: Register, c2: Register, b: Register) -> None:
    u'''ビット列 c1, c2 を制御ビットとして、ビット列 b を反転する。
    '''
    c1.check()
    c2.check()
    b.check_not_locked()
    b.v = (c1.v & c2.v) ^ b.v
```

今回の箱玉系実装では、レジスタ上の右ビットシフト（回転）、左ビットシフト（回転）が必要なので、これも準備しておきます（これは、スワップゲートあるいは制御ノットゲートの組み合わせで書き替え可能です）：

```python
def rshift(b: Register) -> None:
    u'''ビット列 b を右シフト（回転）する。
    '''
    b.check_not_locked()
    size = b.NBITS
    mask = (1 << size) - 1
    b.v = (((b.v << (size - 1)) & mask) |
           (b.v & mask) >> 1)

def lshift(b: Register) -> None:
    u'''ビット列 b を左シフト（回転）する。
    '''
    b.check_not_locked()
    size = b.NBITS
    mask = (1 << size) - 1
    b.v = (((b.v << 1) & mask) |
           (b.v & mask) >> (size - 1))
```

## 箱玉漸化式を使った箱玉系の実現

ある瞬間での箱玉系の状態を $X(t)$ と、次の瞬間の箱玉系を $X(t+1)$ とします。

参照した論文によると箱玉系の状態を $X(t)$ から $X(t+1)$ に進めることは、初期値 $A_0 = B_0 = X$ とした変換

```math
\begin{align}
A_{n+1} &= A_n \lor B_n \\
B_{n+1} &= T_1 ( A_n \land B_n)
\end{align}

```

を $B_n = 0$ になるまで繰り返すことで実現できるそうです（ここで $T_1$ は右方向あるいは左方向へのシフト）。
$X(t+1)$ は $B_n = 0$ のときの $A_n \oplus X$ として得られます。

上記変換は、一旦 $B_n = 0$ に到達するとそれ以上を繰り返しても $A$, $B$ は変化しません。なので箱玉漸化式と呼ばれます。

この漸化式を可逆ゲート（量子ゲート）でどう実装すればよいかも載っているのでそのまま上で実装した可逆ゲートの組み合わせで書いたものが以下です：


```python
def bbs_recurrence_relation_r(a: Register, b: Register, c: Register) -> None:
    u'''右方向に動く箱玉系の漸化式
    '''
    with Lock(a):
        CCX(a, b, c)
        CX(c, b)
        rshift(c)
        CX(a, b)
```

右シフトを左シフトに置き換えると、左方向に動く箱玉系の漸化式になります。
この漸化式を用いて箱玉系 $X(t)$ の次の瞬間の状態 $X(t+1)$ を得る手続きはこんな感じになります：

```python
    # x に現在時刻 t の箱玉系の状態が入っている
    with Lock(x):
        a = x
        b = Register('B')  # allocate
        CX(x, b)
        ts = []
        for i in range(cycle):
            t = Register('temp{0}'.format(i))  # allocate
            bbs_recurrence_relation_r(a, b, t)
            (t, a, b) = (a, b, t)  # rename
            ts.append(t)
        # 時刻 t+1 での x の値を x2 に得る
        CX(x, x2)
        CX(a, x2)
```

古典可逆計算なので毎回 b の値を観測して値が 0 になるまで繰り返す実装にしてもよいのですが、一応量子計算に移行しやすいように決め打ちで（ここでは cycle として与えている）十分な回数繰り返すことにしています。

これで箱玉系は実現できるのですが、このままだと `ts` に一時変数が溜まったままになります。不可逆計算ならゴミとしてそのまま捨ててしまってよいのですが、可逆計算としてはよろしくありません。

`ts` 上の変数を 0 に純化する方法も論文に書いてあって、純化も込みで１ステップの手続きを書くと以下のようになります：

```python
def bbs_recurrence_relation_r_rev(
        a: Register, b: Register, c: Register) -> None:
    u'''bbs_recurrence_relation_r の逆操作
    '''
    with Lock(a):
        CX(a, b)
        lshift(c)
        CX(c, b)
        CCX(a, b, c)

def box_and_ball_step_r(x: Register, x2: Register,
                        cycle: int = 20) -> None:
    u'''右方向に動く箱玉系のステップ
    x  : 時刻 t での箱玉系の状態（入力）
    x2 : 時刻 t+1 での箱玉系の状態（出力）
    '''
    with Lock(x):
        a = x
        b = Register('B')  # allocate
        CX(x, b)
        ts = []
        for i in range(cycle):
            t = Register('temp{0}'.format(i))  # allocate
            bbs_recurrence_relation_r(a, b, t)
            (t, a, b) = (a, b, t)  # rename
            ts.append(t)
        # 時刻 t+1 での x の値を x2 に得る
        CX(x, x2)
        CX(a, x2)
        while ts:
            t = ts[-1]
            ts = ts[:-1]
            (a, b, t) = (t, a, b)  # rename
            bbs_recurrence_relation_r_rev(a, b, t)
            t.deallocate()
        CX(x, b)
        b.deallocate()
```

![q_bbm2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/412308/0a8ab38e-34f6-7ab2-7d72-3ea84e3b0533.png)


## 履歴の消去

上記で実装した `box_and_ball_step_r` は、現在の箱玉系の状態 $X(t)$ を元に $X(t+1)$ を得るものでした。これを単純に繰り返し実行すると、初期状態から各ステップの箱玉系の状態全ての履歴を得ることになります。
量子計算のようにビットリソースが限られた環境でこれをそのまま実現するのは厳しいので、前回ステップの情報を消去して最新の情報だけ残すことを考えます。

最新の情報だけにするには、$X(t) \rightarrow X(t+1)$ の各ステップ毎に元の $X(t)$ の情報を消去する必要があります。可逆計算なので、$X(t+1) \rightarrow X(t)$ という変換も定義可能なはずで、これを使えば以下の手順で履歴を消しながら処理を続けることができます。

1. 次のステップの $X$ を得る（$X(t) \rightarrow X(t+1)$）。
2. 得られた $X(t+1)$ を元に $X(t)$ を得る（$X(t+1) \rightarrow X(t)$）。
3. 元の $X(t)$ の情報に対して 2 で得られた $X(t)$ の情報を CNOT でぶつけて消す。
4. ぶつけるのに使った $X(t)$ の情報は $X(t) \rightarrow X(t+1)$ の逆変換で消去する。

箱玉系の場合、時間を遡るのと進行方向を逆にするのは同じです。

というわけで、`box_and_ball_step_r` の逆向き版 `box_and_ball_step_l` を作ると、履歴を消去する（ステップを繰り返しても必要なビット数が増え続けない）箱玉系シミュレーションが可能になります。

```python
# 1 ステップに付き 20 回漸化式を適用する
cycle=20
x = Register('X(t)')  # allocate
x2 = Register('X(t+1)')  # allocate
# 初期値
X(x, 31)
X(x, 30)
X(x, 29)
X(x, 28)
X(x, 18)
X(x, 17)
X(x, 12)
for n in range(100):
    # X(t) の状態を表示する
    print('{0:3}: {1}'.format(n, str(x)))
    # X(t) から X(t+1) を作る
    with Lock(x):  # x は入力として与えるので更新されない
        box_and_ball_step_r(x, x2, cycle=cycle)
    # X(t+1) から X(t) を作り直すことで、x の値を 0 に戻す
    with Lock(x2):  # x2 は入力として与えるので更新されない
        box_and_ball_step_l(x2, x, cycle=cycle)
    # 確認: x の値は 0 に戻っている
    x.must_be_zero()
    # 次のステップに備えて、x と x2 の値を入れ替え。
    # (x, x2) = (x2, x) でもよい。この場合は名前の付け替えになる。
    CX(x, x2)
    CX(x2, x)
    CX(x, x2)
    # 確認: 途中で確保した変数はすべて 0 に戻して返却している
    assert(REGCOUNT == 2)
```

実行例:

```
$ python box_and_ball_system_qq.py
  0: 11110000000001100001000000000000(X(t))
  1: 00001111000000011000100000000000(X(t))
  2: 00000000111100000110010000000000(X(t))
  3: 00000000000011110001101000000000(X(t))
  4: 00000000000000001110010111000000(X(t))
  5: 00000000000000000001101000111100(X(t))
  6: 11000000000000000000010110000011(X(t))
  7: 00111100000000000000001001100000(X(t))
  8: 00000011110000000000000100011000(X(t))
  9: 00000000001111000000000010000110(X(t))
 10: 10000000000000111100000001000001(X(t))
 11: 01100000000000000011110000100000(X(t))
 12: 00011000000000000000001111010000(X(t))
 13: 00000110000000000000000000101111(X(t))
 14: 11110001100000000000000000010000(X(t))
 15: 00001110011100000000000000001000(X(t))
 16: 00000001100011110000000000000100(X(t))
 17: 00000000011000001111000000000010(X(t))
 18: 00000000000110000000111100000001(X(t))
 19: 10000000000001100000000011110000(X(t))
 20: 01000000000000011000000000001111(X(t))
 21: 10111100000000000110000000000000(X(t))
 22: 01000011110000000001100000000000(X(t))
 23: 00100000001111000000011000000000(X(t))
 24: 00010000000000111100000110000000(X(t))
 25: 00001000000000000011110001100000(X(t))
 26: 00000100000000000000001110011100(X(t))
 27: 11000010000000000000000001100011(X(t))
 28: 00111101000000000000000000011000(X(t))
 29: 00000010111100000000000000000110(X(t))
 30: 10000001000011110000000000000001(X(t))
 (略)
```

古典可逆計算なので随時ビットの値は観測可能ということで、毎ステップ毎に箱玉系の状態を表示しています。

量子計算の場合には、各ステップの状態の観測をするたびに最初からやり直す形で実行します。

## まとめ

量子回路向けの箱玉漸化式について書かれた論文を見て、古典可逆計算を実装しました。
セル・オートマトンというと近隣セルの状態を元に次の状態を決めるルールを定義した表を作って、その表を引く、というやり方で実装するものだと思っていたのですが、漸化式を使うことで箱玉系のようなものを可逆（量子）回路として実装できるのは面白いと思いました。

## ソースコード

```python
# -*- coding:utf-8; mode:python -*-
from types import TracebackType
from typing import Optional, Tuple, Type

# 現在使用中のレジスタの数
REGCOUNT = 0

class Register(object):
    u'''ビット列をひとつのレジスタとして扱うクラス。
    '''
    NBITS = 32  # レジスタのビット幅は固定
    def __init__(self, title: str) -> None:
        u'''ビット列を確保する。
        全ビットの値は 0。
        '''
        self.title = title
        self.v : int = 0  # このレジスタが保持するビット列
        self.lock_count = 0  # ロック回数
        self.is_deallocated = False
        global REGCOUNT
        REGCOUNT += 1
    def deallocate(self) -> None:
        u'''当該ビット列を返却する。
        全ビットの値は 0 でなければならない。
        '''
        self.check()
        assert(self.v == 0)
        self.is_deallocated = True
        global REGCOUNT
        REGCOUNT -= 1
    def __str__(self) -> str:
        u'''ビット列の内容を表示する。
        '''
        return '{0}({1})'.format(
            'DEALLOCATED' if self.is_deallocated else
            ''.join(reversed([('0' if (self.v & (1 << i)) == 0 else '1')
                              for i in range(Register.NBITS)])),
            self.title)
    def must_be_zero(self) -> None:
        u'''検証用: 当該ビット列の全ビットの値が 0 であることを確認する。
        '''
        self.check()
        assert(self.v == 0)
    def set_locked(self, b: bool) -> None:
        u'''検証用: 引数に True を与えた場合、このレジスタの値の変更を禁止する。
        '''
        if b:
            self.lock_count += 1
        else:
            self.lock_count -= 1
    def is_locked(self) -> bool:
        u'''検証用: このレジスタの値の変更が禁止されているなら True を返す。
        '''
        return self.lock_count > 0
    def check(self) -> None:
        u'''検証用: 解放済みでないか確認する。
        '''
        if self.is_deallocated:
            raise RuntimeError('Register {0} already DEALOOCATED'.format(
                self.title))
    def check_not_locked(self) -> None:
        u'''検証用: 解放済みでなく、行進禁止でもないことを確認する。
        '''
        self.check()
        if self.is_locked():
            raise RuntimeError('Register {0} is locked'.format(self.title))
class Lock(object):
    u'''検証用: with 区間でレジスタを更新禁止状態にする。
    呼び出した関数でレジスタが破壊されていないことを確認するために用いる。
    '''
    def __init__(self, r: Register) -> None:
        u'''更新禁止対象のレジスタを引数に与える。
        '''
        self.r = r
    def __enter__(self) -> 'Lock':
        self.r.set_locked(True)
        return self
    def __exit__(self,
                 exc_type: Optional[Type[BaseException]],
                 exc_value: Optional[BaseException],
                 traceback: Optional[TracebackType]) -> None:
        self.r.set_locked(False)
def X(b: Register, index: int) -> None:
    u'''ビット列 b の index 番目の位置のビットを反転する。
    '''
    b.check_not_locked()
    b.v = b.v ^ (1 << index)
def CX(c: Register, b: Register) -> None:
    u'''ビット列 c を制御ビットとして、ビット列 b を反転する。
    '''
    c.check()
    b.check_not_locked()
    b.v = c.v ^ b.v
def CCX(c1: Register, c2: Register, b: Register) -> None:
    u'''ビット列 c1, c2 を制御ビットとして、ビット列 b を反転する。
    '''
    c1.check()
    c2.check()
    b.check_not_locked()
    b.v = (c1.v & c2.v) ^ b.v

def rshift(b: Register) -> None:
    u'''ビット列 b を右シフト（回転）する。
    '''
    b.check_not_locked()
    size = b.NBITS
    mask = (1 << size) - 1
    b.v = (((b.v << (size - 1)) & mask) |
           (b.v & mask) >> 1)

def lshift(b: Register) -> None:
    u'''ビット列 b を左シフト（回転）する。
    '''
    b.check_not_locked()
    size = b.NBITS
    mask = (1 << size) - 1
    b.v = (((b.v << 1) & mask) |
           (b.v & mask) >> (size - 1))

def bbs_recurrence_relation_r(a: Register, b: Register, c: Register) -> None:
    u'''右方向に動く箱玉系の漸化式
    '''
    with Lock(a):
        CCX(a, b, c)
        CX(c, b)
        rshift(c)
        CX(a, b)

def bbs_recurrence_relation_r_rev(
        a: Register, b: Register, c: Register) -> None:
    u'''bbs_recurrence_relation_r の逆操作
    '''
    with Lock(a):
        CX(a, b)
        lshift(c)
        CX(c, b)
        CCX(a, b, c)

def bbs_recurrence_relation_l(
        a: Register, b: Register, c: Register) -> None:
    u'''左方向に動く箱玉系の漸化式
    '''
    with Lock(a):
        CCX(a, b, c)
        CX(c, b)
        lshift(c)
        CX(a, b)

def bbs_recurrence_relation_l_rev(
        a: Register, b: Register, c: Register) -> None:
    u'''bbs_recurrence_relation_l の逆操作
    '''
    with Lock(a):
        CX(a, b)
        rshift(c)
        CX(c, b)
        CCX(a, b, c)

def box_and_ball_step_r(x: Register, x2: Register,
                        cycle: int = 20) -> None:
    u'''右方向に動く箱玉系のステップ
    x  : 時刻 t での箱玉系の状態（入力）
    x2 : 時刻 t+1 での箱玉系の状態（出力）
    '''
    with Lock(x):
        a = x
        b = Register('B')  # allocate
        CX(x, b)
        ts = []
        for i in range(cycle):
            t = Register('temp{0}'.format(i))  # allocate
            bbs_recurrence_relation_r(a, b, t)
            (t, a, b) = (a, b, t)  # rename
            ts.append(t)
        # 時刻 t+1 での x の値を x2 に得る
        CX(x, x2)
        CX(a, x2)
        while ts:
            t = ts[-1]
            ts = ts[:-1]
            (a, b, t) = (t, a, b)  # rename
            bbs_recurrence_relation_r_rev(a, b, t)
            t.deallocate()
        CX(x, b)
        b.deallocate()
def box_and_ball_step_l(x: Register, x2: Register,
                        cycle: int = 20) -> None:
    u'''左方向に動く箱玉系のステップ
    x  : 時刻 t での箱玉系の状態（入力）
    x2 : 時刻 t+1 での箱玉系の状態（出力）
    '''
    with Lock(x):
        a = x
        b = Register('B')  # allocate
        CX(x, b)
        ts = []
        for i in range(cycle):
            t = Register('temp{0}'.format(i))  # allocate
            bbs_recurrence_relation_l(a, b, t)
            (t, a, b) = (a, b, t)  # rename
            ts.append(t)
        # 時刻 t+1 での x の値を x2 に得る
        CX(x, x2)
        CX(a, x2)
        while ts:
            t = ts[-1]
            ts = ts[:-1]
            (a, b, t) = (t, a, b)  # rename
            bbs_recurrence_relation_l_rev(a, b, t)
            t.deallocate()
        CX(x, b)
        b.deallocate()

# 1 ステップに付き 20 回漸化式を適用する
cycle=20
x = Register('X(t)')  # allocate
x2 = Register('X(t+1)')  # allocate
# 初期値
X(x, 31)
X(x, 30)
X(x, 29)
X(x, 28)
X(x, 18)
X(x, 17)
X(x, 12)
for n in range(100):
    # X(t) の状態を表示する
    print('{0:3}: {1}'.format(n, str(x)))
    # X(t) から X(t+1) を作る
    with Lock(x):  # x は入力として与えるので更新されない
        box_and_ball_step_r(x, x2, cycle=cycle)
    # X(t+1) から X(t) を作り直すことで、x の値を 0 に戻す
    with Lock(x2):  # x2 は入力として与えるので更新されない
        box_and_ball_step_l(x2, x, cycle=cycle)
    # 確認: x の値は 0 に戻っている
    x.must_be_zero()
    # 次のステップに備えて、x と x2 の値を入れ替え。
    # (x, x2) = (x2, x) でもよい。この場合は名前の付け替えになる。
    CX(x, x2)
    CX(x2, x)
    CX(x, x2)
    # 確認: 途中で確保した変数はすべて 0 に戻して返却している
    assert(REGCOUNT == 2)
```