FPGA開発日記

FPGAというより、コンピュータアーキテクチャかもね! カテゴリ別記事インデックス <a href="https://msyksphinz.github.io/github_pages/">

Googleの量子コンピュータフレームワークCirqを使ってみる(3. Cirqのチュートリアルを試す)

量子コンピュータの勉強を続けており、Cirqの資料をひたすら読み進めている。

以下の資料を読みながら勉強中である。

  • Cir: Tutorial

Tutorial — Cirq 0.1 documentation

Cirqのチュートリアルは以下に日本語訳と解説がなされており、非常にありがたい。

qiita.com

波動方程式を作ってシミュレーションを行うチュートリアルが紹介されているのだが、量子力学における状態ベクトルの数学的な要素について分かっていないのでいろいろ資料を読みながら進めている。 これは一度しっかり量子力学の本を買った方が良いのだろうか。。。?

f:id:msyksphinz:20180926011313p:plain
図. Cirqで作成される量子回路の概念図 (https://cirq.readthedocs.io/en/latest/tutorial.html#background-variational-quantum-algorithms より抜粋)

固有値 - Wikipedia

上記のCirqの例では、イジングモデルを作成し、まずは量子ビット間の関係を示す$jr$, $jc$をランダムに決める。 その値に従って量子ビットに対してXゲートを適用している。

つまり、ランダムなので最初の量子回路の形は毎回異なるらしい。なんだ。

  • 1回目
transverse fields: [[-1, 1, -1], [-1, -1, -1], [-1, -1, -1]]
row j fields: [[-1, -1, -1], [1, 1, -1]]
column j fields: [[-1, 1], [-1, -1], [1, 1]]
(0, 0): ───X^0.1───X───────@───────X───────────────────────────────────────────────────────────────────────M('x')───
                           │                                                                               │
(0, 1): ───X^0.1───Z^0.2───┼───────X───@───────X───────────────────────────────────────────@───────────────M────────
                           │           │                                                   │               │
(0, 2): ───X^0.1───────────┼───────────┼───────X───@───────X───────────────────────────────@^0.3───────────M────────
                           │           │           │                                                       │
(1, 0): ───X^0.1───X───────@^0.3───X───┼───────────┼───────────@───────────────────────────────────────────M────────
                                       │           │           │                                           │
(1, 1): ───X^0.1───────────────────X───@^0.3───X───┼───────────┼───────@───────────────────────────────────M────────
                                                   │           │       │                                   │
(1, 2): ───X^0.1───────────────────────────────X───@^0.3───X───┼───────┼───────X───@───────X───────────────M────────
                                                               │       │           │                       │
(2, 0): ───X^0.1───────────────────────────────────────────────@^0.3───┼───────────┼───────@───────────────M────────
                                                                       │           │       │               │
(2, 1): ───X^0.1───────────────────────────────────────────────────────@^0.3───────┼───────@^0.3───@───────M────────
                                                                                   │               │       │
(2, 2): ───X^0.1───────────────────────────────────────────────────────────────X───@^0.3───X───────@^0.3───M────────
  • 2回目
transverse fields: [[1, 1, -1], [1, 1, -1], [1, -1, -1]]
row j fields: [[1, -1, 1], [1, 1, -1]]
column j fields: [[1, -1], [-1, 1], [-1, 1]]
(0, 0): ───X^0.1───Z^0.2───@───────────────────────────────────────────────────────────@───────────────M('x')───
                           │                                                           │               │
(0, 1): ───X^0.1───Z^0.2───┼───────X───@───────X───────────────────────────────────────@^0.3───────────M────────
                           │           │                                                               │
(0, 2): ───X^0.1───────────┼───────────┼───────────@───────────────────────────────────────────────────M────────
                           │           │           │                                                   │
(1, 0): ───X^0.1───Z^0.2───@^0.3───────┼───────────┼───────@───────────────────────────────────────────M────────
                                       │           │       │                                           │
(1, 1): ───X^0.1───Z^0.2───X───────────@^0.3───X───┼───────┼───────@───────────────────────────@───────M────────
                                                   │       │       │                           │       │
(1, 2): ───X^0.1───────────────────────────────────@^0.3───┼───────┼───────X───@───────X───────@^0.3───M────────
                                                           │       │           │                       │
(2, 0): ───X^0.1───Z^0.2───────────────────────────────────@^0.3───┼───────────┼───────────────────────M────────
                                                                   │           │                       │
(2, 1): ───X^0.1───────────────────────────────────────────────────@^0.3───────┼───────────────@───────M────────
                                                                               │               │       │
(2, 2): ───X^0.1───────────────────────────────────────────────────X───────────@^0.3───X───────@^0.3───M────────

最終的にシミュレーションを行っているのだが、これはまだ何をやっているのか良く分かっていない。。。

CirqでのQubitsの表現

CirqでのQubitsの表現はどのようになっているのだろう? Cirqの実装を探ってみた。on()という関数で量子ビットに演算が適用されるのだが、そのあと具体的にどのように演算が適用されるのかが分からない。。。

  • Cirq/cirq/devices/grid_qubit.py
class GridQubit(QubitId):
    """A qubit on a 2d square lattice.

    GridQubits use row-major ordering:

        GridQubit(0, 0) < GridQubit(0, 1) < GridQubit(1, 0) < GridQubit(1, 1)
    """

    def __init__(self, row, col):
        self.row = row
        self.col = col
...
  • Cirq/cirq/ops/common_gates.py
class CNotGate(eigen_gate.EigenGate,
               gate_features.TextDiagrammable,
               gate_features.CompositeGate,
               gate_features.TwoQubitGate,
               gate_features.QasmConvertibleGate):
    """When applying CNOT (controlled-not) to QuBits, you can either use
    positional arguments CNOT(q1, q2), where q2 is toggled when q1 is on,
    or named arguments CNOT(control=q1, target=q2).
    (Mixing the two is not permitted.)"""
...
    def on(self, *args: raw_types.QubitId,
           **kwargs: raw_types.QubitId) -> gate_operation.GateOperation:
        if not kwargs:
            return super().on(*args)
        if not args and set(kwargs.keys()) == {'control', 'target'}:
            return super().on(kwargs['control'], kwargs['target'])
        raise ValueError(
            "Expected two positional argument or else 'target' AND 'control' "
            "keyword arguments. But got args={!r}, kwargs={!r}.".format(
                args, kwargs))

TinyEMU Buildrootを試行する

TinyEMUはQEMUを開発したFabrice Bellard氏が開発したRISC-Vシミュレータだ。

さらに彼はRISC-VのLinuxを起動するためのbuildrootも用意したらしい。RISC-VのMLにアナウンスがあった。

TinyEMU

とりあえず超高速のシミュレータでRISC-V Linuxを起動することができるみたいなので、試行してみた。

Buildroot / ツールチェインのダウンロード

まずはBellard氏が作ったツールチェインとイメージをすべてダウンロードする。

curl -L https://bellard.org/tinyemu/buildroot-riscv-2018-09-22.tar.gz | tar xz
curl -L https://bellard.org/tinyemu/diskimage-linux-riscv-2018-09-22.tar.gz | tar xz
curl -L https://bellard.org/tinyemu/tinyemu-2018-09-15.tar.gz | tar xz

次に、buildrootを構成する。今回は64-bit版で構成した。

cd buildroot-riscv-2018-09-22/
cp configs/riscv64_defconfig .config
make menuconfig
make

正常に終了すると、bulidrootのイメージが作成される。

TinyEMUのビルド

次に、RISC-VシミュレータであるTinyEMUをビルドする。

cd temu-2018-09-15
make

Linuxの実行

最後にLinuxを実行する。先ほど展開したdiskimage-linux-riscv-2018-09-22ディレクトリに移り、設定ファイルを変更する。 buildroot-riscv64.cfg内のファイルパスに問題があれば変更する。今回は変更せずに実行できる。

  • buildroot-riscv64.cfg
/* VM configuration file */
{
    version: 1,
    machine: "riscv64",
    memory_size: 256,
    bios: "bbl64.bin",
    cmdline: "console=hvc0 root=/dev/vda rw",
    drive0: { file: "../buildroot-riscv-2018-09-22/output/images/rootfs.ext2" },
    eth0: { driver: "user" },
}

これでシミュレータを実行した。簡単に実行できる。さすがだなあ。

../temu-2018-09-15/temu buildroot-riscv64.cfg
...
Jan  1 00:00:00 localhost user.info kernel: [    0.199194] devtmpfs: mounted
Jan  1 00:00:00 localhost user.info kernel: [    0.199851] Freeing unused kernel memory: 80K
Jan  1 00:00:00 localhost user.warn kernel: [    0.199943] This architecture does not have kernel memory protection.
Jan  1 00:00:00 localhost user.info kernel: [    0.233068] EXT4-fs (vda): re-mounted. Opts: block_validity,barrier,user_xattr,errors=remount-ro
Starting network...
ifup: can't open '/var/run/ifstate': No such file or directory
Starting dropbear sshd: OK
Jan  1 00:00:00 localhost daemon.info init: starting pid 74, tty '/dev/console': '/sbin/getty -L  console 0 vt100 '
Jan  1 00:00:00 localhost authpriv.info dropbear[73]: Running in background

Welcome to Buildroot
localhost login:
f:id:msyksphinz:20180923220808g:plain

Googleの量子コンピュータフレームワークCirqを使ってみる(2. 加算器を作ってみる)

量子プログラミングについて少しずつ分かってきたので、もう少しステップを進めていきたい。

自作の量子コンピュータ・シミュレータの拡張を考えるためにも、GoogleMicrosoftが発表している量子コンピュータフレームワークについて勉強し、そこでいろいろな量子プログラミングの手法を勉強すべきだろう。

という訳で、まずはGoogle量子コンピュータフレームワークであるcirqを使ってみることにした。 試してみたのは、もちろん全加算器(1-bit)からだ。

Google Cirqのチュートリアル

Google Cirqのチュートリアルは以下で公開されている。

Tutorial — Cirq 0.1 documentation

これによると、基本的な考え方は自作シミュレータで培ったものと同じで、少し違うのはグリッドというものに拡張されていることだ。 グリッドといっても、2次元であるだけで別に1次元で考えても大して変わらない。

まずは、このグリッドを使って4つの量子ビットを定義する。

qbits = [cirq.GridQubit(0, i) for i in range(4)]

次に、量子回路を作って全加算器を構成する。単純だ。

circuit = cirq.Circuit()
circuit.append(cirq.CCX.on(qbits[0], qbits[1], qbits[3]))
circuit.append(cirq.CNOT.on(qbits[0], qbits[1]))
circuit.append(cirq.CCX.on(qbits[1], qbits[2], qbits[3]))
circuit.append(cirq.CNOT.on(qbits[1], qbits[2]))
circuit.append(cirq.CNOT.on(qbits[0], qbits[1]))
print(circuit)

以下のようなグラフが出力された。正しく構成できていると思う。

f:id:msyksphinz:20180923014256p:plain
図. cirqで構成した全加算器(1-bit)

全加算器のシミュレーションをする

全加算器のシミュレーションをするためには、最初の量子ビットの状態(qbits[0]qbits[1]の値)を変更したいのだが、どうもその方法が分からない。

仕方がないので、量子ビットの先頭に反転ゲートを挿入して、それっぽく動作させることにした。 つまり、a+b = 1+0 を実行したければ、aが格納されている量子ビットqbits[0]を最初に反転させるわけだ。

...
circuit = cirq.Circuit()
# circuit.append(cirq.X(qbits[0]))   # Set initial value of a
# circuit.append(cirq.X(qbits[1]))   # Set initial value of b
circuit.append(cirq.CCX.on(qbits[0], qbits[1], qbits[3]))
circuit.append(cirq.CNOT.on(qbits[0], qbits[1]))
circuit.append(cirq.CCX.on(qbits[1], qbits[2], qbits[3]))

次にシミュレーションを行う。最後に結果を観測するために観測器を取り付けた。measureがそれに当たる。

simulator = cirq.google.XmonSimulator()
circuit.append(cirq.measure(qbits[2], key='s'))
circuit.append(cirq.measure(qbits[3], key='c'))
print(circuit)

results = simulator.run(circuit, repetitions=100, qubit_order=qbits)
print (results.histogram(key='s'))
print (results.histogram(key='c'))

実行すると以下のようになった。

$ python full-adder.py
[GridQubit(0, 0), GridQubit(0, 1), GridQubit(0, 2), GridQubit(0, 3)]
(0, 0): ───@───@───────────@───
           │   │           │
(0, 1): ───@───X───@───@───X───
           │       │   │
(0, 2): ───┼───────@───X───────
           │       │
(0, 3): ───X───────X───────────
(0, 0): ───@───@───────────@─────────────────────
           │   │           │
(0, 1): ───@───X───@───@───X─────────────────────
           │       │   │
(0, 2): ───┼───────@───X───────M('s')────────────
           │       │
(0, 3): ───X───────X────────────────────M('c')───
Counter({0: 100})
Counter({0: 100})

最後に、qbits[2]の値(=s)が0、qbits[3]の値(=c)が0となっている。 厳密には、100回シミュレーションしてそれぞれのビットが0となっている回数が100回(つまり全部)となった。

たとえば、上記のコメントを外してa+b=1+1を実行してみる。

[GridQubit(0, 0), GridQubit(0, 1), GridQubit(0, 2), GridQubit(0, 3)]
(0, 0): ───X───────@───@───────────@───
                   │   │           │
(0, 1): ───────X───@───X───@───@───X───
                   │       │   │
(0, 2): ───────────┼───────@───X───────
                   │       │
(0, 3): ───────────X───────X───────────
(0, 0): ───X───────@───@───────────@─────────────────────
                   │   │           │
(0, 1): ───────X───@───X───@───@───X─────────────────────
                   │       │   │
(0, 2): ───────────┼───────@───X───────M('s')────────────
                   │       │
(0, 3): ───────────X───────X────────────────────M('c')───
Counter({0: 100})
Counter({1: 100})

この場合は、s=0, c=1となっているのキャリーが発生している。正しく動作したようだ。

ちなみに、3ビットの加算器だと以下のようになる。まだテストしてないけど。

f:id:msyksphinz:20180923015903p:plain
図. 3ビットの量子加算回路

ゼロから作る量子コンピュータ・シミュレータ(3. 量子回路で3ビット加算器を作ろう)

前回は半加算器が完成した。という訳で、これを拡張させて全加算器を作成し、さらにこれを数珠つなぎにして3ビットの加算回路を作ろう。

量子回路で全加算器を作ろう

量子回路を使って全加算器を作るにはどうしたらよいのか、少し検索するとすぐ回路を見つけた。 以下のように構成すれば全加算器を作ることができるらしい。

f:id:msyksphinz:20180922164723p:plain
図.  a+b+c_\text{in} = s, c_\text{out} の全加算器

これは私が開発している量子コンピュータシミュレータを使えば、以下のように表現できる。

void FullAdder(Qbits *qbits, int a_idx, int b_idx, int c0_idx, int c1_idx)
{
  qbits->CCNot(a_idx, b_idx, c1_idx);
  qbits->CNot (a_idx, b_idx);
  qbits->CCNot(b_idx, c0_idx, c1_idx);
  qbits->CNot (b_idx, c0_idx);
}

これだとbの値は最終的に崩れてしまうのだけれども、使わないので良しとする。

3ビットの加算器を作ろう

3ビットの加算器を作るためには、この全加算器を3つ並べれば良いわけだ。 もっとも単純に作るためには以下のようにする。 ちなみに、各計算毎にキャリーを保持するために量子ビットを用意している。 なので、3ビットの加算器を作るために用意したビットは、 - 3ビットで表現される入力a - 3ビットで表現される入力b - 4ビットの、答え兼carry-bitを保持するための答えc

という訳で、3+3+4ビットで10ビット必要である。この時点で内部的には 1024\times 1024の行列計算が行われるわけで、こりゃ確かに大変だ。。。

f:id:msyksphinz:20180922165310p:plain
図. 3ビット量子加算回路

一度全加算器を作ったので、これを並べるだけである。 量子ビットは10ビット用意して、{s[3:0], b[2:0], a[2:0]}とした。なので、最後に6ビットシフトしてs[3:0]のみ出力している。

      Qbits qbits (b << 3 | a, 10);
      FullAdder (&qbits, 0, 3, 6, 7);
      FullAdder (&qbits, 1, 4, 7, 8);
      FullAdder (&qbits, 2, 5, 8, 9);
      int res = qbits.Measure();
      res = res >> 6;
      printf ("%x\n", res);

加算回路の実行

とりあえず全部回してみる。目視だが間違いがあればすぐに分かるはずだ。

  for (int a = 0; a < 8; a++) {
    for (int b = 0; b < 8; b++) {
      printf ("%x + %x = ", a, b);

      Qbits qbits (b << 3 | a, 10);
      FullAdder (&qbits, 0, 3, 6, 7);
      FullAdder (&qbits, 1, 4, 7, 8);
      FullAdder (&qbits, 2, 5, 8, 9);
      int res = qbits.Measure();
      res = res >> 6;
      printf ("%x\n", res);
    }
  }

実行すると以下のようになった。ちゃんと動いているな。

$ ./qbit_simulator
0 + 0 = 0
0 + 1 = 1
0 + 2 = 2
0 + 3 = 3
0 + 4 = 4
0 + 5 = 5
0 + 6 = 6
0 + 7 = 7
1 + 0 = 1
...
7 + 4 = b
7 + 5 = c
7 + 6 = d
7 + 7 = e

RISC-V Day Tokyo 2018 / Fukuoka で登壇します

f:id:msyksphinz:20180921223849p:plain

RISC-V Day Tokyo 2018が10/18(木) に開催されます。私はRISC-Vを支える技術周辺、エコシステム周りを発表する予定です。 15分しかしゃべらないのですが、興味のある方はどうぞ。

tmt.knect365.com

やっぱりメインはRISC-V Readerの日本語訳の配布、そしてDavid Ditzel氏が来日して開発中のRISC-Vプロセッサについて登壇される予定です。 またRISC-Vの開発の中心であるKrste Asanovic氏、Andrew Waterman氏も来日されます。いろいろお話が聞けるチャンスです。

ちなみに、1週間前でもないのに今頃告知をしているかというと、今日09/21まで、参加費が安くなるEarly Accessの期間中だからです。 投稿されてから1時間くらいしか残されていないけど、参加したい人は早く参加登録を!!今なら$50で参加可能です。

なお、10/20(土)にはIEEE MICRO51のワークショップとしてRISC-V Dayが開催される予定です。 こちらは西日本在住で東京に来られない方、IEEEの会員の方は是非どうぞ。MICROとの併設なので参加費めっちゃ高いですけど。

f:id:msyksphinz:20180921225219p:plain

ゼロから作る量子コンピュータ・シミュレータ(2. 複数量子ビットをモデル化して加算回路を作ろう)

前回は1量子ビットC++でモデル化したのだが、それだけではほとんど何も表現出来ないので複数量子ビットに拡張したい。

複数量子ビットが表現できるようになれば、まずは例として加算器を動かしてみよう。 半加算器と全加算器が実現できればかなり面白い。

量子計算に大規模なコンピュータが必要な理由

最初は複数量子ビットを表現するのに量子ビットの数に線形比例するサイズの配列を用意すればよいのかと思っていたが、 考え方を完全に勘違いしていた。

そもそも量子ビットの表現には\alpha, \betaの2種類が必要というところから始まる。 量子ビットが0であるということは\begin{pmatrix}
\alpha \\
\beta
\end{pmatrix} = \begin{pmatrix}
1 \\
0
\end{pmatrix}
であり、一方で量子ビットが1であるということは\begin{pmatrix}
\alpha \\
\beta
\end{pmatrix} = \begin{pmatrix}
0 \\
1
\end{pmatrix} であるということを示している。 これはつまり量子ビットの1ビットを表現するために空間としては2通り必要と言う事である。 「0である場合」の量子ゲートの処理と「1である場合」の量子ゲートの処理が一度に表現できる。

つまり、上記の1量子ビットに対して2\times 2の行列積を適用すると言う事は「量子の状態はわからないけれども、0にときはこっちの計算、1のときはこっちの計算を適用するよ」ということを同時に表現しているに等しい。

例えば、量子を向きを反転させるゲートだと$$\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}\begin{pmatrix} \alpha \\ \beta \end{pmatrix}$$ という行列を適用するが、量子ビット\begin{pmatrix}
\alpha \\
\beta
\end{pmatrix}

  • 上側の要素は量子が0であることを示す。
  • 下側の要素は量子が1であることを示す。

というようにデコードできるとする。 上記のゲートを適用すると\begin{pmatrix}
\beta\\
\alpha
\end{pmatrix} となり、量子が0である確率と1である確率は反転する。

このように考えると、量子状態の表現では、量子ビットの数が増えるとその量子の状態を示す数だけ行列の大きさが増えると言う事がわかる。 つまり、

  • 2量子ビットの場合 : 00, 01, 10, 11 の4状態 →行列のサイズは4
  • 3量子ビットの場合 : 000, 001, 010, 011, 100, 101, 110, 111 の8状態 → 行列のサイズは8

となり、量子ビットの数がn_qbitsであるとすると、シミュレータ上で確保しなければならない行列のサイズは 2^\text{n_qbits}となり指数関数的に増大する。 これが例えば50量子ビットになると2^{50}となり膨大な空間が必要なため、メモリを節約しつつコンピュータ上で計算する方法が必要不可欠となる訳だ。(といってもかなり疎行列なので圧縮する方法は沢山ありそうだが)

複数量子ビットC++で表現しよう

量子ビットを表現する原理が少しわかってきたので、次にこれをどのようにプログラムで表現するか考えた。 基本的には、nビットの量子ビットを定義する場合、2^{n}長の配列を定義するのがよさそうだ。 配列の各要素には、量子がその状態をとる確率が入るという理解で良い気がする。

理解があいまいなのだが、複数量子ビットを表現する場合は行列の値にはどのような制限があるのだろうか? 1量子ビットだと、ブロッホ球上で\alpha^{2}+\beta^{2}=1という空間内に収まっている必要がある (重ね合わせの状態だと\alpha=\dfrac{1}{\sqrt{2}}, \beta=\dfrac{1}{\sqrt{2}}とか)。

これが複数量子ビットになるとどうなるんだ? 基本的にはブロッホ球が複数あるとして、それを1つの行列で表現したときにどのように制約を付け加えればよいのか良く分からない。

2018/09/30 追記。複数量子ビットの場合、2^{n}長の配列の各要素を\psi_iと表現すると、制約は\sum \psi_{i}^2 = 1 となるらしい。 コメントで教えていただいた。ありがとうございます。

複数量子ビットを表現するためにQbitsクラスを定義した。

class Qbits {
  std::complex<double> *m_comp_amp;
  const size_t m_size;

コンストラクタでは量子ビットの初期値と、量子ビットのビット長を指定する。

  Qbits (uint8_t val, const size_t size)
      : m_size(pow(2, size))
  {
    m_comp_amp = new std::complex<double>[m_size];
    for (size_t i = 0; i < size; i++) {
      m_comp_amp[i] = 0.0 + 0.0i;
    }
    m_comp_amp[val] = 1.0 + 0.0i;
  }

例えば3ビットの量子ビットを定義して、初期値が4である場合は配列の4番目が1.0で、それ以外は0.0とする。

複数量子ビットに採用する量子ゲートを作ろう

さて、複数量子ビットが定義できたのでさっそく何かしたい。 まずは半加算器を作る工程を考えて、必要なゲートを定義することにする。

半加算器は一般的に以下のような量子回路で定義されるらしい。

f:id:msyksphinz:20180920230051p:plain
図. 半加算器

ここでCCNOTはトリフォゲートと呼ばれる。詳しくはWikipediaを参照のこと。 要するにある量子ビットを反転させるときに、その量子ビット以外の量子ビットがすべて1ならば反転し、そうでないときは反転しない。

トフォリゲート - Wikipedia

f:id:msyksphinz:20180920230146p:plain
図. トフォリゲートの行列表現 (https://ja.wikipedia.org/wiki/トフォリゲートより抜粋)

一方でCNOTはこれの2ビット版である。2量子ビット(q1, q0)において、q1が1ならばq0が反転し、q1が0ならばq0は反転しない。

行列式では下記のように表現できるが、3ビット以上の量子ビットにも対応させたいことと、 反転するターゲットの量子ビットと残りの2量子ビットを指定する必要がある。

このため、CCNot()では3引数を取り、

  • 反転のターゲットとなる量子ビットのインデックス: v_idx
  • 反転の判断をするための量子ビットのインデックス : c0_idx, c1_idx

とした。CNOTはその特殊版である(c0_idx == c1_idx)。

f:id:msyksphinz:20180920231536p:plain
図. 3量子ビットにおけるCNOT(Q1, Q0)
f:id:msyksphinz:20180920231739p:plain
図. 3量子ビットにおけるCNOT(Q2, Q0)
f:id:msyksphinz:20180920231728p:plain
  void CNot (int c_idx, int v_idx) {
    CCNot (c_idx, c_idx, v_idx);
  }
  void CCNot (int c0_idx, int c1_idx, int v_idx) {
    std::complex<double> *new_comp_amp = new std::complex<double>[m_size];

    uint8_t c_bitmask = (1 << c0_idx) | (1 << c1_idx);
    uint8_t v_bitmask = 1 << v_idx;
    for (int idx = 0; idx < m_size; idx++) {
      int act_idx = idx;
      if ((idx & c_bitmask) == c_bitmask) {
        act_idx = idx ^ v_bitmask;
      }
      new_comp_amp[act_idx] = m_comp_amp[idx];
    }
    for (int idx = 0; idx < m_size; idx++) {
      m_comp_amp[idx] = new_comp_amp[idx];
    }
  }

たぶんもっと簡単に実装できると思う。プログラミング力が無い。。。 要するにbitmaskとなっているビットがすべて1であれば、ターゲットのビット位置を反転する。

勘違いしがちだが値を反転するのではなくて、配列の要素を入れ替える。 例えば、3ビットの量子ビット[xyz]=[110]で、反転するターゲットの量子ビットzとする。 判断を行うxyの2ビットが両方1なので、zは反転するのだが、あくまで8要素(3ビットの量子ビットなので) [110](=6)[111](=7)の要素を入れ替えるだけなのだ。 Wikipediaに乗っている8\times 8の配列はまさにそれを示している。

量子ビットを使って半加算器を作ろう

というわけで、さっそく作ったCNot()CCNot()を使って半加算器を構成した。

  for (int val = 0; val < 4; val++) {
    Qbits qbits(val, 3);
    qbits.CCNot(1, 0, 2); // 量子ビット0,1の値によって、量子ビット2を反転させる
    qbits.CNot(0, 1);     // 量子ビット0の値によって、量子ビット1を反転させる
    printf("qbit val = %d\n", qbits.Measure());
  }

3量子ビットを定義した。

  • ビット0 : 入力a
  • ビット1 : 入力b, かつ計算後にはsum
  • ビット2 : 計算後にcarryが入る

とする。入力を0から3まで変えて計算すると、それぞれ以下のようになった。

qbit val = 0
qbit val = 3
qbit val = 2
qbit val = 5

そっけない結果だが、

  •  0+0 : sum=0, carry=0([000]->[000])
  •  0+1 : sum=1, carry=0([001]->[011])
  •  1+0 : sum=1, carry=0([010]->[010])
  •  1+1 : sum=0, carry=1([011]->[101])

というビットが正しく表現できている。半加算器の完成だ。

ゼロから作る量子コンピュータ・シミュレータ(1. 1量子ビットの表現)

次世代計算機講座に参加して、量子コンピュータについて俄然興味がわいてきた。

手っ取り早く量子コンピュータを試すためには、最近になっていろいろ手段が出てきている。

MicrosoftからはQ#、GoogleからはフレームワークとしてCirqがリリースされている。 これらの環境を使うのもいいけれども、あまり情報もないし折角なので自分で作ってみたい。

という訳で量子コンピュータシミュレータを1から自分で作るにはどうすればよいのか調べてみる。

参考にしたのは、次世代計算機講座で講演なされた京都大学の藤井啓祐先生の資料と、同じく藤井先生がCQ出版のInterfaceに寄稿された量子コンピュータC言語モデルだ。 Interfaceに分かりやすい形で量子コンピュータのモデルが解説されているのは非常にありがたい。

という訳で、まずはInterfaceの記事のコードをどんどんコピペしていく。

1量子ビットの表現

1量子ビット複素数の変数として\alpha, \betaの組で表現されるので、複素型として2変数を用意する。

class Qbit {
  double _Complex m_comp_amp[2];

1量子ビットに対する操作

\alpha, \beta がどの程度傾いているかにより、0に近いか1に近いかが決まる。 そのランダム性は観測時に影響を与える。 ここでは、 0.0<= val <1.0 の乱数を生成して、その値と\alpha, \betaの大小で最終的な0/1を決める。

  int Measure () {
    double rand_val = rand() / static_cast<double>(RAND_MAX);
    double measure_val = pow(cabs(m_comp_amp[0]), 2);
    if (rand_val < pow(cabs(m_comp_amp[0]), 2)) {
      m_comp_amp[0] = 1;
      m_comp_amp[1] = 0;
      return 0;
    } else {
      m_comp_amp[0] = 0;
      m_comp_amp[1] = 1;
      return 1;
    }
  }

各種ゲートの適用

とりあえずは、アダマールゲートと回転ゲートを作ってみる。 アダマールゲートは以下のように表現される。

$$ \dfrac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} $$

次に、回転ゲートは以下のように表現される。

$$ \dfrac{1}{\sqrt{2}}\begin{pmatrix} 1 & 0 \\ 0 & e^{i\theta} \end{pmatrix} $$

というわけで、それぞれQbitクラスのメソッドとして以下のように表現した(というかかなりコピペした)。

  void Hadamard () {
    m_comp_amp[0] = (m_comp_amp[0] + m_comp_amp[1]) / sqrt(2);
    m_comp_amp[1] = (m_comp_amp[0] - m_comp_amp[1]) / sqrt(2);
  }

  void Rotate (double theta) {
    m_comp_amp[0] = m_comp_amp[0];
    m_comp_amp[1] = (cos(theta) + 1.0i * sin(theta)) * m_comp_amp[1];
  }

回転ゲートの動作確認

という訳で、簡単な1量子ビットとそれに適用できるゲートは作ったので、藤井先生の例に出てくる1量子ビット演算を実行してみる。

f:id:msyksphinz:20180919003653p:plain

ここで、回転ゲートの回転量を \dfrac{2\pi}{10}ずつ細切れにして回転させて、1が出現する確率を測定する。

出現する確率は、1000回同様の量子計算を繰り返して測定する。つまり、以下のようになる。

  for (int trial = 0; trial <= 10; trial++) {
    int count_one = 0;
    for (int count = 0; count < 1000; count++) {
      Qbit q_bit;
      q_bit.SetZero();
      q_bit.Hadamard();
      q_bit.Rotate(2 * M_PI * 0.1 * trial);
      q_bit.Hadamard();

      if (q_bit.Measure() == 1) {
        count_one++;
      }
    }
    printf ("Trial = %d : count_one = %d\n", trial, count_one);
  }

上記のプログラムを実行すると、以下のようになった。

$ ./qbit_simulator
Trial = 0 : count_one = 277
Trial = 1 : count_one = 340
Trial = 2 : count_one = 515
Trial = 3 : count_one = 737
Trial = 4 : count_one = 918
Trial = 5 : count_one = 982
Trial = 6 : count_one = 914
Trial = 7 : count_one = 748
Trial = 8 : count_one = 495
Trial = 9 : count_one = 331
Trial = 10 : count_one = 278

グラフを書くと、以下のようになった。

f:id:msyksphinz:20180919003956p:plain

綺麗な三角関数が描けている。正しく計算できているようだ。

とりあえずはここまで。