Haskellで10を作るプログラムを書いてみたので動画で公開してみた

最近Rui Ueyamaさんがコーディング動画をアップされているのを見て、私も動画を撮りたくなりました。題材をしばらく考えていたんですが、10を作るプログラムを書いてみることにしました。

www.youtube.com

後から見直すと色々ミスっていて、わりと焦っていることがわかります。なにかの癖で適当に bc -l とかやったのだけど、音声をあてる時は関係ないオプションだと勘違いしてしまいました。確かにglobされていたのはよくなかったけど、 echo '5 / (5 / (5 + 5))' | bc -l とかで考えてみると -l も必要なんですよね。2つの問題が起きていて混乱してしまった…

もともとはProject Euler 93を昔解いたことがあったので、こういう系の問題と分数の扱い方とか括弧の付け方みたいなところはイメージありました。ただ、Project Eulerの問題は式木を表示する必要がなかったので、式の型を作ってなかったんですよね。それで今回改めてコードを書いていたらなんか意外と色んなところでハマってしまいました。

最近のHaskellは少し雲行きが怪しいなという印象がありますね。ここ10年で出てきた新しい言語より勢いがないのは否定できないし、若者があまり興味を示さないのも分かる気がします。私も最近はRustばかり触っているのですが… なんかHaskellの楽しさを直接伝えるものがあまりないんですよね。ブログエントリーのレベルも二分化されてしまっている印象があります。

私は昔、Project EulerCodeforcesの問題を順番に解いていた時期があって、Haskellを書いていて楽しかったなぁと思い出すのはその頃の思い出なのです。まぁ数学のパズルが好きじゃないとささらないのかもしれませんが…

Haskellerの皆さんはどういうところに一番関心がありますか?Haskellを書いていて楽しい時はどういう瞬間でしょうか。若い世代に興味を持ってもらうにはどういうものがあると良いのでしょうか。

追記: 括弧の省略アルゴリズム

括弧を省略しようとして動画では失敗してしまいました。改めて落ち着いて括弧が必要なケースを調べてみると、以下のようになりました。

  • 左辺の括弧が必要なケース: 演算子が乗算か除算で、かつ左辺が加算か減算かのどちらか
  • 右辺の括弧が必要なケース: 演算子が減算か乗算で、かつ右辺が加算か減算、もしくは演算子が除算で右辺が二項演算子

というわけで、次の実装で大丈夫だと思います。

isAddOrSub :: Expr -> Bool
isAddOrSub (NumInt _) = False
isAddOrSub (BinOp (Op op _) _ _) = op == "+" || op == "-"

isBinOp :: Expr -> Bool
isBinOp (NumInt _) = False
isBinOp (BinOp _ _ _) = True

instance Show Expr where
  show (NumInt n) = show n
  show (BinOp op@(Op opstr _) lhs rhs) = lparen (show lhs) ++ " " ++ show op ++ " " ++ rparen (show rhs)
    where lparen | (opstr == "*" || opstr == "/") && isAddOrSub lhs = \str -> "(" ++ str ++ ")"
                 | otherwise = id
          rparen | (opstr == "-" || opstr == "*") && isAddOrSub rhs || opstr == "/" && isBinOp rhs = \str -> "(" ++ str ++ ")"
                 | otherwise = id

後は自明な左辺右辺の入れ替えの除去くらいですかね。いやはやめんどくさそうだ…

Haskellの普及について色々と考えていたんですが、中高生にプログラマーってすごい、なりたいって思わせるのにライブコーディング動画って有用なんじゃないかと思いました。どうせ若い頃に興味を持つ言語はころころ変わるので、あまりそこにこだわってもよくなくて、むしろもう少し若い世代にプログラミングにどう興味をもってもらうかというところに力を入れるとソフトウェア産業も発展していくのではないでしょうか… うまい問題設定のライブコーディング動画って子供にとって結構刺さると思うんですよねぇ。ここまで考えたところでテトリス一時間で実装動画を思い出しました。これは強烈でしたねぇ。

gitのファイル変更日時をファイルのアクセス日時に設定

普段使っているファイラーはファイルのアクセス日時でソートされるように設定しています。大きめのリポジトリをcloneしてコードを読む時に、意外とファイルの最終変更日時が参考になったりします。仕事で使うリポジトリや、定期的にpullしているなら、徐々に変更のないファイルはファイラーの下の方に移動していく (上の方からアクセス日時の降順として) のですが、cloneしたばかりだとこうは行きません。

要はgitリポジトリ内の各ファイルのアクセス日時を、そのファイルのgit履歴上での最終変更日時に戻したいという気持ちになるわけです。そうするとファイラー上でもいい感じにファイルがソートされるのです。

#!/usr/bin/env bash

if command -v gdate >/dev/null 2>&1; then
  DATE=gdate
else
  DATE=date
fi

while IFS= read -r -d '' file; do
  logtime=$(git log -1 --format=%ci "$file")
  mtime=$("$DATE" -d "$logtime" +%Y%m%d%H%M.%S)
  echo "$(printf "%-80s" "$file")"$'\t'"$logtime"$'\t'"$mtime"
  touch -t "$mtime" "$file"
done < <(git ls-files -z)

実行してみるとこういう感じ。

 $ git-touch
.ctags                                                                              2015-08-28 10:47:55 +0900   201508281047.55
.tmux.conf                                                                          2016-02-05 08:43:46 +0900   201602050843.46
.vimrc                                                                              2017-05-01 18:38:35 +0900   201705011838.35
.vimshrc                                                                            2017-03-10 20:45:22 +0900   201703102045.22
.zshrc                                                                              2017-04-20 00:32:39 +0900   201704200032.39
 $ stat -x .zshrc
  File: ".zshrc"
  Size: 10518        FileType: Regular File
  Mode: (0644/-rw-r--r--)         Uid: (  501/ itchyny)  Gid: (   20/   staff)
Device: 1,4   Inode: 203408150    Links: 1
Access: Thu Apr 20 00:32:39 2017
Modify: Thu Apr 20 00:32:39 2017
Change: Sun May  7 01:03:05 2017

atimeとmtimeがgitでのファイル更新時刻になりました。やったー!

git logはタイムゾーン込みで保存しているようなので、touchする前にローカルタイムゾーンに変換しないといけない。最初はsedタイムゾーン取り除いていたけれど、流石に乱暴だったのでdateコマンドを挟むようにしました。

できたーと思ってふと検索してみるといっぱいでてきた… まぁみんなやるよねぇ…

macOSの (というかFreeBSDの) dateコマンドの-dは日付を指定するオプションじゃないのが罠っぽい。date -jfであれこれしたほうが本当はいいけれど、coreutilsのgdateに逃げた。あとスペース含むファイル名の扱いにはご注意を。

dateコマンドのところはPythonでやろうかと思ったけれど、実装してみたら意外とオーバーヘッドが大きかった。まぁオーバーヘッドと言っても二倍時間がかかるくらいなんですが、リポジトリ内のファイルを全部なめて更新していくので、できれば速いほうが嬉しい。速度を取るかポータビリティを取るかはいつも難しい。PythonPythonで2.7と3系両方で動くスクリプトを書ける自信はない。ところで jq fromdateiso8601 を使う手も考えたのですが (ポータビリティ的には良いアイディア)、タイムゾーン部分をパースできなくて (Zしか対応してない)、なーにがISO 8601じゃと思いました。ああ、いい話だなぁ。それじゃーね!ばいばい!

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;
}