量子コンピュータの華、ショアのアルゴリズムを 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) となる。

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 の系。

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

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

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 をかけ、測定を行う

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

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