モジュロ加算器実装の試み

目的
----

[加減算](https://qiita.com/tadashi9e/items/44656ac35d97ba05ab9d)ができるようになったので、モジュロ加算器を実装しました。
ぱっとソースコードを見てどの程度処理内容を想像できるが重要だと思いますので、ここでは敢えてコメントは取り除いています。具体的な処理内容は後で説明します。

```python
class TestProc_MODULOADD(UO_Proc):
    def __init__(self) -> None:
        UO_Proc.__init__(self, 'MODULOADD')
    def __enter__(self) -> 'TestProc_MODULOADD':
        UO_Proc.__enter__(self)
        Integer.NBITS = 3
        with UO_Proc('MODULOADD.init'):
            a = Integer(0)
            a.hadamard(2)
            b = Integer(0)
            b.hadamard(2)
            n = Integer(4)
            n.hadamard(2)
            b0 = Integer(0)
            b0 << b
        with UO_Proc('MODULOADD.main'):
            b += a
            b -= n
            flag = QBit()
            QBit.cx(b.carry(), flag)
            n0 = Integer(0)
            Integer.cswap(flag, n, n0)
            b += n0
            Integer.cswap(flag, n, n0)
            n0.deallocate()
        with UO_Proc('MODULOADD_freeflag'):
            b -= a
            QBit.x(b.carry())
            QBit.cx(b.carry(), flag)
            QBit.x(b.carry())
            b += a
            flag.deallocate()
        self.a = a
        self.b0 = b0
        self.b = b
        self.n = n
        self.flag = flag
        return self

with TestProc_MODULOADD() as p:
    print('---- circuit ----')
    print(str(p))
    print('a =' + str(p.a))
    print('b0=' + str(p.b0))
    print('n =' + str(p.n))
    print('b =' + str(p.b))
    print('---- simulator ----')
    circuit = p.synthesis(blueqat.Circuit())
    print(str(circuit))
    print('---- result ----')
    for _ in range(100):
        result = str(circuit.m[:].run(shots=1))
        a = p.a.parse(result)
        b0 = p.b0.parse(result)
        n = p.n.parse(result)
        b = p.b.parse(result)
        if ((n == 0 and a + b0 == b) or
            ((a + b0) - b) % n == 0):
            print(str(result) + ' OK {0}+{1} mod {2} = {3}'.format(
                a, b0, n, b))
        else:
            print(str(result) + ' NG {0}+{1} mod {2} = {3}'.format(
                a, b0, n, b))
            sys.exit(1)
        if p.flag.parse(result) != 0:
            print('FLAG IS DIRTY')
```

高位合成っぽいことをしたいので、いかにも量子計算をやっていますという感じの見た目にならないことを目指しているのですが、量子ビット操作を顕に行わざるを得ない部分があって、そこは QBit クラスに関する部分として表現しています。

これを実行すると、以下のような結果が得られます：

```
---- circuit ----
MODULOADD{
 MODULOADD.init{
  Integer.hadamard{
   H[0],
   H[1]
  },
  Integer.hadamard{
   H[3],
   H[4]
  },
  Integer.init{
   X[8]
  },
  Integer.hadamard{
   H[6],
   H[7]
  },
  Integer.xor{
   CX[3,9],
   CX[4,10],
   CX[5,11]
  }
 },
 MODULOADD.main{
  Integer.add{
   Carry[12,0,3,13]{ ccx[0,3,13], cx[0,3], ccx[12,3,13] },
   Carry[13,1,4,14]{ ccx[1,4,14], cx[1,4], ccx[13,4,14] },
   Carry[14,2,5,15]{ ccx[2,5,15], cx[2,5], ccx[14,5,15] },
   CX[2,5],
   Sum[14,2,5]{ cx[2,5], cx[14,5] },
   RCarry[13,1,4,14]{ ccx[13,4,14], cx[1,4], ccx[1,4,14] },
   Sum[13,1,4]{ cx[1,4], cx[13,4] },
   RCarry[12,0,3,13]{ ccx[12,3,13], cx[0,3], ccx[0,3,13] },
   Sum[12,0,3]{ cx[0,3], cx[12,3] }
  },
  Integer.sub{
   RSum[12,6,3]{ cx[12,3], cx[6,3] },
   Carry[12,6,3,13]{ ccx[6,3,13], cx[6,3], ccx[12,3,13] },
   RSum[13,7,4]{ cx[13,4], cx[7,4] },
   Carry[13,7,4,14]{ ccx[7,4,14], cx[7,4], ccx[13,4,14] },
   RSum[14,8,5]{ cx[14,5], cx[8,5] },
   CX[8,5],
   RCarry[14,8,5,15]{ ccx[14,5,15], cx[8,5], ccx[8,5,15] },
   RCarry[13,7,4,14]{ ccx[13,4,14], cx[7,4], ccx[7,4,14] },
   RCarry[12,6,3,13]{ ccx[12,3,13], cx[6,3], ccx[6,3,13] }
  },
  CX[15,12],
  Integer.cswap{
   CSWAP[12,6,13],
   CSWAP[12,7,14],
   CSWAP[12,8,16]
  },
  Integer.add{
   Carry[17,13,3,18]{ ccx[13,3,18], cx[13,3], ccx[17,3,18] },
   Carry[18,14,4,19]{ ccx[14,4,19], cx[14,4], ccx[18,4,19] },
   Carry[19,16,5,15]{ ccx[16,5,15], cx[16,5], ccx[19,5,15] },
   CX[16,5],
   Sum[19,16,5]{ cx[16,5], cx[19,5] },
   RCarry[18,14,4,19]{ ccx[18,4,19], cx[14,4], ccx[14,4,19] },
   Sum[18,14,4]{ cx[14,4], cx[18,4] },
   RCarry[17,13,3,18]{ ccx[17,3,18], cx[13,3], ccx[13,3,18] },
   Sum[17,13,3]{ cx[13,3], cx[17,3] }
  },
  Integer.cswap{
   CSWAP[12,6,13],
   CSWAP[12,7,14],
   CSWAP[12,8,16]
  }
 },
 MODULOADD_freeflag{
  Integer.sub{
   RSum[13,0,3]{ cx[13,3], cx[0,3] },
   Carry[13,0,3,14]{ ccx[0,3,14], cx[0,3], ccx[13,3,14] },
   RSum[14,1,4]{ cx[14,4], cx[1,4] },
   Carry[14,1,4,16]{ ccx[1,4,16], cx[1,4], ccx[14,4,16] },
   RSum[16,2,5]{ cx[16,5], cx[2,5] },
   CX[2,5],
   RCarry[16,2,5,15]{ ccx[16,5,15], cx[2,5], ccx[2,5,15] },
   RCarry[14,1,4,16]{ ccx[14,4,16], cx[1,4], ccx[1,4,16] },
   RCarry[13,0,3,14]{ ccx[13,3,14], cx[0,3], ccx[0,3,14] }
  },
  X[15],
  CX[15,12],
  X[15],
  Integer.add{
   Carry[13,0,3,14]{ ccx[0,3,14], cx[0,3], ccx[13,3,14] },
   Carry[14,1,4,16]{ ccx[1,4,16], cx[1,4], ccx[14,4,16] },
   Carry[16,2,5,15]{ ccx[2,5,15], cx[2,5], ccx[16,5,15] },
   CX[2,5],
   Sum[16,2,5]{ cx[2,5], cx[16,5] },
   RCarry[14,1,4,16]{ ccx[14,4,16], cx[1,4], ccx[1,4,16] },
   Sum[14,1,4]{ cx[1,4], cx[14,4] },
   RCarry[13,0,3,14]{ ccx[13,3,14], cx[0,3], ccx[0,3,14] },
   Sum[13,0,3]{ cx[0,3], cx[13,3] }
  }
 }
}
a =Integer[0,1,2]
b0=Integer[9,10,11]
n =Integer[6,7,8]
b =Integer[3,4,5,15]
---- simulator ----
Circuit(20).h[0].h[1].h[3].h[4].x[8].h[6].h[7].cx[3, 9].cx[4, 10].cx[5, 11].ccx[0, 3, 13].cx[0, 3].ccx[12, 3, 13].ccx[1, 4, 14].cx[1, 4].ccx[13, 4, 14].ccx[2, 5, 15].cx[2, 5].ccx[14, 5, 15].cx[2, 5].cx[2, 5].cx[14, 5].ccx[13, 4, 14].cx[1, 4].ccx[1, 4, 14].cx[1, 4].cx[13, 4].ccx[12, 3, 13].cx[0, 3].ccx[0, 3, 13].cx[0, 3].cx[12, 3].cx[12, 3].cx[6, 3].ccx[6, 3, 13].cx[6, 3].ccx[12, 3, 13].cx[13, 4].cx[7, 4].ccx[7, 4, 14].cx[7, 4].ccx[13, 4, 14].cx[14, 5].cx[8, 5].cx[8, 5].ccx[14, 5, 15].cx[8, 5].ccx[8, 5, 15].ccx[13, 4, 14].cx[7, 4].ccx[7, 4, 14].ccx[12, 3, 13].cx[6, 3].ccx[6, 3, 13].cx[15, 12].ccx[12, 6, 13].ccx[12, 13, 6].ccx[12, 6, 13].ccx[12, 7, 14].ccx[12, 14, 7].ccx[12, 7, 14].ccx[12, 8, 16].ccx[12, 16, 8].ccx[12, 8, 16].ccx[13, 3, 18].cx[13, 3].ccx[17, 3, 18].ccx[14, 4, 19].cx[14, 4].ccx[18, 4, 19].ccx[16, 5, 15].cx[16, 5].ccx[19, 5, 15].cx[16, 5].cx[16, 5].cx[19, 5].ccx[18, 4, 19].cx[14, 4].ccx[14, 4, 19].cx[14, 4].cx[18, 4].ccx[17, 3, 18].cx[13, 3].ccx[13, 3, 18].cx[13, 3].cx[17, 3].ccx[12, 6, 13].ccx[12, 13, 6].ccx[12, 6, 13].ccx[12, 7, 14].ccx[12, 14, 7].ccx[12, 7, 14].ccx[12, 8, 16].ccx[12, 16, 8].ccx[12, 8, 16].cx[13, 3].cx[0, 3].ccx[0, 3, 14].cx[0, 3].ccx[13, 3, 14].cx[14, 4].cx[1, 4].ccx[1, 4, 16].cx[1, 4].ccx[14, 4, 16].cx[16, 5].cx[2, 5].cx[2, 5].ccx[16, 5, 15].cx[2, 5].ccx[2, 5, 15].ccx[14, 4, 16].cx[1, 4].ccx[1, 4, 16].ccx[13, 3, 14].cx[0, 3].ccx[0, 3, 14].x[15].cx[15, 12].x[15].ccx[0, 3, 14].cx[0, 3].ccx[13, 3, 14].ccx[1, 4, 16].cx[1, 4].ccx[14, 4, 16].ccx[2, 5, 15].cx[2, 5].ccx[16, 5, 15].cx[2, 5].cx[2, 5].cx[16, 5].ccx[14, 4, 16].cx[1, 4].ccx[1, 4, 16].cx[1, 4].cx[14, 4].ccx[13, 3, 14].cx[0, 3].ccx[0, 3, 14].cx[0, 3].cx[13, 3]
---- result ----
Counter({'10011001101000000000': 1}) OK 1+2 mod 6 = 3
Counter({'11010000101000000000': 1}) OK 3+2 mod 4 = 1
Counter({'00010010110000000000': 1}) OK 0+1 mod 5 = 1
Counter({'01001000100000000000': 1}) OK 2+0 mod 4 = 2
Counter({'00000000100000000000': 1}) OK 0+0 mod 4 = 0
Counter({'00011000111000000000': 1}) OK 0+3 mod 4 = 3
Counter({'10011001101000000000': 1}) OK 1+2 mod 6 = 3
Counter({'01011010110000000000': 1}) OK 2+1 mod 5 = 3
Counter({'01011000110000000000': 1}) OK 2+1 mod 4 = 3
Counter({'01001011100000000000': 1}) OK 2+0 mod 7 = 2
Counter({'00010011110000000000': 1}) OK 0+1 mod 7 = 1
Counter({'10011011101000000000': 1}) OK 1+2 mod 7 = 3
Counter({'10011000101000000000': 1}) OK 1+2 mod 4 = 3
Counter({'01001001100000000000': 1}) OK 2+0 mod 6 = 2
Counter({'10010000100000000000': 1}) OK 1+0 mod 4 = 1
(略)
```

前回との違い
------------

前回、加減算回路実装を書くときには、「一連の量子計算の手続きのブロック」を UO_Procedure インスタンスで明示的に作成して、明示的に演算操作に引数渡ししようとしていたのですが、'+=' のように引数で明示的に渡すことができない状況がありました。
このため、UO_Procedure インスタンスを渡したり渡さなかったりと一貫性がなくなってしまい分かりづらくなっていました。

今回の実装では、「現在の手続きの流れ」をグローバル変数に置いて、ユニタリ演算は常にそこに追加していく形にし、引数渡しを不要にしました。
新しい手続きの流れを作ったり、あるいは現在の手続きをコミットして直前の手続きの流れに戻ったりする操作は with ステートメントを使って自動化しています。

モジュロ加算器の実装
-------------------

モジュロ加算器は $A + B \bmod N$ を計算する手続きです。$\bmod N$ を計算するには 対象の値が負にならない程度に $N$ の減算を繰り返せばよいのですが、量子計算でのモジュロ加算器は

$A + B$ が $N$ 未満なら $A + B$ を返す
$A + B$ が $N$ 以上なら $A + B - N$ を返す

という手続きのようです（このあたり、理解が間違っているかも）。

このような雑な理解のもとでモジュロ加算器を記述したのが前掲のソースです。

```python
class TestProc_MODULOADD(UO_Proc):
    u'''モジュロ加算器のテスト実装
    '''
    def __init__(self) -> None:
        UO_Proc.__init__(self, 'MODULOADD')
```

まず、回路ブロック宣言部分です。名前 'MODULOADD' を付けて宣言しています。この名前は量子回路のダンプのときに用いられるだけなので適当なもので構いません。

```python
    def __enter__(self) -> 'TestProc_MODULOADD':
        UO_Proc.__enter__(self)
        Integer.NBITS = 3
```

with ステートメントでインスタンス化されたときに __enter__ が呼ばれて回路の合成の開始です。

ここで、一つの整数は（キャリービットを除いて）３量子ビットで表現することにします。４量子ビットでは現実的的な時間内でシミュレーションできませんでした。

```python
        with UO_Proc('MODULOADD.init'):
            a = Integer(0)
            a.hadamard(2)
            b = Integer(0)
            b.hadamard(2)
            n = Integer(4)
            n.hadamard(2)
            b0 = Integer(0)
            b0 << b
```

まずは、３つの整数 $A$, $B$, $N$ を変数 a, b, c に用意する処理です。これらの値に対して、$A + B \bmod N$ を得る量子回路を記述することがクラス全体のゴールです。

$A$, $B$, $N$ はいずれも正の整数で、アルゴリズムの都合上 $A<N, B<N$ とします。

ここではテストなのでいろいろな数で試してみたいので、アダマールゲートを用いてランダムに値を生成することにしました。
a, b は下位 2 ビットにアダマールゲートを作用させることで 3 以下の整数としています。
n は初期値 4 としたうえで下位 2 ビットにアダマールゲートを作用させることで 4〜7 の整数としています。

変数 b は計算途中の値、計算結果によって上書きされてしまうので、元の値 $B$ を変数 b0 にコピーしています。ここで `<<` はシフト演算でなく XOR (C での `^=`) の表現として用いています。

```python
        with UO_Proc('MODULOADD.main'):
            # 1
            b += a
            # 2
            b -= n
            # 3-1
            flag = QBit()
            # 3-2
            QBit.cx(b.carry(), flag)
            # 3-3
            n0 = Integer(0)
            # 3-4
            Integer.cswap(flag, n, n0)
            b += n0
            # 4
            Integer.cswap(flag, n, n0)
            n0.deallocate()
```

モジュロ加算器の核となる部分です。

1. まずは $A+B$ を計算し、変数 b に格納します。

2. 投機的に b から $N$ を減算します。b の値は $A + B - N$ となります。

3. b が負になった場合には $N$ の減算を取り消して（$N$ を加算して）、b の値を $A+B$ に戻します。b が正になった場合にはそのままにします。

古典計算機だと `if b < 0: b += n` と書けば済む話なのですが、量子計算機の場合には制御ゲートを用いて記述する必要があります。

具体的には、b が負になっている場合に 1 となる b のキャリービットを制御ビットとして、キャリービットが立っていたら b に $N$ を加える、ということをします。

3-1. 制御用フラグビットを別に用意しておく。
3-2. b のキャリービットの値を制御用フラグビットにコピーしておく。
3-3 $N$ あるいは $0$ を格納する変数を用意する。この変数の初期値は $0$。
3-4 制御スワップゲートを用いて、制御用フラグビットが 1 なら（b の値が負だったなら） n0 の値を $N$ に置き換える。そうでないなら n0 の値は $0$ のまま。
3-5 b に n0 を加える。

ということをします。

4 わざわざキャリービットを直接見ないで別に制御用フラグビットを持ってきて用いるのは、n0 の値を $0$ に戻して量子ビットを再利用可能にするためです。b += n0 が終わるとキャリービットは 0 に変化してしまうので、n0 に $N$ を格納したことを別のビットに記録しておかなければ元に戻すことができません。制御スワップで再び n と n0 を交換して、n の値が $N$、n0 の値が $0$ の状態に戻します。
deallocate で $0$ に戻った n0 を構成する量子ビットを再利用可能とします。

ここまでで $A + B \bmod N$ の計算は終わったのですが、制御用フラグビットが汚れたままです。汚れたままだと再利用できないので、これも 0 に戻します。

```python
        with UO_Proc('MODULOADD_freeflag'):
            b -= a
            QBit.x(b.carry())
            QBit.cx(b.carry(), flag)
            QBit.x(b.carry())
            # b から減算した a を加算しなおして、b の値を元に戻す。
            b += a
            flag.deallocate()
```

b の値は $A + B$ と $N$ の大小関係に従って２つのパターンがあります。

1. b = $A + B - N$ の場合 ($A + B \ge N$ だった、即ち制御フラグビットが 0 だった場合)
2. b = $A + B$ の場合 ($A + B < N$ だった、即ち制御フラグビットが 1 だった場合)

ここで b から $A$ を引くと上記の場合それぞれについて b の正負は以下のように分かれます：

- b = $B - N$ となる場合 ($A + B \ge N$ だった、即ち制御フラグビットが 0 だった場合) → $B < N$ なので b < 0
- b = $B$ となる場合 ($A + B < N$ だった、即ち制御フラグビットが 1 だった場合) → b > 0

つまり、b - $A$ が正だったら制御フラグを反転する、という方法で制御フラグを 0 に戻せることになります。制御フラグを 0 に戻したら b から $A$ を減算していたのを取り消します（$A$ を加算する）。

ここまででモジュロ加算器の記述は終わりです。

```python
        self.a = a
        self.b0 = b0
        self.b = b
        self.n = n
        self.flag = flag
        return self
```

残りの self に格納している部分は、どの量子ビットをどの変数に見立てたかを記録して、シミュレーション実行結果を表示可能にするためのものです。

まとめ
------

前回の実装に比べると with ステートメントを使ったことで見やすくなったと思います。

しかし、モジュロ加算器を書いてみると「可逆回路として実装するには量子ビットを意識してアルゴリズムを考えなければならない」という現実が見えてしまったため、加減算を手軽に書けるようになった、という以上の意義が見いだせていません。

任意のプログラムを可逆回路に変換できるものすごい高位合成ができるようになったらなったで、Python で実行される部分と量子回路になる部分との見分けが難しくなり可読性が却って下がるかもしれません。
自由にオペレータや構文規則を追加できればあまりそういう心配をしなくてよいのですが、Python の枠内では難しそうです。

ソースコード
------------

```python
import blueqat
import sys

from typing import List, Optional, Set

# ----------------------------------------------------------------------
class QubitAllocator(object):
    u'''量子ビット管理
    '''
    def __init__(self) -> None:
        self.allocated: Set[int] = set()
    def reset(self) -> None:
        u'''割付を初期化する。系全体をリセットする場合に呼び出す。
        '''
        self.allocated = set()
    def allocate1(self) -> int:
        u'''量子ビット一つを割り付け、そのインデックスを返す。
        '''
        if not self.allocated:
            self.allocated.add(0)
            return 0
        for i in range(0, max(self.allocated)):
            if i not in self.allocated:
                self.allocated.add(i)
                return i
        i = max(self.allocated) + 1
        self.allocated.add(i)
        return i
    def deallocate1(self, index: int) -> None:
        u'''量子ビット一つの割り付けを解除し、再割り当て可能にする。
        |0> に初期化された状態で返却されることを期待している。
        '''
        self.allocated.remove(index)

ALLOCATOR = QubitAllocator()

def reset_all() -> None:
    ALLOCATOR.reset()
# ----------------------------------------------------------------------
# 量子ビット操作
# ----------------------------------------------------------------------
class UnitaryOperation(object):
    u'''ユニタリ操作
    '''
    def __init__(self) -> None:
        pass
    def reverse(self) -> None:
        u'''登録されたユニタリ操作を逆順にする
        '''
        pass
    def synthesis(self, bq_circuit : blueqat.Circuit) -> blueqat.Circuit:
        u'''量子回路シミュレータに回路を書き込む
        '''
        return bq_circuit
    def __str__(self) -> str:
        return self.string(0)
    def string(self, depth : int) -> str:
        return ''
class UO_Procedure(UnitaryOperation):
    u'''ユニタリ操作リスト。
    複数ユニタリ操作を一連のユニタリ操作としてまとめたもの。
    '''
    def __init__(self, title: str) -> None:
        UnitaryOperation.__init__(self)
        self.title = title
        self.ops : List[UnitaryOperation] = []
    def append_op(self, op : UnitaryOperation) -> None:
        u'''ユニタリ操作をリストに追加登録する
        '''
        self.ops.append(op)
    def reverse(self) -> None:
        u'''登録されたユニタリ操作をエルミート共役（逆順）にする
        '''
        for op in reversed(self.ops):
            op.reverse()
        self.ops.reverse()
    def synthesis(self, bq_circuit : blueqat.Circuit) -> blueqat.Circuit:
        u'''量子回路シミュレータに回路を書き込む
        '''
        for op in self.ops:
            op.synthesis(bq_circuit)
        return bq_circuit
    def string(self, depth: int) -> str:
        s = self.title + '{\n'
        s += ',\n'.join(
            [' ' * (depth + 1) + op.string(depth + 1) for op in self.ops])
        s += '\n'
        s += ' ' * depth + '}'
        return s
class UO_H(UnitaryOperation):
    u'''アダマール演算
    '''
    def __init__(self, q : int) -> None:
        UnitaryOperation.__init__(self)
        self.q = q
        append_uo(self)
    def reverse(self) -> None:
        u'''反転しても同じ
        '''
        pass
    def synthesis(self, bq_circuit : blueqat.Circuit) -> blueqat.Circuit:
        bq_circuit.h[self.q]
        return bq_circuit
    def string(self, depth: int) -> str:
        return 'H[{0}]'.format(str(self.q))
class UO_X(UnitaryOperation):
    u'''パウリ X ゲート
    '''
    def __init__(self, q: int) -> None:
        UnitaryOperation.__init__(self)
        self.q = q
        append_uo(self)
    def reverse(self) -> None:
        u'''反転しても同じ
        '''
        pass
    def synthesis(self, bq_circuit : blueqat.Circuit) -> blueqat.Circuit:
        bq_circuit.x[self.q]
        return bq_circuit
    def string(self, depth: int) -> str:
        return 'X[{0}]'.format(str(self.q))
class UO_CX(UnitaryOperation):
    u'''制御 NOT ゲート
    '''
    def __init__(self, c: int, x: int) -> None:
        UnitaryOperation.__init__(self)
        self.c = c
        self.x = x
        append_uo(self)
    def reverse(self) -> None:
        u'''反転しても同じ
        '''
        pass
    def synthesis(self, bq_circuit: blueqat.Circuit) -> blueqat.Circuit:
        bq_circuit.cx[self.c, self.x]
        return bq_circuit
    def string(self, depth: int) -> str:
        return 'CX[{0},{1}]'.format(str(self.c), str(self.x))
class UO_CSWAP(UnitaryOperation):
    u'''制御 SWAP ゲート
    '''
    def __init__(self, c: int, a: int, b: int) -> None:
        UnitaryOperation.__init__(self)
        self.c = c
        self.a = a
        self.b = b
        append_uo(self)
    def reverse(self) -> None:
        u'''反転しても同じ
        '''
        pass
    def synthesis(self, bq_circuit: blueqat.Circuit) -> blueqat.Circuit:
        (c, a, b) = (self.c, self.a, self.b)
        bq_circuit.ccx[c, a, b].ccx[c, b, a].ccx[c, a, b]
    def string(self, depth: int) -> str:
        return 'CSWAP[{0},{1},{2}]'.format(
            str(self.c), str(self.a), str(self.b))
class QBit(object):
    u'''量子ビット操作を表現するクラス。
    '''
    def __init__(self, index: Optional[int] = None) -> None:
        u'''与えられた量子ビット位置の操作インスタンスを作成する。
        位置が与えられなかった場合、新しい量子ビット位置を確保する。
        '''
        self.index : int = (index if index is not None else
                            ALLOCATOR.allocate1())
    def deallocate(self) -> int:
        u'''このインスタンスが指す量子ビットを返却する。
        '''
        ALLOCATOR.deallocate1(self.index)
        return self.index
    @staticmethod
    def h(q: 'QBit') -> None:
        UO_H(q.index)
    @staticmethod
    def x(q: 'QBit') -> None:
        UO_X(q.index)
    @staticmethod
    def cx(c: 'QBit', q: 'QBit') -> None:
        UO_CX(c.index, q.index)
    @staticmethod
    def cswap(c: 'QBit', a: 'QBit', b: 'QBit') -> None:
        UO_CSWAP(c.index, a.index, b.index)
    def parse(self, s: str) -> int:
        u'''blueqat の結果出力を元に、当該量子ビットが 0 か 1 かを返す。
        '''
        i = s.find("'")
        s = s[i + 1:]
        return 0 if s[self.index] == '0' else 1
    def __str__(self) -> str:
        return str(self.index)
# ----------------------------------------------------------------------
# ユニタリ操作のグループ化
# ----------------------------------------------------------------------
u'''
ユニタリ操作リストスタックトップ。
ユニタリ操作は暗黙のうちにスタックトップにある操作リストに追加される
'''
PROCEDURE_STACK_TOP : Optional['UO_Proc'] = None
class UO_Proc(object):
    u'''ユニタリ操作リストスタック
    '''
    def __init__(self, title: str) -> None:
        self._procedure = UO_Procedure(title)
    def __enter__(self) -> 'UO_Proc':
        global PROCEDURE_STACK_TOP
        self._parent = PROCEDURE_STACK_TOP
        PROCEDURE_STACK_TOP = self
        return self
    def __exit__(self, type, value, traceback) -> None:
        global PROCEDURE_STACK_TOP
        PROCEDURE_STACK_TOP = self._parent
        if PROCEDURE_STACK_TOP:
            PROCEDURE_STACK_TOP.get_procedure().append_op(
                self._procedure)
    def __str__(self) -> str:
        return str(self._procedure)
    def get_procedure(self) -> 'UO_Procedure':
        return self._procedure
    def synthesis(self, bq_circuit: blueqat.Circuit) -> blueqat.Circuit:
        return self._procedure.synthesis(bq_circuit)
def dump_uo_stack() -> None:
    if PROCEDURE_STACK_TOP is None:
        raise RuntimeError('procedure not started')
    p : Optional[UO_Proc] = PROCEDURE_STACK_TOP
    while p is not None:
        print(str(p.get_procedure()))
        p = p._parent
def append_uo(uo : 'UnitaryOperation') -> None:
    u'''ユニタリ操作リストスタック先頭のユニタリ操作リストに追加する。
    '''
    if PROCEDURE_STACK_TOP is None:
        raise RuntimeError('procedure not started')
    PROCEDURE_STACK_TOP.get_procedure().append_op(uo)
def reverse_current_proc() -> None:
    u'''ユニタリ操作リストスタック先頭のユニタリ操作リストを
    エルミート共役に変換する。
    '''
    if PROCEDURE_STACK_TOP is None:
        raise RuntimeError('procedure not started')
    PROCEDURE_STACK_TOP.get_procedure().reverse()
# ----------------------------------------------------------------------
# 整数操作
# ----------------------------------------------------------------------
class Carry(UnitaryOperation):
    u'''キャリー伝搬付き加算
    '''
    def __init__(self,
                 a: QBit, b: QBit, c: QBit, d: QBit,
                 is_reversed: bool=False) -> None:
        UnitaryOperation.__init__(self)
        self.a = a
        self.b = b
        self.c = c
        self.d = d
        append_uo(self)
        self.is_reversed = is_reversed
    def reverse(self) -> None:
        u'''逆演算
        '''
        self.is_reversed = not self.is_reversed
    def synthesis(self, bq_circuit: blueqat.Circuit) -> blueqat.Circuit:
        u'''量子計算シミュレータに回路を書き込む。
        '''
        (a, b, c, d) = (self.a.index, self.b.index, self.c.index, self.d.index)
        if not self.is_reversed:
            bq_circuit.ccx[b, c, d].cx[b, c].ccx[a, c, d]
        else:
            bq_circuit.ccx[a, c, d].cx[b, c].ccx[b, c, d]
        return bq_circuit
    def string(self, depth: int) -> str:
        return (
            'Carry[{0},{1},{2},{3}]{{ ccx[{1},{2},{3}], cx[{1},{2}], ccx[{0},{2},{3}] }}'.format(str(self.a), str(self.b), str(self.c), str(self.d))
            if not self.is_reversed else
            'RCarry[{0},{1},{2},{3}]{{ ccx[{0},{2},{3}], cx[{1},{2}], ccx[{1},{2},{3}] }}'.format(str(self.a), str(self.b), str(self.c), str(self.d)))
class Sum(UnitaryOperation):
    def __init__(self,
                 a: QBit, b: QBit, c: QBit,
                 is_reversed: bool = False) -> None:
        UnitaryOperation.__init__(self)
        self.a = a
        self.b = b
        self.c = c
        append_uo(self)
        self.is_reversed = is_reversed
    def reverse(self) -> None:
        u'''逆演算
        '''
        self.is_reversed = not self.is_reversed
    def synthesis(self, bq_circuit: blueqat.Circuit) -> blueqat.Circuit:
        (a, b, c) = (self.a.index, self.b.index, self.c.index)
        if not self.is_reversed:
            bq_circuit.cx[b, c].cx[a, c]
        else:
            bq_circuit.cx[a, c].cx[b, c]
        return bq_circuit
    def string(self, depth: int) -> str:
        return (
            'Sum[{0},{1},{2}]{{ cx[{1},{2}], cx[{0},{2}] }}'.format(
                str(self.a), str(self.b), str(self.c))
            if not self.is_reversed else
            'RSum[{0},{1},{2}]{{ cx[{0},{2}], cx[{1},{2}] }}'.format(
                str(self.a), str(self.b), str(self.c)))
class Integer(object):
    NBITS = 4
    def __init__(self, n: int) -> None:
        nbits = Integer.NBITS
        self.qbits = [QBit() for _ in range(nbits)]
        # 補助ビットは必要になるまで確保を遅延する
        self.cs : Optional[List[QBit]] = None
        # キャリービットも必要になるまで確保を遅延する
        self._carry : Optional[QBit] = None
        if n != 0:
            # 初期値を与える
            with UO_Proc('Integer.init') as p:
                for i in range(0, nbits):
                    if n & (1 << i) == 0:
                        continue
                    QBit.x(self.qbits[i])
    def deallocate(self) -> None:
        u'''Integer が保持している整数値が 0 であるとわかっている場合にのみ
        呼び出すこと。
        '''
        [c.deallocate() for c in self.qbits]
    def _deallocate(self) -> None:
        u'''補助ビットは加算・減算が終わったら即座に返却する
        （補助ビットは加減算が終わった時点で |0> になっている）。
        '''
        if self.cs:
            [c.deallocate() for c in self.cs]
            self.cs = None
    def _cs(self, i: int) -> QBit:
        u'''補助ビット・キャリービットの獲得
        '''
        if i >= len(self.qbits):
            return self.carry()
        if not self.cs:
            self.cs = [QBit() for _ in range(len(self.qbits))]
        return self.cs[i]
    def carry(self) -> QBit:
        u'''キャリービットの獲得
        '''
        if self._carry is None:
            self._carry = QBit()
        return self._carry
    def __iadd__(self, other: 'Integer') -> 'Integer':
        u'''引数に与えられた Integer を加算する。
        '''
        self._synthesis_iadd('Integer.add', other)
        self._deallocate()
        return self
    def __isub__(self, other: 'Integer') -> 'Integer':
        u'''引数に与えられた Integer で減算する。加算の逆操作を行う。
        '''
        self._synthesis_iadd('Integer.sub', other,
                             is_reversed = True)
        self._deallocate()
        return self
    def _synthesis_iadd(self, title: str, other: 'Integer',
                        is_reversed : bool = False) -> None:
        u'''加算回路を合成する。
        '''
        if len(other.qbits) != len(self.qbits):
            raise RuntimeError('invalid Integer bit width')
        with UO_Proc(title):
            for i in range(len(self.qbits)):
                a = other.qbits[i]
                b = self.qbits[i]
                c = self._cs(i)
                c2 = self._cs(i + 1)
                Carry(c, a, b, c2)
            QBit.cx(other.qbits[-1], self.qbits[-1])
            for i in range(len(self.qbits) - 1, 0, -1):
                a = other.qbits[i]
                b = self.qbits[i]
                c = self._cs(i)
                a1 = other.qbits[i - 1]
                b1 = self.qbits[i - 1]
                c1 = self._cs(i - 1)
                Sum(c, a, b)
                Carry(c1, a1, b1, c, is_reversed=True)
            Sum(self._cs(0), other.qbits[0], self.qbits[0])
            if is_reversed:
                reverse_current_proc()
    def __lshift__(self, orig: 'Integer') -> 'Integer':
        u'''引数に与えられた Integer との XOR をとる。
        '''
        if len(orig.qbits) != len(self.qbits):
            raise RuntimeError('invalid Integer bit width')
        with UO_Proc('Integer.xor'):
            for c, x in zip(orig.qbits, self.qbits):
                QBit.cx(c, x)
            # orig にキャリービットがあれば、キャリービットも XOR をとる
            if orig._carry is not None:
                QBit.cx(orig._carry, self.carry())
        return self
    def hadamard(self, n: int = -1) -> 'Integer':
        u'''下位 n ビットだけアダマール演算を作用させる。
        '''
        with UO_Proc('Integer.hadamard'):
            if n < 0:
                n = len(self.qbits)
            for i in range(0, n):
                QBit.h(self.qbits[i])
        return self
    @staticmethod
    def cswap(c: QBit, a: 'Integer', b: 'Integer') -> None:
        u'''制御付きで、２つの整数 a と b を交換する。
        '''
        with UO_Proc('Integer.cswap'):
            if len(a.qbits) != len(b.qbits):
                raise RuntimeError('invalid Integer bit width')
            for a1, b1 in zip(a.qbits, b.qbits):
                QBit.cswap(c, a1, b1)
            if a._carry is not None:
                QBit.cswap(c, a._carry, b.carry())
            elif b._carry is not None:
                QBit.cswap(c, a.carry(), b._carry)
    def parse(self, s: str) -> int:
        u'''blueqat の結果出力を元に、当該量子ビット群が表す正整数値を返す。
        '''
        i = s.find("'")
        s = s[i + 1:]
        n = 0
        for i in range(0, len(self.qbits)):
            if s[self.qbits[i].index] == '1':
                n += (1 << i)
        if self._carry is not None and s[self._carry.index] == '1':
            n += (1 << len(self.qbits))
        return n
    def parse_signed(self, s: str) -> int:
        u'''blueqat の結果出力を元に、当該量子ビット群が表す整数値を返す。
        キャリービットが存在し、かつその値が |1> の場合は負数とみなす。
        '''
        n = self.parse(s)
        if n > (1 << len(self.qbits)):
            return n - (1 << (len(self.qbits) + 1))
        return n
    def __str__(self) -> str:
        args = ','.join([str(q) for q in self.qbits])
        if self._carry is not None:
            args += ',' + str(self._carry)
        return 'Integer[{0}]'.format(args)
# ----------------------------------------------------------------------
# 加算のテスト
# ----------------------------------------------------------------------
class TestProc_ADD(UO_Proc):
    def __init__(self) -> None:
        UO_Proc.__init__(self, 'TestProc_ADD')
    def __enter__(self) -> 'TestProc_ADD':
        UO_Proc.__enter__(self)
        # 整数は 4 量子ビットで表現する。
        Integer.NBITS = 4
        a = Integer(0)
        b = Integer(0)
        a.hadamard()
        b.hadamard()
        b0 = Integer(0)
        b0 << b
        b += a
        self.a = a
        self.b0 = b0
        self.b = b
        return self
with TestProc_ADD() as p_add:
    print('---- circuit ----')
    print(str(p_add))
    print('a =' + str(p_add.a))
    print('b0=' + str(p_add.b0))
    print('b =' + str(p_add.b))
    print('---- simulator ----')
    circuit = p_add.synthesis(blueqat.Circuit())
    print(str(circuit))
    print('---- result ----')
    for _ in range(10):
        result = str(circuit.m[:].run(shots=1))
        a = p_add.a.parse(result)
        b0 = p_add.b0.parse(result)
        b = p_add.b.parse(result)
        if a + b0 == b:
            print(str(result) + ' OK {0}+{1}={2}'.format(a, b0, b))
        else:
            print(str(result) + ' NG {0}+{1}={2}'.format(a, b0, b))
            sys.exit(1)
reset_all()

# ----------------------------------------------------------------------
# 減算のテスト
# ----------------------------------------------------------------------
class TestProc_SUB(UO_Proc):
    def __init__(self) -> None:
        UO_Proc.__init__(self, 'TestProc_SUB')
    def __enter__(self) -> 'TestProc_SUB':
        UO_Proc.__enter__(self)
        # 整数は 4 量子ビットで表現する。
        Integer.NBITS = 4
        a = Integer(0)
        b = Integer(0)
        a.hadamard()
        b.hadamard()
        b0 = Integer(0)
        b0 << b
        b -= a
        self.a = a
        self.b0 = b0
        self.b = b
        return self
with TestProc_SUB() as p_sub:
    print('---- circuit ----')
    print(str(p_sub))
    print('a =' + str(p_sub.a))
    print('b0=' + str(p_sub.b0))
    print('b =' + str(p_sub.b))
    print('---- simulator ----')
    circuit = p_sub.synthesis(blueqat.Circuit())
    print(str(circuit))
    print('---- result ----')
    for _ in range(10):
        result = str(circuit.m[:].run(shots=1))
        a = p_sub.a.parse(result)
        b0 = p_sub.b0.parse(result)
        b = p_sub.b.parse_signed(result)
        if b0 -a == b:
            print(str(result) + ' OK {0}-{1}={2}'.format(b0, a, b))
        else:
            print(str(result) + ' NG {0}-{1}={2}'.format(b0, a, b))
            sys.exit(1)
reset_all()
# ----------------------------------------------------------------------
# モジュロ加算のテスト
# ----------------------------------------------------------------------
class TestProc_MODULOADD(UO_Proc):
    u'''モジュロ加算器のテスト実装
    '''
    def __init__(self) -> None:
        UO_Proc.__init__(self, 'MODULOADD')
    def __enter__(self) -> 'TestProc_MODULOADD':
        u'''回路の作成
        '''
        UO_Proc.__enter__(self)
        # 整数は 3 量子ビットで表現する。
        Integer.NBITS = 3
        with UO_Proc('MODULOADD.init'):
            # 加減算に用いる２つの整数 a, b を用意する。
            # 下位 2 bit をランダムに決定する（3 以下の数になる）
            a = Integer(0)
            a.hadamard(2)
            b = Integer(0)
            b.hadamard(2)
            # n は 4 以上にしておく
            n = Integer(4)
            n.hadamard(2)
            # テスト結果確認のために、b の値を bo にコピーしておく。
            # ここで << はビットシフトでなく XOR を使った b0 への値コピー。
            b0 = Integer(0)
            b0 << b
        with UO_Proc('MODULOADD.main'):
            # 足し算 (A + B) を行う。
            b += a
            # A + B から投機的に N を減算してみる
            b -= n
            # A + B - N が負なら +N して A + B に戻す
            flag = QBit()
            QBit.cx(b.carry(), flag)
            n0 = Integer(0)
            Integer.cswap(flag, n, n0)
            b += n0
            # n0 の値を 0 に戻して解放
            Integer.cswap(flag, n, n0)
            n0.deallocate()
        with UO_Proc('MODULOADD_freeflag'):
            # flag も解放したい。
            # flag を |1> にしたかチェックして元に戻す
            #
            # flag が |0> の場合：
            # b = A + B - N になっている。
            # ここからから A を引くと、b = B - N なので負（B < N なので）。
            #
            # flag が |1> の場合：
            # b = A + B になっている。
            # ここからから A を引くと、b = B なので正。
            #
            # よって、b -= a が正の場合に flag を反転すれば flag は |0> に戻る。
            # b -= a が負の場合にキャリーは |1> なので、キャリーの値を反転した上で
            # flag に制御ノットを作用させる。
            # b の値が変わらないように、キャリーの値は再反転して元に戻す。
            b -= a
            QBit.x(b.carry())
            QBit.cx(b.carry(), flag)
            QBit.x(b.carry())
            # b から減算した a を加算しなおして、b の値を元に戻す。
            b += a
            flag.deallocate()
        # どの量子ビットを割り付けたか確認できるように self に記録しておく。
        self.a = a
        self.b0 = b0
        self.b = b
        self.n = n
        self.flag = flag
        return self
with TestProc_MODULOADD() as p:
    print('---- circuit ----')
    print(str(p))
    print('a =' + str(p.a))
    print('b0=' + str(p.b0))
    print('n =' + str(p.n))
    print('b =' + str(p.b))
    print('---- simulator ----')
    circuit = p.synthesis(blueqat.Circuit())
    print(str(circuit))
    print('---- result ----')
    for _ in range(100):
        result = str(circuit.m[:].run(shots=1))
        a = p.a.parse(result)
        b0 = p.b0.parse(result)
        n = p.n.parse(result)
        b = p.b.parse(result)
        if ((n == 0 and a + b0 == b) or
            ((a + b0) - b) % n == 0):
            print(str(result) + ' OK {0}+{1} mod {2} = {3}'.format(
                a, b0, n, b))
        else:
            print(str(result) + ' NG {0}+{1} mod {2} = {3}'.format(
                a, b0, n, b))
            sys.exit(1)
        if p.flag.parse(result) != 0:
            print('FLAG IS DIRTY')
```