定数倍高速化の技術

https://www.youtube.com/live/j1KGQHlHJ5E

競技プログラミングにおいて，定数倍高速化は TLE を回避するための重要な技術である．定数倍高速化により，想定より遅いオーダーの解法を通したり，作問時にそれを防ぐことができたりする．もちろん，競技プログラミング以外の一般のプログラミングにおいても，定数倍高速化はプログラムの質に関わる重要な技術である．
本講義では，

- 計算の仕組み
- 命令の速さ
- コンパイラによる最適化
- アセンブラ
- メモリアクセスの仕組み
- ベクトル化

などのトピックを定数倍高速化の観点から扱い，C++ での様々な定数倍高速化のテクニックを教える．プログラムの定数倍を予想・改善できるようになることが，本講義の目標である．

#### 前提知識 (なくても得るものはあると思います)

- C++ のプログラミング経験
- 競技プログラミングの経験

Transcript

1. 定数倍高速化の技術 2024/07/18 @ traP アルゴリズム班集会 @tatyam_prime 1
2. JOI2018/2019 春合宿 Day1 の出来事 2
3. 試験 (Examination) 3
4. 試験 (Examination) 2 次元平面に 個の点がある． N 4
5. 試験 (Examination) 整数 が 与えられるので，この範囲に ある点の個数を求める を 回 A, B, C Q 5
6. 試験 (Examination) 2 次元平面に 個の点がある． 番目の点は にある． クエリに答えよ． 整数 が与えられる． かつ かつ を満たす 点の個数は？ N i (X ​ , Y ​ ) i i Q A, B, C X ​ ≥ i A Y ​ ≥ i B X ​ + i Y ​ ≥ i C N, Q ≤ 105 6
7. 想定解 7
8. 想定解 定理 辺の多角形領域 以外の部分 は， 個の平行四辺形領域に 分割できる K K 8
9. 想定解 辺の多角形領域 3 9
10. 想定解 辺の多角形領域以外の部分 3 10
11. 想定解 辺の多角形領域以外の 部分は， 個の 平行四辺形領域に分割できる 3 3 11
12. 想定解 辺の多角形領域 3 12
13. 想定解 辺の多角形領域以外の部分 3 13
14. 想定解 辺の多角形領域以外の 部分は， 個の 平行四辺形領域に分割できる 3 3 14
15. この平行四辺形領域にある 点の数は？ を 回 Q 15
16. 座標変換 に より長方形領域に変換 y ← x + y 16
17. この長方形領域にある点の数は？ を 回 Q 17
18. この長方形領域にある点の数は？ を 回 → 平面走査 で 時間 Q O((N + Q) log(N + Q)) 18
19. 平面走査とは… 次元クエリ (変更なし) がたくさんあるとき， ある次元でクエリをソートしておくことで， 次元クエリ (変更あり) で処理するテクニック d d − 1 19
20. 座標の昇順に点を 追加していき，クエリを処理する 点 の処理 に到達したら， に 点を 個追加する クエリ の処理 に到達したら の 点の個数を数える y (X, Y ) y = Y x = X 1 (A, C) y = C x ≥ A 20
21. 私の解法 2 次元クエリをそのまま処理する 領域木で 時間 なら間に合う！ O(N log N + Q(log N) ) 2 N, Q ≤ 105 21
22. 私の解法 2 次元クエリをそのまま処理する 領域木で 時間 なら間に合う！ 領域木 : 座標の区間でセグメント木を作り， セグメント木の各ノードが， 「 座標がその範囲に含まれる点の 座標のリスト」をソートしたものを持つ O(N log N + Q(log N) ) 2 N, Q ≤ 105 x x y 22
23. QCFium の解法 1. 時間の 愚直解を書く while(Q--) { int ans = 0; for(int i = 0; i < N; i++) { if(S[i] >= A && T[i] >= B && S[i] + T[i] >= C) { ans++; } } cout << ans << '\n'; } Θ(NQ) 23
24. QCFium の解法 1. 時間の 愚直解を書く 2. コードの最初に 謎の呪文を追加 #pragma GCC target("avx2") #pragma GCC optimize("O3") #pragma GCC optimize("unroll-loops") // (中略) while(Q--) { int ans = 0; for(int i = 0; i < N; i++) { if(S[i] >= A && T[i] >= B && S[i] + T[i] >= C) { ans++; } } cout << ans << '\n'; } } Θ(NQ) 24
25. QCFium の解法 1. 時間の 愚直解を書く 2. コードの最初に 謎の呪文を追加 3. AC 参考資料 : Speeding Up for Naive Algorithm - E869120 #pragma GCC target("avx2") #pragma GCC optimize("O3") #pragma GCC optimize("unroll-loops") // (中略) while(Q--) { int ans = 0; for(int i = 0; i < N; i++) { if(S[i] >= A && T[i] >= B && S[i] + T[i] >= C) { ans++; } } cout << ans << '\n'; } } Θ(NQ) 25
26. 定数倍高速化ってすごい 定数倍高速化の技術を習得すると… が想定の問題を で通す が想定の問題を で通す こういうことをされない問題を作る 定数倍で落ちたときの原因を見つける ことができるようになる O(N log N) Θ(N(log N) ) 2 O(N(log N) ) 2 Θ(N ) 2 26
27. 定数倍高速化の技術を習得しよう プログラムの計算量だけではなく，定数倍を 予測できるようになろう 定数倍高速化をできるようになろう 27
28. 定数倍高速化の極意 28
29. 定数倍高速化の極意 ボトルネックを知る 29
30. 計算量でボトルネックを知る プログラムに 時間の部分 時間の部分 時間の部分 があるとき，たいてい改善すべきは 時間の部分 (定数倍で差がついていと逆転することも) Θ(N) Θ(N log N) Θ(N(log N) ) 2 Θ(N(log N) ) 2 30
31. [高度なテク] 定数倍でボトルネックを知る 四則演算，ビット演算など色々な演算の速度 コンパイラによる最適化 メモリアクセスの速度 体感のデータ構造の速度 メモリ確保の速度 入出力の速度 などからプログラムの定数倍を知る 31
32. 計測でボトルネックを知る C++ の場合 std::chrono::system_clock::now() で時刻を取得 時刻の差分を std::chrono::duration_cast<精度>() で間隔に 変換 間隔を duration.count() で整数に変換 32
33. #include <chrono> using namespace std::chrono; int main() { auto t0 = system_clock::now(); // なにか処理 auto t1 = system_clock::now(); // 間隔に変換 auto d01 = duration_cast<milliseconds>(t1 - t0); // 整数に変換して出力 cout << d01.count() << " ms" << endl; } 33
34. 色々な演算の速度を知ろう の前に… CPU のしくみを知ろう 34
35. 参考資料 (CPU のつくりかた) ディスクリート半導体の基礎 第1章 半導体の基礎 – TOSHIBA ディスクリート半導体の基礎 第3章 トランジスタ – TOSHIBA コンピュータ講座 応用編 第1回 - 富士通 ASCII.jp：半導体プロセスまるわかり　トランジスタの配線と 形成 35
36. 注 : 雑な理解で雑に説明します 半導体 Q. 半導体の主成分となる元素といえば？ 36
37. A. ケイ素 37
38. 4 価のケイ素は，ダイヤモンドと同じ形の結合で結晶をつくる 38
39. ここに 5 価のリンを少量混ぜると，電子が余る 余った電子が動き回る → n 型半導体 39
40. ここに 3 価のホウ素を少量混ぜると，電子が足りなくなる 正孔 (= 電子の足りない部分, ホール) が動き回る → p 型半導体 40
41. ダイオード n 型半導体と p 型半導体をくっつけると… 41
42. ダイオード 電子と正孔が打ち消しあって空乏層（絶縁体と同じ状態）ができる 42
43. ダイオード こっち向きに電圧をかけると，電子・正孔が増え， 結合部で打ち消しあって電流が流れる 43
44. ダイオード 逆向きに電圧をかけると，電子・正孔が減り，電流が流れない 44
45. トランジスタ (MOSFET) p p Source Drain Gate n 電流の⼀⽅通⾏ p 型半導体と n 型半導体をうまく組み合わせると， 電気を通す・通さないを制御できるようになる！ 45
46. トランジスタ (MOSFET) n p p Source Drain Gate n 電流の⼀⽅通⾏ ー ー ー ー S - G 間に電圧をかけてゲートに負の電荷をためると， 近くの電子が追い出されて正孔ができ，電気を通すようになる 46
47. 論理回路 トランジスタを組み合わせると，論理回路ができる 例えば，NOT 回路は以下のように作れる p p n n n p IN OUT GND 0V 電源 1V 47
48. n p p n n n p 0V 1V GND 0V ー ー ー ー 電源 1V p p n n n n p 1V 0V GND 0V + + + + 電源 1V 48
49. 豆知識 ↓ NOT はこのように表記され，無から有を生み出しているように 　見えるが，ちゃんと電源にも GND にもつながっている 49
50. 論理回路を組み合わせると 機能を持った回路を作れる！ 半加算器 から 1 bit を受け取り， 繰り上がりを に， 繰り上がりを除いた和を に 出力する A, B C S 50
51. 64 bit 加算器 たくさん繋げると 64 bit 加算器に なる 本当に直列につなぐと遅延が大きすぎるので， 分割統治して繰り上がりがある場合とない場合の 両方を計算し，最後に選ぶテクニックが存在する 51
52. D-フリップフロップ回路 電源が与えられている間… C が ON になるたび，D からの入力を保存し， 次に C が ON になるまで Q に出力し続ける 52
53. レジスタ これを 64 個繋げれば，64 bit 整数を 保持できる C が ON になるたび D から 与えられた 64 bit 整数を保存し， 次に C が ON になるまで Q に出力 これを レジスタ と呼ぶ． 53
54. スイッチ ＋ レジスタ 加算器 1 を出⼒ 太⽮印は 64 bit 計算のしくみ 64 bit のレジスタと 加算器でこんな回路を 作ってみる 54
55. スイッチ ＋ レジスタ 加算器 1 を出⼒ 太⽮印は 64 bit 計算のしくみ スイッチを ON にすると, レジスタは後ろからの 入力を保存し… 55
56. スイッチ ＋ レジスタ 加算器 1 を出⼒ 太⽮印は 64 bit 計算のしくみ スイッチを ON にすると, レジスタは後ろからの 入力を保存し，前に出力 56
57. スイッチ ＋ レジスタ 加算器 1 を出⼒ 太⽮印は 64 bit 計算のしくみ 加算器が足し算を計算 57
58. スイッチ ＋ レジスタ 加算器 1 を出⼒ 太⽮印は 64 bit 計算のしくみ 最後にスイッチを OFF に 戻して，1 クロックが 終了 58
59. スイッチ ＋ レジスタ 加算器 1 を出⼒ 太⽮印は 64 bit 計算のしくみ スイッチを ON にすると, レジスタは後ろからの 入力を保存し，前に出力 59
60. スイッチ ＋ レジスタ 加算器 1 を出⼒ 太⽮印は 64 bit 計算のしくみ 加算器が足し算を計算 60
61. スイッチ ＋ レジスタ 加算器 1 を出⼒ 太⽮印は 64 bit 計算のしくみ 最後にスイッチを OFF に 戻して，1 クロックが 終了 61
62. スイッチ ＋ レジスタ 加算器 1 を出⼒ 太⽮印は 64 bit 計算のしくみ スイッチが ON になった 回数がレジスタに 記録されている！ 62
63. ＋ － × ÷ MUX 命令 命令に対応した 計算結果 万能演算装置 足し算以外も 計算できるようにしよう！ → 引き算や掛け算， 割り算などの演算回路を 用意して，命令によって 使い分ける 63
64. スイッチ レジスタ 1 を出⼒ 太⽮印は 64 bit 万能演算装置 + × − ÷ 命令 ⾜し算 万能演算装置 色々な演算に対応した 万能演算装置を入れ， 命令を与えて計算する ように 64
65. スイッチ レジスタ 万能演算装置 太⽮印は 64 bit + × − ÷ 命令 ⾜し算 MUX レジスタ 番号 計算のしくみ レジスタをたくさん増やした スイッチを ON にするたび， 「レジスタ A の値と レジスタ B の値を掛けて， レジスタ C に保存」 のような命令を処理する 65
66. クロック 発振器 レジスタ 万能演算装置 太⽮印は 64 bit + × − ÷ 命令 ⾜し算 MUX レジスタ 番号 計算のしくみ スイッチをクロック発振器 （一定の周期で 0 / 1 を繰り返す, 音叉のようなもの） に変更して，計算が自動で 進むように 66
67. クロック 発振器 レジスタ 万能演算装置 太⽮印は 64 bit + × − ÷ 命令 ⾜し算 MUX レジスタ 番号 計算のしくみ 1. クロック発振器が ON になり， レジスタが後ろの入力を 保存して前に出力 2. 演算装置が計算を実行し，計算 結果がレジスタの後ろまで伝播 3. クロック発振器が OFF になる を 1 クロックとして，これを高速に 繰り返す → 計算機の完成！ 67
68. クロック 発振器 レジスタ 万能演算装置 太⽮印は 64 bit + × − ÷ 命令 ⾜し算 MUX レジスタ 番号 計算のしくみ 1. クロック発振器が ON になり， レジスタが後ろの入力を保存して前に 出力 2. 演算装置が計算を実行し，計算結果が レジスタの後ろまで伝播 3. クロック発振器が OFF になる もしクロック発振器が再び ON になるまでに 計算が終わらないと…？ 68
69. クロック 発振器 レジスタ 万能演算装置 太⽮印は 64 bit + × − ÷ 命令 ⾜し算 MUX レジスタ 番号 1. クロック発振器が ON になり， レジスタが後ろの入力を保存して前に 出力 2. 演算装置が計算を実行し，計算結果が レジスタの後ろまで伝播 3. クロック発振器が OFF になる もしクロック発振器が再び ON になるまでに 計算が終わらないと…？ → 間違った計算結果がレジスタに 保存されてしまう！ 69
70. 間違った計算結果が保存しないために → 最も時間のかかる計算 (割り算とか) に 1 クロックの長さを 合わせる 70
71. 間違った計算結果が保存しないために → 最も時間のかかる計算 (割り算とか) に 1 クロックの長さを 合わせる → 1 命令を複数のクロックにまたがって実行できるようにする！ 71
72. 参考情報 x86_64 (AVX-512 対応) の (論理) レジスタの数 整数レジスタ : 64 bit × 16 個 ベクトルレジスタ : 512 bit × 32 個 小数の計算はこっち レジスタで足りない分はメモリに置く必要がある 72
73. 参考情報 最近の CPU の速度 クロック周波数 : 2 ~ 4 GHz 3 GHz のとき，1 クロックに光は 10 cm しか進めない！ この時間内にすべての計算が終わるように設計され，回路が 詰め込まれている 73
74. 色々な演算の速度を知ろう 74
75. すごいサイト https://www.uops.info/table.html 75
77. あるいはすごい PDF https://www.agner.org/optimize/instruction_tables.pdf 77
78. レイテンシ (Latency) 1 回の命令に何クロックかかるか？ スループット (Throughput) 同じ命令を並列にたくさん与えたら 1 つあたり 何クロックかかるか？ よくあるのスループットの逆数であることに注意 78
79. クイズ プログラムの定数倍を予想できるようになろう！ AtCoder での実行時間は？ using u32 = uint32_t; int main() { u32 N = 1e9; u32 s = 0; for(u32 i = 1; i <= N; i++) s += N / i; cout << s << endl; } 79
80. ヒント s += N / i を 回 加算の部分は並列にできないが, 除算の部分は並列にできる N と s の 2 変数しかないので メモリは使わず，レジスタだけで 完結 AtCoder の CPU は 3.5 GHz using u32 = uint32_t; int main() { u32 N = 1e9; u32 s = 0; for(u32 i = 1; i <= N; i++) { s += N / i; } cout << s << endl; } 109 80
81. 計算してみよう レジスタの 32 bit 除算の行を見る 並列にたくさんできるのでスループットを見る → 6 クロック これを 回で クロック AtCoder の CPU は 3.5 GHz → 1 秒に クロック実行 くらい 109 6 × 109 3.5 × 109 (6 × 10 clk)/(3.5 × 9 10 clk/s) = 9 1.714 s 81
82. 答え 82
83. 色々な演算の速度を知ろう CPU : Icelake，64 bit 演算，メモリアクセスのない場合 命令 意味 Lat TP MOV 代入 1 0.25 ADD 加算 1 0.25 AND bit ごとの AND 1 0.25 SHR 右シフト 1 0.5 source : https://www.uops.info/table.html 83
84. 命令 意味 Lat TP POPCNT 立っている bit 数 3 1 ADDSD 浮動小数点加算 4 0.5 MULSD 浮動小数点乗算 4 0.5 IMUL 乗算 (符号つき) 4 1 DIVSD 浮動小数点除算 13-14 4 IDIV 除算 (符号つき) 15 10 SQRTSD 15-16 4-6 除算は遅い！ ​ x 84
85. 各演算の速度をざっくりと 足し算，引き算，ビット演算はコスト 1 掛け算はコスト 4 割り算はコスト 15 85
86. 最近 64 bit 除算は速くなった 割り算はコスト 15 と言っても…？ Skylake 現在の Codeforces は 除算が遅い！ 命令 意味 Lat TP IDIV 除算 (符号つき) 42-95 24-90 DIVSD 浮動小数点除算 13-14 4 Icelake (2019 年 〜) AtCoder では 2023 年の 言語アップデートから速い 命令 意味 Lat TP IDIV 除算 (符号つき) 15 10 DIVSD 浮動小数点除算 13-14 4 86
87. 素数判定するコードが…？ u64 n = (u64)4e16 + 63; const u64 sq = sqrtl(n); for(u64 i = 2; i <= sq; i++) if(n % i == 0) return 0; cout << "prime!" << endl; 旧ジャッジ 新ジャッジ 87
88. 整数除算が遅い場合の対策 53 bit 除算で十分なら double で除算して 切り捨てる！ 64 bit 除算で十分なら long double で 除算して切り捨てる！ 命令 意味 Lat TP IDIVL 32 bit 除算 26 6 IDIVQ 64 bit 除算 42-95 24-90 DIVSD 64 bit 浮動小数点除算 13-14 4 FDIV 80 bit 浮動小数点除算 14-16 4-5 88
89. ところで， 除算が並列にできるってなんだ？ s += N / i の部分が ボトルネック 加算の部分は並列にできないが, 除算の部分は並列にできる N と s の 2 変数しかないので メモリは使わず，レジスタだけで 完結 AtCoder の CPU は 3.5 GHz using u32 = uint32_t; int main() { u32 N = 1e9; u32 s = 0; for(u32 i = 1; i <= N; i++) { s += N / i; } cout << s << endl; } 89
90. 除算が並列にできるってなんだ？ プログラムは順番に実行されるから，並列にはできないのでは 90
91. 現代の CPU はいろいろなことをやっている… パイプライン処理 命令の実行をいくつかのフェーズに分け， 「流れ作業」で 命令列を処理することで並列化する！ 91
92. だいたいこんな感じ (本当はもっとすごい) 92
93. 問題発生！ 書き込みが反映される前に後ろの命令で読み込みをしている！ 93
94. 問題発生！ 分岐があるときにどこから命令を取ってくれば良いかわからない！ 94
95. 実際のパイプライン処理 参考資料 : Sunny Cove - Microarchitectures - Intel - WikiChip Sunny cove block diagram.png by Chipwikia; CC BY-SA 4.0 95
96. 命令をメモリから 読み込み 命令を構文解析し, μOP (より単純な 命令) に分解 実装されていない 命令を実装済みの 命令列に置き換え などしている 96
97. 分岐予測 分岐があるとき，どこから命令を 取ってくれば良いかわからない → とりあえず片方実行してみて， 間違っていたらやり直そう！ やり直すと何クロックも遅くなってしまう → どっちに分岐するかを予想しよう！ 97
98. if(rand() % 2) { // ... } else { // ... } 1/2 の確率で分岐予測に失敗し， やり直しに何クロックもかかる → 遅い！ for(int i = 0; i < 100; i++) { // ... // ... // ... } for 文の条件式はたいてい true → 分岐予測が基本的に成功し， 　 1 クロック程度で実行できる 98
99. 余談 「分岐予測して実行を進める」がセキュリティーホールになる話， おもしろいです 本当にわかる Spectre と Meltdown – Hirotaka Kawata 99
100. Register Renaming 命令を μOP に分解した後，論理レジスタに物理レジスタを割り当てる 実際のレジスタの数はプログラムが使えるレジスタの数よりかなり多い 論理 物理 整数レジスタ 16 280 ベクトルレジスタ 32 224 ※ Sunny Cove の値 (現在の AtCoder はこれ)．物理レジスタの数は CPU により異なります． 100
101. Register Renaming 命令を μOP に分解した後，論理レジスタに物理レジスタを 割り当てる 「同じレジスタに連続して書き込む」のような依存関係は， 書き込みごとに異なるレジスタを割り当てれば解消 101
102. スケジューラ できる計算から順にどんどん計算していってしまう 計算順序が変わることもある！ (Out-of-Order 実行) 102
103. スケジューラ できる計算から順にどんどん計算していってしまう (計算順序が変わることもある) ALU の回路が 4 つ → 1 クロックに加算が 4 回できる！ 103
104. コンパイラによる最適化を知ろう 104
105. クイズ AtCoder での実行時間は？ constexpr u64 MOD = 998244353; int main() { u64 s = 1; for(u64 i = 1; i < MOD; i++) { s *= i; s %= MOD; } cout << s << endl; } 105
106. ヒント 掛け算と割り算を 回 さっきは並列になったが， 今回は並列にならない 乗算は 4 クロック 除算は 15 クロック AtCoder の CPU は 3.5 GHz 出力は MOD - 1 (ウィルソンの定理) constexpr u64 MOD = 998244353; int main() { u64 s = 1; for(u64 i = 1; i < MOD; i++) { s *= i; s %= MOD; } cout << s << endl; } 109 106
107. 計算してみよう constexpr u64 MOD = 998244353; int main() { u64 s = 1; for(u64 i = 1; i < MOD; i++) { s *= i; s %= MOD; } cout << s << endl; } ​ = 3.5 × 10 clk/s 9 19 × 0.998 × 10 clk 9 5.419 s 107
108. 答え 108
109. なぜ予測より速くなった？ コンパイル結果を見てみよう！ 109
110. すごいサイト https://godbolt.org/ 110
111. すごいサイト https://godbolt.org/ 左にソースコードを入れると，右にコンパイル結果の命令列を 表示してくれる どの部分がどの命令にコンパイルされたかがわかる AtCoder とだいたい同じコンパイル結果がほしいときは， GCC 12.2 で -O2 -std=c++23 -march=icelake-server 111
112. コンパイル結果を見る s %= MOD; に対応する命令が 青くなっている movq : 64 bit 代入 mulq : 64 bit 乗算 shrq : 64 bit 右シフト subq : 64 bit 減算 除算命令が消えている！ 112
113. 除算の最適化 割る数がコンパイル時定数のとき，除算を掛け算 + 右シフトで 行うことがある ざっくりこんな気持ち 整数部 64 bit，小数部 64 bit の小数で掛け算を行う で割る代わりに， を 整数部 64 bit，小数部 64 bit の 小数で表しておき，これを掛けて整数部を取る 参考資料 : コンパイラによる整数除算最適化の証明 | nu50218 blog x 1/x 113
114. アセンブラを読もう！ アセンブラ : 右側にみえているこの言語 機械語を人間に読めるように 文字情報で表したもの 命令と引数が順に並んでいる 114
115. アセンブラを読もう！ 左のコードを -O2 で最適化すると，右のアセンブラになる using u64 = uint64_t; u64 div(u64 x) { return x / 3; } div(unsigned long): movabsq $-6148914691236517205, %rax mulq %rdi movq %rdx, %rax shrq %rax ret 115
116. アセンブラを読もう！ 命令 意味 div(unsigned long): 関数の 1 つ目の引数 x は %rdi に入る movabsq $-6148914691236517205, %rax %rax に 0xAAAAAAAAAAAAAAAB を代入 mulq %rdi %rax * %rdi の上位 64 bit を %rdx に, 下位 64 bit を %rax に 代入 movq %rdx, %rax %rdx を %rax に代入 shrq %rax %rax を 1 bit 右シフト ret %rax を返り値として関数を終了 116
117. つまり… u64 div(u64 x) { return x / 3; } は u64 div(u64 x) { return (u128)x * 0xAAAAAAAAAAAAAAAB >> 65; } に最適化された 117
118. 大前提 コンパイル時に最適化オプション（ -O2 など）を付けないと 最適化は行われない プログラムを高速に動かしたいときは，最適化オプションを 付けよう！ g++ -O2 main.cpp 118
119. 豆知識 : 命令の suffix 例 末尾 意味 bit 長 addb b Byte 8 bit addw w Word 16 bit addl l Long Word 32 bit addq q Quad Word 64 bit 昔ワードサイズが 16 bit だったときの命名を継承している Double Word で D のことも…？ 119
120. この命令なに？というときは すごいサイト https://www.felixcloutier.com/x86/ 120
121. 豆知識 : 乗算命令 C++ で 64 bit 整数の乗算をすると，結果は mod で返ってくる． u64 a = ULLONG_MAX; assert(a * a == 1); 264 121
122. 豆知識 : 乗算命令 C++ で 64 bit 整数の乗算をすると，結果は mod で返ってくる． 一方，64 bit 乗算命令は… 命令 意味 mulq %rdi %rax * %rdi の上位 64 bit を %rdx に， 下位 64 bit を %rax に 代入 128 bit で計算されている！？ 264 122
123. ご存知でしたか？ – https://www.felixcloutier.com/x86/mul 64 bit 乗算をすると結果は 128 bit で計算されている！ 我々は大抵，そのうち上位 64 bit を捨ててしまう 123
124. ご存知でしたか？ – https://www.felixcloutier.com/x86/div 64 bit 除算も 128 bit ÷ 64 bit ができる！ (商が 64 bit に収まらないと RE ) 124
125. 最適化の例 2 冪の掛け算・割り算はビット演算に最適化される 計算 最適化された命令 (例) x / 32 shrq $5, %rax x % 32 andl $31, %eax x * 33 salq $5, %rax addq %rdi, %rax 符号なし整数にした方がこの最適化には有利かも 125
126. 最適化の例 関数呼び出しを削除して，関数の中身を埋め込む (Inline 化) u64 square(u64 x) { return x * x; } u64 cube(u64 x) { return square(x) * x; } square(unsigned long): imulq %rdi, %rdi movq %rdi, %rax ret cube(unsigned long): movq %rdi, %rax imulq %rdi, %rax imulq %rdi, %rax ret square を呼び出している（ call 命令が入る）はずだが，消えている 126
127. メモリアクセスの速度を知ろう！ 127
128. クイズ AtCoder での実行時間は？ const u64 siz = 1 << 27; u8 A[siz]; int main() { u64 i = 1; rep(t, 100'000'000) { A[i] += t; i *= 5; i %= siz; } } 128
129. ヒント A の大きさは バイト 回， バイトの範囲に ランダムアクセス アクセス位置は先読み可能 → 並列にできる i は 5 のべき乗 を 通るが，この周期は const u64 siz = 1 << 27; u8 A[siz]; int main() { u64 i = 1; rep(t, 100'000'000) { A[i] += t; i *= 5; i %= siz; } } 227 108 227 mod 227 225 129
130. 答え 130
131. ランダムアクセスは メモリが律速 131
132. メモリアクセスの速度 メモリアクセスのレイテンシはだいたいこれくらい 大きさ (目安) レイテンシ (目安) レジスタ 2 KB 1 クロック メモリ 8 GB ~ ≥100 クロック メモリ制限はたいてい 1 GB だけど… メモリアクセスは遅い！ 132
133. メモリアクセスは遅い！ → 読み書きしたデータを一時的に保存して， 再利用できるようにしよう！ 133
134. メモリアクセスの速度 メモリアクセスのレイテンシはだいたいこれくらい 大きさ (目安) レイテンシ (目安) レジスタ 2 KB 1 クロック L1 データキャッシュ 48 KB 5 クロック L2 キャッシュ 1.25 MB 13 クロック L3 キャッシュ 54 MB 42 クロック メモリ 8 GB ≥100 クロック 注 : 大きさ / レイテンシは CPU により異なります 参考資料 : Intel Xeon Platinum 8375C – CPUWorld 134
135. 実測 : ランダムアクセスの速度 さっきのコードの siz を変化させたときの実行時間 (縦軸は アクセスあたりの平均クロック数) 1 135
136. 48 KB まで L1 キャッシュが効いて，平均 2.3 クロックくらい 1.25 MB まで L2 キャッシュが効いて，平均 3.5 クロックくらい 10 MB くらいまで L3 キャッシュが効いて，平均 7 クロックくらい それ以上はどんどん遅くなっていく 136
137. 0x00000000 0x00000040 0x00000080 0x000000C0 0x00000100 0x00000140 0x00000180 0x000001C0 0x00000200 0x00000240 0x00000280 0x000002C0 0x00000300 0x00000340 キャッシュメモリのしくみ メモリ・キャッシュメモリは 64 バイト単位 で読み書きする この 64 バイト単位を ライン と呼ぶ 137
138. キャッシュメモリのしくみ 64 バイト単位で読み書きする → シーケンシャルアクセス（for 文で順番にアクセス）すると… vector<int> A(N); for(int i = 0; i < N; i++) { A[i] = i; } int は バイトなので，16 回中 15 回は必ずキャッシュにある！ 4 138
139. シーケンシャルアクセスは速い！ 139
140. キャッシュメモリのしくみ (L1 キャッシュの場合) L1 キャッシュの大きさ : 注 : 大きさ は CPU により異なります ラインからなるグループが 64 個ある 64 × 12 × 64 Bytes 12 140
141. キャッシュメモリのしくみ (L1 キャッシュの場合) 0x12307080 の値が欲しい！ → 64 で割って，ライン番号 0x48C1C2 → ライン番号 mod 64 は 2 なので，グループ 2 の中を探す 0x12307080 0x12307080 141
142. キャッシュメモリのしくみ (L1 キャッシュの場合) 0x12307080 から始まるラインを保存したい！ → 64 で割って，ライン番号 0x48C1C2 → ライン番号 mod 64 は 2 なので，グループ 2 の中にしまう 0x12307080 0x12307080 142
143. キャッシュメモリのしくみ (L1 キャッシュの場合) ラインを保存したいがすでにいっぱいだった！ → 最近使ったのが最も遅いラインを削除する （削除戦略にもいろいろある） 0x12327080 0x12337080 0x12347080 0x12357080 0x12367080 0x12377080 0x12387080 0x12397080 0x123A7080 0x123B7080 0x123C7080 0x12307080 0x12317080 143
144. キャッシュメモリのしくみ (L2 キャッシュの場合) L2 キャッシュの大きさ : 注 : 大きさ は CPU により異なります ラインからなるグループが 個ある 1024 × 20 × 64 Bytes 20 1024 144
145. クイズ A は B の何倍速い？ （3 重ループ以外の部分の実行時間も含む） mt19937 rnd; u64 A[1000][1000]; rep(i, 1000) rep(j, 1000) A[i][j] = rnd(); A. rep(k, 1000) rep(i, 1000) rep(j, 1000) { chmin(A[i][j], A[i][k] + A[k][j]); } B. rep(k, 1000) rep(i, 1000) rep(j, 1000) { chmin(A[j][i], A[j][k] + A[k][i]); } 145
146. ヒント のワーシャルフロイド法 の入力生成部分は無視できる A はシーケンシャルアクセス B は配列が転置されていて， A[j][i] , A[j][k] が シーケンシャルアクセスにならない A. rep(k, 1000) rep(i, 1000) rep(j, 1000) chmin(A[i][j], A[i][k] + A[k][j]); B. rep(k, 1000) rep(i, 1000) rep(j, 1000) chmin(A[j][i], A[j][k] + A[k][i]); n = 103 O(n ) 2 146
147. 実測 #include <cstdint> #include <random> using namespace std; using u64 = uint64_t; #define rep(i, a) for(u64 i = 0; i < a; i++) void chmin(u64& a, u64 b) { if(a > b) a = b; } mt19937 rnd; u64 A[1000][1000]; int main() { rep(i, 1000) rep(j, 1000) A[i][j] = rnd(); rep(k, 1000) rep(i, 1000) rep(j, 1000) { chmin(A[i][j], A[i][k] + A[k][j]); } } 147
148. 答え 148
149. 答え 配列を転置させるだけで， 倍の差がついた！ ​ = 783 ms 1678 ms 2.143 149
150. シーケンシャルアクセスになるよう 添字の順番に気をつけよう！ 150
151. 実は，簡単にもっと遅くすることができる 151
152. 実は，簡単にもっと遅くすることができる 配列の大きさを 2 冪に揃えるだけ u64 A[1000][1000]; → u64 A[1024][1024]; が左シフトに変わって高速化するのでは ×1000 152
153. u64 A[1000][1000]; のとき， rep(j, 1000) A[j][k] は ライン間隔のメモリアクセスを 回 → 各ラインは L2 キャッシュの異なるグループに入る u64 A[1024][1024]; のとき， rep(j, 1000) A[j][k] は ライン間隔のメモリアクセスを 回 → L2 キャッシュのグループ 個のうち 8 個しか使えない！ → ラインしか保存できない → L2 キャッシュが無効化！ 8000 Byte = 125 1000 8192 Byte = 128 1000 1024 160 153
154. 実測 u64 A[1024][1024]; にするだけでさらに 倍 の差がついた！ ​ = 1678 ms 4970 ms 2.962 154
155. シーケンシャルアクセスなら問題は起こらない！ シーケンシャルアクセスになるよう添字の順番に気を つけよう！ 155
156. どうしてもシーケンシャルアクセスにできないとき 行列の転置 u64 A[3000][3000]; rep(i, 3000) rep(j, i) { swap(A[i][j], A[j][i]); } A[i][j] か A[j][i] のどちらかはシーケンシャルアクセスに できない！ 156
157. A[0][0] A[0][7] A[7][0] A[7][7] これは シーケンシャルアクセス 157
158. A[0][0] A[0][7] A[7][0] A[7][7] これは シーケンシャルアクセス ではない 158
159. A[0][0] A[0][7] A[7][0] A[7][7] じゃあこれは…？ 159
160. A[0][0] A[0][7] A[7][0] A[7][7] 縦方向のアクセス → シーケンシャルアクセスでは ない…？ 短い期間内に横方向に 連続してアクセス → 同じラインだから キャッシュが効く！ 160
161. A[0][0] A[0][7] A[7][0] A[7][7] ブロック化 A[i][j] の の範囲と の範囲を いくつかのブロックに分け， ブロック内のアクセスを キャッシュに乗せるテクニック i j 161
162. A[0][0] A[0][7] A[7][0] A[7][7] ブロック化 キャッシュに乗るくらいの 大きさのブロックに 切り分け，ブロックごとに 転置するだけで転置が 高速化される！ 162
163. Cache-Oblivious どんな大きさのキャッシュメモリも効率的に使って高速化できる アルゴリズムが存在する… 163
164. 1 2 3 4 Cache-Oblivious 1. 全体を 4 個のブロックに 分割し，ブロックごとに 転置を行う 164
165. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Cache-Oblivious 1. 全体を 4 個のブロックに 分割し，ブロックごとに 転置を行う 2. 各ブロックを再帰的に 4 個のブロックに分割し， これを繰り返す 165
166. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Cache-Oblivious 1. 全体を 4 個のブロックに分割し， ブロックごとに転置を行う 2. 各ブロックを再帰的に 4 個のブロックに 分割し，これを繰り返す 3. ブロックの大きさが確実にキャッシュに 乗るところで再帰を止めてループで処理 → 「ブロックごとに処理」を様々な 大きさでできるので，どんな大きさの キャッシュも効率的に使える！ 166
167. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Cache-Oblivious ちょっと順序を工夫すると ヒルベルト曲線 ができる → ヒルベルト曲線の順で 処理すると Cache- Oblivious！ (ヒルベルト曲線の計算が大変) 167
168. メモリ確保は遅い 168
169. クイズ AtCoder での実行時間は？ vector A(1e7, vector<int>{1}); ヒント 1 が 1 個入った vector を 個作っている 107 169
170. 答え 170
171. メモリ確保は遅い！ メモリ確保は とされている (？) が，100 クロックくらいは かかってしまう → std::vector のメモリ再確保や平衡二分探索木における 　 new など，逐一メモリ確保する操作が遅い！ 参考資料 malloc の旅 (glibc 編) – Motohiro KOSAKI 動画版 mallocの動作を追いかける – @kaityo256 O(1) 171
172. メモリ確保は遅い！ C/C++ では malloc 関数を使って，OS から ざっくり割り当ててもらったメモリを細分化して渡している メモリ確保は とされている (？) が，100 クロックくらいは かかってしまう 逆操作の free も同様に遅い O(1) 172
173. 対策① std::vector の個数を減らす 小さな std::vector をたくさん作ったり，たくさんの std::vector に push_back したりするとメモリ確保が たくさん起こって遅い 長さが固定ならスタック領域を使う std::array にしたり， 多次元の std::vector を flatten したりして std::vector の 個数を減らすと良い 173
174. 対策② メモリ解放しない (= メモリリーク) メモリ解放せずプログラムを終了すれば，メモリ解放の分の 実行時間を減らせる OS が割り当てたメモリは後で OS がいい感じにしてくれる 174
175. 対策③ メモリ確保を 1 回にまとめる 十分な量のメモリを最初に確保しておき，それを前から 切り分けていく メモリ解放はしない 確保したい量の総和が計算できるなら，計算するのも良い 175
176. ポインタ木は遅い 同じ でも， ポインタ木（ std::set など）の は重い セグメント木の は軽い のはなぜか…？ log log log 176
177. ポインタ木は遅い ポインタ木 では，親ノードのメモリにアクセスしないと 子ノードの位置がわからない → メモリアクセスが直列 セグメント木 では，親ノードのメモリにアクセスしなくても 子ノードの位置がわかる → 並列にメモリアクセスできる！ 177
178. お待たせしました ベクトル化 やっていきます 178
179. ベクトル化 並列化可能な繰り返し処理を，専用命令 (SIMD 命令) を用いてまと めて処理する SIMD (Single Instruction Multiple Data) 命令 1 回の命令で，複数のデータに対し，並列に同じ処理を行う命令 179
180. 例 u32 A[8]; rep(i, 8) A[i] = 1; ↓ コンパイル 命令 意味 movl $1, %eax %eax に を代入 vpbroadcastd %eax, %ymm0 %ymm0 レジスタ (256 bit) を 8 分割， それぞれに %eax を代入 vmovdqa %ymm0, A(%rip) メモリ上で， A のある位置から 始まる 256 bit に %ymm0 を代入 1 180
181. 命令 意味 movl $1, %eax %eax に を代入 vpbroadcastd %eax, %ymm0 %ymm0 レジスタ (256 bit) を 8 分割， それぞれに %eax を代入 vmovdqa %ymm0, A(%rip) メモリ上で， A のある位置から 始まる 256 bit に %ymm0 を代入 1 命令で 8 要素に同時に代入ができた！ 1 181
182. SIMD 命令の歴史 MMX : 64 bit SSE シリーズ : 128 bit AVX シリーズ : 256 bit AVX-512 : 512 bit 現在 AtCoder, QOJ (UCup), CodeChef 等で使える 182
183. SIMD 命令一覧 参考資料 : インテル® C++ コンパイラー 17.0 デベロッパー・ ガイドおよびリファレンス 183
184. SIMD 命令一覧 Icelake, 512 bit，メモリアクセスのない場合 命令 意味 Lat. TP VPADDB 8 bit 整数 64 個を同時に加算 1 0.5 VPADDW 16 bit 整数 32 個を同時に加算 1 0.5 VPADDD 32 bit 整数 16 個を同時に加算 1 0.5 VPADDQ 64 bit 整数 8 個を同時に加算 1 0.5 184
185. 命令 意味 Lat. TP VPADDD 32 bit 整数 16 個を同時に加算 1 0.5 32 bit 整数 16 個同時に加算が 1 クロックで！？ 185
186. 命令 意味 Lat. TP VPADDB 8 bit 整数 64 個を同時に加算 1 0.5 VPADDW 16 bit 整数 32 個を同時に加算 1 0.5 VPADDD 32 bit 整数 16 個を同時に加算 1 0.5 VPADDQ 64 bit 整数 8 個を同時に加算 1 0.5 キャッシュに乗らないとすぐメモリ律速になります (それでも十分速い) 精度を小さくするほど速い！ 186
187. Icelake, 512 bit，メモリアクセスのない場合 命令 意味 Lat TP VPADDQ i64 の加算 * 8 1 0.5 VPAND 512 bit の AND 1 0.5 VPSRLVQ i64 の右シフト * 8 1 1 VPMAXSQ i64 の max * 8 3 1 VPOPCNTQ u64 の popcnt * 8 3 1 VPCMPQ i64 の 比較 * 8 5 1 187
188. 命令 意味 Lat TP VADDPD f64 の加算 * 8 4 1 VMULPD f64 の乗算 * 8 4 1 VCVTQQ2PD i64 → f64 * 8 4 1 VPBROADCASTQ i64 を 8 個に複製 ≤6 1 VPMULLQ i64 の乗算 (mod ) * 8 15 3 VDIVPD f64 の除算 * 8 ≤23 ≤16 VSQRTPD ≤32 ≤18 整数除算の命令が存在しない！ 264 ​ x 188
189. ベクトルレジスタ 浮動小数点数の計算や SIMD 命令で使うのがベクトルレジスタ レジスタ 長さ 意味 XMM レジスタ 128 bit ZMM レジスタの先頭 1/4 YMM レジスタ 256 bit ZMM レジスタの先頭 1/2 ZMM レジスタ 512 bit 189
190. ジャッジで使われている CPU を確認するには… Bash で cat /proc/cpuinfo を実行 Bash が使えない場合は，https://gist.github.com/t-mat/3769328 のようなコードを実行 190
191. ベクトル化されたコードを書く インテル® C++ コンパイラー 17.0 デベロッパー・ガイドおよびリファレンス とかを いろいろ見ながら #include <immintrin.hpp> u32 A[1024]; int main() { __m512i T0 = _mm512_loadu_si512((__m512i*)A); // (略) } と頑張って書くのは大変なので…？ 191
192. コンパイラに任せよう 1. 並列化しやすいようにプログラムを書く A[i] += B[i]; のように， シーケンシャルアクセスで同じ操作をするコードを書く AVX-512 の命令一覧を見て，それにコンパイルできそうな コードを書く Compiler Explorer を見ながら， 思った通りにベクトル化されていることを確認する 192
193. コンパイラに任せよう 1. 並列化しやすいようにプログラムを書く 2. 最適化オプション -ftree-vectorize をつけてコンパイルする 193
194. コンパイラに任せよう 1. 並列化しやすいようにプログラムを書く 2. 最適化オプション -ftree-vectorize をつけてコンパイルする 競プロではコンパイルオプションが変更できないので， プログラム中に最適化オプションを書く！ #pragma GCC optimize("tree-vectorize") ↑ これ以下に書かれたコードに -ftree-vectorize を適用する 194
195. 最適化オプションを適用 #pragma GCC target("arch=icelake-server") ↑ これ以下に書かれたコードに -march=icelake-server を適用， AVX-512 を使って良いことを教える (AtCoder は -march=native が指定されているので必要ない) 195
196. 最適化オプションを適用 #pragma GCC optimize("Ofast") ↑ これ以下に書かれたコードに -Ofast を適用， -Ofast には 　 -ftree-vectorize が含まれる． ( -Ofast は浮動小数の結果が変わって危険とか言われたりするが， そういう最適化は競プロでは大歓迎) 196
197. ベクトル化の注意点 ジャッジが対応している命令より新しい命令を 実行してしまうと， Illigal Instruction で RE になる → ジャッジが使いたい命令に対応していることを確認しよう！ 197
198. ベクトル化の注意点 重い SIMD 命令を使用するとクロック周波数が低下する – Xeon Platinum 8180 の場合 : WikiChip 競プロでは 1 コアしか使わないので，クロック低下は 5 〜 10% 程度 198
199. ベクトル化の注意点 重い SIMD 命令を使用するとクロック周波数が低下する → AtCoder では -march=native がついているが， -march=native でコンパイルすると，256 bit 幅の方が 512 bit 幅より速いと判断されて 512 bit 幅を使ってくれない！ 1 コア (競プロ用途なら) なら 512 bit 幅の方が速いので， -mperfer-vector-width=512 を指定して 512 bit 幅を 使ってもらう 199
200. 定数倍高速化の極地 キャッシュとベクトル化をうまく使ってチューニングすると驚異的な 速度になる Intel の x86-simd-sort std::sort より最大 10 倍高速なクイックソート Nyaan さんの AVX2 FFT ACL より 2 倍高速な FFT QCFium さんの quick_floyd_warshall 普通に書くより 8 倍高速な Floyd–Warshall 法 200
201. ベクトル化ライブラリ Google の highway std::simd; 将来標準ライブラリに入るかも…？ 201
202. 高速なデータ構造を使う 区間和のセグ木の代わりに Fenwick Tree Dijkstra 法の priority_queue の代わりに Radix Heap std::set<int> の代わりに Fast Set (word-size tree) 202
203. まとめ CPU・メモリの仕組み，アルゴリズムを理解して，プログラムの 定数倍を予測・改善できるようになろう！ 203