可逆論理で状態可逆分割セル・オートマトンシミュレーション

## 状態可逆分割セル・オートマトン

[箱玉系](https://qiita.com/tadashi9e/items/3f1a3adde974524b0d69)の次は二次元可逆 CA を実装してみます。

二次元の可逆 CA としてはビリヤードボールモデルが有名ですが、ここでは少し毛色の異なる CA [^cite] を実装してみます。

[^cite]: ナチュラルコンピューティングシリーズ第５巻可逆計算 第５章 状態可逆分割 CA モデル 1

![example.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/412308/344a6798-ff61-aa37-0578-96bf17ec4b83.gif)

この CA では、一つのセルは４つの部分からなり、４つの部分それぞれが２状態を持ちます。

図にするとこんな感じです。

╋━━━╋
┃　□　┃
┃□　□┃
┃　□　┃
╋━━━╋

４つの部分すべてが '0' の状態

╋━━━╋
┃　■　┃
┃■　■┃
┃　■　┃
╋━━━╋

４つの部分すべてが '1' の状態

セルの状態遷移規則は、隣接する４つのセルの隣接部分の状態に応じて以下のように変化するというものです：


　┃　□　┃
━╋━━━╋━　　　　╋━━━╋
　┃　　　┃　　　　　┃　□　┃
□┃　　　┃□　→　　┃□　□┃
　┃　　　┃　　　　　┃　□　┃
━╋━━━╋━　　　　╋━━━╋
　┃　□　┃
周りから■が来ていなければ□



　┃　■　┃
━╋━━━╋━　　　　╋━━━╋
　┃　　　┃　　　　　┃　□　┃
□┃　　　┃□　→　　┃□　□┃
　┃　　　┃　　　　　┃　■　┃
━╋━━━╋━　　　　╋━━━╋
　┃　□　┃
どちらかの方向から■が来ていれば、その方向に■を進める


　┃　■　┃
━╋━━━╋━　　　　╋━━━╋
　┃　　　┃　　　　　┃　■　┃
□┃　　　┃■　→　　┃□　■┃
　┃　　　┃　　　　　┃　□　┃
━╋━━━╋━　　　　╋━━━╋
　┃　□　┃
９０度の角度で２つの■が来ていれば、２つとも来ていた方向に打ち返す


　┃　■　┃
━╋━━━╋━　　　　╋━━━╋
　┃　　　┃　　　　　┃　□　┃
□┃　　　┃□　→　　┃■　■┃
　┃　　　┃　　　　　┃　□　┃
━╋━━━╋━　　　　╋━━━╋
　┃　■　┃
正面衝突する角度で２つの■が来ていれば、２つとも９０度向きを替えて打ち返す

　┃　■　┃
━╋━━━╋━　　　　╋━━━╋
　┃　　　┃　　　　　┃　■　┃
■┃　　　┃■　→　　┃■　■┃
　┃　　　┃　　　　　┃　□　┃
━╋━━━╋━　　　　╋━━━╋
　┃　□　┃
３つの■が来ていれば、３つとも来ていた方向に打ち返す。

　┃　■　┃
━╋━━━╋━　　　　╋━━━╋
　┃　　　┃　　　　　┃　■　┃
■┃　　　┃■　→　　┃■　■┃
　┃　　　┃　　　　　┃　■　┃
━╋━━━╋━　　　　╋━━━╋
　┃　■　┃
４つの■が来ていれば、４つとも来ていた方向に打ち返す。

ビリヤードボールモデルと同様に、これらの規則の CA 上で可逆論理を実装することも可能なようなのです。

ビリヤードボールモデルは、セル境界をステップ毎に切り替える形でルールを適用していきますが、状態可逆分割セル・オートマトンはセルを分割することでセル境界切り替えを不要としています。

## 実装

### 抽象化

最終的には、すべての操作はビットとその可逆操作で行うのですが、どうそれらを組み合わせれば CA を実現できるか考えるために、それらの要素の組み合わせを抽象化して考えることにします。

抽象度の高いものから低いものに向かって、あるいは全体から部分に向かって並べると以下のようなものです。

- 可逆並列計算機（可逆計算機を可逆通信網で接続したもの）
- 可逆計算機
- 可逆 CPU・可逆メモリ
- レジスタと可逆レジスタ間操作
- ビットとその可逆ゲート操作

一番下にあるのは一ビットずつあるビットとその可逆ゲート操作（ノットゲート、制御ノットゲート）で、その一つ上にあるのがビットを一つ以上並べて名前を付けたレジスタです。

レジスタの可逆操作であるノットゲート、制御ノットゲートもビットに対するそれと同様に定義します。

- レジスタのノットゲート: レジスタの全ビットを反転する。

- レジスタの制御ノットゲート：あるレジスタの i 番目のビットを制御ビットとして、別のレジスタ上の i 番目のビットに制御ノット処理する、ということを当該レジスタ上の全ビットについて行う。

- レジスタのスワップゲート：あるレジスタと別のレジスタとの間で全ビットの値を交換する。

これを組合せたものとして、可逆 CPU・可逆メモリを考えることができます。可逆 CPU・可逆メモリがあれば可逆計算機が構成可能で、可逆計算機を通信網で接続したものとして可逆並列計算機も考えることができます。
通信網も可逆ゲート操作で実現するので、上のリストでは可逆通信網と銘打っています。

### セル・オートマトン専用 SIMD 計算機

ここでは、CA を可逆並列計算機として実装することにします。前に作っていた[ビットシリアル SIMD シミュレータ](https://qiita.com/tadashi9e/items/32a97f845a8f74bb9caa)やその上で動かした[ライフゲーム](https://qiita.com/tadashi9e/items/0d696fea8364db29db6c)の可逆論理版といったところです（やっと伏線回収できた...）。

セルには４つの部分があってそれぞれ２状態を持つので、それに対応させて各 CPU には４つの 1bit レジスタ N / E / W/ S があり、東西南北方向に隣接する４つの CPU と直結した通信路で通信可能ということにしておきます。
各 CPU は通信を行うための受信バッファも必要です。東西南北から受信したデータを保持する必要があるので、N / E / W / S レジスタと同様に受信バッファも N / E / W / S の 1 bit 受信バッファ４つがあるものとします。
セル・オートマトンなので、全セル（＝CPU）は一斉に同じ処理を行えばよいので、全体として SIMD 型の並列計算機を構成します。

この並列計算機上では以下のようにセル・オートマトンのステップを進めることができます：

1. 通信フェーズ

  受信バッファの中身が空（=0）の状態で、以下を行います：
  a. 全 CPU は北方向のセル（CPU）に、自身の北部分の状態を送信します。
  b. 全 CPU は東方向のセル（CPU）に、自身の東部分の状態を送信します。
  c. 全 CPU は西方向のセル（CPU）に、自身の西部分の状態を送信します。
  d. 全 CPU は南方向のセル（CPU）に、自身の南部分の状態を送信します。

2. 決定フェーズ

  各 CPU は受信バッファに集まった隣接セル状態を元に、次ステップのセル状態を決定します。

3. 受信バッファクリアフェーズ

  各 CPU は受信バッファの中身を空（=0）に戻します。

4. 決定された状態の観測と表示

  決定されたセル状態を表示します。

### 通信フェーズ

各 CPU は東西南北方向の CPU と繋がっていて、それぞれの方向の CPU からの受信データを受け取る受信バッファも持っています。

SIMD なので、全 CPU は一斉に同じ方向にデータ送信を行います。全 CPU 一斉なのでデータ送信と受信は同時で不可分です。

次のフェーズで受信バッファ上のデータを元に、送信元だったレジスタの値を新しいステップのセル状態の値に更新したいので、0 クリアも同時にできれば好都合です。
なので、送信先バッファの状態は 0 の状態であることを前提として、送信を送信元と送信先のスワップ操作として定義します。

### 決定フェーズ

決定フェーズでは、受信バッファに集められた近隣セルの値を元に、レジスタ N / E / W / S の値を次ステップのセル状態の値に更新します。

例えば、

　┃　■　┃
━╋━━━╋━　　　　╋━━━╋
　┃　　　┃　　　　　┃　■　┃
■┃　　　┃■　→　　┃■　■┃
　┃　　　┃　　　　　┃　■　┃
━╋━━━╋━　　　　╋━━━╋
　┃　■　┃

は、４つの受信バッファ全ての値が 1 ならば、セル状態レジスタ N / E / W / S 全てを 1 にすることに対応します。
これは以下のような可逆回路で実現可能です。

![0000.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/412308/e755905a-3298-92f0-34a8-3b37fe99a755.png)

もうすこしだけ面白みのある例だと

　┃　■　┃
━╋━━━╋━　　　　╋━━━╋
　┃　　　┃　　　　　┃　□　┃
□┃　　　┃□　→　　┃■　■┃
　┃　　　┃　　　　　┃　□　┃
━╋━━━╋━　　　　╋━━━╋
　┃　■　┃

は以下のような可逆回路で実現できます：

![1010.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/412308/593bc9d1-584d-797c-4670-53fdf9706a06.png)

さらに一般化すると、受信バッファをマッチさせたいパターンに応じて各ビットを反転させたうえで４入力制御ノットゲートに与え、４入力制御ノットゲートが検知したらレジスタ N/E/W/S に値を設定する、という形になります：

![general.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/412308/15353217-9d0a-48ac-4010-0ed936cb9e0a.png)


受信バッファはパターンマッチに使うために反転で書き換えても元に戻すこと、受信バッファがマッチするパターンがたかだか一つしかないことに留意すれば、このパターンマッチを元にしたレジスタ N/E/W/S への出力設定は全パターンについて繰り返すことができます。

#### グレイコードを応用した決定フェーズの効率化

入力パターンは N/E/W/S の 4 bit をまとめて二進数で表せば簡単に全パターン列挙できます。

とはいえ、単純に入力パターン 0000 から 1111 まで順にパターンマッチを繰り返すのは若干効率がよくありません。

例えば、入力パターン 0001 を検査してから、次の入力パターン 0010 を検査する手順を以下のように書くことができます：

1. X(受信 N) X(受信 E) X(受信 W) で受信バッファ N/E/W を反転
2. ４入力制御ノットゲートで受信バッファをチェックしてマッチすればレジスタ N/E/W/S に出力設定
3. X(受信 N) X(受信 E) X(受信 W) で受信バッファ N/E/W を反転して元に戻す
4. X(受信 N) X(受信 E) X(受信 S) で受信バッファ N/E/S を反転
5. ４入力制御ノットゲートで受信バッファをチェックしてマッチすればレジスタ N/E/W/S に出力設定
6. X(受信 N) X(受信 E) X(受信 S) で受信バッファ N/E/S を反転

のように受信バッファの無駄な反転を繰り返すことになります（3 と 4 で受信バッファ N/E は反転２回で結局元に戻している）。
無駄な反転を省くプログラムを書くのは難しくありませんが、そうしたとしても入力パターン 0111 の次に 1000 を検査する、といった状況ではより多くのビット反転が必要になります。

そこで、

~~~
0000
0001
0010
0011
0100
0101
0110
0111
1000
1001
1010
1011
1100
1101
1110
1111
~~~

という一般的な二進数の並び順に検査する代わりに

~~~
0000
0001
0011
0010
0110
0111
0101
0100
1100
1101
1111
1110
1010
1011
1001
1000
~~~

という順に検査すれば、前後のパターンは 1bit しか違わないので検査の度に 1bit 反転するのみで済みます。
この順序はグレイコードと呼ばれるものです。グレイコードは前後で隣接する符号が 1bit しか違わない（ハミング距離 1 となる）ように作られた二進表現です。
上記グレイコードそのままだと最初にいきなり全部のビットを反転して 0000 から検査する必要があるので、論理反転した

~~~
1111
1110
1100
1101
1001
1000
1010
1011
0011
0010
0000
0001
0101
0100
0110
0111
~~~

を検査順序とします（但し、今回の CA ルールでは入力パターン 0000 に対しては出力パターンも 0000 なので検査不要です）。

### 受信バッファクリアフェーズ

次のステップに進む前に、受信バッファを 0 クリアしておかなければなりません。
幸いにも、今回の CA のルールは逆転しても同じなので

`受信バッファ検査→状態レジスタ設定`

で行ったことをそのまま

`状態レジスタ検査→受信バッファ設定`

という向きで適用するだけで受信バッファを 0 クリアできます。

### 決定された状態の観測と表示

量子計算だと観測するたびに初期状態設定からやり直す必要がありますが、古典可逆計算の場合にはそういう制約は不要です。

キャラクタ表示だと描画できるセルの数が少数に制限されてしまうため、pygame などでセルの状態を GUI 表示するのが便利です。上のアニメーション gif は pygame で作りました。

## 実行する可逆回路

32×32 の二次元空間についてプログラムが生成した可逆回路を掲載しておきます（但し、送受信部分は単純な繰り返しな上に分量が多くなるので割愛します）。

ここで `X`(ノットゲート), `CX`(制御ノットゲート), `SWAP`(交換), `CCCCX`(４入力制御ノットゲート) は可逆ゲートで、続く [] 内は可逆ゲート操作対象

- `N` : 状態レジスタ N
- `E` : 状態レジスタ E
- `W` : 状態レジスタ W
- `S` : 状態レジスタ S
- `recv_N`: 受信バッファ N
- `recv_E` : 受信バッファ E
- `recv_W` : 受信バッファ W
- `recv_S` : 受信バッファ S
- `flag` : フラグ

で、どのセルのレジスタなのかは `N:1` のようにコロンとインデックス（セル位置 x,y を一次元の値に変換したもの。具体的には x + y * 32）を表す数字で表記しています。

本来、全セル個別に

~~~
CCCCX[N:0,E:0,W:0,S:0,flag:0]
CCCCX[N:1,E:1,W:1,S:1,flag:1]
CCCCX[N:2,E:2,W:2,S:2,flag:2]
 ...
CCCCX[N:1023,E:1023,W:1023,S:1023,flag:1023]
~~~

などと書き並べるべきかもしれませんが、あまりに大量になるうえに無駄に読みづらいため

~~~
CCCCX[N,E,W,S,flag]
~~~

のように圧縮した表記としています。

```
# N 方向への送信（スワップ）
SWAP[N:0,recv_S:992]
SWAP[N:1,recv_S:993]
SWAP[N:2,recv_S:994]
SWAP[N:3,recv_S:995]
SWAP[N:4,recv_S:996]
SWAP[N:5,recv_S:997]
SWAP[N:6,recv_S:998]
SWAP[N:7,recv_S:999]
SWAP[N:8,recv_S:1000]
SWAP[N:9,recv_S:1001]
(略)
# E 方向への送信（スワップ）
SWAP[E:0,recv_W:1]
SWAP[E:1,recv_W:2]
SWAP[E:2,recv_W:3]
SWAP[E:3,recv_W:4]
SWAP[E:4,recv_W:5]
SWAP[E:5,recv_W:6]
SWAP[E:6,recv_W:7]
SWAP[E:7,recv_W:8]
SWAP[E:8,recv_W:9]
(略)
# W 方向への送信（スワップ）
SWAP[W:0,recv_E:31]
SWAP[W:1,recv_E:0]
SWAP[W:2,recv_E:1]
SWAP[W:3,recv_E:2]
SWAP[W:4,recv_E:3]
SWAP[W:5,recv_E:4]
SWAP[W:6,recv_E:5]
SWAP[W:7,recv_E:6]
SWAP[W:8,recv_E:7]
SWAP[W:9,recv_E:8]
(略)
# S 方向への送信（スワップ）
SWAP[S:0,recv_N:32]
SWAP[S:1,recv_N:33]
SWAP[S:2,recv_N:34]
SWAP[S:3,recv_N:35]
SWAP[S:4,recv_N:36]
SWAP[S:5,recv_N:37]
SWAP[S:6,recv_N:38]
SWAP[S:7,recv_N:39]
SWAP[S:8,recv_N:40]
SWAP[S:9,recv_N:41]
(略)
# --- 受信バッファパターンマッチによる次ステップ状態決定手続き ---
# Rule: 1111 -> 1111
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
CX[flag,E]
CX[flag,N]
CX[flag,W]
CX[flag,S]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
# Rule: 1110 -> 1110
X[recv_E]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
CX[flag,N]
CX[flag,W]
CX[flag,S]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
# Rule: 1100 -> 1100
X[recv_N]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
CX[flag,W]
CX[flag,S]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
# Rule: 1101 -> 1101
X[recv_E]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
CX[flag,E]
CX[flag,W]
CX[flag,S]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
# Rule: 1001 -> 1001
X[recv_W]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
CX[flag,E]
CX[flag,S]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
# Rule: 1000 -> 0010
X[recv_E]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
CX[flag,N]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
# Rule: 1010 -> 0101
X[recv_N]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
CX[flag,E]
CX[flag,W]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
# Rule: 1011 -> 1011
X[recv_E]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
CX[flag,E]
CX[flag,N]
CX[flag,S]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
# Rule: 0011 -> 0011
X[recv_S]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
CX[flag,E]
CX[flag,N]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
# Rule: 0010 -> 1000
X[recv_E]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
CX[flag,S]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
# Rule: 0001 -> 0100
X[recv_E]
X[recv_N]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
CX[flag,W]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
# Rule: 0101 -> 1010
X[recv_W]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
CX[flag,N]
CX[flag,S]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
# Rule: 0100 -> 0001
X[recv_E]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
CX[flag,E]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
# Rule: 0110 -> 0110
X[recv_N]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
CX[flag,N]
CX[flag,W]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
# Rule: 0111 -> 0111
X[recv_E]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
CX[flag,E]
CX[flag,N]
CX[flag,W]
CCCCX[recv_N,recv_E,recv_W,recv_S,flag]
X[recv_S]
# --- 次ステップ状態を元に受信バッファをクリアする手続き ---
# Rule: 1111 -> 1111
CCCCX[N,E,W,S,flag]
CX[flag,recv_E]
CX[flag,recv_N]
CX[flag,recv_W]
CX[flag,recv_S]
CCCCX[N,E,W,S,flag]
# Rule: 1110 -> 1110
X[E]
CCCCX[N,E,W,S,flag]
CX[flag,recv_N]
CX[flag,recv_W]
CX[flag,recv_S]
CCCCX[N,E,W,S,flag]
# Rule: 1100 -> 1100
X[N]
CCCCX[N,E,W,S,flag]
CX[flag,recv_W]
CX[flag,recv_S]
CCCCX[N,E,W,S,flag]
# Rule: 1101 -> 1101
X[E]
CCCCX[N,E,W,S,flag]
CX[flag,recv_E]
CX[flag,recv_W]
CX[flag,recv_S]
CCCCX[N,E,W,S,flag]
# Rule: 1001 -> 1001
X[W]
CCCCX[N,E,W,S,flag]
CX[flag,recv_E]
CX[flag,recv_S]
CCCCX[N,E,W,S,flag]
# Rule: 1000 -> 0010
X[E]
CCCCX[N,E,W,S,flag]
CX[flag,recv_N]
CCCCX[N,E,W,S,flag]
# Rule: 1010 -> 0101
X[N]
CCCCX[N,E,W,S,flag]
CX[flag,recv_E]
CX[flag,recv_W]
CCCCX[N,E,W,S,flag]
# Rule: 1011 -> 1011
X[E]
CCCCX[N,E,W,S,flag]
CX[flag,recv_E]
CX[flag,recv_N]
CX[flag,recv_S]
CCCCX[N,E,W,S,flag]
# Rule: 0011 -> 0011
X[S]
CCCCX[N,E,W,S,flag]
CX[flag,recv_E]
CX[flag,recv_N]
CCCCX[N,E,W,S,flag]
# Rule: 0010 -> 1000
X[E]
CCCCX[N,E,W,S,flag]
CX[flag,recv_S]
CCCCX[N,E,W,S,flag]
# Rule: 0001 -> 0100
X[E]
X[N]
CCCCX[N,E,W,S,flag]
CX[flag,recv_W]
CCCCX[N,E,W,S,flag]
# Rule: 0101 -> 1010
X[W]
CCCCX[N,E,W,S,flag]
CX[flag,recv_N]
CX[flag,recv_S]
CCCCX[N,E,W,S,flag]
# Rule: 0100 -> 0001
X[E]
CCCCX[N,E,W,S,flag]
CX[flag,recv_E]
CCCCX[N,E,W,S,flag]
# Rule: 0110 -> 0110
X[N]
CCCCX[N,E,W,S,flag]
CX[flag,recv_N]
CX[flag,recv_W]
CCCCX[N,E,W,S,flag]
# Rule: 0111 -> 0111
X[E]
CCCCX[N,E,W,S,flag]
CX[flag,recv_E]
CX[flag,recv_N]
CX[flag,recv_W]
CCCCX[N,E,W,S,flag]
X[S]
```

## ソース

`from reversible 〜` の部分は未公開（準備中）の可逆計算ライブラリなのでこのままで動かすことはできませんが、なんとなくやっていることはわかると思います。

~~~python
# -*- coding:utf-8; mode:python -*-
import pygame
from reversible import (CIRCUIT, RecordBlock)
from reversible.classical import (Bit, Register, ConserveBlock)
from reversible.classical.gate import (x, cx, ccccx, swap)

from types import TracebackType
from typing import List, Optional, Tuple, Type

WIDTH=32
HEIGHT=32
CELL_SIZE=16
GRID_COLOR=(0, 128, 128)
ALIVE_COLOR=(255, 255, 255)
DEAD_COLOR=(0, 0, 0)
FPS=100

Register.NBITS = WIDTH * HEIGHT

######################################################################
# Rules
######################################################################
        
class Rule(object):
    def __init__(self, from_pattern: int, to_pattern: int) -> None:
        self.from_pattern = from_pattern
        self.to_pattern = to_pattern
    def apply(self,
              current_pattern: int,
              n: Register, e: Register, w: Register, s: Register,
              n2: Register, e2: Register, w2: Register, s2: Register,
              flag: Register) -> int:
        CIRCUIT('# Rule: {0:04b} -> {1:04b}'.format(
            self.from_pattern, self.to_pattern))
        p = current_pattern ^ self.from_pattern
        if (p & 0b0001) != 0:
            Register.x(e)
        if (p & 0b0010) != 0:
            Register.x(n)
        if (p & 0b0100) != 0:
            Register.x(w)
        if (p & 0b1000) != 0:
            Register.x(s)
        with ConserveBlock(flag):
            Register.ccccx(n, e, w, s, flag)
            if (self.to_pattern & 0b0001) != 0:
                Register.cx(flag, e2)
            if (self.to_pattern & 0b0010) != 0:
                Register.cx(flag, n2)
            if (self.to_pattern & 0b0100) != 0:
                Register.cx(flag, w2)
            if (self.to_pattern & 0b1000) != 0:
                Register.cx(flag, s2)
            Register.ccccx(n, e, w, s, flag)
        flag.must_be_zero()
        return self.from_pattern
def clear_pattern(
        current_pattern: int,
        n: Register, e: Register, w: Register, s: Register) -> None:
    if (current_pattern & 0b0001) == 0:
        Register.x(e)
    if (current_pattern & 0b0010) == 0:
        Register.x(n)
    if (current_pattern & 0b0100) == 0:
        Register.x(w)
    if (current_pattern & 0b1000) == 0:
        Register.x(s)
def gray_to_int(i: int) -> int:
    u'''converts Gray code number to int number
    '''
    m = i
    while m != 0:
        m >>= 1
        i = i ^ m
    return i
def negative_4bit_gray_order(rules: List[Rule]) -> List[Rule]:
    u'''re-order rules
    '''
    result = sorted(
        rules,
        key = lambda rule: gray_to_int(0x0f & ~rule.from_pattern))
    print('\n'.join(['{0:04b}'.format(rule.from_pattern) for rule in result]))
    return result
        
# Rule table (ordered by nagative logic Gray code)
#   S, W, N, E  ->  S, W, N, E
RULES = negative_4bit_gray_order(
    [Rule(0b0001, 0b0100),
     Rule(0b0010, 0b1000),
     Rule(0b0011, 0b0011),
     Rule(0b0100, 0b0001),
     Rule(0b0101, 0b1010),
     Rule(0b0110, 0b0110),
     Rule(0b0111, 0b0111),
     Rule(0b1000, 0b0010),
     Rule(0b1001, 0b1001),
     Rule(0b1010, 0b0101),
     Rule(0b1011, 0b1011),
     Rule(0b1100, 0b1100),
     Rule(0b1101, 0b1101),
     Rule(0b1110, 0b1110),
     Rule(0b1111, 0b1111)])

def index_of_cell(x: int, y: int) -> int:
    if x < 0:
        x += WIDTH
    elif x > WIDTH - 1:
        x -= WIDTH
    if y < 0:
        y += HEIGHT
    elif y > HEIGHT - 1:
        y -= HEIGHT
    return x + y * WIDTH

class Cells(object):
    def __init__(self) -> None:
        self.reg_n = Register('N')
        self.reg_e = Register('E')
        self.reg_w = Register('W')
        self.reg_s = Register('S')
        self.recv_n = Register('recv_N')
        self.recv_e = Register('recv_E')
        self.recv_w = Register('recv_W')
        self.recv_s = Register('recv_S')
        self.flag = Register('flag')
    def send_n(self) -> None:
        CIRCUIT('# send_n')
        self.reg_n.check()
        self.recv_s.must_be_zero()
        for y in range(HEIGHT):
            for x in range(WIDTH):
                Bit.swap(self.reg_n.at(index_of_cell(x, y)),
                         self.recv_s.at(index_of_cell(x, y - 1)))
        self.reg_n.must_be_zero()
    def send_e(self) -> None:
        CIRCUIT('# send_e')
        self.reg_e.check()
        self.recv_w.must_be_zero()
        for y in range(HEIGHT):
            for x in range(WIDTH):
                Bit.swap(self.reg_e.at(index_of_cell(x, y)),
                         self.recv_w.at(index_of_cell(x + 1, y)))
        self.reg_e.must_be_zero()
    def send_w(self) -> None:
        CIRCUIT('# send_w')
        self.reg_w.check()
        self.recv_e.must_be_zero()
        for y in range(HEIGHT):
            for x in range(WIDTH):
                Bit.swap(self.reg_w.at(index_of_cell(x, y)),
                         self.recv_e.at(index_of_cell(x - 1, y)))
        self.reg_w.must_be_zero()
    def send_s(self) -> None:
        CIRCUIT('# send_s')
        self.reg_s.check()
        self.recv_n.must_be_zero()
        for y in range(HEIGHT):
            for x in range(WIDTH):
                Bit.swap(self.reg_s.at(index_of_cell(x, y)),
                         self.recv_n.at(index_of_cell(x, y + 1)))
        self.reg_s.must_be_zero()
    def next_state(self) -> None:
        with ConserveBlock(self.recv_n, self.recv_e, self.recv_w, self.recv_s):
            current_pattern = 0b1111
            for rule in RULES:
                current_pattern = rule.apply(
                    current_pattern,
                    self.recv_n, self.recv_e, self.recv_w, self.recv_s,
                    self.reg_n, self.reg_e, self.reg_w, self.reg_s,
                    self.flag)
            clear_pattern(current_pattern,
                  self.recv_n, self.recv_e, self.recv_w, self.recv_s)
    def purify_recv(self) -> None:
        with ConserveBlock(self.reg_n, self.reg_e, self.reg_w, self.reg_s):
            current_pattern = 0b1111
            for rule in RULES:
                current_pattern = rule.apply(
                    current_pattern,
                    self.reg_n, self.reg_e, self.reg_w, self.reg_s,
                    self.recv_n, self.recv_e, self.recv_w, self.recv_s,
                    self.flag)
            clear_pattern(current_pattern,
                  self.reg_n, self.reg_e, self.reg_w, self.reg_s)
        self.recv_n.must_be_zero()
        self.recv_e.must_be_zero()
        self.recv_w.must_be_zero()
        self.recv_s.must_be_zero()
def reversible_ca_step(cells: Cells) -> None:
    cells.send_n()
    cells.send_e()
    cells.send_w()
    cells.send_s()
    cells.next_state()
    cells.purify_recv()
######################################################################
# display CA
######################################################################
class Display(object):
    def __init__(self, cells: Cells) -> None:
        pygame.init()
        self.screen = pygame.display.set_mode((WIDTH * CELL_SIZE, HEIGHT * CELL_SIZE))
        pygame.display.set_caption('CA')
        self.clock = pygame.time.Clock()
        # display grid
        for y in range(1, HEIGHT):
            py = y * CELL_SIZE
            pygame.draw.line(self.screen, GRID_COLOR,
                             (0, py), (WIDTH * CELL_SIZE, py))
        for x in range(1, WIDTH):
            px = x * CELL_SIZE
            pygame.draw.line(self.screen, GRID_COLOR,
                             (px, 0), (px, HEIGHT * CELL_SIZE))
        # display initial CA states
        self.update(cells)
    def update(self, cells: Cells) -> None:
        u'''display CA status
        '''
        for y in range(WIDTH):
            for x in range(HEIGHT):
                px = x * CELL_SIZE
                py = y * CELL_SIZE
                cell_mask = (1 << index_of_cell(x, y))
                pygame.draw.rect(
                    self.screen,
                    (ALIVE_COLOR if (cells.reg_n.get() & cell_mask) != 0 else
                     DEAD_COLOR),
                    pygame.Rect(px + CELL_SIZE * 3 / 8,
                                py + CELL_SIZE     / 8,
                                CELL_SIZE / 4, CELL_SIZE / 4)) 
                pygame.draw.rect(
                    self.screen,
                    (ALIVE_COLOR if (cells.reg_e.get() & cell_mask) != 0 else
                     DEAD_COLOR),
                    pygame.Rect(px + CELL_SIZE * 5 / 8,
                                py + CELL_SIZE * 3 / 8,
                                CELL_SIZE / 4, CELL_SIZE / 4)) 
                pygame.draw.rect(
                    self.screen,
                    (ALIVE_COLOR if (cells.reg_w.get() & cell_mask) != 0 else
                     DEAD_COLOR),
                    pygame.Rect(px + CELL_SIZE     / 8,
                                py + CELL_SIZE * 3 / 8,
                                CELL_SIZE / 4, CELL_SIZE / 4)) 
                pygame.draw.rect(
                    self.screen,
                    (ALIVE_COLOR if (cells.reg_s.get() & cell_mask) != 0 else
                     DEAD_COLOR),
                    pygame.Rect(px + CELL_SIZE * 3 / 8,
                                py + CELL_SIZE * 5 / 8,
                                CELL_SIZE / 4, CELL_SIZE / 4))
        pygame.display.flip()
        self.clock.tick(FPS)

######################################################################
# initialize CA
######################################################################
def put_box(cells: Cells, x: int, y: int) -> None:
    Bit.x(cells.reg_e.at(index_of_cell(x, y)))
    Bit.x(cells.reg_s.at(index_of_cell(x, y)))
    Bit.x(cells.reg_w.at(index_of_cell(x + 1, y)))
    Bit.x(cells.reg_s.at(index_of_cell(x + 1, y)))
    Bit.x(cells.reg_n.at(index_of_cell(x, y + 1)))
    Bit.x(cells.reg_e.at(index_of_cell(x, y + 1)))
    Bit.x(cells.reg_n.at(index_of_cell(x + 1, y + 1)))
    Bit.x(cells.reg_w.at(index_of_cell(x + 1, y + 1)))
    
def initialize(cells: Cells) -> None:
    x0 = 16
    y0 = 16
    put_box(cells, x0 + 10, y0 + 10)
    put_box(cells, x0 + 9, y0 + 11)
    put_box(cells, x0 + 10, y0 - 10)
    put_box(cells, x0 + 9, y0 - 11)
    put_box(cells, x0 - 10, y0 - 10)
    put_box(cells, x0 - 9, y0 - 11)
    put_box(cells, x0 - 10, y0 + 10)
    put_box(cells, x0 - 9, y0 + 11)
    
    Bit.x(cells.reg_e.at(index_of_cell(x0 , y0 + 10)))
    Bit.x(cells.reg_e.at(index_of_cell(x0 - 2, y0 + 10)))

######################################################################
# main function
######################################################################
def main() -> None:
    cells = Cells()
    initialize(cells)
    display = Display(cells)
    first_time = True
    while True:
        with RecordBlock() as circuit:
            reversible_ca_step(cells)
            if first_time:
                circuit.show()
        first_time = False
        display.update(cells)
        assert(Register.allocated_count() == 9)

if __name__ == '__main__':
    main()
~~~

## まとめ

量子計算を一旦諦めて、古典可逆回路とすることで結構な規模の回路が書けるようになりました。
量子回路をいきなり書くとデバッグが大変ですが、古典可逆回路だとあちこちで不変チェックやら０チェックやらのアサートチェックを入れて確認できるのでデバッグもしやすいです（もしかしたら量子回路シミュレータでそういう機能をもつものがあるかも）。
将来量子計算機で信頼性のある大量の量子ビットを使い放題！という状況になってもなにかしらに使えそうな気分になってきました。古典回路を量子計算機で実行するメリットは一見なさげですが、それでも低消費電力という点で意義が見出される時代がやってくるかもしれません。