乗算回路の高位合成の試み

目的
----

[加算器を `+=` で書けるようにした](https://qiita.com/tadashi9e/items/44656ac35d97ba05ab9d)ので、乗算も `*` で書けるようにしました。

アルゴリズム
----------

二進数での掛け算を筆算で書くとこんな感じです：

|  |  |  |  |  |
|:-:|:-:|:-:|:-:|:-:|
|  |  | 1| 1| 0|
|  |×| 1| 0| 1|
|  |─|─|─|─|

ここでは、上に書いた数を被乗算数、下に書いた数を乗算数と呼ぶことにします。

まず、乗算数の最下位ビットが 1 なので、$110 \times 1 = 110$

|  |  |  |  |  |
|:-:|:-:|:-:|:-:|:-:|
|  |  | 1| 1| 0|
|  |×| 1| 0| 1|
|  |─|─|─|─|
|  |  | 1| 1| ***0*** |
| :| :| :| :|  |
|─|─|─|─|─|
| ?| ?| ?| ?| ***0*** |

これで、掛け算結果の最下位ビットが $0$ に確定しました。

乗算数の次のビットは 0 なので、$110 \times 0 = 000$ なので $110$ の下に
$000$ を書いて両者の和を考えます。$110$ の最下位ビットは既に結果として確定しているので、それを右シフトで取り除いた $11$ と $000$ の和を考えます。和の最下位ビットが掛け算結果として確定します。

|  |  |  |  |  |
|:-:|:-:|:-:|:-:|:-:|
|  |  | 1| 1| 0|
|  |×| 1| 0| 1|
|  |─|─|─|─|
|  |  | 1| ***1*** | 0|
|  | 0| 0| ***0*** |  |
| :| :| :|  |  |
|─|─|─|─|─|
| ?| ?| ?| ***1*** | 0|

最後に残った乗算数のビット 1 について $110 \times 1 = 110$ なので、先程の $11 + 000 = 011$ を右シフトした値 $01$ に $110$ を足して $111$ を得ます。$111$ と先程「確定した」と書いた２つのビット $10$ とをつなげると、掛け算全体の結果 $11110$ が得られます。

|  |  |  |  |  |
|:-:|:-:|:-:|:-:|:-:|
|  |  | 1| 1| 0|
|  |×| 1| 0| 1|
|  |─|─|─|─|
|  |  | ***1***| 1| 0|
|  | ***0***| ***0***| 0|  |
| ***1***| ***1***| ***0***|  |  |
|─|─|─|─|─|
| ***1***| ***1***| ***1***| 1| 0|

このような手続きをそのまま量子回路として書けば動きそうな気がします:

1. 引数として整数 a と b が与えられます。処理を簡単にするために、両者は同じビット数で表現されているものとします。
2. ワークエリアとなる整数 t、a0 を値０で準備します。
3. 変数 b の各ビットについて、以下を行います:
 1. b の i 番目のビットが 1 なら t に b を加算します。0 なら t はそのままにします。
これは、具体的には以下の手続きでできます：
<dl>
  <dt>ステップ 1</dt><dd>b の i 番目のビットをコントロールビットとして a の値と a0 の値を交換。</dd>
  <dt>ステップ 2</dt><dd>`t += a0` を実行。</dd>
  <dt>ステップ 3</dt><dd>もう一度、b の i 番目のビットをコントロールビットとして a の値と a0 の値を交換</dd>
</dl>
 2. t の最下位ビットをメモしたうえで、最下位ビット以外のビットを元に新たに t を作り直します（右シフト）。このとき、最上位ビットには新しい量子ビットを割り付けます。
4. ここまでで乗算結果が得られました。a0 は値 $0$ なので量子ビットを解放します。3.2 でメモしておいて各最下位ビットの値のリストと、最終的な t の値のリストを繋げたものを乗算結果の整数として返します。結果として返す整数のビット数は、元の整数のビット数の二倍です。

実装
----

考えたとおりに動くか検証のために書いたテストコードが以下です：

```python
class Multiplier(UO_Proc):
    def __init__(self, a: int, b: int, nbits: int = Integer.NBITS) -> None:
        UO_Proc.__init__(self, 'Multiplier')
        self._a = a
        self._b = b
        self.nbits = nbits
    def __enter__(self) -> 'Multiplier':
        UO_Proc.__enter__(self)
        # 入力量子ビットの準備
        a = Integer(self._a, nbits=self.nbits)
        b = Integer(self._b, nbits=self.nbits)
        # a あるいは 0 が入る領域。
        a0 = Integer(0, nbits=a.nbits())
        # 乗算結果を保持する量子ビットのリスト
        results: List[QBit] = []
        t = Integer(0, nbits=a.nbits())
        for i in range(a.nbits()):
            # b の特定 bit との積を a0 に格納する
            Integer.cswap(b.qbits[i], a, a0)
            t += a0
            # a0 の値を 0 に戻す
            Integer.cswap(b.qbits[i], a, a0)
            # t の最下位ビットを results に追加し、t をビットシフト。
            results.append(t.rshift())
        a0.deallocate()
        self.a = a
        self.b = b
        # 最終的に、results および t に乗算結果が入っている。
        self.c = Integer(results + t.qbits, nbits=self.nbits * 2)
        self.c._carry = t.carry()
        return self
for a in range(16):
    for b in range(16):
        with Multiplier(a, b) as p_m:
            print('---- circuit ----')
            circuit = p_m.synthesis(qiskit.QuantumCircuit(
                get_qiskit_register()))
            print(str(p_m))
            circuit.measure_all()
            print('---- result ----')
            sv_sim = qiskit.Aer.get_backend('aer_simulator')
            circuit.save_statevector()
            print(str(circuit))
            job = sv_sim.run(qiskit.assemble(circuit, shots=1))
            result = str(job.result().get_counts())
            print('a={0}[{1}]'.format(p_m.a.parse(result), str(p_m.a)))
            print('b={0}[{1}]'.format(p_m.b.parse(result), str(p_m.b)))
            print('c={0}[{1}]'.format(p_m.c.parse(result), str(p_m.c)))
            if p_m.c.parse(result) == a * b:
                print('{0} OK {1}*{2}={3}'.format(result,
                                                  p_m.a.parse(result),
                                                  p_m.b.parse(result),
                                                  p_m.c.parse(result)))
            else:
                print('{0} NG {1}*{2}={3}'.format(result,
                                                  p_m.a.parse(result),
                                                  p_m.b.parse(result),
                                                  p_m.c.parse(result)))
                sys.exit(1)
        reset_all()
```
Integer クラスのメンバを直接参照しているのが気持ちよくありませんが、とりあえずテストコードなのでこれでよしとします。

ここで、rshift は Integer を表現するビットを右シフトしたうえで、右シフトの結果として Integer 上からなくなるビットを返すメソッドです。これを Integer クラスに追加するなどの改造により上記のコードを動作可能にしました。

テストコードを動かしてみて期待通り動作したので、Integer クラスに __mul__ オペレータとして追加しました。これで、`c = b * a` などと書くと乗算器が合成できます。

```python
    def __mul__(self, other: 'Integer') -> 'Integer':
        u'''引数に与えられた値との乗算を行う。
        '''
        if len(self.qbits) != len(other.qbits):
            raise RuntimeError('invalid Integer bit width')
        with UO_Proc('Integer.mul'):
            # a0 には 0 あるいは other の値が入る
            a0 = Integer(0, nbits=len(self.qbits))
            # 乗算結果（の下位ビット群）を保持する量子ビットのリスト
            results: List[QBit] = []
            # ワークエリアとして用いる量子ビット
            # 最終的に乗算結果の上位ビット群の値となる
            t = Integer(0, nbits=len(self.qbits))
            for i in range(len(self.qbits)):
                # t = other * self.qbits[i]
                Integer.cswap(self.qbits[i], other, a0)
                t += a0
                Integer.cswap(self.qbits[i], other, a0)
                # t の最下位ビットを results に追加し、t をビットシフト。
                results.append(t.rshift())
            a0.deallocate()
            return Integer(results + t.qbits, len(self.qbits) * 2)
```

もう一度、上記のテストコードを __mul__ オペレータを使ったものに置き換えて動作確認しました。


生成例
------

[GitHub](https://github.com/tadashi9e/benchmark_quantum_simulators) にランダムな値の３量子ビットの掛け算を行うベンチマークコードを追加しました。

以下のような回路が得られます：

```
---- circuit ----
TestMul{
 Integer.hadamard{
  H[0],
  H[1],
  H[2]
 },
 Integer.hadamard{
  H[3],
  H[4],
  H[5]
 },
 Integer.mul{
  Integer.cswap{
   CSWAP[0,3,6],
   CSWAP[0,4,7],
   CSWAP[0,5,8]
  },
  Integer.add{
   Carry[12,6,9,13]{ ccx[6,9,13], cx[6,9], ccx[12,9,13] },
   Carry[13,7,10,14]{ ccx[7,10,14], cx[7,10], ccx[13,10,14] },
   Carry[14,8,11,15]{ ccx[8,11,15], cx[8,11], ccx[14,11,15] },
   CX[8,11],
   Sum[14,8,11]{ cx[8,11], cx[14,11] },
   RCarry[13,7,10,14]{ ccx[13,10,14], cx[7,10], ccx[7,10,14] },
   Sum[13,7,10]{ cx[7,10], cx[13,10] },
   RCarry[12,6,9,13]{ ccx[12,9,13], cx[6,9], ccx[6,9,13] },
   Sum[12,6,9]{ cx[6,9], cx[12,9] }
  },
  Integer.cswap{
   CSWAP[0,3,6],
   CSWAP[0,4,7],
   CSWAP[0,5,8]
  },
  Integer.cswap{
   CSWAP[1,3,6],
   CSWAP[1,4,7],
   CSWAP[1,5,8]
  },
  Integer.add{
   Carry[12,6,10,13]{ ccx[6,10,13], cx[6,10], ccx[12,10,13] },
   Carry[13,7,11,14]{ ccx[7,11,14], cx[7,11], ccx[13,11,14] },
   Carry[14,8,15,16]{ ccx[8,15,16], cx[8,15], ccx[14,15,16] },
   CX[8,15],
   Sum[14,8,15]{ cx[8,15], cx[14,15] },
   RCarry[13,7,11,14]{ ccx[13,11,14], cx[7,11], ccx[7,11,14] },
   Sum[13,7,11]{ cx[7,11], cx[13,11] },
   RCarry[12,6,10,13]{ ccx[12,10,13], cx[6,10], ccx[6,10,13] },
   Sum[12,6,10]{ cx[6,10], cx[12,10] }
  },
  Integer.cswap{
   CSWAP[1,3,6],
   CSWAP[1,4,7],
   CSWAP[1,5,8]
  },
  Integer.cswap{
   CSWAP[2,3,6],
   CSWAP[2,4,7],
   CSWAP[2,5,8]
  },
  Integer.add{
   Carry[12,6,11,13]{ ccx[6,11,13], cx[6,11], ccx[12,11,13] },
   Carry[13,7,15,14]{ ccx[7,15,14], cx[7,15], ccx[13,15,14] },
   Carry[14,8,16,17]{ ccx[8,16,17], cx[8,16], ccx[14,16,17] },
   CX[8,16],
   Sum[14,8,16]{ cx[8,16], cx[14,16] },
   RCarry[13,7,15,14]{ ccx[13,15,14], cx[7,15], ccx[7,15,14] },
   Sum[13,7,15]{ cx[7,15], cx[13,15] },
   RCarry[12,6,11,13]{ ccx[12,11,13], cx[6,11], ccx[6,11,13] },
   Sum[12,6,11]{ cx[6,11], cx[12,11] }
  },
  Integer.cswap{
   CSWAP[2,3,6],
   CSWAP[2,4,7],
   CSWAP[2,5,8]
  }
 }
}
a =Integer[0,1,2]
b =Integer[3,4,5]
c =Integer[9,10,11,15,16,17]
---- simulator ----
Circuit(18).h[0].h[1].h[2].h[3].h[4].h[5].ccx[0, 3, 6].ccx[0, 6, 3].ccx[0, 3, 6].ccx[0, 4, 7].ccx[0, 7, 4].ccx[0, 4, 7].ccx[0, 5, 8].ccx[0, 8, 5].ccx[0, 5, 8].ccx[6, 9, 13].cx[6, 9].ccx[12, 9, 13].ccx[7, 10, 14].cx[7, 10].ccx[13, 10, 14].ccx[8, 11, 15].cx[8, 11].ccx[14, 11, 15].cx[8, 11].cx[8, 11].cx[14, 11].ccx[13, 10, 14].cx[7, 10].ccx[7, 10, 14].cx[7, 10].cx[13, 10].ccx[12, 9, 13].cx[6, 9].ccx[6, 9, 13].cx[6, 9].cx[12, 9].ccx[0, 3, 6].ccx[0, 6, 3].ccx[0, 3, 6].ccx[0, 4, 7].ccx[0, 7, 4].ccx[0, 4, 7].ccx[0, 5, 8].ccx[0, 8, 5].ccx[0, 5, 8].ccx[1, 3, 6].ccx[1, 6, 3].ccx[1, 3, 6].ccx[1, 4, 7].ccx[1, 7, 4].ccx[1, 4, 7].ccx[1, 5, 8].ccx[1, 8, 5].ccx[1, 5, 8].ccx[6, 10, 13].cx[6, 10].ccx[12, 10, 13].ccx[7, 11, 14].cx[7, 11].ccx[13, 11, 14].ccx[8, 15, 16].cx[8, 15].ccx[14, 15, 16].cx[8, 15].cx[8, 15].cx[14, 15].ccx[13, 11, 14].cx[7, 11].ccx[7, 11, 14].cx[7, 11].cx[13, 11].ccx[12, 10, 13].cx[6, 10].ccx[6, 10, 13].cx[6, 10].cx[12, 10].ccx[1, 3, 6].ccx[1, 6, 3].ccx[1, 3, 6].ccx[1, 4, 7].ccx[1, 7, 4].ccx[1, 4, 7].ccx[1, 5, 8].ccx[1, 8, 5].ccx[1, 5, 8].ccx[2, 3, 6].ccx[2, 6, 3].ccx[2, 3, 6].ccx[2, 4, 7].ccx[2, 7, 4].ccx[2, 4, 7].ccx[2, 5, 8].ccx[2, 8, 5].ccx[2, 5, 8].ccx[6, 11, 13].cx[6, 11].ccx[12, 11, 13].ccx[7, 15, 14].cx[7, 15].ccx[13, 15, 14].ccx[8, 16, 17].cx[8, 16].ccx[14, 16, 17].cx[8, 16].cx[8, 16].cx[14, 16].ccx[13, 15, 14].cx[7, 15].ccx[7, 15, 14].cx[7, 15].cx[13, 15].ccx[12, 11, 13].cx[6, 11].ccx[6, 11, 13].cx[6, 11].cx[12, 11].ccx[2, 3, 6].ccx[2, 6, 3].ccx[2, 3, 6].ccx[2, 4, 7].ccx[2, 7, 4].ccx[2, 4, 7].ccx[2, 5, 8].ccx[2, 8, 5].ccx[2, 5, 8]
```