FPGA開発日記

カテゴリ別記事インデックス https://msyksphinz.github.io/github_pages , English Version https://fpgadevdiary.hatenadiary.com/

Rocket-Chip実装のリグレッションスイートを試す

RISC-V

全く資料が無いのだが、RISC-V実装であるRocket-Chipでは、リグレッションスイートが格納されているregression/というディレクトリが存在する。

make emulator-regression-tests SUITE=UnittestSuite

まずはAMBAのテストを実行してみるためにいくつかテストを削減する。 AMBAUnitTestConfigのみを実行する。

diff --git a/regression/Makefile b/regression/Makefile
index 7953daa6..43b1e85c 100644
--- a/regression/Makefile
+++ b/regression/Makefile
@@ -58,7 +58,8 @@ endif

 ifeq ($(SUITE),UnittestSuite)
 PROJECT=freechips.rocketchip.unittest
-CONFIGS=AMBAUnitTestConfig TLSimpleUnitTestConfig TLWidthUnitTestConfig
+CONFIGS=AMBAUnitTestConfig
 endif

この辺について、Chiselでどのようにテストを構成すればいいのか良く分かっていないのだが、TestHarnessにAMBAの実装を当てはめて、テストのコンフィグレーションを挿入するらしい。

  • rocket-chip/src/main/scala/unittest/Configs.scala
class AMBAUnitTestConfig extends Config(new WithAMBAUnitTests ++ new WithTestDuration(10) ++ new BaseSubsystemConfig)

WithAMBAUnitTestsの実装は以下のようになっている。

class WithAMBAUnitTests extends Config((site, here, up) => {
  case UnitTests => (q: Parameters) => {
    implicit val p = q
    val txns = 100 * site(TestDurationMultiplier)
    val timeout = 50000 * site(TestDurationMultiplier)
    Seq(
      Module(new AHBBridgeTest(true, txns=8*txns, timeout=timeout)),
      Module(new AHBNativeTest(true, txns=6*txns, timeout=timeout)),
      Module(new AHBNativeTest(false,txns=6*txns, timeout=timeout)),
      Module(new APBBridgeTest(true, txns=6*txns, timeout=timeout)),
      Module(new APBBridgeTest(false,txns=6*txns, timeout=timeout)),
      Module(new AXI4LiteFuzzRAMTest(txns=6*txns, timeout=timeout)),
      Module(new AXI4FullFuzzRAMTest(txns=3*txns, timeout=timeout)),
      Module(new AXI4BridgeTest(     txns=3*txns, timeout=timeout)),
      Module(new AXI4XbarTest(       txns=1*txns, timeout=timeout)),
      Module(new AXI4RAMAsyncCrossingTest(txns=3*txns, timeout=timeout))) }
})

例えば AXI4BridgeTestについて実装を見てみた。

  • rocket-chip/src/main/scala/amba/axi4/Test.scala
class AXI4BridgeTest(txns: Int = 5000, timeout: Int = 500000)(implicit p: Parameters) extends UnitTest(timeout) {
  val dut = Module(LazyModule(new AXI4FuzzBridge(txns)).module)
  io.finished := dut.io.finished
}
...
class AXI4FuzzBridge(txns: Int)(implicit p: Parameters) extends LazyModule
{
  val master = LazyModule(new AXI4FuzzMaster(txns))
  val slave  = LazyModule(new AXI4FuzzSlave)

  slave.node := master.node

  lazy val module = new LazyModuleImp(this) with UnitTestModule {
    io.finished := master.module.io.finished
  }
}

以下をmakeを実行するとgenerated-srcディレクトリに生成されたVerilogが格納された。

cd regression
make emulator-regression-tests SUITE=UnittestSuite
cd ../emulator/generated-src
$ ls -1
freechips.rocketchip.unittest.AMBAUnitTestConfig
freechips.rocketchip.unittest.AMBAUnitTestConfig.anno.json
freechips.rocketchip.unittest.AMBAUnitTestConfig.behav_srams.v
freechips.rocketchip.unittest.AMBAUnitTestConfig.conf
freechips.rocketchip.unittest.AMBAUnitTestConfig.d
freechips.rocketchip.unittest.AMBAUnitTestConfig.fir
freechips.rocketchip.unittest.AMBAUnitTestConfig.plusArgs
freechips.rocketchip.unittest.AMBAUnitTestConfig.v

これからどのようにしてテストスイートを実行すればよいのだろう？

色々解析した結果、全てのテストが同時に動作してテストを実施しているらしい。 つまり、上記の - AHBBridgeTest - AHBNativeTest - AHBNativeTest - APBBridgeTest - APBBridgeTest - AXI4LiteFuzzRAMTest - AXI4FullFuzzRAMTest - AXI4BridgeTest - AXI4XbarTest - AXI4RAMAsyncCrossingTest

それぞれが勝手にMasterとSlaveを作ってテストを実行しているっていうことか。

波形を開いてみるとこのようになった。

f:id:msyksphinz:20181008145559p:plain

それぞれで独立してテストを実施している。

f:id:msyksphinz:20181008145947p:plain

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

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

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

  • Cir: Tutorial

Tutorial — Cirq 0.1 documentation

一度チュートリアルを読んだのだが、意味が分からずに諦めてしまったのであった。リトライする。

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

固有値 - Wikipedia

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

つまり、ランダムなので最初の量子回路の形は毎回異なるらしい。量子ビットの位相回転のパラメータを自由に設定しているということになる。

チュートリアルでは、3つの量子ビットの回転角度をそれぞれ$\alpha, \beta, \gamma$としている。 $\alpha=0.1, \beta=0.2, \gamma=0.3$として設定したときの量子回路は以下のようになった。

(0, 0): ───X^0.1───────────@───────────────────────────────────────────────────────────────────────────────────────────M('x')───
                           │                                                                                           │
(0, 1): ───X^0.1───────────┼───────X───@───────X───────────────────────────────@───────────────────────────────────────M────────
                           │           │                                       │                                       │
(0, 2): ───X^0.1───────────┼───────────┼───────────@───────────────────────────@^0.3───────────────────────────────────M────────
                           │           │           │                                                                   │
(1, 0): ───X^0.1───Z^0.2───@^0.3───────┼───────────┼───────@───────────────────────────────────@───────────────────────M────────
                                       │           │       │                                   │                       │
(1, 1): ───X^0.1───────────X───────────@^0.3───X───┼───────┼───────X───@───────X───────────────@^0.3───@───────────────M────────
                                                   │       │           │                               │               │
(1, 2): ───X^0.1───────────────────────────────────@^0.3───┼───────────┼───────────────@───────────────@^0.3───────────M────────
                                                           │           │               │                               │
(2, 0): ───X^0.1───Z^0.2───────────────────────────────────@^0.3───────┼───────────────┼───────────────@───────────────M────────
                                                                       │               │               │               │
(2, 1): ───X^0.1───Z^0.2───────────────────────────────────────────X───@^0.3───X───────┼───────────────@^0.3───@───────M────────
                                                                                       │                       │       │
(2, 2): ───X^0.1───Z^0.2───────────────────────────────────────────────────────────────@^0.3───────────────────@^0.3───M────────

この時にシミュレーションを行った結果、最終的に9つの量子ビットはどのような状態になりうるだろうか。 それをシミュレーションしたのが以下の結果である。 100回シミュレーションを繰り返し、量子ビットの観測結果をヒストグラムで表示してある。

...
circuit.append(one_step(h, jr, jc, 0.1, 0.2, 0.3))
circuit.append(cirq.measure(*qubits, key='x'))
print(circuit)
results = simulator.run(circuit, repetitions=100, qubit_order=qubits)
...

実行結果。

$ python3 cirq_tutorial.py
Counter({0: 83, 4: 4, 32: 3, 256: 2, 64: 2, 128: 2, 8: 1, 2: 1, 144: 1, 320: 1})

この時の量子状態において、エネルギー関数を計算する。 このためにenergy_func()という関数を定義する。

def energy_func(length, h, jr, jc):
    def energy(measurements):
        meas_list_of_lists = [measurements[i:i + length] for i in range(length)]
        pm_meas = 1 - 2 * np.array(meas_list_of_lists).astype(np.int32)
        tot_energy = np.sum(pm_meas * h)
        for i, jr_row in enumerate(jr):
            for j, jr_ij in enumerate(jr_row):
                tot_energy += jr_ij * pm_meas[i, j] * pm_meas[i + 1, j]
        for i, jc_row in enumerate(jc):
            for j, jc_ij in enumerate(jc_row):
                tot_energy += jc_ij * pm_meas[i, j] * pm_meas[i, j + 1]
        return tot_energy
    return energy
def obj_func(result):
    energy_hist = result.histogram(key='x', fold_func=energy_func(3, h, jr, jc))
    return np.sum(k * v for k,v in energy_hist.items()) / result.repetitions

この時のエネルギー関数は以下のようになった。

$ python3 cirq_tutorial.py
[cirq.GridQubit(0, 0), cirq.GridQubit(0, 1), cirq.GridQubit(0, 2), cirq.GridQubit(1, 0), cirq.GridQubit(1, 1), cirq.GridQubit(1, 2), cirq.GridQubit(2, 0), cirq.GridQubit(2, 1), cirq.GridQubit(2, 2)]
Counter({0: 81, 128: 5, 2: 3, 64: 3, 16: 2, 1: 1, 320: 1, 32: 1, 264: 1, 8: 1, 4: 1})
Value of the objective function -0.72

最後に、この3つのパラメータを0.0から1.0まで変化させていき、その時の最小値を求める。このためにはrun_sweep()関数を用いて複数の条件を一気にシミュレーションする。

circuit = cirq.Circuit()
alpha = cirq.Symbol('alpha')
beta = cirq.Symbol('beta')
gamma = cirq.Symbol('gamma')
circuit.append(one_step(h, jr, jc, 0.1, 0.2, 0.3))
circuit.append(cirq.measure(*qubits, key='x'))
results = simulator.run(circuit, repetitions=100, qubit_order=qubits)
print(results.histogram(key='x'))
print('Value of the objective function {}'.format(obj_func(results)))

$ python3 cirq_tutorial.py
...
OrderedDict([('alpha', 1.0), ('beta', 1.0), ('gamma', 0.3333333333333333)]) -4.64
OrderedDict([('alpha', 1.0), ('beta', 1.0), ('gamma', 0.4444444444444444)]) -4.58
OrderedDict([('alpha', 1.0), ('beta', 1.0), ('gamma', 0.5555555555555556)]) -4.82
OrderedDict([('alpha', 1.0), ('beta', 1.0), ('gamma', 0.6666666666666666)]) -4.68
OrderedDict([('alpha', 1.0), ('beta', 1.0), ('gamma', 0.7777777777777778)]) -4.44
OrderedDict([('alpha', 1.0), ('beta', 1.0), ('gamma', 0.8888888888888888)]) -4.7
OrderedDict([('alpha', 1.0), ('beta', 1.0), ('gamma', 1.0)]) -4.5
Minimum objective value is -4.96.

"GAFA 四騎士が創り変えた世界"を読んだ

原題は"The Four"という。面白そうなので読んでみた。

f:id:msyksphinz:20181005225632p:plain

ご存知の通りGAFAというのは

という世界を支配する4社のことであるが、この4社がどのようにして優位性を持つに至ったのか、どのようにして世界を支配するようになったのか、逆に言うと既存の企業はどのようにしてGAFAに自分のシェアを奪われたのか、いろいろな目線で解決してある。

GAFAはどのようにして世界を支配したのか

上記の写真は日本語版には掲載されておらず(かなり挑発的だから？)、英語の原本にのみ掲載されているのだが、4社はそれぞれ人間の働きかけてくる部位が異なる。

  • Google : 人間の知識欲に働きかける
  • Facebook : 人間の心に働きかける
  • Apple : 人間の性器に働きかける
  • (上記の写真ではAmazonが抜けているが、別の資料ではAmazonは人間の胃(つまり物欲、食欲)に働きかけるとあった)。

読みながら少し笑ってしまったのだが、特にAppleに対する言い方がヒドい(まあ私もAppleはそこまで好きではないが)。

Appleの製品は高級ブランドと一緒なのだそうだ。つまりAppleの製品を身に着けていれば、高級品を購入するだけの財力があり生き残れる可能性が高い。 異性に対して魅力的に映るための道具となるそうだ。

GoogleやFacebookは広告により成果を伸ばした。どちらも広告により勢力を伸ばしていったのだが、2社には違いがある。Googleは個人情報を(建前上)マーケティングなどに活用することはできないが、一方でFacebookは自分の行動と個人情報がしっかり関連づいている。 Facebookの方がより広告の作用ととしては大きくなる。

一方でGoogleは検索という万能ツールと、膨大なデータを使って未来を予見できるようになった。 また、Googleで検索しなければコンテンツが見つからない、という自社の圧倒的な優位性を利用し、他社の勢力を弱めることに成功した(本文ではニューヨークタイムズが自社の記事をGoogleにより無料でアクセスできるようにしたせいで自社の価値を大きく落としたという実例が挙げられている)。

Amazonはこれまでの小売業からシェアを完全に奪ってしまった。 Amazonは、人間の本能的な物欲(収集癖)を刺激するようにできている。購入するためのステップをなるべく減らし、そのためのインフラを徹底手系に効率化する。 Amazon Goなどの取り組みにより、なるべく面倒なことをせずに購入することができるようにステップを簡略化し、消費者の物欲を刺激している。

後半では、Next GAFAになれる可能性のある会社についていくつか紹介がなされている。どの会社も勢いがあるが、それぞれに欠点もある(特にUber社の扱いがヒドイ)。

次のGAFAになれずとも、GAFAの農奴にならないためにはどうするべきか

  • 大学に行くこと。良い大学に行くことは、素晴らしい環境とレベルの高い仲間を手に入れられる。また、学歴は現在でも非常に重要な要素である。
  • 何かを成し遂げること。何かを成し遂げるためには根性が必要だ。これは社会で生き残るために必要である。
  • 自分のキャリアをよく見せること。Linkedinなどでも良いし、もし自分の名前を検索してよくない結果が出れば、自分をもっとよく見せるための努力をすべきだ。

SiFiveはクラウド上でRISC-Vチップを設計する

RISC-V
f:id:msyksphinz:20181004233231p:plain
図. RISC-Vチップ設計を手掛けるSiFiveのシリコンクラウドサービスについて (https://content.riscv.org/wp-content/uploads/2018/07/1330-18.07.18-Huzefa-Cutlerywala-Open-Silicon.pdf より抜粋)

本日のTSMCのアナウンス記事、面白いと思ったのでメモ。

前々から少し聞いていたけれども、SiFiveはRISC-Vチップの設計環境をローカルの自前環境からクラウドへ移している。 TSMCのオープンエコシステムOIP VDE(Oven Intonation Platform Virtual Design Environment)を使って64-bit RISC-VマルチコアチップFreedom Unleashed 540をテープアウトしたちう話だ。 これはTSMCのこのサービスを使った最初のテープアウト事例となる。

www.digitimes.com

TSMCも公式サイトでアナウンスを出した。

TSMC Announces OIP Ecosystem Enabled in the Cloud

このTSMCのアライアンスにはEDAベンダのSynopsys, Cadenceだけでなく、クラウド環境を提供するAmazon, Microsoftも加入している。

MicrosoftとCadenceはSiFiveと手を組み、TSMCのOIP VDEを使った最初のフルデザインSoCをテープアウトした。このチップには64-bitマルチコアRISC-V CPUであるFreedom Unleashed 540が含まれている。 シノプシスは7nmプロセスを使用してPCI Express 5.0のDesignWare PHY IPをクラウドソリューションとしてテープアウトした。

一方でRISC-Vに続いてArmのこの流れに乗っていく。同様にArmもクラウド上でIPをライセンスし、Armエコシステムを構築するようだ。

これらの流れにより、これまで非常に高価だったEDAツールのコストや、設計に必要な手間が削減される方向に向かえば良いのだが。

クラウド化により参入障壁が下がることで、より半導体業界が発展することに期待したい。

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]