量子コンピュータで因数分解 (sympy)

量子コンピュータの華、ショアのアルゴリズムを sympy を使って試してみる。

# 古典的に解くよ

以下のような手順で、整数 N の因数分解が可能

- Nより小さい適当な数字 a を持ってくる
- a の位数 r を求める
    - 位数 r とは、aを何乗するとaに戻ってくるか(mod N)の数値
- N, a, r を元にちょっとした計算をして、二つの整数を作る
- その2つの数字と、N の最大公約数（GCD）を計算すると
- N の因数がでてくる

具体的には、以下のコードを見たほうが早かろう。15 を因数分解するために適当に 8という数字を選んで、位数を計算すると 4 であることがわかる。この数字をもとに 65, 63 という数字を得る。この数字と 15 の最大公約数は 5, 3 となり、これが答え ( 3x5 = 15) となる。

```py
import random
random.seed(0)

def pick_a(number):
  return random.randrange(2, number - 1)

def find_order(n, a):
  order, work = 0, a
  while True:
    order += 1
    work = (work * a) % n
    if work == a:
        return order

def gcd(a, b):
  if a < b:
    a, b = b, a
  return b if a % b == 0 else gcd(b, a % b)

def generate_factor(n, a, r):
  return gcd(n, a**int(r/2)+1), gcd(n, a**int(r/2)-1)

n = 15
a = pick_a(n) # ->8
r = find_order(n, a) # ->4
answer = generate_factor(n, a, r) # -> 5, 3
```


# 量子的に解くよ
上記の位数を求める部分 `find_order(n,a)` が、ものすごい時間がかかる（らしい）ので、この部分を量子回路で実装する。以下の手順で位数が求められることが知られている。

- `|a>|x>` を用意
    - `|a>` は a の値をエンコードした状態
    - `|x>` は作業用で、N を格納できるサイズのqubit からなる状態の linear superposition
- CModゲートを使ってこの二つをエンタングルさせて
- `|x>` 側で量子フーリエ変換し、観測をする
    - 非常に限られた値しか測定にかからない
    - 具体的には `|16/r>, |2*16/r>, |3*16/r>, ....`
    - よって `r` の値を推定できる。

n = 15, a = 8 を例に、sympy で実装してみよう。**目標は N = 15, a = 8 の際の、位数 r = 4 を求めること**である。まずは、始状態を構築しよう。n =15 なので、4-qubit あれば十分格納できる。a 格納用、 x 格納用の合わせて 8-qubit の系。

```py
# |a> = |1000> = |8>, |x> = H * |0000>
state0 = qapply(H(0) * H(1) * H(2) * H(3) * Qubit('10000000'))
```

これに CMod と QFT16 をかけてあげればよいのだが、本日時点(2018April)では sympyの **CMod実装が動作しない**。GateAPIの変更に追従していないからのようだ。仕方ないので自分でそれっぽいものを作る。

```py
class MyCMod(Gate):
  def _apply_operator_Qubit(self, qubits, **options):
    assert len(qubits) == 8
    N = 15
    left = qubits.qubit_values[0:4] # a is encoded, left-hand
    right = qubits.qubit_values[4:8]  # x is encoded, right-hand
    x = IntQubit(*right).as_int()
    a = IntQubit(*left).as_int()
    assert a == 8
    perm = int(a**x % N)

    out = []
    for i in reversed(range(4)):
      out.append(perm >> i & 1)
    out = out + list(right)
    return Qubit(*out)
```

改めて、始状態にCMod, QFT16 をかけ、測定を行う

```py
state0 = qapply(H(0) * H(1) * H(2) * H(3) * Qubit('10000000'))
state1 = qapply(QFT(0,4) * MyCMod() * state0)
for state, prob in measure_partial(state1, (0,1,2,3)):
  print ( prob, state )
```

すると、測定結果は、4パターンしか出てこない。

```
1/4 |00010000>/2 + |00100000>/2 + |01000000>/2 + |10000000>/2
1/4 |00010100>/2 - I*|00100100>/2 - |01000100>/2 + I*|10000100>/2
1/4 |00011000>/2 - |00101000>/2 + |01001000>/2 - |10001000>/2
1/4 |00011100>/2 + I*|00101100>/2 - |01001100>/2 - I*|10001100>/2
```

すなわち、 `|0>, |4>, |8>, |12>` である。これらが `|16 * k/r>` の候補であるので、位数`r = 4`であることが分かる。（注: 1回目の測定で |0> が出た場合などは即座に r はわからないが、何度か測定すればわかる)

以上、量子回路で位数を求めることができた。


# 深く学ぶために

## 考慮1: 一般の位数 r について
実は、上記の処方は、位数r が 2のべき乗の時にしかうまく動作しない。例では位数が 4 であり、確かに2のべき乗であった。一般の r の場合、測定結果 ≈ `|N *k/r >`  から `r` を類推するのにちょっとしたテクニック（連分数）を使うが、本稿がこれ以上長くなるのは嫌なので、参考資料を参考されたい。

## 考慮2: 失敗例
このアルゴリズムは、確率的にしかうまくいかない。何回かやるとうまくいく、というタイプである。そこで、どんな時に失敗するのかを見てみよう。

- 観測結果が |0> である
何度か計測して、0以外の観測結果を得よう

- 位数が奇数になった場合
`a**int(r/2)+1` が整数にならないので、これ以上計算が進められない。別の a を選んで計算をやり直そう

- 因数がトリビアルの場合
つまり、因数分解結果が 1 * N となることがある。この場合は、別の a を選んで計算をやり直そう

## 考慮3: 物理的な実装
今回使った `MyCMod`, `QFT16` はいずれも具体的な物理的な実装を意識していない。`CNOT` や `H, Z` などの基本的な素子でどのようにこれらの量子ゲートが表現されるかは、参考文献を参考のこと。


参考文献: https://www.appi.keio.ac.jp/Itoh_group/abe/pdf/qc6.pdf
https://www.appi.keio.ac.jp/Itoh_group/abe/pdf/qc7.pdf