コンパイル時計算完全に理解したリクルートテクノロジーズの竹迫です。
この記事はRecruit Engineers Advent Calendar 2019の4日目(12/4)のエントリーです。先月、社内のTGIFで飛び入り発表したスライドを記事にまとめなおしました。
adventar.org
皆さんは素数を数えることが好きでしょうか？私は大好きです。
問題：MAX以下の素数の個数を数えよ。
与えられた数以下の素数の個数を予測する「ガウスの素数階段」を作る簡単な問題です。素数かどうかを判別するにはどんな数字でも割り切れないことを確認する必要があるため、ナイーブに素数判別を計算するとどうしても時間がかかってしまいます。（※厳密に素数の個数を予測する魔法の「リーマンの素数公式」はとても面白いのですが、今回の記事では解説対象外です。）
引用元：http://tsujimotter.hatenablog.com/entry/riemann-prime-number-formula
（ガウスの素数階段を改良した）リーマンの素数公式に興味のある方は上記URLに日曜数学者のtsujimotterさんによる可視化アプリもあるので（ぜひ解説も含めて）お楽しみください。
C++でコンパイル時計算というテクニックを使うと、コンパイル時に計算結果が定数として出力され、プログラム実行時に計算コストが一切かからないというメリットがあります。
C++03のテンプレートメタプログラミングでMAX以下の素数の個数を計算してみます。
#include <iostream>
template<int N, int n=N-1, int m=N%n> struct isPrime {
enum { ok = isPrime<N, n-1>::ok };
};
template<int N, int n> struct isPrime<N, n, 0> {
enum { ok = 0 };
};
template<int N> struct isPrime<N, 1, 0> {
enum { ok = 1 };
};
template<int N> struct countPrime {
enum { sum = isPrime<N>::ok + countPrime<N-1>::sum };
};
template<> struct countPrime<1> {
enum { sum = 0 };
};
int main() {
std::cout << countPrime<MAX>::sum << std::endl;
}
複雑にならないよう以前herumiさんに教えてもらったenumを使うテクニックを使っています。
このプログラム prime.cpp をコンパイルして、整数MAX=100までの素数の個数を数えてみます。
$ g++ -DMAX=100 prime.cpp && ./a.out
25
100までの素数の個数25が出力されました。
objdumpでディスアセンブルしてみる
本当にコンパイル時に計算されているのかどうか、実行時に計算しているのではないかという疑惑を払拭するため、コンパイルされた結果のバイナリをobjdumpコマンドでディスアセンブルして中身を見てみます。
$ objdump -d a.out
＜抜粋＞
000000000000088a <main>:
88a: 55 push %rbp
88b: 48 89 e5 mov %rsp,%rbp
88e: be 19 00 00 00 mov $0x19,%esi
893: 48 8d 3d 86 07 20 00 lea 0x200786(%rip),%rdi # 201020
89a: e8 c1 fe ff ff callq 760 <_ZNSolsEi@plt>
89f: 48 89 c2 mov %rax,%rdx
8a2: 48 8b 05 27 07 20 00 mov 0x200727(%rip),%rax # 200fd0
8a9: 48 89 c6 mov %rax,%rsi
8ac: 48 89 d7 mov %rdx,%rdi
8af: e8 8c fe ff ff callq 740 <_ZNSolsEPFRSoS_E@plt>
8b4: b8 00 00 00 00 mov $0x0,%eax
8b9: 5d pop %rbp
8ba: c3 retq
ちなみにobjdumpのデフォルトはGAS形式のAT&T表記で出力されるため、私にとって読みやすいIntel表記で出力するにはobjdumpに-M intelオプションを指定します。
$ objdump -M intel -d a.out
＜抜粋＞
000000000000088a <main>:
88a: 55 push rbp
88b: 48 89 e5 mov rbp,rsp
88e: be 19 00 00 00 mov esi,0x19
893: 48 8d 3d 86 07 20 00 lea rdi,[rip+0x200786] # 201020
89a: e8 c1 fe ff ff call 760 <_ZNSolsEi@plt>
89f: 48 89 c2 mov rdx,rax
8a2: 48 8b 05 27 07 20 00 mov rax,QWORD PTR [rip+0x200727] # 200fd0
8a9: 48 89 c6 mov rsi,rax
8ac: 48 89 d7 mov rdi,rdx
8af: e8 8c fe ff ff call 740 <_ZNSolsEPFRSoS_E@plt>
8b4: b8 00 00 00 00 mov eax,0x0
8b9: 5d pop rbp
8ba: c3 ret
これを見ると「mov esi,0x19」のコードがあり、100までの素数の個数（25=0x19）がコンパイル時に定数となってバイナリに埋め込まれていることがわかります。コンパイル時計算すごいですね。
線形再帰の問題
ここで、調子に乗ってMAX=1000までの素数の個数を計算してみます。
$ g++ -DMAX=1000 prime.cpp
prime.cpp: In instantiation of ‘struct isPrime<997, 101, 88>’:
prime.cpp:5:8: recursively required from ‘struct isPrime<997, 995, 2>’
prime.cpp:15:31: recursively required from ‘struct countPrime<999>’
prime.cpp:15:31: required from ‘struct countPrime<1000>’
prime.cpp:60:31: required from here
prime.cpp:5:8: fatal error: template instantiation depth exceeds maximum of 900
(use -ftemplate-depth= to increase the maximum)
enum { ok = isPrime<N, n-1>::ok };
^
compilation terminated.
おっと。コンパイル時にfatal errorが出て止まりました！
よくエラーメッセージを見ると「template instantiation depth exceeds maximum of 900」と書いてあり、ご丁寧にも(use -ftemplate-depth= to increase the maximum)と深さを増やすオプション指定をコンパイラが教えてくれています。
$ g++ -ftemplate-depth=1000 -DMAX=1000 prime.cpp
$ ./a.out
168
このようにC++03の古典的なテンプレートメタプログラミングにおいては、線形再帰だと探索の深さが問題になるので、わざわざ二分再帰の形に書き直して深さではなく幅を広く探索するようにプログラムを書きなおすテクニックが必要でした。テンプレートメタプログラミングが魔術と言われる所以はこういうところにもありました。
C++14 constexpr を使おう
C++11の言語仕様からconstexprが追加され、コンパイル時定数式として評価できる関数を定義できるようになりました。ただ、C+11ではconstexpr関数の中で事実上return文一個しか書くことができず、変数を使った処理やif文やforループの繰り返しを書くことが困難で、変数は引数を使って表現し直して、繰り返しのループは再帰関数呼び出しに置き換えたりする必要があり、実際のところかなり面倒でした。
そこで、C++14の仕様ではconstexpr関数の制限が大きく緩和され、変数やif文の条件分岐や繰り返しのforループなどが使えるようになりました。これで普通の人でもカジュアルにconstexprを書くことができます。
わかりやすくなるよう、いったんコンパイル速度は気にせずナイーブに、C++14 constexpr関数を使ってコンパイル時計算で素数の個数を数えるプログラムを書いてみました。
#include <iostream>
template<size_t N>
struct constexpr_Prime {
int debug = 0x33333333;
int isPrime[N+1];
constexpr constexpr_Prime() : isPrime() {
int ok = 0;
isPrime[0] = 0;
isPrime[1] = 0;
for(int i = 2; i <= N; i++) {
for (int j = 2; j <= i/2; j++) {
if (i % j == 0) { --ok; break; }
}
isPrime[i] = ++ok;
}
}
};
int main() {
constexpr auto cp = constexpr_Prime<MAX>();
std::cout << cp.isPrime[MAX] << std::endl;
}
ここではMAX以下の素数の個数をすべてコンパイル時にisPrimeの定数配列に計算して格納しています。実行時に使うときは、constexpr auto cp = constexpr_Prime<MAX>();で実体化してから
for (int i = 0; i <= MAX; i++) {
std::cout << i << ":" << cp.isPrime[i] << std::endl;
}
のようにして呼び出します。実行時は定数配列を添え字で参照するだけなので超高速です。
$ g++ -DMAX=100 prime.cpp && ./a.out
25
実際にobjdumpコマンドで定数配列が生成されているかどうかを確認してみます。
$ g++ -DMAX=100 prime.cpp && ./a.out
25
$ objdump -z -S -j .rodata a.out
a.out: file format elf64-x86-64
Disassembly of section .rodata:
0000000000000a60 <_IO_stdin_used>:
a60: 01 00 02 00 00 00 00 00 00 00 00 00 00 00 00 00 ................
a70: 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 ................
0000000000000a80 <_ZStL19piecewise_construct>:
a80: 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 ................
a90: 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 ................
aa0: 33 33 33 33 00 00 00 00 00 00 00 00 01 00 00 00 3333............
ab0: 02 00 00 00 02 00 00 00 03 00 00 00 03 00 00 00 ................
ac0: 04 00 00 00 04 00 00 00 04 00 00 00 04 00 00 00 ................
ad0: 05 00 00 00 05 00 00 00 06 00 00 00 06 00 00 00 ................
ae0: 06 00 00 00 06 00 00 00 07 00 00 00 07 00 00 00 ................
af0: 08 00 00 00 08 00 00 00 08 00 00 00 08 00 00 00 ................
b00: 09 00 00 00 09 00 00 00 09 00 00 00 09 00 00 00 ................
b10: 09 00 00 00 09 00 00 00 0a 00 00 00 0a 00 00 00 ................
b20: 0b 00 00 00 0b 00 00 00 0b 00 00 00 0b 00 00 00 ................
b30: 0b 00 00 00 0b 00 00 00 0c 00 00 00 0c 00 00 00 ................
b40: 0c 00 00 00 0c 00 00 00 0d 00 00 00 0d 00 00 00 ................
b50: 0e 00 00 00 0e 00 00 00 0e 00 00 00 0e 00 00 00 ................
b60: 0f 00 00 00 0f 00 00 00 0f 00 00 00 0f 00 00 00 ................
b70: 0f 00 00 00 0f 00 00 00 10 00 00 00 10 00 00 00 ................
b80: 10 00 00 00 10 00 00 00 10 00 00 00 10 00 00 00 ................
b90: 11 00 00 00 11 00 00 00 12 00 00 00 12 00 00 00 ................
ba0: 12 00 00 00 12 00 00 00 12 00 00 00 12 00 00 00 ................
bb0: 13 00 00 00 13 00 00 00 13 00 00 00 13 00 00 00 ................
bc0: 14 00 00 00 14 00 00 00 15 00 00 00 15 00 00 00 ................
bd0: 15 00 00 00 15 00 00 00 15 00 00 00 15 00 00 00 ................
be0: 16 00 00 00 16 00 00 00 16 00 00 00 16 00 00 00 ................
bf0: 17 00 00 00 17 00 00 00 17 00 00 00 17 00 00 00 ................
c00: 17 00 00 00 17 00 00 00 18 00 00 00 18 00 00 00 ................
c10: 18 00 00 00 18 00 00 00 18 00 00 00 18 00 00 00 ................
c20: 18 00 00 00 18 00 00 00 19 00 00 00 19 00 00 00 ................
c30: 19 00 00 00 19 00 00 00 ........
_ZStL19piecewise_constructのシンボル上に100個までの素数の個数が定数展開されていることがわかります。
ここまで自力で出来たので「C++完全に理解した」Tシャツを着ることができる人権が得られました。
引用元：https://twitter.com/shapoco/status/958302684880109568
最近のC++20では、constexprにtry-catchブロックが書けるようになったり、インラインアセンブリ(!)が書けるようになったりしています。C++の進化はすごいですね。
C++の言語仕様とコンパイラの対応状況などは以下のページでまとめられています。
cpprefjp.github.io
ちなみに、ワタシハC++チョットデキルTシャツはこの世界に存在してはならないと個人的には思います。
実行時間は？
さきほどのコンパイル時に素数を計算する時間をシェルスクリプトを書いて計測してみました。
n="10 100 1000 5000 10000 20000 30000 40000 50000"
for i in $n
do
echo -n "$i, "
/usr/bin/time -f "%U" g++ -DMAX=$i prime.cpp
done
実行結果
10, 0.20
100, 0.21
1000, 0.32
5000, 1.65
10000, 5.29
20000, 20.65
30000, 45.87
40000, 81.28
50000, 129.12
グラフにすると以下のような形に
素数の数え上げにアトキンのふるいのアルゴリズムを実装するなどすれば、もう少しまともに早くなりますが。
コンパイル時計算は、プログラム実行時を超高速に処理するために、事前コンパイルをものすごく頑張るというポリシーです。
さいごに：素数表を持ち歩こう
私は歩いているときに視界に入ってきた数字が素数かどうかをパッと見で判別するのが趣味で、暗黒通信団の素数表をいつもカバンの中に入れて持ち歩いています。
引用元：https://www.amazon.co.jp/dp/4873101565/
本体価格は税別で357円とお買い求めしやすい素数となっております。（いや、3で割り切れるので違いますね）
計算機資源を使わない最高のキャッシュとして、紙を使った素数表、一家に一冊いかがでしょうか？
おわり
