FPGA開発日記

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

RISC-Vブート時のDevice TreeをParseする仕組み

f:id:msyksphinz:20181004231514p:plain

SiFiveをはじめとするRISC-Vの多くの実装では、プログラムの起動時にDevice Treeを読み取りに行くケースが多い。 事実として、RISC-Vの命令セットシミュレータであるSpike(riscv-isa-sim)は内部にDevice Treeが格納されており、ブートする前に読み取ることができる。 このDevice Treeによってブート時にどのような処理がなされているのか、解析してみることにした。

Spikeは実行時に以下のオプションを追加することでDevice Treeの情報を出力することができる。

$ spike --dump-dts [実行ファイル]

ちなみにこの実行ファイルは必要ないと思うのだが、今の仕様上必要なので入れておく。

/dts-v1/;

/ {
  #address-cells = <2>;
  #size-cells = <2>;
  compatible = "ucbbar,spike-bare-dev";
  model = "ucbbar,spike-bare";
  cpus {
    #address-cells = <1>;
    #size-cells = <0>;
    timebase-frequency = <10000000>;
    CPU0: cpu@0 {
      device_type = "cpu";
      reg = <0>;
      status = "okay";
      compatible = "riscv";
      riscv,isa = "rv64imafdc";
      mmu-type = "riscv,sv48";
      clock-frequency = <1000000000>;
      CPU0_intc: interrupt-controller {
        #interrupt-cells = <1>;
        interrupt-controller;
        compatible = "riscv,cpu-intc";
      };
    };
  };
  memory@80000000 {
    device_type = "memory";
    reg = <0x0 0x80000000 0x0 0x80000000>;
  };
  soc {
    #address-cells = <2>;
    #size-cells = <2>;
    compatible = "ucbbar,spike-bare-soc", "simple-bus";
    ranges;
    clint@2000000 {
      compatible = "riscv,clint0";
      interrupts-extended = <&CPU0_intc 3 &CPU0_intc 7 >;
      reg = <0x0 0x2000000 0x0 0xc0000>;
    };
  };
  htif {
    compatible = "ucb,htif0";
  };
};

上記のDevice Tree情報はDevice Tree Compilerによりコンパイルされ、メモリ上に格納される。これをブート時に読み取るわけだ。

読み取るコードはriscv-pkから確認することができる。

  • riscv-pk/machine/minit.c
void init_first_hart(uintptr_t hartid, uintptr_t dtb)
{
  // Confirm console as early as possible
  query_uart(dtb);
  query_htif(dtb);

  hart_init();
  hls_init(0); // this might get called again from parse_config_string

  // Find the power button early as well so die() works
  query_finisher(dtb);

  query_mem(dtb);
  query_harts(dtb);
  query_clint(dtb);
  query_plic(dtb);

  wake_harts();

  plic_init();
  hart_plic_init();
  //prci_test();
  memory_init();
  boot_loader(dtb);
}

ここでは、dtsにも出てきているhtifの場所を解析するシーケンスについてみていく。query_htif()は以下のようになっている。

  • riscv-pk/machine/uart.c
void query_htif(uintptr_t fdt)
{
  struct fdt_cb cb;
  struct htif_scan scan;

  memset(&cb, 0, sizeof(cb));
  cb.open = htif_open;
  cb.prop = htif_prop;
  cb.done = htif_done;
  cb.extra = &scan;

  fdt_scan(fdt, &cb);
}

ここでcbというのはコールバック関数のことだが

  • cb.open() : htifの情報を読み取るための下準備(初期化)
  • cb.prop() : Propertyを読み取るコード(本体)
  • cb.done() : htifの情報を読み終わった後の最終的な処理

に分けられている。cb.extraには、読み取った情報を格納しておくための領域として使用される。

fdt_scan()はDevice Tree BlobをParseするためのルーチンである。 これはDevice Treeのフォーマットに基づいてdtbをParseしていく。 そして、Propertyのツリーにたどり着くとcb.prop() (= htif_prop())を読み出してParseする。 (実際のParseについてはfdt_scan_helper()に実装があるが、ここでは本題ではないので割愛する)。

  • riscv-pk/machine/fdt.c
static uint32_t *fdt_scan_helper(
  uint32_t *lex,
  const char *strings,
  struct fdt_scan_node *node,
  const struct fdt_cb *cb)
{
  struct fdt_scan_node child;
  struct fdt_scan_prop prop;
  int last = 0;
...
      case FDT_PROP: {
        assert (!last);
...
        lex += 3 + (prop.len+3)/4;
        cb->prop(&prop, cb->extra);
        break;
      }

htif_prop()で実行されているのは、Propertyの内容を確認して、HTIFのデバイス情報と合致すればHTIFは存在するという情報を掴むだけである。

  • riscv-pk/machine/htif.c
static void htif_prop(const struct fdt_scan_prop *prop, void *extra)
{
  struct htif_scan *scan = (struct htif_scan *)extra;
  if (!strcmp(prop->name, "compatible") && !strcmp((const char*)prop->value, "ucb,htif0")) {
    scan->compat = 1;
  }
}

全てのParseが終了すると、htif_done()が実行される。もし(上記のhtif_prop()により)HTIFに関するデバイス情報を発見していれば グローバル変数htifに1が設定される。 これは実際のプログラム実行中にHTIFを使用できるか否かを確認するために使用される。

  • riscv-pk/machine/htif.c
static spinlock_t htif_lock = SPINLOCK_INIT;
uintptr_t htif;
...

static void htif_done(const struct fdt_scan_node *node, void *extra)
{
  struct htif_scan *scan = (struct htif_scan *)extra;
  if (!scan->compat) return;

  htif = 1;
}

ここで示したのはHTIFの簡単な例だが、それ以外の例えばUARTやメモリ領域の設定についても同様にDevice Treeの探索が実行され、アドレス情報などが設定用の変数に書き込まれて使用される。

ゼロから作る量子コンピュータ・シミュレータ(6. Deutsch-Jozsaアルゴリズムを実装する)

量子コンピュータ・シミュレータの実装、骨格となるところは完成しており、あとはいろいろなゲートを挿入するだけとなっている。

かなり進んできたので、少し気分転換にDeutsch-Jazsaのアルゴリズムを実装してみよう。 Deutsch-Jazsaのアルゴリズムについて詳細はここでは説明しないが、量子コンピュータにより古典コンピュータよりも高速に解くことができるといわれているアルゴリズムだ。

詳細なアルゴリズムは以下を参考にした。

http://quantum.classcat.com/category/quantum-algorithms/

この中では、関数 f(x)が与えられたときに、0か1から構成されるnビットの入力が与えられて、

  • 全ての入力ビットが0か1のどちらかである = constant
  • 入力ビットの半分が0、残り半分が1である = balanced

とし、そのどちらであるかどうかを判定するアルゴリズムだ。何に使うんだ?といえばそこまでだが、量子コンピュータの威力を見ることができる簡単な例なので追ってみることにする。

アルゴリズムの原理はさておき、以下のようにすれば f(x)の判定を行うことができるらしい。

  1. n量子ビットを0で初期化する。
  2. n量子ビットにHadamardゲートを適用する。
  3. n量子ビットに対象となる関数f(x)を適用する。
  4. n量子ビットにHadamardゲートを適用する。
  5. n量子ビットを測定する。

上記のウェブサイトで登場するのは  f(x)=x _ {0} \oplus x _ {1} x _ {2} という関数だ。これはBalancedとなる。これを自作量子コンピュータ・シミュレータで実装してみる。

真ん中に入っているのが関数f(x)を表現したオラクル回路U_{f}である。これに対して、前後に全ビットに対してHadamardゲートを適用する。

cat ../test/test_dj_balanced.py
qubits = [qbit.make_qbit() for i in range(3)]

for q in qubits:
    qbit.h(q)

qbit.h(qubits[2])
qbit.cnot(qubits[1], qubits[2])
qbit.h(qubits[2])
qbit.z(qubits[0])

for q in qubits:
    qbit.h(q)

result = [0 for i in range(8)]

for i in range(0, 1000):
    val = qbit.val(qubits[2]) << 2 | qbit.val(qubits[1]) << 1 | qbit.val(qubits[0])
    result[val] = result[val] + 1

print(result)

実行結果は以下のようになった。量子ビットが000になる確率は0であり、Balancedであることが分かる。

$ ./qbit_simulator ../test/test_dj_balanced.py
[0, 265, 0, 244, 0, 255, 0, 236]

一方でConstantの関数を入力してみる。というか今度は U_{f}=Iなだけである(たぶん)。

qubits = [qbit.make_qbit() for i in range(3)]

for q in qubits:
    qbit.h(q)

for q in qubits:
    qbit.h(q)

result = [0 for i in range(8)]

for i in range(0, 1000):
    val = qbit.val(qubits[2]) << 2 | qbit.val(qubits[1]) << 1 | qbit.val(qubits[0])
    result[val] = result[val] + 1

print(result)

実行結果は以下のようになった。量子ビットが000になる確率は100%であり、Constantであることが分かる。

$ ./qbit_simulator ../test/test_dj_constant.py
[1000, 0, 0, 0, 0, 0, 0, 0]
f:id:msyksphinz:20181003000801p:plain

ゼロから作る量子コンピュータ・シミュレータ(5. Pythonインタフェースから全加算器を構成してシミュレート)

量子コンピュータ・シミュレータの実装、骨格となるところは完成しており、あとはいろいろなゲートを挿入するだけとなっている。

まずはアダマールゲート(Hadamard Gate)を実装してみることにする。

Hadamard Gate(Hゲート)の役割

Hadamard Gateを 1量子ビットに適用するためには、以下の行列を適用する。

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

これだけなら非常に簡単なのだが、これを複数量子ビットに対して適用する場合にはどうしたら良いのだろう?

例えば、以下のように2量子ビットを定義して1ビットにだけHゲートを適用する場合はどのように計算すればよい?

qbit_a = qbit.make_qbit()
qbit_b = qbit.make_qbit()
qbit.h(qbit_a)
# qbit.h(qbit_b)

この場合には、2量子ビットが定義されているのでとりあえず行列の空間としては4要素があり得る。 1ビット目の量子ビットにHゲートが適用されるので、インデックスとしては00⇔01 / 10⇔11 が互いに影響する。 つまり、量子状態ベクトルの要素0と要素1、要素2と要素3の確率が互いに影響される。

一方で、量子ビットの2ビット目に対してHゲートを適用する場合には、インデックスとして00⇔10 / 01⇔11が互いに影響する。 つまり、量子状態ベクトルの要素0と要素2、要素1と要素3の確率が互いに影響される。

f:id:msyksphinz:20181001232715p:plain

これらを実現するためには、以下のように行列式を構成すればよいと思う。

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

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

自作量子コンピュータ・シミュレータへのHadamard Gateの実装

という訳で自作量子コンピュータ・シミュレータに実装してみた。とりあえず正しく実装できてるか確認したかったので非常に単純な実装にしてある。

行列の計算とかそのままだ。ライブラリとかちゃんと使え、って感じだが。。。

void Qbits::Hadamard(int idx)
{
  std::vector<std::vector<double>> matrix;
  matrix.resize(m_qbit_size);
  for (size_t i = 0; i < m_qbit_size; i++) {
    matrix[i].resize(m_qbit_size);
    for (size_t j = 0; j < m_qbit_size; j++) {
      matrix[i][j] = 0.0;
    }
  }
  uint64_t c_bitmask = (1 << idx);
  for (size_t idx = 0; idx < m_qbit_size; idx++) {
    if ((idx & c_bitmask) == c_bitmask) {
      matrix[idx][idx^c_bitmask] =  1.0 / sqrtf(2);
      matrix[idx][idx          ] = -1.0 / sqrtf(2);
    } else {
      matrix[idx][idx^c_bitmask] = 1.0 / sqrtf(2);
      matrix[idx][idx          ] = 1.0 / sqrtf(2);
    }
  }

  std::complex<double> *new_comp_amp = new std::complex<double>[m_qbit_size];
  for (size_t v_idx = 0; v_idx < m_qbit_size; v_idx++) {
    new_comp_amp[v_idx] = 0.0;
    for (size_t h_idx = 0; h_idx < m_qbit_size; h_idx++) {
      new_comp_amp[v_idx] += matrix[v_idx][h_idx] * m_comp_amp[h_idx];
    }
  }

  for (size_t idx = 0; idx < m_qbit_size; idx++) {
    m_comp_amp[idx] = new_comp_amp[idx];
  }
}

量子ビットへ適用される行列を表示してみると、以下のようになる。正しく構成できているようだ。

0.707107 0.707107 0.000000 0.000000
0.707107 -0.707107 0.000000 0.000000
0.000000 0.000000 0.707107 0.707107
0.000000 0.000000 0.707107 -0.707107

Hadamard Gateのテスト

という訳で、さっそくテストプログラムを動かしてみる。以下は量子ビットを2つ定義し、1つ目にのみHゲートを適用したあと結果を表示する。 Hゲートは量子ビットを重ね合わせの状態にするので、0と1の状態は1000回観測するとおおよそ半分になるはずだ。

qbit_a = qbit.make_qbit()
qbit_b = qbit.make_qbit()
qbit.h(qbit_a)

result_a = [0, 0]
result_b = [0, 0]

for i in range(0, 1000):
    val = qbit.val(qbit_a)
    result_a[val] = result_a[val] + 1
    val = qbit.val(qbit_b)
    result_b[val] = result_b[val] + 1

print(result_a)
print(result_b)

実行してみると、以下のようになった。qbit_aは重ね合わせの状態なので確率は0と1で半分ずつ。qbit_bは何もゲートを適用していないのでひたすら0のみである。正しく実装できているようだ。

[477, 523]
[1000, 0]

Googleの量子コンピュータフレームワークCirqを使ってみる(5. 資料を見ながらCirqの計算フローを追ってみる)

Cirqの計算プロセスの続きである。Cirqがどのようにして量子状態の計算を行っているのかを見ていく。

とりあえず、前回から少しコードを変えた。アダマールゲートを使用するように変更する。アダマールゲートというのは、量子を重ね合わせの状態に持っていくゲートのことで、具体的には以下のような行列を計算する。1

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

ただし、これはもう少し一般化できて、XZ軸の周りで回転させるときは下記のように表される(以下は https://www.appi.keio.ac.jp/Itoh_group/abe/pdf/ap2016_4.pdf より抜粋)。

f:id:msyksphinz:20181001011942p:plain
図. 一般的なアダマールゲートの表現(https://www.appi.keio.ac.jp/Itoh_group/abe/pdf/ap2016_4.pdf より抜粋)

さて、下記のようにCirqで記述した量子回路はどのようにして計算されるのだろうか。

import cirq

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

circuit = cirq.Circuit()
circuit.append(cirq.H.on(qbits[0]))

simulator = cirq.google.XmonSimulator()

results = simulator.run(circuit, repetitions=1000, qubit_order=qbits)

まずは回路の構成をprintを入れて見ていくと、以下のようになった。ExpWGateExpZGateというもので構成されているようだ。そんな気がする。

(Pdb) c
> /home/msyksphinz/work/Cirq/cirq/google/sim/xmon_simulator.py(535)_simulator_iterator()
-> if isinstance(gate, xmon_gates.ExpZGate):
(Pdb) p circuit
cirq.Circuit(moments=[
    cirq.Moment(operations=[
        ExpWGate(half_turns=-0.5, axis_half_turns=0.5).on(cirq.GridQubit(0, 0)),
    ]),
    cirq.Moment(operations=[
        ExpZGate(half_turns=1.0).on(cirq.GridQubit(0, 0)),
    ]),
])

PDBなどでPythonのコードを追いかけていくと、実際のシミュレーションはcirq/google/sim/xmon_stepper.pyで実装されている。 ゲートを適用しているのは、以下のコードのようだ。

  • cirq/google/sim/xmon_stepper.py
class ThreadlessPool(object):
    """A Pool that does not use any processes or threads.

    Only supports map, close, and join, the later two being trivial.
    No enforcement of closing or joining is done, so map can be called
    repeatedly.
    """

    # noinspection PyMethodMayBeStatic
    def map(self, func, iterable, chunksize=None):
        assert chunksize is None, 'Chunking not supported by SimplePool'
        return [func(x) for x in iterable]

それぞれfuncにはどのようなものが入っているのかを調査すると、以下のようになった。

func(x) = <function _w_within_shard at 0x7fbbbe0a62f0>
func(x) = <function _norm_squared at 0x7fbbbe0a6510>
func(x) = <function _renorm at 0x7fbbbe0a6598>
func(x) = <function _clear_scratch at 0x7fbbbe0a6048>
func(x) = <function _apply_scratch_as_phase at 0x7fbbbe0a6268>
func(x) = <function _clear_scratch at 0x7fbbbe0a6048>
func(x) = <function _single_qubit_accumulate_into_scratch at 0x7fbbbe0a60d0>
func(x) = <function _apply_scratch_as_phase at 0x7fbbbe0a6268>
func(x) = <function _clear_scratch at 0x7fbbbe0a6048>
func(x) = <function _apply_scratch_as_phase at 0x7fbbbe0a6268>
func(x) = <function _state_shard at 0x7fbbbe0a4c80>

例えば、Hゲートを実行しているのは_w_within_shardのような気がする。

  • `cirq/google/sim/xmon_stepper.py``
def _w_within_shard(args: Dict[str, Any]):
    """Applies a W gate when the gate acts only within a shard."""
    index = args['index']
    half_turns = args['half_turns']
    axis_half_turns = args['axis_half_turns']
    state = _state_shard(args)
    pm_vect = _pm_vects(args)[index]
    num_shard_qubits = args['num_shard_qubits']
    shard_size = 2 ** num_shard_qubits

    reshape_tuple = (2 ** (num_shard_qubits - 1 - index), 2, 2 ** index)
    perm_state = np.reshape(
        np.reshape(state, reshape_tuple)[:, ::-1, :], shard_size)
    cos = np.cos(-0.5 * np.pi * half_turns)
    sin = np.sin(-0.5 * np.pi * half_turns)

    cos_axis = np.cos(np.pi * axis_half_turns)
    sin_axis = np.sin(np.pi * axis_half_turns)

    new_state = cos * state + 1j * sin * perm_state * (
        cos_axis - 1j * sin_axis * pm_vect)
    np.copyto(state, new_state)

各種関数にprintを挟んで実行してくと、以下のようになった。

func(x) = <function _w_within_shard at 0x7fbbbe0a62f0>
<perm_state = [0.+0.j 1.+0.j], pm_state = [ 1 -1]
<              cos_axis - 1j * sin_axis * pm_vect  = [6.123234e-17-1.j 6.123234e-17+1.j]>
<perm_state * (cos_axis - 1j * sin_axis * pm_vect) = [0.000000e+00+0.j 6.123234e-17+1.j]>
half_turns = -0.5, state = [1.+0.j 0.+0.j], new_state = [ 0.70710677+0.00000000e+00j -0.70710677+4.32978021e-17j]
func(x) = <function _norm_squared at 0x7fbbbe0a6510>
func(x) = <function _renorm at 0x7fbbbe0a6598>
state = [ 0.70710677+0.0000000e+00j -0.70710677+4.3297803e-17j]
new_state = [ 0.7071068+0.000000e+00j -0.7071068+4.329781e-17j]
func(x) = <function _clear_scratch at 0x7fbbbe0a6048>
func(x) = <function _apply_scratch_as_phase at 0x7fbbbe0a6268>
state = [ 0.7071068+0.000000e+00j -0.7071068+4.329781e-17j], _scratch_shard(args) = [0.+0.j 0.+0.j]
new_state = [ 0.7071068+0.000000e+00j -0.7071068+4.329781e-17j]
func(x) = <function _clear_scratch at 0x7fbbbe0a6048>
func(x) = <function _single_qubit_accumulate_into_scratch at 0x7fbbbe0a60d0>
func(x) = <function _apply_scratch_as_phase at 0x7fbbbe0a6268>
state = [ 0.7071068+0.000000e+00j -0.7071068+4.329781e-17j], _scratch_shard(args) = [-1.+0.j  1.+0.j]
new_state = [-3.090862e-08-0.7071068j  3.090862e-08-0.7071068j]
func(x) = <function _clear_scratch at 0x7fbbbe0a6048>
func(x) = <function _apply_scratch_as_phase at 0x7fbbbe0a6268>
state = [-3.090862e-08-0.7071068j  3.090862e-08-0.7071068j], _scratch_shard(args) = [0.+0.j 0.+0.j]
new_state = [-3.090862e-08-0.7071068j  3.090862e-08-0.7071068j]
func(x) = <function _state_shard at 0x7fbbbe0a4c80>
Counter({1: 502, 0: 498})

最初のアダマールゲートの適用で、[1.+0.j 0.+0.j][ 0.70710677+0.00000000e+00j -0.70710677+4.32978021e-17j] に変換されていることが分かる。 そのあとで_apply_scratch_as_phaseにより[ 0.7071068+0.000000e+00j -0.7071068+4.329781e-17j] から [-3.090862e-08-0.7071068j 3.090862e-08-0.7071068j] になっているが良く分からない。

つまりこのあたりを見ていけばCirqの量子計算の核の部分を見ることができそうだ。

ちなみに、Cirqはやはり量子ビットの数を大きくするとその 2^{n}の行列を用意するようだ。 17ビットを超えるとそれ以上は_within_shardは呼ばれなくなった。

perm_state = [0.+0.j 1.+0.j]
pm_vect    = [ 1 -1]
state      = [1.+0.j 0.+0.j]
new_state  = [ 0.70710677+0.00000000e+00j -0.70710677+4.32978021e-17j]
perm_state = [0.+0.j 0.+0.j 1.+0.j 0.+0.j]
pm_vect    = [ 1  1 -1 -1]
state      = [1.+0.j 0.+0.j 0.+0.j 0.+0.j]
new_state  = [ 0.70710677+0.00000000e+00j  0.        +0.00000000e+00j
 -0.70710677+4.32978021e-17j  0.        +0.00000000e+00j]
perm_state = [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 1.+0.j 0.+0.j 0.+0.j 0.+0.j
 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
 0.+0.j 0.+0.j]
pm_vect    = [ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 -1 -1 -1 -1 -1 -1 -1 -1
 -1 -1 -1 -1 -1 -1 -1 -1]
state      = [1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
 0.+0.j 0.+0.j]
new_state  = [ 0.70710677+0.00000000e+00j  0.        +0.00000000e+00j
  0.        +0.00000000e+00j  0.        +0.00000000e+00j
  0.        +0.00000000e+00j  0.        +0.00000000e+00j
  0.        +0.00000000e+00j  0.        +0.00000000e+00j
  0.        +0.00000000e+00j  0.        +0.00000000e+00j
  0.        +0.00000000e+00j  0.        +0.00000000e+00j
  0.        +0.00000000e+00j  0.        +0.00000000e+00j
  0.        +0.00000000e+00j  0.        +0.00000000e+00j
 -0.70710677+4.32978021e-17j  0.        +0.00000000e+00j
  0.        +0.00000000e+00j  0.        +0.00000000e+00j
  0.        +0.00000000e+00j  0.        +0.00000000e+00j
  0.        +0.00000000e+00j  0.        +0.00000000e+00j
  0.        +0.00000000e+00j  0.        +0.00000000e+00j
  0.        +0.00000000e+00j  0.        +0.00000000e+00j
  0.        +0.00000000e+00j  0.        +0.00000000e+00j
  0.        +0.00000000e+00j  0.        +0.00000000e+00j]

ただし量子ビットを18ビットよりも大きくすると_w_within_shardは呼ばれなくなった。 別のアルゴリズムを使っているのだろうか。

Googleの量子コンピュータフレームワークCirqを使ってみる(4. Cirqの中身を追ってみる)

Google量子コンピュータフレームワークであるCirqを使っており、簡単な回路なら記述できるようになった。

次に疑問点として沸いてくるのは、Cirqはどのようにして量子状態のシミュレーションを行っているのかというところだ。 Pythonのコードを見るのは苦手だけれども、どうにかしてこれを追ってみることにした。

Pythonパッケージを変更して試行してい見る。

まずはCirqでシミュレーションを実行するにあたり以下の簡単な量子回路を用意した。 非常に簡単な回路だ。Xゲートを量子ビットに適用しているだけである。

  • cnot_test.py
import cirq

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

circuit = cirq.Circuit()
circuit.append(cirq.X.on(qbits[0]))

simulator = cirq.google.XmonSimulator()
circuit.append(cirq.measure(qbits[0], key='c'))

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

量子回路は以下のようになる。

(0, 0): ───X───M('c')───

いちいちパッケージをビルドしてインストールするのは面倒なので良い方法はないかと悩んでいたのだが、Twitterpython setup.py developでリンクを張られた状態でパッケージを使用できることが分かった。これは便利。

python setup.py develop

シミュレーションの様子を追いかけてみる

まずはcircuit.X.on()circuit.CNOT.on()がどこで実行されているのかを確認する。これはどうもXを実行するためのゲートオブジェクトを返しているだけのようだ。

  • cirq/ops/common_gates.py
CNOT = CNotGate()  # Controlled Not Gate.
  • cirq/ops/common_gates.py
class CNotGate(eigen_gate.EigenGate,
               gate_features.TextDiagrammable,
               gate_features.CompositeGate,
               gate_features.TwoQubitGate,
               gate_features.QasmConvertibleGate):
...
    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))
  • cirq/ops/raw_types.py
class Gate:
...
    def on(self, *qubits: QubitId) -> 'gate_operation.GateOperation':
        """Returns an application of this gate to the given qubits.

        Args:
            *qubits: The collection of qubits to potentially apply the gate to.
        """
        # Avoids circular import.
        from cirq.ops import gate_operation

        if len(qubits) == 0:
            raise ValueError(
                "Applied a gate to an empty set of qubits. Gate: {}".format(
                    repr(self)))
        self.validate_args(qubits)
        return gate_operation.GateOperation(self, list(qubits))
  • cirq/ops/gate_operation.py
class GateOperation(raw_types.Operation,
                    extension.PotentialImplementation[Union[
                        gate_features.BoundedEffect,
                        gate_features.CompositeOperation,
                        gate_features.ExtrapolatableEffect,
                        gate_features.ParameterizableEffect,
                        gate_features.PhaseableEffect,
                        gate_features.ReversibleEffect,
                        gate_features.TextDiagrammable,
                        gate_features.QasmConvertibleOperation,
                    ]]):
    """An application of a gate to a collection of qubits.

    Attributes:
        gate: The applied gate.
        qubits: A sequence of the qubits on which the gate is applied.
    """

    def __init__(self,
                 gate: raw_types.Gate,
                 qubits: Sequence[raw_types.QubitId]) -> None:
        self._gate = gate
        self._qubits = tuple(qubits)

次に、量子ビットを計測したときにどのように実行されるのかを観測した。

  • cirq/google/sim/xmon_stepper.py
class Stepper(object):
...
    def sample_measurements(
            self,
            indices: List[int],
            repetitions: int=1) -> List[List[bool]]:
...
        # Calculate probabilities and reshape to tensor of qubits.
        tensor = np.reshape(np.abs(self.current_state) ** 2,
                            self._num_qubits * [2])

        # Tensor axis order is reverse of index order, so we transpose here.
        tensor = np.transpose(tensor)

このtensorという変数が怪しい気がする。いろいろ量子回路を変えながらtensorの値を出力してみた。

まずはXゲート→CNOTゲートを当てたもの。

circuit = cirq.Circuit()
circuit.append(cirq.X.on(qbits[0]))
circuit.append(cirq.CNOT.on(qbits[0], qbits[1]))

[cirq.GridQubit(0, 0), cirq.GridQubit(0, 1)]
(0, 0): ───X───@───
               │
(0, 1): ───────X───

current_state = [ 2.0323548e-18-8.4182900e-19j  1.1845434e-17+2.8597407e-17j
 -2.2809713e-08+5.5067517e-08j  1.0000001e+00+0.0000000e+00j]

tensor = [[4.8391427e-36 9.5812607e-34]
 [3.5527145e-15 1.0000002e+00]]
Counter({1: 1})

次にCNOTゲートのみ当てたもの。

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

current_state = [-2.9802322e-08+1.0000000e+00j -5.7024274e-09-1.3766877e-08j
 -2.8597404e-17-1.1845433e-17j -8.4182890e-19+2.0323546e-18j]

[cirq.GridQubit(0, 0), cirq.GridQubit(0, 1)]
(0, 0): ───@───
           │
(0, 1): ───X───
tensor = [[1.000000e+00 2.220446e-16]
 [9.581259e-34 4.839141e-36]]
Counter({0: 1})

やはりtensorの値が変わったので、この変数が量子ビットの状態を保持していそうな気がする。

ただし具体的にどうやって量子ゲートの計算をしているのかについてはたどり着かなかった。 もう少し掘り下げてみていきたい。

f:id:msyksphinz:20180923015903p:plain

ゼロから作る量子コンピュータ・シミュレータ(5. Pythonインタフェースから全加算器を構成してシミュレート)

フルスクラッチ量子コンピュータ・シミュレータを作る話、骨格となる部分はおおよそ完成し、量子ビットを取り扱う方法を考える。

量子ビットを作成する方法の検討

量子ビットは、1量子ビットPythonの1変数として宣言できるようにしたい。 例えば、以下のようにPythonで記述すると2つの量子ビットが個別に宣言される。

qbit1 = qbit()
qbit2 = qbit()

で、この状態だと2量子ビット存在するので、状態ベクトルとしては 2^{2}のベクトルを内部に保持する必要がある。 さらに量子ビットを追加するとベクトル長は8にまで拡大するので、そのたびに状態ベクトルを拡大する必要がある。

いろいろ考えて、とりあえず今はグローバルに量子ビットを管理するクラスを1つだけ用意し、量子ビットが追加される毎に状態ベクトルを再定義することにした。 そして、量子ビットそのものを示す変数(上記だとqbit1, qbit2)は単なる整数として表現することにした。

PyMethodDef qbit_methods[] = {
...
  { "make_qbit" , (PyCFunction)MakeQbits, METH_VARARGS, "Make Qbits"                     },
...
PyObject* MakeQbits (PyObject* self, PyObject* args)
{
  int new_qbit_idx = g_qbits.AddQbit();
  PyObject *ret = Py_BuildValue("I", new_qbit_idx);
  return ret;
}

量子ビットの操作

例えば、2つの量子ビットに対してCNOTゲートを適用する場合、以下のように記述するイメージだ。

cnot(qbit1, qbit2)

結局のところはqbit1, qbit2は整数で割り当てられたインデックスなので、インデックスに基づいて内部で量子状態を変えているに過ぎない。

PyMethodDef qbit_methods[] = {
...
  { "cnot"      , (PyCFunction)CNot     , METH_VARARGS, "Get Applied CNot"               },
...
PyObject* CNot (PyObject* self, PyObject* args)
{
  int a_idx, b_idx;
  if (!PyArg_ParseTuple(args, "ii", &a_idx, &b_idx))
    return NULL;
  g_qbits.CNot (a_idx, b_idx);

  return Py_BuildValue("i", 0);
}

g_qbits.CNot (a_idx, b_idx)は、行列ベクトルの操作する対象の量子ビット2つ(結局はインデックス)指定して量子状態を変える操作を実行する。

Pythonインタフェースを使って全加算器を作る

という訳で、Pythonインタフェースを整備したので、これを使って全加算器を作る。

量子ビットとしては、入力ビットとしてqbit_a, qbit_b 答えとしてqbit_c0, qbit_c1を用意した。

  • test/test_fulladder.py
qbit_a  = qbit.make_qbit()
qbit_b  = qbit.make_qbit()
qbit_c0 = qbit.make_qbit()
qbit_c1 = qbit.make_qbit()

for b in range (2):
    for a in range (2):
        if a != qbit.val(qbit_a):
            qbit.x(qbit_a)
        if b != qbit.val(qbit_b):
            qbit.x(qbit_b)

        if qbit.val(qbit_c0) == 1:
            qbit.x(qbit_c0)

        print('a + b = s1,s0 => %d + %d = ' % (qbit.val(qbit_a), qbit.val(qbit_b)), end="")

        qbit.ccnot(qbit_a, qbit_b, qbit_c1);
        qbit.cnot (qbit_a, qbit_b);
        qbit.ccnot(qbit_b, qbit_c0, qbit_c1);
        qbit.cnot (qbit_b, qbit_c0);

        print('%d,%d' % (qbit.val(qbit_c1), qbit.val(qbit_c0)))

実際に全加算器を実行しているのは、後半のccnotから始まる4行に過ぎないのだが、それ以外のところでは量子ビットの初期化を行っている。

a=0, b=0 ... a=1, b=1の順番に4種類の加算を実行しているのだが、その度に量子の状態を初期化している。

さっそく実行してみた。正しく計算できている。問題なさそうだ。

$ ./qbit_simulator ../test/test_fulladder.py
a + b = s1,s0 => 0 + 0 = 0,0
a + b = s1,s0 => 1 + 0 = 0,1
a + b = s1,s0 => 0 + 1 = 0,1
a + b = s1,s0 => 1 + 1 = 1,0
f:id:msyksphinz:20180927230554p:plain

ゼロから作る量子コンピュータ・シミュレータ(4. シミュレータにPythonのインタフェースを追加する)

フルスクラッチ量子コンピュータ・シミュレータを作る話、骨格となる部分はおおよそ完成してきたので次はインタフェースの部分を調整する。

現在は量子回路を生成する部分はC++で記述しており、量子回路を書き換えるたびにリコンパイルが発生して格好悪い。

void Adder3Bit ()
{
      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);
}

そこで、Pythonのインタフェースを接続し(私の十八番)、Pythonで回路を記述するとC++にオフロードされて量子回路が作成されるようにしよう。

Python 3.6をC++と接続する

まずコンパイルの方法で戸惑った。基本的にC++側からPythonを呼び出すのは以下のように記述する。 これだけでPythonコマンドプロンプトを立ち上げることができる。

int main (int argc, char **argv)
{
  wchar_t *_program_name = Py_DecodeLocale(argv[0], NULL);
  wchar_t **_argv = nstrws_array(argc, argv);

  Py_SetProgramName(_program_name);
  Py_Initialize();

  Py_Main(argc, _argv);
  Py_Finalize ();

  return 0;
}

次に、コンパイルの方法だがかなり戸惑った。基本的には、コンパイラに渡すオプションは以下のようにして生成するのが正統派のようだ。

CFLAGS += $(shell python3-config --cflags)
LDFLAGS += $(shell python3-config --ldflags)

これはつまりpython3-config --cflagsというコマンドを実行すればコンパイル時とリンク時に必要なオプションが文字列で返されるので、それをg++に渡すのが最も簡単だということ。これを採用した。

さらに最終的なリンクの方法だが、リンクするオブジェクトより前に-lpythonを書いてはならないらしい。つまり、

$(TARGET): $(OBJS) $(HEADER)
        g++ -o $@ $^ $(LDFLAGS)

ならコンパイルできるが、

$(TARGET): $(OBJS) $(HEADER)
        g++ $(LDFLAGS) -o $@ $^

とするとエラーとなる。なんでだ!

g++ -o -Wall --std=c++14 -L/usr/lib/python3.6/config-3.6m-x86_64-linux-gnu -L/usr/lib -lpython3.6m -lpthread -ldl  -lutil -lm  -Xlinker -export-dynamic -Wl,-O1 -Wl,-Bsymbolic-functions qbit_simulator ../src/qbit.o ../src/qbit_simulator.o ../src/util.o ../src/py_env.o ../src/qbit.hpp ../src/util.hpp ../src/py_env.hpp
/usr/bin/x86_64-linux-gnu-ld: qbit_simulator: undefined reference to symbol 'PyModule_GetState'
/usr/lib/python3.6/config-3.6m-x86_64-linux-gnu/libpython3.6m.so: error adding symbols: DSO missing from command line

C++で作ったQbitsクラスをPython用のWrapperで囲む

C++で開発したQbitsクラスは、複数量子ビットそのものを表現しており、このオブジェクトをPython上でやり取りできるようにした方がよさそうだ。 そこで、Python側ではQbitsを扱う一連の関数群をqbitsパッケージとしてまとめることにする。 このパッケージは、アプリケーションの起動時に読み込むこともできるし、パッケージとしてPython側に埋め込むこともできることにする。

int main (int argc, char **argv)
{
  wchar_t *_program_name = Py_DecodeLocale(argv[0], NULL);
  wchar_t **_argv = nstrws_array(argc, argv);

  /* Add a built-in module, before Py_Initialize */
  PyImport_AppendInittab("qbit", InitPyEnv);

  Py_SetProgramName(_program_name);
  Py_Initialize();

  /* Optionally import the module; alternatively,
     import can be deferred until the embedded script
     imports it. */
  PyImport_ImportModule("qbit");

  PyRun_SimpleString ("import qbit");
  Py_Main(argc, _argv);
  Py_Finalize ();

  return 0;
}

パッケージ内で使える関数をいくつか定義し、それをC++側でも実装していく。

PyMethodDef qbit_methods[] = {
  { "py_add",             (PyCFunction)HelloAdd,           METH_VARARGS, "Example ADD"                    },
  { "qbits" ,             (PyCFunction)MakeQbits,          METH_VARARGS, "Make Qbits"                     },
  { "get" ,               (PyCFunction)GetByIdx,           METH_VARARGS, "Get Qbits by Index"             },
  { "cnot" ,              (PyCFunction)CNot,               METH_VARARGS, "Get Applied CNot"               },
  { "ccnot" ,             (PyCFunction)CCNot,              METH_VARARGS, "Get Applied CCNot"              },
  { "val" ,               (PyCFunction)Measure,            METH_VARARGS, "Get Qbits value"                },
  { NULL,                 NULL,                            0,            NULL                             } /* Sentinel */
};

static int qbit_traverse(PyObject *m, visitproc visit, void *arg) {
  Py_VISIT(GETSTATE(m)->error);
  return 0;
}

static int qbit_clear(PyObject *m) {
  Py_CLEAR(GETSTATE(m)->error);
  return 0;
}

static struct PyModuleDef moduledef = {
  PyModuleDef_HEAD_INIT,
  "qbit",
  NULL,
  sizeof(struct module_state),
  qbit_methods,
  NULL,
  qbit_traverse,
  qbit_clear,
  NULL
};

このqbitsパッケージを生成するためのモジュールを作成する関数は以下で定義した。

PyMODINIT_FUNC InitPyEnv (void)
{
  PyObject *module = PyModule_Create (&moduledef);

  if (module == NULL) {
    return NULL;
  }

  struct module_state *st = GETSTATE(module);

  st->error = PyErr_NewException("qbit.Error", NULL, NULL);
  if (st->error == NULL) {
    Py_DECREF(module);
    return NULL;
  }

  return module;
}

ここまでで一応シミュレータ上で立ち上がったPythonプロンプトて量子ビットを作れるようになった。

次は前に作った量子加算回路をPythonで実現させたい。

f:id:msyksphinz:20180922165310p:plain