読者です 読者をやめる 読者になる 読者になる

VimプラグインのTravis CIテストを複数のVimのバージョンで動かそう

Vimプラグインにテストがあるのはあたりまえ。 そういう空気になってきたのはここ3年くらいのことでしょうか。

私自身、昔はあまりテスト文化に慣れておらず、「Vimプラグインみたいな小さなスクリプトにテストなんているのか?自分のプラグインは普段から使う、バグっていたらすぐ気がつくからテストなんていらないでしょ」と思っていました。 しかし、そういうテストのない自作プラグインがどんどん増えていき、3年4年と経ってしまうと自分のプラグインのコードを触りにくくなってきました。 昔はあまりVimプラグインの書き方に慣れていなかったので、酷いコードが絡み合っているのだけど、普段使う分には普通に便利なプラグイン、しかしリファクタリングしにくいというのがいくつかある状態です。

やっぱり、Vimプラグインもテストを書いたほうがよいですね。

id:thincaさんのvim-themisは、Vimプラグインのテストフレームワークとして広く使われています。 よくできていて、とても気に入っています。 github.com テストの書き方などはvim-themisのドキュメントや、これを使っている他のVimプラグインなどを参照されるとよいでしょう。

ところで、これを読まれているあなたが書かれているVimプラグインは、どのVimのバージョンで動きますか。 Vim 7.4では動きますか?Vim 7.3もサポートしてますか?それともVim 8.0のみのサポートですか?

これはVimプラグイン作者が決めることです。 7.3もサポートしてあげたいと思えばサポートしてあげて下さい。8.0以降しかサポートしないなら、そう書いておきましょう。 これは作者のポリシーなので、8.0以降しかサポートしないのはどうこうという風には思いません。

私はlightline.vimというプラグインを作ってメンテナンスしているのですが、かつてmap()v:keyを使うとVim 7.2で動かないという問い合わせが来たことがありました。 一瞬「そんな古いVimでこれまで動いていたのか」と思いましたが、せっかく自分のプラグインを便利と思って使っていただいている方なので、7.2でも動くようにコードを修正しました。

普段の開発は、新しいバージョンのVimを使っています。 私は毎日Vimをビルドし、最新のバージョンを使うようにしています。 新しいものを使っているとたまにバグにも遭遇しますし、新しい機能もすぐに知ることができます。 しかし、新しい機能ばかり見ていると、古いバージョンで動かない変更をプラグインに入れてしまうかもしれません。 動作確認用に古いバージョンのバイナリーはストックしていますが、手で古いバージョンでのテストを回すのはあまりに面倒です (実はそういうシェルスクリプトは手元にあるのですが…)。

そこで、Travis CIでVimの複数のバージョンでテストを動かそうという発想になるわけです。 lightline.vimも一年くらい前から複数バージョンでのテストを行ってはいましたが、改めて設定を見直してみたところ、次のような感じになりました。

language: generic

sudo: false

install:
  - git clone --depth=1 https://github.com/thinca/vim-themis /tmp/themis
  - (if ! test -d $HOME/vim-$VIM_VERSION/bin; then
       git clone https://github.com/vim/vim $HOME/vim &&
       cd $HOME/vim &&
       git checkout v$VIM_VERSION &&
       ./configure --prefix=$HOME/vim-$VIM_VERSION &&
       make &&
       make install;
     fi)

cache:
  directories:
    - $HOME/vim-$VIM_VERSION

env:
  - VIM_VERSION=8.0.0000
  - VIM_VERSION=7.4
  - VIM_VERSION=7.3

script:
  - export PATH=$HOME/vim-$VIM_VERSION/bin:$PATH
  - vim --version
  - /tmp/themis/bin/themis --reporter spec

VIM_VERSIONの部分は、対応したいバージョンを入れておくと良いでしょう (タグと対応しているのでgit tagで確認)。 かつてはTravisへの知識がなくて、複数のVimのバージョンでのテストを直列にしていたのですが、Environmentを使うと並列になってハッピーになりました。 Travis CIのキャッシュを使うことで、毎回Vimをビルドすることもcloneすることもありません。 Vimリポジトリをgit cloneするだけでも10秒くらいかかっていたので、かなり速くなりました。 あとcache.directoriesの中で環境変数使えるんですね… 便利…

ここまでダラダラと書いてきたけど、あれ実はみんなやってるのかなと思ってググったらaiya000さんの記事が見つかりました。 qiita.com そして自分のidが最後に書いてあった… Twitterでやり取りしたような気がするけどもう覚えてなかった…

まとめ

Vimプラグインもテストを書きましょう。 どのバージョンをサポートするかをきちんと決めましょう。 古いバージョンでもテストを動かしましょう。

みたいな感じ!じゃあね〜!

円周率の16進数表現100億桁目を求めてみた! ― 円周率の世界記録をどのように検証するか ―

あなたは円周率を何桁言えますか。3.14159…という、あの数字です。 円周率の小数部分は無限に続き、循環することもありません。 古来より、数学者は円周率の値を様々な幾何学的な近似や公式を用いて計算してきました。 その桁数は計算機の発明により飛躍的に伸び、収束の速い公式の発見や効率の良いアルゴリズムの発明などによって加速してきました *1

5年前、私がまだ学生だった頃、円周率1億桁の計算に挑んだことがありました。 私にとって高精度計算の初めての挑戦で、様々な試行錯誤で苦労したのをよく覚えています。

itchyny.hatenablog.com

2017年現在、円周率計算の世界記録は22兆桁です。 円周率計算の歴史をご覧いただくとよく分かると思いますが、近年の円周率計算の世界記録からは次のような特徴が読み取れます。

  • 2002年に1兆を超え、最新の記録 (2016年) は22兆桁 (10進数)
  • 企業ではなく個人によって記録が塗り替えられてきている
  • Alexander J. Yee氏によって作られたy-cruncherが多くの記録で使われている
  • ストレージは100TB単位
  • 計算時間は100日〜1年がかり

世界記録が22兆桁である一方で、5年前に私が計算したのはたったの1億桁です。 当時、私のブログのコードを少し改変して10億桁計算された方もいらっしゃいました (ありがとうございます) が、それでも世界記録とは1万倍の差があります。 そこには、かけた計算時間以上の技術の差があります。

  • 数字の桁数が多すぎてメモリーに乗らないので、ファイルに書きつつ計算する必要がある
  • 掛け算1つするにしても、桁を分割してファイルを読み書きしながら行う必要がある
  • そういうファイルに書きつつ計算するソフトウェアを自前で実装したり工夫する必要がある (億桁はGMPで計算できてしまった)
  • 世界記録の桁は他の人が計算したことがない値なので、計算結果を検証する方法を工夫する必要がある (私はオンラインサービスの結果と突き合わせてしまった)

本記事では、4番目の、円周率計算の結果の検証について説明してみようと思います。 円周率計算の記録は、どうやって検証するのでしょうか。 渋いテーマですが、世界記録をどう検証するかには興味がありますし、5年前の自分がスルーしてしまった部分でした。 当時の私は、オンラインで任意の桁を表示できるページで答えを見て済ませてしまっていました。 それができたのは、1億桁程度だからだったのと、たまたまそういうオンラインサービスがあったからです。 「人類未踏の値」の検証の裏には、数学と計算機科学の織りなす美しい公式の発見がありました。

円周率1000000桁表

円周率1000000桁表

円周率の世界記録の検証

円周率の世界記録が正しいことを検証するにはどうすればいいのでしょうか。 求めた桁の値は、その計算をした人しか見たことがない値です。 別の公式を使って同じ桁数計算して、結果を突き合わせるという方法がありますが、これには膨大な時間がかかります。 20兆桁計算するのに100日ほどかかっているので、別のアルゴリズムで全く同じ桁数分計算するにはさらに100日以上かかってしまいます。

ところが円周率計算の歴史を見ていると、最新の世界記録の検証は28時間で済んでいると書かれています。 2014年の記録でも182時間、その前年の記録でも46時間となっています。 全体の掛かった時間が100日という単位ですから、それに比べるとこれらの検証時間は短すぎます。 この時間で、別のアルゴリズムで全ての桁数計算したとは思えません。 もしこの時間で全く同じ桁数の計算を行えるなら、最初からそれを使ってますよね。 これは一体どういうことなのでしょうか。

実は、全ての桁数の検証をしているわけではないのです。 円周率の公式でBBP formulaと呼ばれているものがあります *2

 \displaystyle \pi = \sum_{n = 0}^{\infty} \frac{1}{16^{n}} \left(\frac{4}{8n+1} - \frac{2}{8n+4} - \frac{1}{8n+5} - \frac{1}{8n+6}\right)

この公式を使えば、円周率の16進数表現の任意の桁を、高速に求めることができるのです。 つまり16進数表現の小数点以下 N 桁目の数を、 N-1 桁までの数字を使うことも保存することもなく、直接求めることができるのです。 円周率計算の世界において、この公式が発見されたことは衝撃的でした。 今から当時の衝撃を感じるのは困難なのですが、それ以前はこの公式がなく、愚直に別のアルゴリズムで検証する他なかったことを想像すると、とても重要な公式の発見だと言えるでしょう。

BBP formulaの肝は、 1/16^{n} の分母が2の累乗であることです。 つまり、この公式を使うと、計算機上の数字の表現の生のデータを、直接検証することができます。 実はBBP公式の発見以降、類似の公式がいくつか見つかっています。 そのうちの1つが1997年に発見されたBellard’s formulaです *3

\displaystyle\pi=\frac{1}{2^{6}}\sum_{n=0}^{\infty}\frac{(-1)^{n}}{2^{10n}}\bigg(-\frac{2^{5}}{4n+1}-\frac{1}{4n+3}+\frac{2^{8}}{10n+1} \displaystyle \quad\quad\quad-\frac{2^{6}}{10n+3}-\frac{2^{2}}{10n+5}-\frac{2^{2}}{10n+7}+\frac{1}{10n+9}\bigg)

この公式も分母が2の累乗なので、円周率の16進数表現の検証に用いることができます。 同じ目的の公式が2つ以上あることで、仮に一方の計算に誤りがあったとしても、もう片方でさらに検証を重ねることができます。

以上のことを踏まえて、円周率計算の世界記録がどう計算され、どう検証されるのかをまとめてみました。

  • 円周率の計算を行う。2進数表現のままファイルに書いたり読んだりしながら計算が進みます。兆桁台にもなると、計算中にハードディスクが壊れるといったリスクもありますので、復帰ポイントを随時作りながら、プログラムをいったん止めても再開できるようにしておく必要があります。近年の世界記録はChudnovskyアルゴリズム*4にもとづいています。1989年に発見されたChudnovsky兄弟によって発見された公式で、収束が高速である (一項計算すると14桁進む) ことが知られています。
  • 2進数表現の円周率を、10進数表現に変換します。
  • 2進数表現の検証を行います。円周率の最後の方の桁を、BBP formulaやBellard’s formulaのような任意桁を求める公式によって検証します。イテレーションの早い段階や浅い桁で起きたエラーはすぐに最後の桁まで伝搬するので、この検証でかなり正当性が保証されます。
  • しかし前の桁で間違っている可能性は否定できません。特に、最後の掛け算中に浅い桁で誤りが起きる可能性があります。最終フェーズの掛け算は剰余を使って検証します。つまり十分大きな素数 p に対して、 nm\,\mathrm{mod}\,p = \left((n \,\mathrm{mod}\,p)(m \,\mathrm{mod}\,p)\right) \,\mathrm{mod}\, p を確認します。
  • 10進数表現の検証を行います。2進数表現が正しくても基数変換に誤りがある可能性があるので、ここも検証する必要があります。10進数表現を2進数表現に変換して元に戻れば確実です。素数で割った余りを比較する方法もあります。つまり素数 p と十分大きな数 n に対して   \lfloor 10^{n} \pi \rfloor \,\mathrm{mod}\, p を計算する方法です。10進数は10進数のまま、16進数は16進数のまま10掛けたり剰余を計算し、最後に計算を突き合わせるのです。

整数演算は2進数で行いますが、最後は10進数表現の円周率に変換しないと何兆桁計算したとはいえません。 この基数変換の検証も重要な部分なのですが、このエントリーでは省略します。 なお、10進数版BBP formulaがもし発見できれば、円周率10進数表現の任意桁をすぐに求められるようになり、計算結果の検証どころか、何百兆桁, 何千兆桁目を現実的な時間の中で求めることさえできてしまうという驚異的な結果なのですが、2017年現在ではそういう式は見つかっていません。

円周率計算の検証の更なる詳細はPi - 5 Trillion DigitsAlexander Yeeさんのstackoverflowの回答2016年世界記録のブログなどをご覧ください。 BBP formulaやBellard’s formulaがいかに重要かがよく分かると思います。

BBP formula

BBP formulaを実装してみましょう。 あらためて公式を眺めてみます。

 \displaystyle \pi = \sum_{n = 0}^{\infty} \frac{1}{16^{n}} \left(\frac{4}{8n+1} - \frac{2}{8n+4} - \frac{1}{8n+5} - \frac{1}{8n+6}\right)

求めたいのは円周率の16進数表現の、小数点以下 d 桁目の数字です。 ここで小数部分を返す \{x\} := x - \lfloor x \rfloor (\lfloor x \rfloorガウス記号で x 以下の最大の整数) を定義すると、

 \displaystyle \{16^{d}\pi\} = \{ 4\{16^{d} S_1\} - 2\{16^{d} S_4\} - \{16^{d} S_5\} - \{16^{d} S_6\} \}

と書くことができます。ただし

 \displaystyle S_j := \sum_{i=0}^{\infty} \frac{1}{16^{i} (8k + j)}

と定義しました。 ここで和を2つに分けて

 \displaystyle \{16^{d} S_j\} = \bigg\{ \sum_{i=0}^{d} \frac{16^{d-i}}{8k+j} + \sum_{i=d+1}^{\infty} \frac{16^{d-i}}{8k+j} \bigg\}

と書けますが、一項目の  16^{d-i} は整数で、かつ小数部分を返すのは和のそれぞれの項に分けても同じ (\{\sum_k a_k\} = \{\sum_k \{a_k\}\}) を用いると次のように書くことができます。

 \displaystyle \{16^{d} S_j\} = \bigg\{ \bigg\{ \sum_{i=0}^{d} \frac{16^{d-i}\,\mathrm{mod}\,(8k+j)}{8k+j} \bigg\} + \sum_{i=d+1}^{\infty} \frac{16^{d-i}}{8k+j} \bigg\}

第一項の d-i は何千万何億となるので、 16^{d-i} が大きくなりすぎないように剰余を取るのが大事なのです (いわゆる「冪剰余」です)。 第二項はどんどん0に収束していくので、適当なところまで計算したら和を取ってしまいます。

近年Rustが大流行していて、利用者が爆発的に増加しているという状況の中で、私も早く入門したいという気持ちが強まっています。 円周率計算はちょうどよい題材なのでRustで実装してみました。

github.com

単一ファイルで実行できるように、ここにもおいておきます。

pihex_bbp.rs

use std::env;
use std::thread;
use std::sync::mpsc;

fn main() {
    let place: u64 = env::args().nth(1).and_then(|x| x.parse().ok()).unwrap_or(0);
    print!("{}:", place);
    for i in 0..4 {
        print!(" {}", pihex_bbp(place + 4 * i));
    }
    print!("\n");
}

fn pihex_bbp(d: u64) -> String {
    let (tx, rx) = mpsc::channel();
    for &(j, k) in &[(1, 4.0), (4, -2.0), (5, -1.0), (6, -1.0)] {
        let tx = tx.clone();
        thread::spawn(move || tx.send(k * series_sum(d, j)).unwrap());
    }
    drop(tx);
    let fraction: f64 = rx.iter().sum();
    (0..4)
        .scan(fraction, |x, _| {
            *x = (*x - x.floor()) * 16.0;
            Some(format!("{:x}", x.floor() as u64))
        })
        .fold(String::new(), |s, t| s + &t)
}

fn series_sum(d: u64, j: u64) -> f64 {
    let fraction1: f64 = (0..d + 1)
        .map(|i| (pow_mod(16, d - i, 8 * i + j) as f64) / ((8 * i + j) as f64))
        .fold(0.0, |x, y| (x + y).fract());
    let fraction2: f64 = (d + 1..)
        .map(|i| 16.0_f64.powi(-((i - d) as i32)) / ((8 * i + j) as f64))
        .take_while(|&x| x > 1e-13_f64)
        .sum();
    fraction1 + fraction2
}

fn pow_mod(n: u64, m: u64, d: u64) -> u64 {
    if m == 0 {
        1 % d
    } else if m == 1 {
        n % d
    } else {
        let k = pow_mod(n, m / 2, d);
        if m % 2 == 0 {
            (k * k) % d
        } else {
            (k * k * n) % d
        }
    }
}

実行してみましょう

 $ rustc -O pihex_bbp.rs
 $ ./pihex_bbp 1000000
1000000: 6c65 e52c b459 3500
# ./pihex_bbp 1000000  4.52s user 0.04s system 305% cpu 1.494 total
 $ ./pihex_bbp 10000000
10000000: 7af5 863e fed8 de97
# ./pihex_bbp 10000000  52.90s user 0.30s system 351% cpu 15.154 total
 $ ./pihex_bbp 100000000
100000000: cb84 0e21 926e c5ae
# ./pihex_bbp 100000000  606.44s user 3.26s system 361% cpu 2:48.52 total
 $ sysctl -n machdep.cpu.brand_string
Intel(R) Core(TM) i7-5557U CPU @ 3.10GHz

およそ3分で1億桁目の数字が表示されました。 このコードでは4回pihex_bbpの計算を行っているため、4倍時間がかかっていることに注意してください。 1億桁目付近ののcb84は45秒ほどで計算することができます。

円周率自体の計算は、桁が増えるほど多くのメモリーやディスクを消費します。 一方で、BBP formulaを使ったアルゴリズムでは計算過程で桁が増えていくことはありません。 小数部分をdoubleでひたすら計算していくので、消費メモリーが増えていくこともありませんし、桁数が増えるために基本演算の計算時間が伸びていくといったこともありません。 これが、100日かかる世界記録をたった数十時間で検証できる理由なのです (実際はこの公式以外にも複合的な観点から検証されます。いずれにせよBBP formulaを使ってランダムサンプリングした桁の数字をチェックしたり、剰余を取って検証したりという感じで、計算結果の数字をくまなくチェックすることはありません)。

Bellard’s formula

BBP formulaと並んでよく使われる公式がBellard’s formulaです。 BBP formulaを使って検証したとしても、BBPの実装自体が誤っている可能性も否定できませんので、Bellard’s formulaも実装して結果を突き合わせておきましょう。

公式はこのようになっています。 \displaystyle\pi=\frac{1}{2^{6}}\sum_{n=0}^{\infty}\frac{(-1)^{n}}{2^{10n}}\bigg(-\frac{2^{5}}{4n+1}-\frac{1}{4n+3}+\frac{2^{8}}{10n+1} \displaystyle \quad\quad-\frac{2^{6}}{10n+3}-\frac{2^{2}}{10n+5}-\frac{2^{2}}{10n+7}+\frac{1}{10n+9}\bigg)

ここで

 \displaystyle S_{j,k} := \frac{1}{2^{6}} \sum_{i=0}^{\infty} \frac{(-1)^{i}}{1024^{i} (ij + k)}

とおくと、次の関係がなりたちます。

 \displaystyle \{16^{d} S_{j,k}\} = \bigg\{ \bigg\{ \sum_{i=0}^{\lfloor(2d-3)/5\rfloor} (-1)^{i} \bigg\{ \frac{4^{2d-3-5i}\,\mathrm{mod}\,(ij+k)}{ij+k} \bigg\} \bigg\}  \displaystyle \quad\quad\quad\quad\quad\quad + \sum_{i=\lfloor(2d-3)/5\rfloor+1}^{\infty} \frac{(-1)^{i} 4^{2d-3-5i}}{ij+k} \bigg\}

この式を使って実装してみました。

pihex_bellard.rs

use std::env;
use std::thread;
use std::sync::mpsc;

fn main() {
    let place: u64 = env::args().nth(1).and_then(|x| x.parse().ok()).unwrap_or(0);
    print!("{}:", place);
    for i in 0..4 {
        print!(" {}", pihex_bellard(place + 4 * i));
    }
    print!("\n");
}

fn pihex_bellard(d: u64) -> String {
    let (tx, rx) = mpsc::channel();
    for &(j, k, l) in &[(4, 1, -32.0),
                        (4, 3, -1.0),
                        (10, 1, 256.0),
                        (10, 3, -64.0),
                        (10, 5, -4.0),
                        (10, 7, -4.0),
                        (10, 9, 1.0)] {
        let tx = tx.clone();
        thread::spawn(move || tx.send(l * series_sum(d, j, k)).unwrap());
    }
    drop(tx);
    let fraction: f64 = rx.iter().sum();
    (0..4)
        .scan(fraction, |x, _| {
            *x = (*x - x.floor()) * 16.0;
            Some(format!("{:x}", x.floor() as u64))
        })
        .fold(String::new(), |s, t| s + &t)
}

fn series_sum(d: u64, j: u64, k: u64) -> f64 {
    let fraction1: f64 = (0..(2 * d + 2) / 5)
        .map(|i| {
            ((if i % 2 == 0 { 1.0 } else { -1.0 }) *
             pow_mod(4, 2 * d - 3 - 5 * i, j * i + k) as f64) / ((j * i + k) as f64)
        })
        .fold(0.0, |x, y| (x + y).fract());
    let fraction2: f64 = ((2 * d + 2) / 5..)
        .map(|i| -(-4.0_f64).powi(-((5 * i + 3 - 2 * d) as i32)) / ((j * i + k) as f64))
        .take_while(|&x| x.abs() > 1e-13_f64)
        .sum();
    fraction1 + fraction2
}

fn pow_mod(n: u64, m: u64, d: u64) -> u64 {
    if m == 0 {
        1 % d
    } else if m == 1 {
        n % d
    } else {
        let k = pow_mod(n, m / 2, d);
        if m % 2 == 0 {
            (k * k) % d
        } else {
            (k * k * n) % d
        }
    }
}

実行してみます。

 $ rustc -O pihex_bellard.rs
 $ ./pihex_bellard 1000000
1000000: 6c65 e52c b459 3500
# ./pihex_bellard 1000000  3.40s user 0.02s system 341% cpu 1.003 total
 $ ./pihex_bellard 10000000
10000000: 7af5 863e fed8 de97
# ./pihex_bellard 10000000  39.30s user 0.23s system 362% cpu 10.905 total
 $ ./pihex_bellard 100000000
100000000: cb84 0e21 926e c5ae
# ./pihex_bellard 100000000  447.77s user 2.75s system 355% cpu 2:06.71 total

先程のBBPの実行結果と見比べると、およそ25%速くなっていることがわかります。 Bellard’s formulaはBBP formulaよりも43%速いと言われます *5。 そこまでの記録は出ませんでしたが、確かにBBPより速いようだということは観測できました。

1000万桁目の悲劇

上記のコードは最初からこの形をしていたわけではありませんでした。 最初修正に手間取った部分は以下のコードでした。

    let fraction1: f64 = (0..d + 1)
        .map(|i| (pow_mod(16, d - i, 8 * i + j) as f64) / ((8 * i + j) as f64))
        .fold(0.0, |x, y| (x + y).fract());

この最後の行は元々次のように書いていました。

        .sum();

当初、総和を取っているから.sum()でという短絡的な考えでコードを書いていました。 その結果1000万桁目は次のようになっていました。

 $ pihex 10000000
10000000: 7af5 863f fed8 de97

この結果が間違っているのがわかりますか? 次の結果が正解です。

10000000: 7af5 863e fed8 de97

sumを使うと、863e863fとなってしまいます。 偶然The BBP Algorithm for Pi (David H. Bailey) に計算結果が掲載されており、それと見比べていた時に気がついたのです。 863efeという数字の並びからお察しの通り、誤差が積もってこういう間違いが起きてしまっています。

計算に用いた式をあらためて見てみましょう。

 \displaystyle \{16^{d} S_j\} = \bigg\{ \bigg\{ \sum_{i=0}^{d} \frac{16^{d-i}\,\mathrm{mod}\,(8k+j)}{8k+j} \bigg\} + \sum_{i=d+1}^{\infty} \frac{16^{d-i}}{8k+j} \bigg\}

この第一項の総和に注目してみます。 各項は剰余を取るために0以上1以下の数なのですが、何千万項も足しているために、結果が何百万という数になって、それだけで6桁も損してしまいます。 総和を取ってから小数部分を返すのではなく、各項を足しては小数部分にすることを繰り返すことで、常に和を1以下に保ち、小数点以下の精度が落ちることを防ぐ必要があったのです。

エラーが起きやすい項が1000万桁付近にあったことはとても幸運でした。 計算結果が載っている論文があったことも恵まれていました。 もし仮にこの誤差に気がついていなければ、誤った実装のまま満足してしまったことでしょう。

総和を一気にとるのではなく、足しては小数部分にすることを返すことでかなり誤差は減らすことができました。 しかし厳密に言うと、これでもあらゆる桁に正しく返せるとは言えません。 倍精度浮動小数点の仮数部、つまり52bit以上0fが揃っていたら、どうしても浮動小数点計算誤差の影響を受けてしまいます。 こういう場所が円周率のどの桁にあるかは知りませんが、そういう場所がないとは言いきれない以上、浮動小数点を使う計算で以って絶対に正しく答えを返すという保証はないのです。

BBP formulaやBellard’s formulaを計算機で実装すると、浮動小数点の誤差を必ず受けうることになります。 上記のように小数部分のみで計算進めるのは必要な対処ですが、その対処をしても避けられないケースは存在するのです。 しかし、このアルゴリズムを使う意図を思い出してみましょう。 計算した円周率の最後の方の桁の検証のために、特定の桁を計算するために使います。 計算結果を見て0fが続いているかどうか、すなわち誤差の影響を受けうるかどうかは目で見て分かるので、そこまで大きな問題にはならないのです。

10億桁目の悪夢

上記の浮動小数点計算誤差に起因するバグを修正したことで、かなりのエラーは減りました。 ところが、いざ10億桁目を計算してみると、まためちゃくちゃな答えを返すようになりました。 これくらいの桁にもなると計算自体に数分かかるようになり、デバッグも困難を極めました。

上記のコードで10億桁目を計算した結果、次のようになっていました (時間がかかるのでpihex_{bbp,bellard}を1回にしています)。

 $ ./pihex_bbp 1000000000
1000000000: 4fa1de2c22c1
 $ ./pihex_bellard 1000000000
1000000000: a3222cec6720

両者の結果が異なっている時点で少なくとも一方が間違っていると分かりますし、そもそも正解は5895585a0428...ですから (David H. Bailey, The BBP Algorithm for Pi)、両者ともめちゃくちゃです。

BBP formulaとBellard’s formulaの両方で誤っていることから、共通の関数であるpow_modを調べました。 しばらくしてからk * k * nが64bit整数から溢れていることに気が付き、まず次のように書き直しました。

        (((k * k) % d) * if m % 2 == 0 { 1 } else { n % d }) % d

つまりk * k * nを一気に計算せずに、一回掛けては剰余を取るようにしたのです。 これで10億桁目を計算してみました。

 $ ./pihex_bbp 1000000000
1000000000: 73ccc138ae45
 $ ./pihex_bellard 1000000000
1000000000: 5895585a7ef4

Bellard’s formulaは5895585aまで合ってますが、BBP formulaの結果はまだ間違っています。 これではだめですね。

そこで  k^{2} \, \mathrm{mod}\, d = (d - k)^{2} \, \mathrm{mod}\, d であることを用いて次のように書き換えてみました。

        let l = if k <= d / 2 { k } else { d - k };
        (((l * l) % d) * if m % 2 == 0 { 1 } else { n }) % d

つまり kd-k の小さい方を二乗するというわけです。 実行結果は以下のようになります。

 $ ./pihex_bbp 1000000000
1000000000: 58955859f966
 $ ./pihex_bellard 1000000000
1000000000: 5895585a7ef4

ようやく両者のプログラム共に、10億桁目の正しい値を得られたようです。

pow_modの第三引数は8i + jですから、10億桁目を求める時は80億、つまり  8 \times 10^{9} まで大きくなります。 上記のコードにおけるkはすなわちpow_modの返り値ですから、これも最大で 8 \times 10^{9} くらいになります。 それを二乗すると 6.4 \times 10^{19} となり、unsigned 64bit整数の限界の 2^{64}-1\simeq 1.84 \times 10^{19} を簡単に超えてしまうのです。

 d/2 を超えたら d-k を使うようにすると、掛け算の最大値を4分の1におさえることができます。 そうすると最大でも 1.6 \times 10^{19} つまりunsigned 64bit整数の限界をぎりぎりで下回り、オーバーフローすることなく計算できるのです。

100億桁目への挑戦

ところが100億桁となると、またいともたやすく64bit整数の限界を超えてしまいます。 あらためて状況を整理すると次のようになります。

  • d 10^{10} を超える数で、d 以下の数 k に対して、k^{2} \, \mathrm{mod}\, d を求めたい

d の最大値の二乗が64bit整数の最大値を超えているので、単純な二乗はできません。

これまで64bit整数の限界と戦ってきたわけですが、よく調べてみるとRustでは128bit整数を使うことができます (記事執筆時点では#![feature(i128_type)]が必要ですが、そのうち書かなくてもよくなるでしょう*6 )。 つまり、最初のコード中のu64u128に書き換えるだけで、オーバーフローの苦労なしに10億桁、そして100億桁を突破することができます。

以下の結果が、実際に128bit整数にして実行してみた様子です。

 $ ./pihex_bbp 1000000000
1000000000: 58955859f966
# ./pihex_bbp 1000000000  7120.63s user 20.24s system 375% cpu 31:47.46 total
 $ ./pihex_bbp 10000000000
10000000000: 21c73c6889eb
# ./pihex_bbp 10000000000  139050.65s user 365.69s system 373% cpu 10:21:40.48 total
 $ ./pihex_bellard 1000000000
1000000000: 5895585a7ef8
# ./pihex_bellard 1000000000  2643.48s user 6.98s system 374% cpu 11:47.56 total
 $ ./pihex_bellard 10000000000
10000000000: 21c73c65d3f4
# ./pihex_bellard 10000000000  80295.34s user 182.43s system 384% cpu 5:48:59.45 total

やった! 100億桁目を求めることができました! 実装上はdが大きくてオーバーフローの可能性がある場合にのみ128bitで計算するようにしています。 64bitで計算できる範囲のものを全て128bitで計算してしまうと、計算時間が余計にかかってしまうからです。

計算がうまく行ってしまうと感覚が麻痺してしまいますが、100億桁目の計算に成功したことに大きな喜びを感じています。 100億桁ですよ、100億桁! 公式をほぼそのまま書き下したコードを実行したら、ノートPCの上で現実的な時間内にこれだけの計算ができるというのは、本当に驚きました。 さらに計算機リソースの並列度を上げたり、計算結果をキャッシュできるものはキャッシュするとか、fract()の呼び出しを減らすなど、改善の余地はあると思います。

私が五年前に計算したのは円周率1億桁でした。 Bellard’s formulaを使えば1億桁目のcb84を計算するのに30秒もかかりません。 10進数と16進数で見ている桁数が異なることに注意は必要ですが、十何分も掛けて計算した1億桁の最後の数字を30秒もかけずに検証できるということを肌で感じることができました。 数式を変形したり、コードを書いて試行錯誤する中で、円周率計算世界記録の技術にまた一歩近づけたのではないかなと思います。

おわりに

世界新記録の結果を検証するのは大変なことです。 同じ結果を持っている人がその人しかいないからです。 全く異なる方法で計算すると、元の計算以上の時間がかかってしまいます。 そういう状況で、特定の桁を直接求めることができるBBP formulaはとても重要な存在なのです。 浮動小数点を書き換えていくだけなので、巨大整数をファイルに書く必要もありませんし、メモリーも一定量しか消費しません。

円周率の16進数表現は、Blowfish暗号でも使われます。 Blowfish暗号の実装のためにBBP formulaを用いることもできますが、Blowfishで必要な桁数は決まっているので、ソースコードにべた書きされることが多いようです。 逆にBBPアルゴリズムを実装したら、検証データとしてBlowfish暗号の実装コードを参照することもできます。

浮動小数点演算に依存している以上、誤差は無視できるものではありません。 53bitの誤差が16進数表現の1桁をひっくり返すこともありえます。 ただし、ひっくり返ったかどうかは計算結果を見れば明らかに分かります。 その次の桁を見て000とかfffとなっているかどうかを見ればよいのです。 このアルゴリズムには、任意の桁の数字を正確に求めるよりも、小数点以下深い特定の桁を高速に求めることが要求されます。 間違える桁があることが否定できなくても、間違えたことがすぐに分かるのならそれでいいわけです。

本記事を書くにあたって、Rustを初めて真面目に書いてみました (以前、チュートリアルを少し触ったくらいでした)。 かなり書き心地がよくてびっくりしました。 普段書いている言語に対する不満が、いくつも解決されているように思えます。 今後はRustでもう少し大きなものを書いてみたいと思っています。

円周率計算は奥深い世界です。 円周率を計算する公式、Chudnovsky algorithmやBBP formula、Bellard’s formulaなど、数論や解析学の結果によって何百日もの効率化につながることがあります。 円周率の16進数表現というわりとニッチな分野を調べていたのですが、Blowfish暗号と円周率の関係性を見つけるなど、意外な出会いもありました。

y-cruncherは円周率計算記録の多くで使われていることから分かるように、速度と安定性に評判があるソフトウェアです。 一方でy-cruncherのソースコードは公開されておらず、トータルで30万行のC++やインラインアセンブラで書かれているなど、なかなか謎の多いソフトウェアです。 ここまでのソフトウェアを書き上げた作者の情熱がどういうところにあるのか、とても気になります。

近年の円周率計算の世界記録は、y-cruncherの一強となっています。 これは面白くありませんので、いつかは兆桁台にも挑戦してみたいなと思っています。 みなさんも、ぜひ世界記録に挑戦してみてください。 Fabrice Bellard氏 *7 のComputation of 2700 billion decimal digits of Pi using a Desktop Computer (2009) は、円周率計算に関する素晴らしい資料です。 円周率計算の世界に少しでも興味を持っていただけたら幸いです。

*1:Chronology of computation of π https://en.wikipedia.org/wiki/Chronology_of_computation_of_%CF%80

*2:David H. Bailey, Peter Borwein, and Simon Plouffe, On the Rapid Computation of Various Polylogarithmic Constants, Mathematics of Computation 66, 903-913, 1997.

*3:Fabrice Bellard, A new formula to compute the n'th binary digit of pi, 1997.

*4:David V. Chudnovsky, Gregory V. Chudnovsky, The Computation of Classical Constants, Proceedings of the National Academy of Sciences of the United States of America 86 (21) 8178–8182, 1989.

*5:David H. Bailey. The BBP Algorithm for Pi, 2006.

*6:Tracking issue for 128-bit integer support (RFC 1504) https://github.com/rust-lang/rust/issues/35118

*7:Bellard’s formulaの発見者であり、2009年12月に円周率を2.7兆桁計算して世界記録をとったプログラマです。FFmpegの作者でもあります。

LLVM APIを使ってみよう! 〜 Brainf**kコンパイラをIRBuilderで書き直してみた 〜

先日LLVMの入門記事を書きました。 clangが吐くLLVM IR (Intermediate representation, 中間表現) を頼りに、Brainf**kコンパイラを書いてみました。 itchyny.hatenablog.com この記事で書いたコードでは、直接printfでLLVM IRの命令を出力していました。 このステップを踏むことで、LLVM IRの命令をどう調べればいいかについて身についたと思います。

しかし、この「コンパイラ」は次のような問題がありました。

  • bf2llvmコマンドが出力するのがLLVM IRのために、llillcといったLLVM IRのランタイムやコンパイラが必要となる
  • 最適化を行うにはさらにoptコマンドが必要になり、やはりLLVMツールチェインをインストールしている環境でしか使えない
  • ソースコードから実行ファイルまでのパスの中で、LLVM IRの文字列に変換してパースしている箇所がある
  • 不正なIRコードを吐く可能性も高い
  • 一時変数のインデックスを自前で管理しており、メンテナンス性も悪い

LLVMツールチェインには、LLVM IRのコードをプログラムから吐くためのAPIも含まれています。 今回はLLVM APIを使ってコンパイラを作る練習をしてみます。

LLVM API, Hello world!

まずはとても短いコードを生成するところから始めましょう。 Kaleidoscopeチュートリアルを横目にしながらコードを書いていきます。

次のようなコードを考えてみます。

test.c

int main() {
  return 42;
}

このコードに対応するLLVM IRのコードは次のようになります (clang -S -emit-llvm -O2 test.cの出力を削ったものです)。

test.ll

define i32 @main() {
  ret i32 42
}

このIRを吐くコードを書いてみます。

test.cpp

#include "llvm/IR/IRBuilder.h"
#include "llvm/IR/LLVMContext.h"
#include "llvm/IR/Module.h"

static llvm::LLVMContext TheContext;
static llvm::IRBuilder<> Builder(TheContext);
static std::unique_ptr<llvm::Module> TheModule;

int main() {
  TheModule = llvm::make_unique<llvm::Module>("top", TheContext);
  llvm::Function* mainFunc = llvm::Function::Create(
      llvm::FunctionType::get(llvm::Type::getInt32Ty(TheContext), false),
      llvm::Function::ExternalLinkage, "main", TheModule.get());
  Builder.SetInsertPoint(llvm::BasicBlock::Create(TheContext, "", mainFunc));

  Builder.CreateRet(Builder.getInt32(42));

  TheModule->dump();
}

コンパイルして実行してみます。

 $ g++ `llvm-config --cxxflags --ldflags --libs --system-libs` test.cpp -o test
 $ ./test
; ModuleID = 'top'
source_filename = "top"

define i32 @main() {
  ret i32 42
}
 $ ./test 2>&1 | lli
 $ echo $?
42

TheModule->dump()標準エラー出力にIRを出力するので2>&1しています。 lliを使うとLLVM IR命令列を直接実行して、動作を確認することができます。

簡単なところから説明していきましょう。 次の行はret i32 42を生成しているコードです。

  Builder.CreateRet(Builder.getInt32(42));

IRBuilderの各種メソッドを調べる時は、doxygenで生成されたヘルプが役に立ちます。 CreateRetの引数は1つで、Value*を取ります。 getInt32IRBuilderBaseで定義されています。

次のコードはdefine i32 @main()に対応します。

  llvm::Function* mainFunc = llvm::Function::Create(
      llvm::FunctionType::get(llvm::Type::getInt32Ty(TheContext), false),
      llvm::Function::ExternalLinkage, "main", TheModule.get());

llvm::FunctionType::getでググって出て来るFunctionTypeを見ると、llvm::FunctionType::getの引数は2つまたは3つで、引数の型は省略できることがわかります。 関数を作った後に中身を生成するには、SetInsertPointを使ってコードを吐く場所を指定します。

  Builder.SetInsertPoint(llvm::BasicBlock::Create(TheContext, "", mainFunc));

llvm::BasicBlock::Createはコードブロックを作り、ラベルを生成します。

  Builder.SetInsertPoint(llvm::BasicBlock::Create(TheContext, "entry", mainFunc));

このようにラベル名を指定すると、次のようになります。

define i32 @main() {
entry:
  ret i32 42
}

ラベルを使って遊んでみましょう。 br (branch) 命令を使って作ったラベルにジャンプしてみます。

test.c

  Builder.SetInsertPoint(llvm::BasicBlock::Create(TheContext, "entry", mainFunc));
  llvm::BasicBlock* newBlock = llvm::BasicBlock::Create(TheContext, "new_label", mainFunc);
  Builder.CreateBr(newBlock);

  Builder.CreateRet(Builder.getInt32(42));

  Builder.SetInsertPoint(newBlock);
  Builder.CreateRet(Builder.getInt32(36));

動作を確認します。

 $ g++ `llvm-config --cxxflags --ldflags --libs --system-libs` test.cpp -o test
 $ ./test
; ModuleID = 'top'
source_filename = "top"

define i32 @main() {
entry:
  br label %new_label
  ret i32 42

new_label:                                        ; preds = %entry
  ret i32 36
}
 $ ./test 2>&1 | lli
 $ echo $?
36

llvm::BasicBlock::CreateBuilder.CreateBrの順に呼んでいますが、きちんと下のラベルへのbr命令を作ることができています。 llvm::BasicBlock::Createは、指定された関数の一番最後にブロックを生成するという関数なのです (第四引数を指定して生成場所を制御することもできます)。

次に、簡単な演算を行うコードを生成してみます。

int main() {
  int a, b;
  a = 32;
  b = a + 8;
  return a + b;
}

72で終了するコードです。 これに対応するLLVM IRコードを手で書いてみます。 IRに慣れてきて、これくらいのサイズなら手で書けるようになってきました。

test.ll

define i32 @main() {
  %a = alloca i32
  %b = alloca i32
  store i32 32, i32* %a
  %1 = load i32, i32* %a
  %2 = add i32 %1, 8
  store i32 %2, i32* %b
  %3 = load i32, i32* %a
  %4 = load i32, i32* %b
  %5 = add i32 %3, %4
  ret i32 %5
}

実際このコードは実行できます。

 $ lli test.ll
 $ echo $?
72

このLLVM IRを出力するコードをIRBuilderを使って書いてみます。

test.cpp

#include "llvm/IR/IRBuilder.h"
#include "llvm/IR/LLVMContext.h"
#include "llvm/IR/Module.h"

static llvm::LLVMContext TheContext;
static llvm::IRBuilder<> Builder(TheContext);
static std::unique_ptr<llvm::Module> TheModule;

int main() {
  TheModule = llvm::make_unique<llvm::Module>("top", TheContext);
  llvm::Function* mainFunc = llvm::Function::Create(
      llvm::FunctionType::get(llvm::Type::getInt32Ty(TheContext), false),
      llvm::Function::ExternalLinkage, "main", TheModule.get());
  Builder.SetInsertPoint(llvm::BasicBlock::Create(TheContext, "", mainFunc));

  llvm::Value* a = Builder.CreateAlloca(Builder.getInt32Ty(), nullptr, "a");
  llvm::Value* b = Builder.CreateAlloca(Builder.getInt32Ty(), nullptr, "b");
  Builder.CreateStore(Builder.getInt32(32), a);
  Builder.CreateStore(Builder.CreateAdd(Builder.CreateLoad(a), Builder.getInt32(8)), b);

  Builder.CreateRet(Builder.CreateAdd(Builder.CreateLoad(a), Builder.CreateLoad(b)));

  TheModule->dump();
}

動作を確認します。

 $ g++ `llvm-config --cxxflags --ldflags --libs --system-libs` test.cpp -o test
 $ ./test
; ModuleID = 'top'
source_filename = "top"

define i32 @main() {
  %a = alloca i32
  %b = alloca i32
  store i32 32, i32* %a
  %1 = load i32, i32* %a
  %2 = add i32 %1, 8
  store i32 %2, i32* %b
  %3 = load i32, i32* %a
  %4 = load i32, i32* %b
  %5 = add i32 %3, %4
  ret i32 %5
}
 $ ./test 2>&1 | lli
72

少しずつIRBuilderの触り方が分かってきました。 LLVM IRのどの命令を使えばいいかわかっていれば、IRBuilderのドキュメントの中でページ内検索してコードを書いていくことができます。

  Builder.CreateStore(Builder.CreateAdd(Builder.CreateLoad(a), Builder.getInt32(8)), b);

この一行が次の三行に対応することがわかります。

  %1 = load i32, i32* %a
  %2 = add i32 %1, 8
  store i32 %2, i32* %b

これならわざわざ一時変数のためのインデックスを自前で管理しなくても良いですね。 だいぶコードを書くのが楽になりますね。

きつねさんでもわかるLLVM ~コンパイラを自作するためのガイドブック~

きつねさんでもわかるLLVM ~コンパイラを自作するためのガイドブック~

Brainf**kのコンパイラを書いてみよう!

前回の記事でもBrainf**kのコンパイラを書きましたが、直接LLVM IRを吐くコードだったために一時変数を管理しないといけない、不正なコードを出力してしまいそうになる、ソースコードがダサいと言った問題がありました。 今回はLLVM APIを使って書き直してみましょう。

まずはBrainf**kのデータ領域の確保と開放を行う部分を書いてみます。 Cで書くと

#include <stdlib.h>

int main() {
  char* data = (char*)calloc(30000, sizeof(char));
  char* ptr = data;
  free(data);
}

に対応するコードです。

bf2llvm.cpp

#include "llvm/IR/IRBuilder.h"
#include "llvm/IR/LLVMContext.h"
#include "llvm/IR/Module.h"

static llvm::LLVMContext TheContext;
static llvm::IRBuilder<> Builder(TheContext);
static std::unique_ptr<llvm::Module> TheModule;

int main() {
  TheModule = llvm::make_unique<llvm::Module>("top", TheContext);
  llvm::Function* mainFunc = llvm::Function::Create(
      llvm::FunctionType::get(llvm::Type::getInt32Ty(TheContext), false),
      llvm::Function::ExternalLinkage, "main", TheModule.get());
  Builder.SetInsertPoint(llvm::BasicBlock::Create(TheContext, "", mainFunc));

  llvm::Value* data = Builder.CreateAlloca(Builder.getInt8PtrTy(), nullptr, "data");
  llvm::Value* ptr = Builder.CreateAlloca(Builder.getInt8PtrTy(), nullptr, "ptr");
  llvm::Function* funcCalloc = llvm::cast<llvm::Function>(
      TheModule->getOrInsertFunction("calloc",
        Builder.getInt8PtrTy(),
        Builder.getInt64Ty(), Builder.getInt64Ty(),
        nullptr));
  llvm::Value* data_ptr = Builder.CreateCall(funcCalloc, {Builder.getInt64(30000), Builder.getInt64(1)});
  Builder.CreateStore(data_ptr, data);
  Builder.CreateStore(data_ptr, ptr);

  llvm::Function* funcFree = llvm::cast<llvm::Function>(
      TheModule->getOrInsertFunction("free",
        Builder.getVoidTy(),
        Builder.getInt8PtrTy(),
        nullptr));
  Builder.CreateCall(funcFree, {Builder.CreateLoad(data)});

  Builder.CreateRet(Builder.getInt32(0));

  TheModule->dump();
}

実行して動作を確かめます。

 $ g++ `llvm-config --cxxflags --ldflags --libs --system-libs` bf2llvm.cpp -o bf2llvm
 $ ./bf2llvm
; ModuleID = 'top'
source_filename = "top"

define i32 @main() {
  %data = alloca i8*
  %ptr = alloca i8*
  %1 = call i8* @calloc(i64 30000, i64 1)
  store i8* %1, i8** %data
  store i8* %1, i8** %ptr
  %2 = load i8*, i8** %data
  call void @free(i8* %2)
  ret i32 0
}

declare i8* @calloc(i64, i64)

declare void @free(i8*)

関数callocfreeを呼んでいますね。

  llvm::Function* funcCalloc = llvm::cast<llvm::Function>(
      TheModule->getOrInsertFunction("calloc",
        Builder.getInt8PtrTy(),
        Builder.getInt64Ty(), Builder.getInt64Ty(),
        nullptr));
  llvm::Value* data_ptr = Builder.CreateCall(funcCalloc, {Builder.getInt64(30000), Builder.getInt64(1)});

getOrInsertFunctionでは関数名、返り値の型と引数の型を取ります。 モジュールの中から関数を探して返しますが、もしなければdeclareによって宣言します (リンク時に解決されます)。 そしてCreateCallによってcall命令を作っています。

次はポインタの移動を書いてみます。 Brainf**kでは>, <命令に相当するコードです。 Cだと++ptr, --ptrに対応します。 前回のエントリーでは次のように書いていました。

int idx = 1;
void emit_move_ptr(int diff) {
  printf("  %%%d = load i8*, i8** %%ptr, align 8\n", idx);
  printf("  %%%d = getelementptr inbounds i8, i8* %%%d, i32 %d\n", idx + 1, idx, diff);
  printf("  store i8* %%%d, i8** %%ptr, align 8\n", idx + 1);
  idx += 2;
}

これは次のようになります。

void emit_move_ptr(llvm::Value* ptr, int diff) {
  Builder.CreateStore(
      Builder.CreateInBoundsGEP(
        Builder.getInt8Ty(),
        Builder.CreateLoad(ptr),
        Builder.getInt32(diff)),
      ptr);
}

CreateInBoundsGEPgetelementptr命令を生成します。

次に、メモリーに足したり引いたりする+, -命令に相当するコードを書いてみます。 Cで書くと++*ptr--*ptrとなるコードですね。

void emit_add(llvm::Value* ptr, int diff) {
  llvm::Value* tmp = Builder.CreateLoad(ptr);
  Builder.CreateStore(
      Builder.CreateAdd(
        Builder.CreateLoad(tmp),
        Builder.getInt8(diff)),
      tmp);
}

そして、Brainf**kの., Cではputchar(*ptr)に相当するコードを書いてみます。 前回見たように、sext命令でi8からi32に変換する必要があります。

void emit_put(llvm::Value* ptr) {
  llvm::Function* funcPutChar = llvm::cast<llvm::Function>(
      TheModule->getOrInsertFunction("putchar",
        Builder.getInt32Ty(),
        Builder.getInt32Ty(),
        nullptr));
  Builder.CreateCall(
      funcPutChar,
      Builder.CreateSExt(
        Builder.CreateLoad(Builder.CreateLoad(ptr)),
        Builder.getInt32Ty()));
}

これで>, <, +, -, .の5命令の準備ができたので、標準入力からBrainf**kを受け取ってLLVM IRを吐くようにしてみましょう。

bf2llvm.cpp

// ...
  llvm::Value* data_ptr = Builder.CreateCall(funcCalloc, {Builder.getInt64(30000), Builder.getInt64(1)});
  Builder.CreateStore(data_ptr, data);
  Builder.CreateStore(data_ptr, ptr);

  char c;
  while (std::cin.get(c)) {
    switch (c) {
      case '>': emit_move_ptr(ptr, 1); break;
      case '<': emit_move_ptr(ptr, -1); break;
      case '+': emit_add(ptr, 1); break;
      case '-': emit_add(ptr, -1); break;
      case '.': emit_put(ptr); break;
    }
  }

  llvm::Function* funcFree = llvm::cast<llvm::Function>(
// ...

動かしてみましょう。

 $ g++ `llvm-config --cxxflags --ldflags --libs --system-libs` bf2llvm.cpp -o bf2llvm
 $ echo "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++.+++++++++++++++++++++++++++++.+++++++..+++.>++++++++++++++++++++++++++++++++.<++++++++.--------.+++.------.--------.>+." | ./bf2llvm 2>&1 | lli
Hello world!

正しくIRを吐けてるようですね。

ループに対応しよう

Brainf**kの[C言語while (*ptr) {に、]}に対応します。 LLVM IRではそれぞれ

    br label %while_cond

  while_cond:
    %5 = load i8*, i8** %%ptr
    %6 = load i8, i8* %5
    %7 = icmp ne i8 %6, 0
    br i1 %7, label %while_body, label %while_end

  while_body:

    br label %while_cond

  while_end:

に対応します。 Brainf**kのコードではループは何度も出てくるので、while_cond6 .. while_body6 .. while_end6のように数字をつけてラベルを区別する必要があります。 前回はこの番号をemit_while_startに渡してIRを直接出力していましたが、今回はジャンプする時にllvm::BasicBlock*を使う必要があるため、構造体を作って管理してみました。

struct WhileBlock {
  llvm::BasicBlock* cond_block;
  llvm::BasicBlock* body_block;
  llvm::BasicBlock* end_block;
};

void emit_while_start(llvm::Function* func, llvm::Value* ptr, WhileBlock* while_block, int while_index) {
  while_block->cond_block = llvm::BasicBlock::Create(
      TheContext, std::string("while_cond") + std::to_string(while_index), func);
  while_block->body_block = llvm::BasicBlock::Create(
      TheContext, std::string("while_body") + std::to_string(while_index), func);
  while_block->end_block = llvm::BasicBlock::Create(
      TheContext, std::string("while_end") + std::to_string(while_index), func);
  Builder.CreateBr(while_block->cond_block);
  Builder.SetInsertPoint(while_block->cond_block);
  Builder.CreateCondBr(
      Builder.CreateICmpNE(
        Builder.CreateLoad(Builder.CreateLoad(ptr)),
        Builder.getInt8(0)),
      while_block->body_block,
      while_block->end_block);
  Builder.SetInsertPoint(while_block->body_block);
}

void emit_while_end(WhileBlock* while_block) {
  Builder.CreateBr(while_block->cond_block);
  Builder.SetInsertPoint(while_block->end_block);
}

// ...
  int while_index = 0;
  WhileBlock while_blocks[1000];
  WhileBlock* while_block_ptr = while_blocks;

// ...
      case '[': emit_while_start(mainFunc, ptr, while_block_ptr++, while_index++); break;
      case ']': emit_while_end(--while_block_ptr); break;

llvm::BasicBlock::Createが関数の最後にブロックをつくるというのと、Builder.SetInsertPointでIRを吐くブロックを移動するということを理解しておけば、そんなに難しくないと思います。

 $ g++ `llvm-config --cxxflags --ldflags --libs --system-libs` bf2llvm.cpp -o bf2llvm
 $ echo "+++++++++[>++++++++>+++++++++++>+++++<<<-]>.>++.+++++++..+++.>-.------------.<++++++++.--------.+++.------.--------.>+." | ./bf2llvm 2>&1 | lli
Hello, world!
 $ cat prime.bf
++++++++++[->+++++>++++>+++<<<]>>++++>++>>+>+++<<<<<.>>>>>[<[-]++<[-]+>>>+[<[->>
+<<]<[->>>>+>+<<<<<]>>>>[-<<<<+>>>>]<[->+>-[>+>>]>[+[-<+>]>+>>]<<<<<<]>[-<<<+>>>
]>[-]>[-<<+>+>]<<[->>+<<]+>[<[-]>[-]]<[<<<<<[-]>>>>>[-]]>[-]>[-]>[-]<<<<<<-<[->>
>+>+<<<<]>>>>[-<<<<+>>>>]<<<[->>>>+>+<<<<<]>>>>>[-<<<<<+>>>>>]<<<+[->>[->+>+<<]>
>[-<<+>>]+<[>-<<->[-]]>[<<<+<[-]>>>>[-]]<[-]<<<]>>[-]<[<<->>[-]]<[-]<<+<[->>>+>+
<<<<]>>>>[-<<<<+>>>>]>+++++++++++++<<[->>[->+>+<<]>>[-<<+>>]+<[>-<<->[-]]>[<<<+<
[-]>>>>[-]]<[-]<<<]>>[-]<[<<[-[-]]>>[-]]<[-]<<<+>>]<<<[<.>>>>>++++++++++<<[->+>-
[>+>>]>[+[-<+>]>+>>]<<<<<<]>[-<+>]>[-]>>>++++++++++<[->-[>+>>]>[+[-<+>]>+>>]<<<<
<]>[-]>>[>++++++[-<++++++++>]<.<<+>+>[-]]<[<[->-<]++++++[->++++++++<]>.[-]]>[-]<
<<++++++[-<++++++++>]<.[-]<<<<<[-]]>>+]
 $ cat prime.bf | ./bf2llvm 2>&1 | /usr/local/opt/llvm/bin/lli
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251

大丈夫そうですね。

オブジェクトファイルを生成しよう

ここまでTheModule->dump()を使ってLLVM IRを出力してきました。 しかし、これでは実行するためにlliが必要となり、結局LLVMツールチェインの入ってない環境では動かすことができません。 KaleidoscopeのChapter 8を参考にしながら (まったく同じです。まだ理解が追いついてなくて…)、オブジェクトファイルを生成してみましょう。

TargetTripleとは、アーキテクチャ, ベンダー, プラットフォームのことです。 clang --version | grep Targetで確認することができます。 私の手元ではx86_64-apple-darwin15.6.0でした。

  std::string TargetTriple = llvm::sys::getDefaultTargetTriple();

このTargetTripleに対応するtargetを引きます。

  std::string err;
  const llvm::Target* Target = llvm::TargetRegistry::lookupTarget(TargetTriple, err);
  if (!Target) {
    std::cerr << "Failed to lookup target " + TargetTriple + ": " + err;
    return 1;
  }

CPUを指定します。 この指定によってメモリーalignmentやendiannessなどが変わるようです。

  llvm::TargetOptions opt;
  llvm::TargetMachine* TheTargetMachine = Target->createTargetMachine(
      TargetTriple, "generic", "", opt, llvm::Optional<llvm::Reloc::Model>());

ここで引いたtargetとtarget machineを、私たちのモジュールにセットします。

  TheModule->setTargetTriple(TargetTriple);
  TheModule->setDataLayout(TheTargetMachine->createDataLayout());

オブジェクトファイルを開いて

  std::string Filename = "output.o";
  std::error_code err_code;
  llvm::raw_fd_ostream dest(Filename, err_code, llvm::sys::fs::F_None);
  if (err_code) {
    std::cerr << "Could not open file: " << err_code.message();
    return 1;
  }

書き込んで終わりです。

  llvm::legacy::PassManager pass;
  if (TheTargetMachine->addPassesToEmitFile(pass, dest, llvm::TargetMachine::CGFT_ObjectFile)) {
    std::cerr << "TheTargetMachine can't emit a file of this type\n";
    return 1;
  }
  pass.run(*TheModule);
  dest.flush();

動かしてみましょう。

 $ g++ `llvm-config --cxxflags --ldflags --libs --system-libs` bf2llvm.cpp -o bf2llvm
 $ echo "+++++++++[>++++++++>+++++++++++>+++++<<<-]>.>++.+++++++..+++.>-.------------.<++++++++.--------.+++.------.--------.>+." | ./bf2llvm
Wrote output.o
 $ gcc output.o
 $ ./a.out
Hello, world!
 $ cat fibonacci.bf
>>+++++++++++[-<<++++>+++>>+<]>>+<++<<->>[>>>++++++++++<<[->+>-[>+>>]>[+[-<+>]>+
>>]<<<<<<]>>[-]>>>++++++++++<[->-[>+>>]>[+[-<+>]>+>>]<<<<<]>[-]>>[>++++++[-<++++
++++>]<.<<+>+>[-]]<[<[->-<]++++++[->++++++++<]>.[-]]<<++++++[-<++++++++>]<.[-]<<
[-<+>]<<[->>+>+<<<]>>>[-<<<+>>>]<-[<<<<<.>.>>>>[-]]<[->+>+<<]>>[-<<+>>]<<<<[->>+
<<]>>>[-<<<+>>>]<<-]
 $ cat fibonacci.bf | ./bf2llvm
Wrote output.o
 $ gcc output.o
 $ ./a.out
1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233

きちんと動いてそうですね。

まとめ

Brainf**kのソースコードを受け取り、オブジェクトファイルを生成するコンパイラを作ることができました。 LLVM APIを使うことで、以下のような恩恵がありました。

  • 一時変数のインデックスをいちいち管理しなくてよい
  • オブジェクトファイルを生成することができる
  • IRの命令が抽象化されており、おかしなIRを吐く確率が減る
  • IRを直接書いていると、外部の関数のdeclareを書かなければいけなかったが、getOrInsertFunctionを使うと気にする必要がなくなる

APIを使ってコードを書くことは意外と骨が折れる事で、結局doxygenの生成ドキュメントを行ったり来たりしながら試行錯誤する必要がありました。 意外だったことは、IRのシンタックス的におかしなコードを吐くこともあるということです。 例えば

  llvm::Value* a = Builder.CreateAlloca(Builder.getInt32Ty(), nullptr, "a");
  Builder.CreateAdd(Builder.CreateLoad(a), Builder.getInt8(8));

このコードは次のようなコードを吐きます。

  %1 = load i32, i32* %a
  %2 = add i32 %1, i8 8

実行してみるとエラーになります。

lli: <stdin>:7:20: error: expected value token
  %2 = add i32 %1, i8 8
                   ^

型が違う、ではなくてシンタックスエラーになることもあるんですね。 将来挙動が変わるかもしれませんが。

前回の記事で、ある程度LLVM IRに慣れ親しんでいたために、IRBuilderもかなり楽に触ることができました。 逆に、IRの命令をなにも知らずにIRBuilderでコード生成するのはかなり難しいと思います。 対応するCのコードがあって、対応するLLVM IRコードがあって、それを出力するようなIRBuilderの関数を探します。 Brainf**kのコンパイラ作りを通じてかなりLLVMの世界の雰囲気がわかってきました。 そろそろKaleidoscopeを参考にしながら、実用的なシンタックスを持つ言語のコンパイラを作ってみようと思います。

ソースコード

github.com

bf2llvm.cpp

/*
 * Brainf**k compiler based on LLVM API
 *  $ g++ `llvm-config --cxxflags --ldflags --libs --system-libs` bf2llvm.cpp -o bf2llvm
 *  $ echo "+++++++++[>++++++++>+++++++++++>+++++<<<-]>.>++.+++++++..+++.\
            >-.------------.<++++++++.--------.+++.------.--------.>+." | \
            ./bf2llvm
 *  $ gcc output.o
 *  $ ./a.out
 */
#include <iostream>
#include <string>
#include "llvm/ADT/Optional.h"
#include "llvm/IR/BasicBlock.h"
#include "llvm/IR/IRBuilder.h"
#include "llvm/IR/LegacyPassManager.h"
#include "llvm/IR/Module.h"
#include "llvm/Support/FileSystem.h"
#include "llvm/Support/TargetRegistry.h"
#include "llvm/Support/TargetSelect.h"
#include "llvm/Support/raw_ostream.h"
#include "llvm/Target/TargetMachine.h"
#include "llvm/Target/TargetOptions.h"

static llvm::LLVMContext TheContext;
static llvm::IRBuilder<> Builder(TheContext);
static std::unique_ptr<llvm::Module> TheModule;

void emit_move_ptr(llvm::Value* ptr, int diff) {
  Builder.CreateStore(
      Builder.CreateInBoundsGEP(
        Builder.getInt8Ty(),
        Builder.CreateLoad(ptr),
        Builder.getInt32(diff)),
      ptr);
}

void emit_add(llvm::Value* ptr, int diff) {
  llvm::Value* tmp = Builder.CreateLoad(ptr);
  Builder.CreateStore(
      Builder.CreateAdd(
        Builder.CreateLoad(tmp),
        Builder.getInt8(diff)),
      tmp);
}

void emit_put(llvm::Value* ptr) {
  llvm::Function* funcPutChar = llvm::cast<llvm::Function>(
      TheModule->getOrInsertFunction("putchar",
        Builder.getInt32Ty(),
        Builder.getInt32Ty(),
        nullptr));
  Builder.CreateCall(
      funcPutChar,
      Builder.CreateSExt(
        Builder.CreateLoad(Builder.CreateLoad(ptr)),
        Builder.getInt32Ty()));
}

void emit_get(llvm::Value* ptr) {
  llvm::Function* funcGetChar = llvm::cast<llvm::Function>(
      TheModule->getOrInsertFunction("getchar",
        Builder.getInt32Ty(),
        nullptr));
  Builder.CreateStore(
      Builder.CreateTrunc(
        Builder.CreateCall(funcGetChar),
        Builder.getInt8Ty()),
      Builder.CreateLoad(ptr));
}

struct WhileBlock {
  llvm::BasicBlock* cond_block;
  llvm::BasicBlock* body_block;
  llvm::BasicBlock* end_block;
};

void emit_while_start(llvm::Function* func, llvm::Value* ptr, WhileBlock* while_block, int while_index) {
  while_block->cond_block = llvm::BasicBlock::Create(
      TheContext, std::string("while_cond") + std::to_string(while_index), func);
  while_block->body_block = llvm::BasicBlock::Create(
      TheContext, std::string("while_body") + std::to_string(while_index), func);
  while_block->end_block = llvm::BasicBlock::Create(
      TheContext, std::string("while_end") + std::to_string(while_index), func);
  Builder.CreateBr(while_block->cond_block);
  Builder.SetInsertPoint(while_block->cond_block);
  Builder.CreateCondBr(
      Builder.CreateICmpNE(
        Builder.CreateLoad(Builder.CreateLoad(ptr)),
        Builder.getInt8(0)),
      while_block->body_block,
      while_block->end_block);
  Builder.SetInsertPoint(while_block->body_block);
}

void emit_while_end(WhileBlock* while_block) {
  Builder.CreateBr(while_block->cond_block);
  Builder.SetInsertPoint(while_block->end_block);
}

int main() {
  TheModule = llvm::make_unique<llvm::Module>("top", TheContext);
  llvm::Function* mainFunc = llvm::Function::Create(
      llvm::FunctionType::get(llvm::Type::getInt32Ty(TheContext), false),
      llvm::Function::ExternalLinkage, "main", TheModule.get());
  Builder.SetInsertPoint(llvm::BasicBlock::Create(TheContext, "", mainFunc));

  llvm::Value* data = Builder.CreateAlloca(Builder.getInt8PtrTy(), nullptr, "data");
  llvm::Value* ptr = Builder.CreateAlloca(Builder.getInt8PtrTy(), nullptr, "ptr");
  llvm::Function* funcCalloc = llvm::cast<llvm::Function>(
      TheModule->getOrInsertFunction("calloc",
        Builder.getInt8PtrTy(),
        Builder.getInt64Ty(), Builder.getInt64Ty(),
        nullptr));
  llvm::Value* data_ptr = Builder.CreateCall(funcCalloc, {Builder.getInt64(30000), Builder.getInt64(1)});
  Builder.CreateStore(data_ptr, data);
  Builder.CreateStore(data_ptr, ptr);

  int while_index = 0;
  WhileBlock while_blocks[1000];
  WhileBlock* while_block_ptr = while_blocks;
  char c;
  while (std::cin.get(c)) {
    switch (c) {
      case '>': emit_move_ptr(ptr, 1); break;
      case '<': emit_move_ptr(ptr, -1); break;
      case '+': emit_add(ptr, 1); break;
      case '-': emit_add(ptr, -1); break;
      case '[': emit_while_start(mainFunc, ptr, while_block_ptr++, while_index++); break;
      case ']': if (--while_block_ptr < while_blocks) {
                  std::cerr << "unmatching ]\n";
                  return 1;
                }
                emit_while_end(while_block_ptr); break;
      case '.': emit_put(ptr); break;
      case ',': emit_get(ptr); break;
    }
  }

  llvm::Function* funcFree = llvm::cast<llvm::Function>(
      TheModule->getOrInsertFunction("free",
        Builder.getVoidTy(),
        Builder.getInt8PtrTy(),
        nullptr));
  Builder.CreateCall(funcFree, {Builder.CreateLoad(data)});

  Builder.CreateRet(Builder.getInt32(0));

  llvm::InitializeAllTargetInfos();
  llvm::InitializeAllTargets();
  llvm::InitializeAllTargetMCs();
  llvm::InitializeAllAsmParsers();
  llvm::InitializeAllAsmPrinters();

  std::string TargetTriple = llvm::sys::getDefaultTargetTriple();

  std::string err;
  const llvm::Target* Target = llvm::TargetRegistry::lookupTarget(TargetTriple, err);
  if (!Target) {
    std::cerr << "Failed to lookup target " + TargetTriple + ": " + err;
    return 1;
  }

  llvm::TargetOptions opt;
  llvm::TargetMachine* TheTargetMachine = Target->createTargetMachine(
      TargetTriple, "generic", "", opt, llvm::Optional<llvm::Reloc::Model>());

  TheModule->setTargetTriple(TargetTriple);
  TheModule->setDataLayout(TheTargetMachine->createDataLayout());

  std::string Filename = "output.o";
  std::error_code err_code;
  llvm::raw_fd_ostream dest(Filename, err_code, llvm::sys::fs::F_None);
  if (err_code) {
    std::cerr << "Could not open file: " << err_code.message();
    return 1;
  }

  llvm::legacy::PassManager pass;
  if (TheTargetMachine->addPassesToEmitFile(pass, dest, llvm::TargetMachine::CGFT_ObjectFile)) {
    std::cerr << "TheTargetMachine can't emit a file of this type\n";
    return 1;
  }
  pass.run(*TheModule);
  dest.flush();
  std::cout << "Wrote " << Filename << "\n";

  return 0;
}

LLVMを始めよう! 〜 LLVM IRの基礎はclangが教えてくれた・Brainf**kコンパイラを作ってみよう 〜

コンパイラを作ってみたいと思っていても、アセンブリ言語はよくわからない。 パーサーみたいなコードは書いたことがあるけれど、コード生成の処理はさっぱりだ。 実行ファイルをバイナリエディターで見るとかなにそれ怖い。

そんな私なのですが、LLVMに興味を持ち始めています。 SwiftやRust、あるいはEmscriptenなど、近年注目されている言語やコンパイラ技術の中枢にはLLVMがあります。 アセンブリはよく分からなくてもLLVMを使いこなせるようになれば、マルチプラットフォームで実行ファイルを生成できる言語処理系を作るのではないか。 コンパイラ作ってみたいな、LLVMを使ってみようかなと思っている今日このごろです。

ところが、いざLLVMを勉強しようと思ってもどこから始めればいいかよく分かりませんでした。 マニュアルは巨大で読む気が起きないし、リファレンスを見てもさっぱりです。 雰囲気はわかるんです、インストラクションっていうのはプリミティブな命令のことでしょう、そのリファレンスは便利なんでしょう。 でも素人にはどこから手を付けてよいのかわからない。 LLVMを用いた処理系の入門として名高いKaleidoscopeをやってみようかなと思って始めても、Lexerを作ってParserを作って、「ああ、C++なんて久々に書いたなぁ」なんて言っているうちにLLVM IRにたどり着かずに飽きてしまう。 C++が肌に合わないんだと思ってGo言語のLLVM bindingを使ってコードを書こうとしても、Go bindingのソースコードLLVM本体のコードをあっちこっち飛んで回って疲弊してしまう。 他の人が作った処理系のコードを読んでみようと思ってcloneしてみても、いつの間にか読む暇もなく忘れてしまう。

こんなダメダメな私はいったいどうすればいいのかと思って調べていたのですが、ある日次のようなブログ記事を見つけました。 そうだ、こんな簡単なことを忘れていた。 LLVMを使っているとても身近なコンパイラがあるじゃないか! 私は慌ててこのブログの人と同じことをしてみました (本記事はLLVMがインストールされており、LLVMツールのパスが通っていることを前提とします。私の手元では/usr/local/opt/llvm/bin/でした)。

 $ cat > test.c <<EOF
int main() {
  return 42;
}
EOF
 $ clang -S -emit-llvm -O2 test.c
 $ lli test.ll
 $ echo $?
42

これなら私にもできる! 高ぶる気持ちを抑えて実行ファイル生成まで試してみます。

 $ llc test.ll
 $ gcc test.s -o test
 $ ./test
 $ echo $?
42

実行ファイルを生成し、実行することができました。

clang -S -emit-llvmが生成したファイル、test.llをおそるおそる覗いてみます。

test.ll

; ModuleID = 'test.c'
source_filename = "test.c"
target datalayout = "e-m:o-i64:64-f80:128-n8:16:32:64-S128"
target triple = "x86_64-apple-macosx10.12.0"

; Function Attrs: norecurse nounwind readnone ssp uwtable
define i32 @main() #0 {
  ret i32 42
}

attributes #0 = { norecurse nounwind readnone ssp uwtable "disable-tail-calls"="false" "less-precise-fpmad"="false" "no-frame-pointer-elim"="true" "no-frame-pointer-elim-non-leaf" "no-infs-fp-math"="false" "no-nans-fp-math"="false" "stack-protector-buffer-size"="8" "target-cpu"="penryn" "target-features"="+cx16,+fxsr,+mmx,+sse,+sse2,+sse3,+sse4.1,+ssse3" "unsafe-fp-math"="false" "use-soft-float"="false" }

!llvm.module.flags = !{!0}
!llvm.ident = !{!1}

!0 = !{i32 1, !"PIC Level", i32 2}
!1 = !{!"Apple LLVM version 8.0.0 (clang-800.0.42.1)"}

うーん、これはさっぱりです。 とりあえずコメントなどいらなさそうな物を削ってみて、どれだけ消しても動くか試してみます。

test.ll

source_filename = "test.c"

define i32 @main() #0 {
  ret i32 42
}

attributes #0 = { norecurse nounwind readnone ssp uwtable "disable-tail-calls"="false" "less-precise-fpmad"="false" "no-frame-pointer-elim"="true" "no-frame-pointer-elim-non-leaf" "no-infs-fp-math"="false" "no-nans-fp-math"="false" "stack-protector-buffer-size"="8" "target-cpu"="penryn" "target-features"="+cx16,+fxsr,+mmx,+sse,+sse2,+sse3,+sse4.1,+ssse3" "unsafe-fp-math"="false" "use-soft-float"="false" }
 $ lli test.ll; echo $?
42

動きます。 #0と書いてあるのがよく分からないですが消してみます。

test.ll

define i32 @main() {
  ret i32 42
}
 $ lli test.ll; echo $?
42

動きました。なんだ、source_filenameもなくても動くんですね…

改めて最も小さくなったLLVM IR (Intermediate representation, 中間表現) のコードを眺めてみます。

test.ll

define i32 @main() {
  ret i32 42
}

defineで関数を定義しているのかな、retreturnだろうな、i32は32bit intなんだろなと、容易に想像が付きます。 こんな小さいコードなら、手で書くこともできそうです。 そうです、すでにdefine, ret, i32といったキーワードを理解し始めているのです。

ここまでやったことはとても簡単で、しかもすぐに試せることです。 とても小さなCのコードを書き、clang -S -emit-llvmLLVM IRファイルを生成し、llillcといったコマンドで実行できることを確かめながら、ファイルを削って本当に必要な部分を探しただけです。 LLVMのリファレンスは全く見なくても、Cのどういうコードがどういう命令列になるかはclangが教えてくれるのです。

この方法なら、長いリファレンスや入門チュートリアルを読まなくても、一歩ずつ着実にLLVM IRを学んでいけるのです。

きつねさんでもわかるLLVM ~コンパイラを自作するためのガイドブック~

きつねさんでもわかるLLVM ~コンパイラを自作するためのガイドブック~

基本的な命令を学ぼう

int main()と書くとdefine i32 @main()となることや、return 42ret i32 42となることはわかりました。 しかしこれだけでは何もできないので、いくつか基本的な文を変換してみてLLVM IRがどうなるかを調べてみましょう。

まず、変数の宣言を調べてみます。 次のようなコードを書いてみます。

test.c

int main() {
  char c;
  int i;
  long l;
}

これのLLVM IRを見てみます。-O0で最適化されないように気をつけます。

 $ clang -S -emit-llvm -O0 test.c
 $ cat test.ll

コメントやattributesなどを削除したら次のようになります。

test.ll

define i32 @main() {
  %1 = alloca i8, align 1
  %2 = alloca i32, align 4
  %3 = alloca i64, align 8
  ret i32 0
}

実際このファイルは実行できます。

 $ lli test.ll

一行ずつ見比べると、chari8, inti32, longi64に対応していることがわかります。 alignはアドレスのアラインメントを表しており、変数のアドレスがこの数字の倍数であることが保証されます。

次に代入を調べます。 先程のコードに代入を追加して調べてみます。

test.c

int main() {
  char c;
  int i;
  long l;
  c = 'a';
  i = 72;
  l = 123456789012345;
}

このコードのLLVM IRは次のようになります。

test.ll

define i32 @main() {
  %1 = alloca i8, align 1
  %2 = alloca i32, align 4
  %3 = alloca i64, align 8
  store i8 97, i8* %1, align 1
  store i32 72, i32* %2, align 4
  store i64 123456789012345, i64* %3, align 8
  ret i32 0
}

store命令となりました。 第一引数が代入したい値で、第二引数が代入先のアドレスとなるようです。

少しコードを変えてみます。

test.c

int main() {
  char a, b;
  a = 32;
  b = a;
}

このコードのLLVM IRは次のようになります。

test.ll

define i32 @main() {
  %1 = alloca i8, align 1
  %2 = alloca i8, align 1
  store i8 32, i8* %1, align 1
  %3 = load i8, i8* %1, align 1
  store i8 %3, i8* %2, align 1
  ret i32 0
}

%132を保存し、%1の値をload命令で取り出して%3とし、それを%2storeしています。

足し算や引き算するとどうなるでしょうか。int同士の足し引きを見てみます。

test.c

int main() {
  int a, b, c;
  a = 32;
  b = a + 24;
  c = b - a;
}

test.ll

define i32 @main() {
  ; int a, b, c
  %1 = alloca i32, align 4
  %2 = alloca i32, align 4
  %3 = alloca i32, align 4

  ; a = 32
  store i32 32, i32* %1, align 4

  ; b = a + 24
  %4 = load i32, i32* %1, align 4
  %5 = add nsw i32 %4, 24
  store i32 %5, i32* %2, align 4

  ; c = b - a
  %6 = load i32, i32* %2, align 4
  %7 = load i32, i32* %1, align 4
  %8 = sub nsw i32 %6, %7
  store i32 %8, i32* %3, align 4

  ret i32 0
}

やはり足し算はadd命令, 引き算はsub命令となりました。 nswは理解を後回しして良いでしょう。

後は標準出力が気になります。Hello worldLLVM IRを見てみます。

test.c

#include <stdio.h>
int main() {
  printf("Hello, world!\n");
}

test.ll

@.str = private unnamed_addr constant [15 x i8] c"Hello, world!\0A\00", align 1

define i32 @main() {
  %1 = call i32 (i8*, ...) @printf(i8* getelementptr inbounds ([15 x i8], [15 x i8]* @.str, i32 0, i32 0))
  ret i32 0
}

declare i32 @printf(i8*, ...)

declare命令でprintfを使うことを宣言し、call命令で呼び出しを行っています。 printfではなくてputcharを使ってみます。

test.c

#include <stdio.h>
int main() {
  char c;
  c = 'a';
  putchar(c);
}

test.ll

define i32 @main() {
  %1 = alloca i8, align 1
  store i8 97, i8* %1, align 1
  %2 = load i8, i8* %1, align 1
  %3 = sext i8 %2 to i32
  %4 = call i32 @putchar(i32 %3)
  ret i32 0
}

declare i32 @putchar(i32)

よく見ると、loadしてputcharするだけではなくて、sext命令が入っています。 そういえばputcharの引数はcharではなくてintでしたね。 sextはsigned extendの略で、型の変換を行う命令です。

Brainf**kのコンパイラを作ってみよう

宣言・代入・演算そして関数と、とても基本的なコードがLLVM IRでどのようになるかを観察してきました。 これだけでは実用的なプログラミング言語コンパイラを作るには知識が足りないのですが、プリミティブな言語の処理系なら作れる気がしてきました。

プリミティブな言語…といえばBrainf**kですよね。 私はこの言語が大好きです。 言語処理系のHello worldと言ってもいいのではないでしょうか。 先ほど紹介した記事でも、Brainf**kのコンパイラを書いています。 目的が同じですが、そっくり真似にならないようにちら見程度にして、自分の手で書いてみようと思います。

まずはBrainf**kのデータ領域の用意をします。

test.c

#include <stdlib.h>

int main() {
  char* data = (char*)calloc(30000, sizeof(char));
  char* ptr = data;
  free(data);
}

このLLVM IRを見てみましょう。

test.ll

define i32 @main() {
  %1 = alloca i8*, align 8
  %2 = alloca i8*, align 8

  %3 = call i8* @calloc(i64 30000, i64 1)
  store i8* %3, i8** %1, align 8

  %4 = load i8*, i8** %1, align 8
  store i8* %4, i8** %2, align 8

  %5 = load i8*, i8** %1, align 8
  call void @free(i8* %5)
  ret i32 0
}

declare i8* @calloc(i64, i64)

declare void @free(i8*)

i8*型の変数dataptrを作り、callocの結果を代入しています。

次に、ポインターの移動>とインクリメント+に対応するCのコードを、LLVM IRで見てみます。

test.ll

  ; >   ++ptr
  %5 = load i8*, i8** %2, align 8
  %6 = getelementptr inbounds i8, i8* %5, i32 1
  store i8* %6, i8** %2, align 8

  ; +   ++*ptr
  %7 = load i8*, i8** %2, align 8
  %8 = load i8, i8* %7, align 1
  %9 = add i8 %8, 1
  store i8 %9, i8* %7, align 1

getelementptrという新しい命令が出てきましたね。 意味は後で調べるとして、今はこういうものだと思っておきましょう。 ++*ptrは、ptrの値を%7に入れて、pointer dereferenceして%8に入れて、1足して%9に入れて、それを%7に代入する、と一行ずつ意味を読み取りやすいですね。

あとBrainf**kの.命令に対応するputchar(*ptr);のコードを見ておきます。

test.ll

  ; .  putchar(*ptr)
  %5 = load i8*, i8** %2, align 8
  %6 = load i8, i8* %5, align 1
  %7 = sext i8 %6 to i32
  %8 = call i32 @putchar(i32 %7)

そうでした、sext (signed extend) 命令でi32に変換する必要がありましたね。

Brainf**kの>, +, .に対応するLLVM IRコードが分かりました。<, -のコードも容易に想像がつきます。 ようやく道具が揃ったので、Brainf**kのLLVM IRコンパイラを書いてみます。

bf2llvm.c

#include <stdio.h>
#include <stdlib.h>

void emit_header() {
  printf("define i32 @main() {\n");
  printf("  %%data = alloca i8*, align 8\n");
  printf("  %%ptr = alloca i8*, align 8\n");
  printf("  %%data_ptr = call i8* @calloc(i64 30000, i64 1)\n");
  printf("  store i8* %%data_ptr, i8** %%data, align 8\n");
  printf("  store i8* %%data_ptr, i8** %%ptr, align 8\n");
}

int idx = 1;
void emit_move_ptr(int diff) {
  printf("  %%%d = load i8*, i8** %%ptr, align 8\n", idx);
  printf("  %%%d = getelementptr inbounds i8, i8* %%%d, i32 %d\n", idx + 1, idx, diff);
  printf("  store i8* %%%d, i8** %%ptr, align 8\n", idx + 1);
  idx += 2;
}

void emit_add(int diff) {
  printf("  %%%d = load i8*, i8** %%ptr, align 8\n", idx);
  printf("  %%%d = load i8, i8* %%%d, align 1\n", idx + 1, idx);
  printf("  %%%d = add i8 %%%d, %d\n", idx + 2, idx + 1, diff);
  printf("  store i8 %%%d, i8* %%%d, align 1\n", idx + 2, idx);
  idx += 3;
}

void emit_put() {
  printf("  %%%d = load i8*, i8** %%ptr, align 8\n", idx);
  printf("  %%%d = load i8, i8* %%%d, align 1\n", idx + 1, idx);
  printf("  %%%d = sext i8 %%%d to i32\n", idx + 2, idx + 1);
  printf("  %%%d = call i32 @putchar(i32 %%%d)\n", idx + 3, idx + 2);
  idx += 4;
}

void emit_footer() {
  printf("  %%%d = load i8*, i8** %%data, align 8\n", idx);
  printf("  call void @free(i8* %%%d)\n", idx);
  printf("  ret i32 0\n");
  printf("}\n\n");
  printf("declare i8* @calloc(i64, i64)\n\n");
  printf("declare void @free(i8*)\n\n");
  printf("declare i32 @putchar(i32)\n");
}

int main() {
  char c;
  emit_header();
  while ((c = getchar()) != EOF) {
    switch (c) {
      case '>': emit_move_ptr(1); break;
      case '<': emit_move_ptr(-1); break;
      case '+': emit_add(1); break;
      case '-': emit_add(-1); break;
      case '.': emit_put(); break;
    }
  }
  emit_footer();
  return 0;
}

動かしてみます。

 $ gcc bf2llvm.c -o bf2llvm
 $ echo "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++." | ./bf2llvm
define i32 @main() {
  %data = alloca i8*, align 8
  %ptr = alloca i8*, align 8
  %data_ptr = call i8* @calloc(i64 30000, i64 1)
  store i8* %data_ptr, i8** %data, align 8
  store i8* %data_ptr, i8** %ptr, align 8
  %1 = load i8*, i8** %ptr, align 8
  %2 = load i8, i8* %1, align 1
  %3 = add i8 %2, 1
  store i8 %3, i8* %1, align 1
  ; 略
  store i8 %195, i8* %193, align 1
  %196 = load i8*, i8** %ptr, align 8
  %197 = load i8, i8* %196, align 1
  %198 = sext i8 %197 to i32
  %199 = call i32 @putchar(i32 %198)
  %200 = load i8*, i8** %data, align 8
  call void @free(i8* %200)
  ret i32 0
}

declare i8* @calloc(i64, i64)

declare void @free(i8*)

declare i32 @putchar(i32)
 $ echo "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++." | ./bf2llvm | lli
A

やった! これでは+.しか試してないので他の命令も使ってみます。

 $ echo "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++.+++++++++++++++++++++++++++++.+++++++..+++.>++++++++++++++++++++++++++++++++.<++++++++.--------.+++.------.--------.>+." | ./bf2llvm | lli
Hello world!

動いてそうですね! 最初の峠は超えたようです!

ループに対応しよう

Brainf**kの[, ]に対応しましょう。 まずはCで書いてclangLLVM IRを調べてみます。 [while (*ptr) {に、]}に対応します。 以下のプログラムは、Brainf**kの++[>++<-]>.に対応するコードで、実行したら"\x04"を出力します。

test.c

#include <stdio.h>
#include <stdlib.h>

int main() {
  char* data = (char*)calloc(30000, sizeof(char));;
  char* ptr = data;
  ++*ptr;
  ++*ptr;
  while (*ptr) {
    ++ptr;
    ++*ptr;
    ++*ptr;
    --ptr;
    --*ptr;
  }
  ++ptr;
  putchar(*ptr);
  free(data);
}

LLVM IRは次のようになりました。

test.ll

define i32 @main() {
  %1 = alloca i32, align 4
  %2 = alloca i8*, align 8
  %3 = alloca i8*, align 8
  store i32 0, i32* %1, align 4
  %4 = call i8* @calloc(i64 30000, i64 1)
  ; 略
  store i8 %11, i8* %9, align 1
  br label %12

; <label>:12                                      ; preds = %16, %0
  %13 = load i8*, i8** %3, align 8
  %14 = load i8, i8* %13, align 1
  %15 = icmp ne i8 %14, 0
  br i1 %15, label %16, label %30

; <label>:16                                      ; preds = %12
  %17 = load i8*, i8** %3, align 8
  %18 = getelementptr inbounds i8, i8* %17, i32 1
  store i8* %18, i8** %3, align 8
  ; 略
  store i8 %29, i8* %27, align 1
  br label %12

; <label>:30                                      ; preds = %12
  %31 = load i8*, i8** %3, align 8
  ; 略
  ret i32 %38
}

declare i8* @calloc(i64, i64)

declare i32 @putchar(i32)

declare void @free(i8*)

while文のジャンプにbr (branch) 命令が、条件文に対応する場所はicmp (integer compare) が使われていることがわかります。 ; <label>:12は、見ての通りラベルで、br label %12でそこにジャンプしているのだと想像がつきます。 次のようにラベルを読みやすく書くこともできます。

  br label %while_cond0

while_cond0:
  ; 略
  br i1 %14, label %while_body0, label %while_end0

while_body0:
  ; 略
  br label %while_cond0

while_end0:

while文の出力の仕方が分かったので、コンパイラの続きを書いてみます。

bf2llvm.c

void emit_while_start(int while_index) {
  printf("  br label %%while_cond%d\n", while_index);
  printf("while_cond%d:\n", while_index);
  printf("  %%%d = load i8*, i8** %%ptr, align 8\n", idx);
  printf("  %%%d = load i8, i8* %%%d, align 1\n", idx + 1, idx);
  printf("  %%%d = icmp ne i8 %%%d, 0\n", idx + 2, idx + 1);
  printf("  br i1 %%%d, label %%while_body%d, label %%while_end%d\n", idx + 2, while_index, while_index);
  printf("while_body%d:\n", while_index);
  idx += 3;
}

void emit_while_end(int while_index) {
  printf("  br label %%while_cond%d\n", while_index);
  printf("while_end%d:\n", while_index);
}

/* 略 */
int main() {
  char c;
  int while_index = 0;
  int while_indices[1000];
  int* while_index_ptr = while_indices;
  emit_header();
  while ((c = getchar()) != EOF) {
    switch (c) {
      /* 略 */
      case '[': emit_while_start(*while_index_ptr++ = while_index++); break;
      case ']': emit_while_end(*--while_index_ptr); break;
    }
  }
  emit_footer();
  return 0;
}

試してみましょう。

 $ gcc bf2llvm.c -o bf2llvm
 $ echo "+++++++++[>++++++++>+++++++++++>+++++<<<-]>.>++.+++++++..+++.>-.------------.<++++++++.--------.+++.------.--------.>+." | ./bf2llvm | lli
Hello, world!

動いてそうですね! もう少しややこしいコードを動かしてみます。

 $ cat fibonacci.bf
>>+++++++++++[-<<++++>+++>>+<]>>+<++<<->>[>>>++++++++++<<[->+>-[>+>>]>[+[-<+>]>+
>>]<<<<<<]>>[-]>>>++++++++++<[->-[>+>>]>[+[-<+>]>+>>]<<<<<]>[-]>>[>++++++[-<++++
++++>]<.<<+>+>[-]]<[<[->-<]++++++[->++++++++<]>.[-]]<<++++++[-<++++++++>]<.[-]<<
[-<+>]<<[->>+>+<<<]>>>[-<<<+>>>]<-[<<<<<.>.>>>>[-]]<[->+>+<<]>>[-<<+>>]<<<<[->>+
<<]>>>[-<<<+>>>]<<-]
 $ cat fibonacci.bf | ./bf2llvm | lli
1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233
 $ cat prime.bf
++++++++++[->+++++>++++>+++<<<]>>++++>++>>+>+++<<<<<.>>>>>[<[-]++<[-]+>>>+[<[->>
+<<]<[->>>>+>+<<<<<]>>>>[-<<<<+>>>>]<[->+>-[>+>>]>[+[-<+>]>+>>]<<<<<<]>[-<<<+>>>
]>[-]>[-<<+>+>]<<[->>+<<]+>[<[-]>[-]]<[<<<<<[-]>>>>>[-]]>[-]>[-]>[-]<<<<<<-<[->>
>+>+<<<<]>>>>[-<<<<+>>>>]<<<[->>>>+>+<<<<<]>>>>>[-<<<<<+>>>>>]<<<+[->>[->+>+<<]>
>[-<<+>>]+<[>-<<->[-]]>[<<<+<[-]>>>>[-]]<[-]<<<]>>[-]<[<<->>[-]]<[-]<<+<[->>>+>+
<<<<]>>>>[-<<<<+>>>>]>+++++++++++++<<[->>[->+>+<<]>>[-<<+>>]+<[>-<<->[-]]>[<<<+<
[-]>>>>[-]]<[-]<<<]>>[-]<[<<[-[-]]>>[-]]<[-]<<<+>>]<<<[<.>>>>>++++++++++<<[->+>-
[>+>>]>[+[-<+>]>+>>]<<<<<<]>[-<+>]>[-]>>>++++++++++<[->-[>+>>]>[+[-<+>]>+>>]<<<<
<]>[-]>>[>++++++[-<++++++++>]<.<<+>+>[-]]<[<[->-<]++++++[->++++++++<]>.[-]]>[-]<
<<++++++[-<++++++++>]<.[-]<<<<<[-]]>>+]
 $ cat prime.bf | ./bf2llvm | lli
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251

Erik Bosmanさんという方が書いた、マンデルブロ集合のプログラムも動かすことができます。 出力は読者が確認してください。

基本構文に対応できてそれらを使った小さなコードが動いてしまえば、あらゆるプログラムが動く (とは言えメモリーの限界はありますが) というのはコンパイラの分野でコード書いていて楽しい所です。

最適化

LLVMにはLLVM IRレベルの命令列を最適化する機能があります。 LLVMをバックエンドに持つコンパイル言語は、LLVMの最適化の恩恵を受けられます。

LLVM IR命令列がどのように最適化されるかを調べるには、optコマンドを使うとよいでしょう。

 $ cat hello.bf 
+++++++++[>++++++++>+++++++++++>+++++<<<-]>.>++.+++++++..+++.>-.------------.<++++++++.--------.+++.------.--------.>+.
 $ cat hello.bf | ./bf2llvm | opt -S -O3
; ModuleID = '<stdin>'
source_filename = "<stdin>"

; Function Attrs: nounwind
define i32 @main() local_unnamed_addr #0 {
while_end0:
  %0 = tail call i32 @putchar(i32 72)
  %1 = tail call i32 @putchar(i32 101)
  %2 = tail call i32 @putchar(i32 108)
  %3 = tail call i32 @putchar(i32 108)
  %4 = tail call i32 @putchar(i32 111)
  %5 = tail call i32 @putchar(i32 44)
  %6 = tail call i32 @putchar(i32 32)
  %7 = tail call i32 @putchar(i32 119)
  %8 = tail call i32 @putchar(i32 111)
  %9 = tail call i32 @putchar(i32 114)
  %10 = tail call i32 @putchar(i32 108)
  %11 = tail call i32 @putchar(i32 100)
  %12 = tail call i32 @putchar(i32 33)
  ret i32 0
}

; Function Attrs: nounwind
declare i32 @putchar(i32) local_unnamed_addr #0

attributes #0 = { nounwind }

なんと、IO以外は全て計算されてしまいました! もともと配列を確保してその上で計算していたはずが、最適化によってそれも消えてしまいました。

もう少し複雑な、マンデルブロ集合のスクリプトで実行時間を比較してみましょう。

 $ time sh -c 'cat mandelbrot.b | ./bf2llvm | lli > /dev/null'
        6.34 real         6.27 user         0.04 sys
 $ time sh -c 'cat mandelbrot.b | ./bf2llvm | opt -S -O3 | lli > /dev/null'
        3.55 real         3.50 user         0.05 sys

LLVMの最適化によって約40%も高速化しました (Hello worldほど単純ではないので、配列ごと消えるということはありません)。 コンパイルして実行ファイルを作ってみます。

 $ cat mandelbrot.b | ./bf2llvm | opt -S -O3 > mandelbrot.ll
 $ llc mandelbrot.ll
 $ gcc mandelbrot.s -o mandelbrot
 $ time ./mandelbrot > /dev/null
        0.98 real         0.97 user         0.00 sys

はい。すごく速いですね。 手でBrainf**kのインタープリタを書いて最適化しようとしても、小手先の改善の組み合わせでLLVMレベルの最適化するのは難しいでしょう。 まずは素朴に命令列を吐いて、最適化をおまかせしてしまえば速度が出るとてもうれしいですね。 私のようにアセンブリの分からない人でも、LLVM IRを正しく吐くことさえできれば実行ファイル生成まで面倒見てくれるのはLLVMの素晴らしいところですね。

まとめ

LLVMの最初の一歩はLLVM IRの命令と親しむところにあります。 しかし多くのチュートリアルでは、難しい言語をC++でパースして構文木を作ってIRを吐いてという説明から始まります。 できる人にはできるのでしょうが、私には理解が追いつきませんでした。 LLVM IRを初めて触った人が、いきなりIRBuilderからコード生成できるのでしょうか。 LLVM素人がいきなりIRBuilderのメソッド一覧を眺めて、これを使えばいいなみたいにすぐに分かるものなのでしょうか。

そんな人にとって、clangは良い先生だと思います。 CのコードがどのようなLLVM IRの命令列になるか教えてくれます。 CとLLVM IRを見比べながら、自分でも命令列を吐くコードを書いてみて、徐々に分かっていくものだと思います。 それはIRBuilderを使ったコードに比べたら泥臭いコードかもしれませんが、入門の壁を超える大事な一歩だと思います。

この記事は、実際に私自身が一歩ずつ理解を深めながら、同時並行して書き進めてきました。 私はLLVMについてはド素人で、Kaleidoscopeレベルにすら到達していません。

そんな私だからこそ、この記事を書こうと思った時には今の私に書ける「Kaleidoscopeの手前のLLVM入門記事にしよう」という思いがありました。 LLVMを触ってみたいけど、最初の壁が超えられない、そういった人に届けばいいなと思います。

いくつかの命令についてあまり説明せずに通り過ぎたところがありました。 命令の正確な意味や引数などの理解を深めるためには、言語マニュアルを参照してください。 なお、用語や記述には気をつけて書いていますが、誤りがあればご指摘ください。

LLVMのページを開いた時、命令の説明書をぼーっと眺めている時、このツールは私に使えるものなのだろうかという不安がありました。 今の気持ちは全く違います。 実際に触って試行錯誤したら、コツが分かってきました。 ようやくLLVMの世界の入り口に立った気持ちです。 IRの命令の雰囲気が掴めてきたので、次はIRBuilderを使ったコード生成や、実用的な言語のコンパイラ作成に挑戦してみようと思います。 今後も理解を深め、コンパイラ技術の勉強を続けていこうと思います。

ソースコード

bf2llvm.c

/*
 * Brainf**k -> LLVM IR Compiler
 *  $ gcc bf2llvm.c -o bf2llvm
 *  $ echo "+++++++++[>++++++++>+++++++++++>+++++<<<-]>.>++.+++++++..+++.\
            >-.------------.<++++++++.--------.+++.------.--------.>+." | \
            ./bf2llvm | opt -S -O3 | lli
 */
#include <stdio.h>
#include <stdlib.h>

void emit_header() {
  printf("define i32 @main() {\n");
  printf("  %%data = alloca i8*, align 8\n");
  printf("  %%ptr = alloca i8*, align 8\n");
  printf("  %%data_ptr = call i8* @calloc(i64 30000, i64 1)\n");
  printf("  store i8* %%data_ptr, i8** %%data, align 8\n");
  printf("  store i8* %%data_ptr, i8** %%ptr, align 8\n");
}

int idx = 1;
void emit_move_ptr(int diff) {
  printf("  %%%d = load i8*, i8** %%ptr, align 8\n", idx);
  printf("  %%%d = getelementptr inbounds i8, i8* %%%d, i32 %d\n", idx + 1, idx, diff);
  printf("  store i8* %%%d, i8** %%ptr, align 8\n", idx + 1);
  idx += 2;
}

void emit_add(int diff) {
  printf("  %%%d = load i8*, i8** %%ptr, align 8\n", idx);
  printf("  %%%d = load i8, i8* %%%d, align 1\n", idx + 1, idx);
  printf("  %%%d = add i8 %%%d, %d\n", idx + 2, idx + 1, diff);
  printf("  store i8 %%%d, i8* %%%d, align 1\n", idx + 2, idx);
  idx += 3;
}

void emit_put() {
  printf("  %%%d = load i8*, i8** %%ptr, align 8\n", idx);
  printf("  %%%d = load i8, i8* %%%d, align 1\n", idx + 1, idx);
  printf("  %%%d = sext i8 %%%d to i32\n", idx + 2, idx + 1);
  printf("  %%%d = call i32 @putchar(i32 %%%d)\n", idx + 3, idx + 2);
  idx += 4;
}

void emit_get() {
  printf("  %%%d = call i32 @getchar()\n", idx);
  printf("  %%%d = trunc i32 %%%d to i8\n", idx + 1, idx);
  printf("  %%%d = load i8*, i8** %%ptr, align 8\n", idx + 2);
  printf("  store i8 %%%d, i8* %%%d, align 1\n", idx + 1, idx + 2);
  idx += 3;
}

void emit_while_start(int while_index) {
  printf("  br label %%while_cond%d\n", while_index);
  printf("while_cond%d:\n", while_index);
  printf("  %%%d = load i8*, i8** %%ptr, align 8\n", idx);
  printf("  %%%d = load i8, i8* %%%d, align 1\n", idx + 1, idx);
  printf("  %%%d = icmp ne i8 %%%d, 0\n", idx + 2, idx + 1);
  printf("  br i1 %%%d, label %%while_body%d, label %%while_end%d\n", idx + 2, while_index, while_index);
  printf("while_body%d:\n", while_index);
  idx += 3;
}

void emit_while_end(int while_index) {
  printf("  br label %%while_cond%d\n", while_index);
  printf("while_end%d:\n", while_index);
}

void emit_footer() {
  printf("  %%%d = load i8*, i8** %%data, align 8\n", idx);
  printf("  call void @free(i8* %%%d)\n", idx);
  printf("  ret i32 0\n");
  printf("}\n\n");
  printf("declare i8* @calloc(i64, i64)\n\n");
  printf("declare void @free(i8*)\n\n");
  printf("declare i32 @putchar(i32)\n\n");
  printf("declare i32 @getchar()\n");
}

int main() {
  char c;
  int while_index = 0;
  int while_indices[1000];
  int* while_index_ptr = while_indices;
  emit_header();
  while ((c = getchar()) != EOF) {
    switch (c) {
      case '>': emit_move_ptr(1); break;
      case '<': emit_move_ptr(-1); break;
      case '+': emit_add(1); break;
      case '-': emit_add(-1); break;
      case '[': emit_while_start(*while_index_ptr++ = while_index++); break;
      case ']': emit_while_end(*--while_index_ptr); break;
      case '.': emit_put(); break;
      case ',': emit_get(); break;
    }
  }
  emit_footer();
  return 0;
}

二週間で簡単なインタープリタ言語を実装してみた (日記)

私は昔から言語処理系に興味があり、自分で言語を作ることを夢見てきました。 しかし、いざ言語を実装しようと思って言語処理系に関する本を読んでも何から手を付けていいか分からず、アセンブラもまともに読めないまま、数年が経ってしまいました。 大学時代は情報系ではなかったため、コンパイラの実験がある情報系の学科のカリキュラムを羨ましく思い、情報系の授業の教科書を手にとって見ても読む気が起きず、自分に作れるのは所詮、構文木をちょこっといじって変換するレベルのもの (例えばsjspなど) にとどまっていました。

そんな中、去年のRebuild.fmで、とても感銘を受けた回がありました。 LLVMのlinkerであるLLDを開発されているrui314さんの回です。 rebuild.fm セルフコンパイルできるC言語コンパイラを実装するという話のなかで、インクリメンタルに開発する重要性について話をされています。

qiita.com

この回を聞いてやっと、コンパイラの書き方がわかった気がしました。 コンパイラに関する本を読んでも実装の進め方がわからなかったのは、コンパイルフェーズごとに章に分かれていたからで、正規言語の難しい話を読んでは本を投げ出していたのです。 Rebuild.fmのrui314さんのお話を聞いて、まず最初は小さい言語をとにかくパースからコード生成まで通して動くようにして、そこから徐々に文法を増やしていくというインクリメンタルなアプローチが有力なのだとやっとわかったのです。 ここで最初の言語は、関数もforループも演算子もない、とてもプリミティブな文法だけもつ言語から始めるのです。 コンパイラの教科書を読んでもさっぱりわからなかった実装の進め方が、ようやくわかった瞬間でした。

最近、まつもとゆきひろさんの『言語のしくみ』が出版されたこともあり、言語実装の熱が高まっているように感じます。 私は今年は言語実装の年にしようと思っていて、その第一歩としてインタープリタ言語を書いてみました。 まずやったことは、mrubyやLuaなどの実装を手元に引っ張ってきて、ctagsを打って数日読んでみました。

Luaのファイル構成は分かりやすいので読みやすいです。mrubyはmrbgems/mruby-compiler/coreやsrc/vm.cなどを中心に見ていきました。ざっくり雰囲気を掴むくらいの気持ちで読みました。 あとはひたすら実装を進めていきました。

つまり参考にしたものはmrubyとLuaです。 書籍としては以下の本を挙げておきます (もうちょっとやわらかいインタープリタ実装の本ある気がするけど知らないです)。

大学時代に買ったけど当時はあまり理解できず、最近開いてみてようやく理解し始めたコンパイラの本。 この本だけ読んでてもよく分からないけど、自分でVMを書いてみてから改めてこの本を読んでみたら、納得する部分がたくさんありました (特に関数のコード生成について)。 大学で使われている教科書で少し硬い。

情報系教科書シリーズ コンパイラ

情報系教科書シリーズ コンパイラ

Direct threadingなどmrubyで使われている技術がよく分かる本 (すみません、最初の方をちら見しただけでちゃんとは読んではないです…)。

rui314先生をリスペクトしているので、私も日記形式でお送りしてみます。

1月3日

実家から戻り、インタープリタ言語を作る意識を高める。 ディレクトリを作り、yaccファイルを用意する。

1月7日

mrubyの構文木の実装を参考にして、ASTを全てnode*で実装することにする。

typedef struct node {
  struct node *car, *cdr;
} node;

とりあえずノードを追加するたびにmallocするような形で実装を進める。 文法は四則演算のみ。

1月8日

構文木のためのメモリープールを実装する。 ASTは全てnode*で表すことにしたので、この構造体のみ考えておけば良いのは楽。

四則演算のみの言語に対して、コード生成と実行器を実装する。 実装が簡単なスタックマシンで作る。

今ある文法は

  • 式 (expression) が並んだものがprogram
  • expressionは四則演算が使える
  • リテラルは整数 (long) または浮動小数点数 (double)
enum OP_CODE {
  OP_ADD,
  OP_MINUS,
  OP_TIMES,
  OP_DIVIDE,
  OP_LOAD_LONG,
  OP_LOAD_DOUBLE,
  OP_PRINT_POP,
};

OP_PRINT_POPはちょっとダサい…

1月9日

コードが読みやすいようにいくつか変数名をリネーム。 変数を作り始めるが難しくて散歩に出かける。

夜、変数を実装。新たに追加した文法は

  • 変数に代入できる
  • 四則演算のなかで変数を使える
  • print文
statement         : IDENTIFIER EQ expression
                    {
                      $$ = cons(nint(NODE_ASSIGN), cons($1, $3));
                    }
                  | PRINT expression
                    {
                      $$ = cons(nint(NODE_PRINT), $2);
                    }
                  ;

OP_PRINT_POPを削除してOP_ASSIGN, OP_PRINT, OP_LOAD_IDENTを追加。

 enum OP_CODE {
+  OP_ASSIGN,
+  OP_PRINT,
   OP_ADD,
   OP_MINUS,
   OP_TIMES,
   OP_DIVIDE,
   OP_LOAD_LONG,
   OP_LOAD_DOUBLE,
-  OP_PRINT_POP,
+  OP_LOAD_IDENT,
 };

OP_ASSIGNOP_LOAD_IDENTのoperandは、変数配列のindexが入っているので、実行器の雰囲気はこんな感じ。

+      case OP_ASSIGN:
+        e->variables[GET_ARG_A(e->codes[i])].value = e->stack[--e->stackidx];
+        break;
       case OP_ADD: BINARY_OP(+); break;
@@ -172,6 +214,9 @@ static void execute_codes(env* e) {
         break;
+      case OP_LOAD_IDENT:
+        e->stack[e->stackidx++] = e->variables[offset - GET_ARG_A(e->codes[i])].value;
+        break;

booleanリテラルを実装する。まだ演算はない。

1月11日

簡単なテストケースを書く。

foo = 3 * 4 / 5 + 6.7 * 7
bar = foo * 3.1 + foo / 7.2
print foo + bar

このスクリプトを実行すると207.281666667と表示される。 Python 2のコンソールで確かめながら実行結果を確かめる。

if文を実装する。Vim scriptに倣って、if … endifにする。 if文の実装に伴って、OP_JMP_NOTを作る。 各ステートメントの命令数をカウントしないといけないことに気が付き、codegen関数の返り値をvoidからuint16_tにする。

1月12日

else文を実装する。 無条件にジャンプするOP_JMPを作る。

if true
  a = 10
  print a
else
  b = 11
  print b
endif

このスクリプトは次のようなコードに変換されて実行される。

bool true
jmp_ifnot 5
long 10
let 0
load 0
print
jmp 4
long 11
let 1
load 1
print

for文の中にprogram counterをincrementするところがあるので、jmp_ifnotのoperandが1少ないように見えるがこれで正しい。

elseif文を実装する。 構文木を構築するときに新しいif文にしてしまうことで、コード生成を触らなくても実装できてしまうことに気がつく。

デバッグしやすいように、構文木LISPのような形式で表示できるようにする。

1月13日

while文を実装する。 これまでoperandはuint16_tだと思っていたが、これではジャンプ命令で負の値を指定できないという初歩的なミスに気がつく。 とりあえず暫定対処でOP_JMP_BACKを追加する。

>=, ==, <= など、数値の比較を行う演算子を実装する。

&&, || を実装する。これらは四則演算などの他の二項演算子とは異なり、右辺を評価しないことがある。 if文を作ったときに追加したOP_JMP_NOTは評価値をスタックからpopしてしまうため、仕方なしにOP_JMP_NOT_KEEPを作る。

if文とwhile文が実装されたので、なんとなく「言語らしさ」が出てきてかわいい。

1月14日

operandをuint16_tからint16_tに変更し、OP_JMP_BACKを削除する。 少しリファクタリングする。

最近は仕事が忙しくて進められていない。

1月16日

組み込み関数min(), max()を実装する。 関数の引数のASTを組み立てていて、ようやくnode*の良さが分かってくる。 関数の実装はべた書きでスタックを直接操作していて危ない。 関数呼び出しのoperandには、何番目の関数かというindexと、引数の数を指定するようにする。 ユーザー定義の関数の実装を妄想する。

1月17日

単項演算子 +, - を実装する。 組み込み関数abs()を実装する。

1月18日

endifとendwhileをendにする。RubyLuaっぽくてかわいい。

a = 10
while a >= 0
  print a
  if a > 5
    print 10
  else
    print 0
  end
  a = a - 1
end

ユーザー定義関数をとりあえずパースできるようにする。

+statement         : FUNC IDENTIFIER LPAREN fargs_opt RPAREN sep statements sep END
+                    {
+                      $$ = cons(nint(NODE_FUNCTION), cons($2, cons($4, $7)));
+                    }
+                  | RETURN expression
+                    {
+                      $$ = cons(nint(NODE_RETURN), $2);
+                    }

まだコード生成はどうすれば良いのかさっぱりわからない。

ファイルが大きくなってきたので分割する。

1月19日

simple virtual machineという意味でsvmと名付けていたが、support vector machineと紛らわしいので、minivmにrenameする。

ユーザー定義関数の実装に想いを馳せる。

1月21日

ユーザー定義関数のコード生成と実行器を実装する。 program counterをsave・unsaveするOP_SAVEPC, OP_UNSAVEPC、ローカル変数の領域を確保と解放を行うOP_ALLOC, OP_UNALLOCを実装する。 少し汚いけど、実行時には変数領域を逆さまに使うようにすると、確保した数を持たなくてよいので楽ということに気がつく。 関数の実行は、おおむね次のような流れになる。

  • 引数を順番に評価してスタックに積む
  • 現在のprogram counter (呼び出し前のpc) をスタックに積む
  • 関数呼び出しを行う
  • ローカル変数のための変数領域を確保
  • popして呼び出し前のpcをローカル変数に入れる
  • 引数をpopして後ろからローカル変数に入れていく
  • 関数の中身を評価する
  • return文で、ローカル変数の変数領域を解放し、呼び出し前のpcに戻す

しかし、関数の中からグローバル変数が見れないとか、再帰呼出しができない (関数自体はグローバル変数領域に確保されるが、関数の中身のコード生成時にローカル変数リストを見ているため) という問題に直面し、頭を抱える。昼寝する。

単項演算子 ! (not) を実装する。

1月22日

グローバル変数とローカル変数は全く性質が違うと悟る。 湯淺先生の本の関数のコード生成の節を読む。 ローカル変数に対するオペコード OP_LET_LOCAL, OP_LOAD_LOCAL_IDENT を追加し、グローバル変数とは区別するように実装し、ようやく関数の中からグローバル変数を参照したり、再帰呼び出しできるようになる。

a = 10
b = 20
func foo(x)
  b = 30
  return a * b + x
end

print foo(50)
print a
print b

このスクリプトが350, 10, 20と表示するようになる (昨日の時点ではグローバル変数領域を見れてなかったので、Undefined variable: aだった)。

フィボナッチ数を計算する (効率の悪い再帰呼び出しの) スクリプトが動き、とても興奮する。テストケースに追加する。

OP_UNALLOCOP_UNSAVEPCを削除してOP_RETに統一する。ローカル変数領域を明け渡すとともに、program counterを呼び出し元に戻す。

break文とcontinue文を実装する。 continueは簡単なのでまあいいとして、break文のジャンプ先がわからずに小一時間悩む。 ワンパスで生成までやっているので、while文の中身を生成しているときは飛び先のpcが分からない。 仕方なく、while文の先頭にwhile文の最後まで飛ぶジャンプ命令を追加する。

スタックの一番上をコピーして積むOP_DUPを追加して、 &&, ||を作ったときに追加したOP_JMP_IF_KEEPOP_JMP_NOT_KEEPを削除する。

パフォーマンス改善のために、即値を足したり引いたりするOP_IADDOP_IMINUSを追加する。 足し算は可換なので即値が左のときも適用できるはずだけど後回し。

そこそこコードを書けるようになったので、軽くスクリプトを書いてみたり速度比較をしてみる。 以下のスクリプトが動く。

func fib(n)
  if n <= 1
    return 1
  end
  return fib(n - 1) + fib(n - 2)
end

n = 0
while n < 38
  print fib(n)
  n = n + 1
end

このスクリプトを9.05sで実行できた。 これと全く同じことを行うスクリプトruby 2.1.9で9.10s, mruby 1.2.0で14.78s, Lua 5.2.4で12.83s, Python 2.7.12で26.30sかかった。 スクリプト言語としてそこそこの速度が出ている事がわかる。

なお、生成された中間コードは次のような感じになった。

long 3
let 0
jmp 20
ret 2
alloc 2
let_local 1
let_local 0
load_local 0
long 1
<=
jmp_ifnot 2
long 1
jmp -10
load_local 0
iminus 1
call 0 1
load_local 0
iminus 2
call 0 1
+
jmp -18
long 0
jmp -20
long 0
let 1
jmp 1
jmp 11
load 1
long 38
<
jmp_ifnot 7
load 1
call 0 1
print
load 1
iadd 1
let 1
jmp -11

こうやって改めて見てみると、やはり関数の最初にretがあったりwhile文の最初のjmpとかして意味不明かもしれない。 returnやbreakなど、生成時に本来飛びたいジャンプ先がわからずにこうなってしまった。 ワンパスで生成するのにこだわらなかったらもう少し素直なコードになるはずだけど、今回は深追いしないことにする。

これからやること

言語処理系としてやらないといけないことは、まだまだたくさんあります。 文字列や配列、辞書なども作りたいし、型とか関数引数の数のチェックなどをスキップしているので、誤った文法のスクリプトを流すと簡単にセグフォしてしまいます… 流石にそういうのは直したい。 構文エラー時のメッセージも雑すぎるのでなんとかしたいと思っています。 今の実装だとおそらくクロージャが作れないので、もう少しスコープをきちんと扱えるようになってから考えたいです。

それっぽく動くものはできたけど色々と気に入らない部分はあるので、今回の実装はいったん捨てて、また1から書いていくと思います。 参考程度にGitHubリポジトリを置いておきます。

まとめ

二週間で初めてのインタープリタ言語を実装できたので、今年中にあと23個くらい言語処理系を書くと思います。 今回はスタックマシンで作りましたが、レジスタマシンのコード生成はどうすればいいのかよくわからないし、並行処理のコード生成とか検討もつきません。 LLVMバックエンドを使ってみたり、アセンブラを吐いてみたり、GCを実装してみたり、あるいはgoroutineスケジューラの実装も読んでみたりと、やりたい事が山積みです。 Scalaのwebアプリケーションを運用している以上はJVMについて詳しくなりたいという思いもあります。

その最初のステップとして、簡単なインタープリタ言語を実装してみました。 とりあえず最初の一歩を踏み出せてよかったです。 小さい構文から始めて、少しずつ構文を足していくというやり方は、言語処理系を作る上で有用だと思います。 よく考えてみればそりゃそうかという方法ですが、rui314さんのお話を始めて聞いたときは目から鱗が落ちる思いがしました。

ASTにmruby方式 (LISPのようにcar cdrで表す) を採用したことは、良い面と悪い面がありました。 メモリー管理は1つの構造体について集中しておけばいいので楽です。 関数の引数部分のように、まさにリスト構造になっている部分は配列にするよりも、consで表すほうが楽だと思います。 悪い面としては、コード生成時にcarやcdrが何を表しているのかよく分からなくなるという問題があります。 パーサーの対応する部分を見れば分かるのですが、生成部分のコードだけを見ていたら何も分かりません。 次に書く言語では、この方式は使わずに実装してみたいと思います。

言語実装の楽しみは、バグっていたら全く動かないし、きちんと実装できたら、その文法で書けるあらゆるコード、あらゆるアルゴリズムが動くようになるというところにあると思います。あるところまで書いたら可能性が一気に広がる楽しさというのは昔Tweetしたことがあります。 この思いは今でも変わっていません。

言語処理系の分野は広大で、いくらでも学ぶ楽しみが広がっています。 いくらやっても学ぶことがあり、しかも実装してみて動くとすごく楽しいというのは、プログラマの趣味としてはうってつけではないでしょうか (趣味レベルに作る程度だと、本業で研究されているレベルまではいくらやっても到達しないままかもしれませんが…)。 実装が公開されている言語もいくつもあり、素晴らしい教材がたくさん見つかります。 第一歩として、スタックマシンのインタープリタ言語を実装してみました。 レジスタマシンやネイティブコード生成、LLVMバックエンドなど、色々な方式を試しながら、言語処理系作りを楽しみたいと思います。