cachegrindを使って高速化してみる

[細かすぎて伝わらないエラトステネスの篩の高速化](http://qiita.com/Ishotihadus/items/73e107271275611f05f2)に触発されてプログラムを書いていた時にキャッシュを意識しないと高速化できなかったので、cachegrindを使ってその影響を調べてみました。

## cachegrindとは
バイナリ書き換えを使用したデバッグツール群Valgrindの中の、キャッシュプロファイラプログラム。```valgrind --tool=cachegrind```として実行する。
ぷろんご。

## 環境
* Core i7-4771@3.5GHz L3cache:8MB
* bash on Windows
* clang 3.8.0-2ubuntu4  -O3 -march=native
* valgrind-3.11.0

## キャッシュを汚さないようになるべく小さなデータ構造を使うと高速化する

### まずは普通に実行時間を計測してみる
```base.c
#include <stdio.h>
#define MAXVAL 200000000

int main() {
    int i, j, count = 0;
    static int flags[MAXVAL + 1] = {};
    for (i = 2; i <= MAXVAL; i++) {
        if (flags[i] == 0) {
            ++count;
            for (j = i * 2; j <= MAXVAL; j += i) {
                flags[j] = 1;
            }
        }
    }
    printf("%d\n", count);
    return 0;
}
```

flagsの型を```long long```, ```int```, ```short```, ```char```のそれぞれに変更して、速度の変化を確認します。計測はtimeコマンドを使い、最低5回計測した平均をとることにします。

|型       |bit幅|実行時間|
|:--------|----:|:-----:|
|long long|64bit|   5.8s|
|int      |32bit|   4.7s|
|short    |16bit|   4.0s|
|char     | 8bit|   3.5s|

bit幅が短い方が実行時間が短くなっています。フラグに必要な領域は1bitで済むので、以下のようにしてみます。

```base_bit.c
#include <stdio.h>
#define MAXVAL 200000000

int main() {
    int i, j, count = 0;
    static long long flags[(MAXVAL + 1) / 64] = {};
    for (i = 2; i <= MAXVAL; i++) {
        if ((flags[i / 64u] & (1ull << i % 64u)) == 0) {
            ++count;
            for (j = i * 2; j <= MAXVAL; j += i) {
                flags[j / 64u] |= 1ull << j % 64u;
            }
        }
    }
    printf("%d\n", count);
    return 0;
}
```

実行結果を以下にまとめてみます。

|型       |bit幅|実行時間|
|:--------|----:|:-----:|
|long long|64bit|   5.8s|
|int      |32bit|   4.7s|
|short    |16bit|   4.0s|
|char     | 8bit|   3.5s|
|*bit*    | 1bit|   2.3s|

### 余談

* -march=nativeによる違いは、レジスタ割り当ての違いを除けば、```inc esi```が```add esi, 1```に変わるだけでした。inc命令はキャリーフラグを変更しない、つまり前の結果に依存する(パーシャルレジスタストールが起こる)ため、遅くなります。intelの最適化マニュアルにもinc/decはadd/subにしろと書いてあります。

* flagsは、main関数中のstatic変数としても、global領域の変数としても、あるいは```= {}```をつけてもつけなくても、速度は変わりませんでした。

* base_bit.cにおいて```i / 64u```や```i % 64u```としている部分があります。C99以降では、割り算の商は0側へ丸められます。このため、64で割るといった場合でも、単純に算術シフト命令を使うわけにはいかず、以下のような命令が生成されてしまいます。

```
mov edx, r9d  ;; r9d は i が入っている。
sar edx, 31   ;; edx = i<0 ? -1 : 0;　と等価。
shr edx, 26   ;; ここまでで、edx = i<0 ? 63 : 0;  と等価。  
add edx, r9d  ;; ここまでで、edx = i<0 ? i+63 : i; と等価。
sar edx, 6    ;; edx >>= 6; により、edx = i / 64; を計算した。
```

また、余りに関しても全く同様の問題により、単純に下6bitのマスクをとるというわけにいかなくなっています。実際には```i```に負の値が入ることはないので、単純に算術/論理シフトや論理積演算をしてもらえれば同じ計算結果になるはずですが、コンパイラには見抜けなかったようです。解決方法としては、

1. シフト演算子```>>```とビット積演算子```&```を使う
2. ```i```を```unsigned```で宣言する
3. ```i```を```unsigned```にキャストする
4. リテラル```64```(これは```int```型である)を```64u```として```unsigned```型のリテラルとする(```unsigned```を含む二項演算はもう片方も```unsigned```に変換されます)

などが考えられます。これらはどれも一長一短です[^1]が、ソースコードの変更がなるべく少なくなるようにここでは4.を選択しました。
実際どのくらい遅くなるかというと、全部```64```と書いた場合が2.68s、全部```64u```と書いた場合が2.36sと、10%以上変わってくることが分かります。

[^1]: シフト演算子は、「見慣れない」とか「最適化はコンパイラに任せろ」との意見があり、また、ループ変数を```unsigned```で宣言することは賛否両論です(たとえば、[Google C++スタイルガイド](https://ttsuki.github.io/styleguide/cppguide.ja.html#Integer_Types))。整数昇格は多くの落とし穴を生んできました。

### cachegrindでキャッシュヒット率を見る

|型       |bit幅|実行時間|LLd misses(read) |LLd misses(write)|LLd misses(sum)|
|:--------|----:|:-----:|-----------------|----------------|----------------|
|long long|64bit|   5.8s|※               |※              |※              |
|int      |32bit|   4.7s| 6.3%(12.5M/200M)|68.7%(430M/627M)|53.6%(443M/827M)|
|short    |16bit|   4.0s| 3.1%( 6.3M/200M)|60.1%(376M/627M)|46.3%(383M/827M)|
|char     | 8bit|   3.5s| 1.6%( 3.2M/200M)|52.3%(328M/627M)|40.0%(331M/827M)|
|*bit*    | 1bit|   2.3s|17.2%(176M/1027M)| 0.0%( 818/200M)|14.4%(176M/1227M)|
※配列のサイズが大きすぎてmmapに失敗

やはりbit数の小さい型を使うにつれてキャッシュミスが減っていることが分かります。特にreadは```int```→```short```→```char```ときれいにミス率が半減しています。一方で*```bit```*は傾向が大きく違いますが、複合代入命令を使っているところでは、参照回数・キャッシュミス回数はread側に計上され(書き込む前に読み込みが発生するため)、write側には計上されないこと、bt命令をなぜかwrite回数に計上していることが原因であると考えられます[^2]。このようにread/writeや参照回数を見ていてもあまり意味のある結果が得られないため、以下では、ミス回数を中心に見ていきたいと思います。

[^2]: base_bit.cの```long long```を```char```や```short```に置き換えた場合、機械語としては複合代入命令ではなく二つの命令に分割されました。この時、参照回数は200Mに627Mを加えた827Mとなったため、この予測はおそらくあっていると思われます。逆に、```int```に置き換えた場合は複合代入命令にコンパイルされ、参照回数は200Mのままでした。次に、コンパイル時にアセンブリ言語で出力し、bt命令の部分にもう一度bt命令を挿入したところ参照回数が400Mになったので、write参照回数200Mの原因はbt命令です(cachegrindのバグ？)。read回数は外周のforループの200M＋内周の複合代入命令627Mで合わせて827Mのはずですが、1027Mになる理由はわかりませんでした。

## if文のコストとキャッシュミスのコスト、どっちが高い？

> 知っている合成数 × 今見ている素数 のパターンは、既にフラグが立っているはずなので、これはフラグを起こさなくてよい。

を用いた最適化をする場合、if文のコストがかかります。逆にそれをしない場合、flagがすでに立っているところのflagを立てようとしてキャッシュミスするという残念なことになる可能性があります。このトレードオフについて、調べてみます。

```if.c
#include <stdio.h>
#include <math.h>
#define MAXVAL 200000000

int main() {
    int i, j, f, s, count = 2;
    int max = (int)sqrt(MAXVAL) + 1;
    static int flags[MAXVAL + 1] = {};
    for (i = 5, f = 4; i <= max; i += (f = 6 - f)) {
        if (!flags[i]) {
            ++count;
            s = MAXVAL / i;
            for (j = s - !(s & 1); j >= i; j -= 2) {
                if (!flags[j]) flags[j * i] = 1;
            }
        }
    }
    for (; i <= MAXVAL; i += (f = 6 - f))
        if (!flags[i])
            ++count;
    printf("%d\n", count);
    return 0;
}
```

ifなしのコードは、14行目(最内for文の中)のif文をなくすだけです。

|型       |bit幅|ifあり  |          |ifなし |          |
|:--------|----:|:-----:|---------:|:-----:|---------:|
|         |     |実行時間|LLd misses|実行時間|LLd misses|
|long long|64bit|  1.15s|※        |  1.83s|※        |
|int      |32bit|  0.88s|    102.4M|  1.58s|    166.0M|
|short    |16bit|  0.68s|     76.4M|  1.31s|    142.7M|
|char     | 8bit|  0.56s|     58.0M|  1.12s|    119.8M|
|*bit*    | 1bit|  0.57s|     23.4M|  0.64s|     58.3M|


ifなしの方がキャッシュミスが発生しやすく、それにより無駄な時間がとられていることが分かります。flagを1bitにした場合は、この影響が小さくif文のオーバーヘッドの方が大きいと考えていましたが、この計測によればまだifありの方が速いようです。

```if_no2x.c
#include <stdio.h>
#include <math.h>
#define MAXVAL 200000000
#define BITS (sizeof(int) * 8)
#define FLAGS_NUM (MAXVAL / BITS + 1)

static int flags[FLAGS_NUM] = {};

inline static void upflag(int i) {
    i = i / 2 - 1;
    flags[i / BITS] |= 1 << (i % BITS);
}

inline static int getflag(int i) {
    i = i / 2 - 1;
    return (flags[i / BITS] >> (i % BITS)) & 1;
}

int main() {
    int i, j, f, s, count = 2;
    int max = (int)sqrt(MAXVAL) + 1;
    for (i = 5, f = 4; i <= max; i += (f = 6 - f)) {
        if (!getflag(i)) {
            ++count;
            s = MAXVAL / i;
            for (j = s - !(s & 1); j >= i; j -= 2) {
                if (!getflag(j))
                    upflag(j * i);
            }
        }
    }
    for (; i <= MAXVAL; i += (f = 6 - f))
        if (!getflag(i))
            ++count;
    printf("%d\n", count);
    return 0;
}
```

ごくわずかに修正して、2の倍数用のbitを削りました。ここまで来るとさすがにキャッシュヒット率が上がるのでif文をなくした方が速度が上がるようです。

|型       |bit幅|ifあり  |          |ifなし |          |
|:--------|----:|:-----:|---------:|:-----:|---------:|
|         |     |実行時間|LLd misses|実行時間|LLd misses|
|long long|64bit|  1.15s|※        |  1.83s|※        |
|int      |32bit|  0.88s|    102.4M|  1.58s|    166.0M|
|short    |16bit|  0.68s|     76.4M|  1.31s|    142.7M|
|char     | 8bit|  0.56s|     58.0M|  1.12s|    119.8M|
|*bit*    | 1bit|  0.57s|     23.4M|  0.64s|     58.3M|
|*halfbit*|0.5bit| 0.60s|     12.2M|  0.55s|     34.1M|


### 余談

* 手元のノートパソコンでは1bitの場合はif文を削除した方が速く、Ishotihadus氏のMacでも同様の結果が得られていたのですが、今回の環境ではif文があったほうが良いようです。
* 平方根はコンパイル時に計算されている(14144という数値が埋め込まれている(<=maxのところが<14144にコンパイルされる) )のでそもそも速いとかはありませんでした。そういえば```-lm```をつけてコンパイルしていないような……。

## キャッシュをヒットさせるためにブロック化する
[fujitanozomu氏がコメントに記載](http://qiita.com/Ishotihadus/items/73e107271275611f05f2#comment-c362a769dde750d8cd1d)した、ブロック化を行ったコードについてもcachegrindで確認してみましょう。
このコードは、$\sqrt N$までifなしで篩った後、それらの倍数の篩いをブロック化して行っているようです。フラグはcharで持っていて、bit演算を使っていないようです。

|ブロックサイズ|実行時間|LLd misses|D1 misses|備考|
|------------:|:-----:|---------:|--------:|----|
|    200000000|  3.85s|   434.10M|  437.85M|実質ブロック化なし|
|     10000000|  1.05s|   101.63M|  231.42M|
|      8388608|  0.89s|     6.32M|  222.01M|L3キャッシュは8MB(ただし共有)|
|      7340032|  0.80s|     6.34M|  226.57M|
|      6291456|  0.74s|     6.28M|  222.01M|
|      4194304|  0.58s|     6.28M|  221.96M|
|      2097152|  0.57s|     6.27M|  221.71M|
|      1048576|  0.56s|     6.26M|  219.77M|
|       262144|  0.46s|     6.25M|  211.41M|
|        65536|  0.40s|     6.25M|  159.83M|
|        40960|  0.39s|     6.25M|  107.38M|
|        34048|  0.26s|     6.25M|   39.27M|
|        32768|  0.23s|     6.25M|   13.86M|L1dキャッシュは32KB|
|        16384|  0.27s|     6.25M|    8.64M|
|         8192|  0.30s|     6.25M|    6.33M|

ブロックサイズをちょうどL1キャッシュのサイズ32KBに設定するとD1 missesがかなり減って高速になっています。それよりブロックサイズを減らすと分岐回数の増加などにより、徐々に遅くなるようです。

## まとめ
高速化をしたいときはアセンブリを見たりプロファイラを見たりするものですが、キャッシュの威力はとても強いのでキャッシュプロファイラも使いましょう。