Lock in $30 Savings on PRO—Offer Ends Soon! ⏳

定数倍高速化の技術

tatyam
July 16, 2024

 定数倍高速化の技術

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

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

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

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

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

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

tatyam

July 16, 2024
Tweet

More Decks by tatyam

Other Decks in Technology

Transcript

  1. 試験 (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
  2. 私の解法 2 次元クエリをそのまま処理する 領域木で 時間 なら間に合う! 領域木 : 座標の区間でセグメント木を作り, セグメント木の各ノードが,

    「 座標がその範囲に含まれる点の 座標のリスト」をソートしたものを持つ O(N log N + Q(log N) ) 2 N, Q ≤ 105 x x y 22
  3. 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
  4. 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
  5. 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
  6. #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
  7. 参考資料 (CPU のつくりかた) ディスクリート半導体の基礎 第1章 半導体の基礎 – TOSHIBA ディスクリート半導体の基礎 第3章

    トランジスタ – TOSHIBA コンピュータ講座 応用編 第1回 - 富士通 ASCII.jp:半導体プロセスまるわかり トランジスタの配線と 形成 35
  8. トランジスタ (MOSFET) p p Source Drain Gate n 電流の⼀⽅通⾏ p

    型半導体と n 型半導体をうまく組み合わせると, 電気を通す・通さないを制御できるようになる! 45
  9. トランジスタ (MOSFET) n p p Source Drain Gate n 電流の⼀⽅通⾏

    ー ー ー ー S - G 間に電圧をかけてゲートに負の電荷をためると, 近くの電子が追い出されて正孔ができ,電気を通すようになる 46
  10. 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
  11. レジスタ これを 64 個繋げれば,64 bit 整数を 保持できる C が ON

    になるたび D から 与えられた 64 bit 整数を保存し, 次に C が ON になるまで Q に出力 これを レジスタ と呼ぶ. 53
  12. スイッチ + レジスタ 加算器 1 を出⼒ 太⽮印は 64 bit 計算のしくみ

    64 bit のレジスタと 加算器でこんな回路を 作ってみる 54
  13. スイッチ + レジスタ 加算器 1 を出⼒ 太⽮印は 64 bit 計算のしくみ

    スイッチを ON にすると, レジスタは後ろからの 入力を保存し… 55
  14. スイッチ + レジスタ 加算器 1 を出⼒ 太⽮印は 64 bit 計算のしくみ

    スイッチを ON にすると, レジスタは後ろからの 入力を保存し,前に出力 56
  15. スイッチ + レジスタ 加算器 1 を出⼒ 太⽮印は 64 bit 計算のしくみ

    最後にスイッチを OFF に 戻して,1 クロックが 終了 58
  16. スイッチ + レジスタ 加算器 1 を出⼒ 太⽮印は 64 bit 計算のしくみ

    スイッチを ON にすると, レジスタは後ろからの 入力を保存し,前に出力 59
  17. スイッチ + レジスタ 加算器 1 を出⼒ 太⽮印は 64 bit 計算のしくみ

    最後にスイッチを OFF に 戻して,1 クロックが 終了 61
  18. スイッチ + レジスタ 加算器 1 を出⼒ 太⽮印は 64 bit 計算のしくみ

    スイッチが ON になった 回数がレジスタに 記録されている! 62
  19. + - × ÷ MUX 命令 命令に対応した 計算結果 万能演算装置 足し算以外も

    計算できるようにしよう! → 引き算や掛け算, 割り算などの演算回路を 用意して,命令によって 使い分ける 63
  20. スイッチ レジスタ 1 を出⼒ 太⽮印は 64 bit 万能演算装置 + ×

    − ÷ 命令 ⾜し算 万能演算装置 色々な演算に対応した 万能演算装置を入れ, 命令を与えて計算する ように 64
  21. スイッチ レジスタ 万能演算装置 太⽮印は 64 bit + × − ÷

    命令 ⾜し算 MUX レジスタ 番号 計算のしくみ レジスタをたくさん増やした スイッチを ON にするたび, 「レジスタ A の値と レジスタ B の値を掛けて, レジスタ C に保存」 のような命令を処理する 65
  22. クロック 発振器 レジスタ 万能演算装置 太⽮印は 64 bit + × −

    ÷ 命令 ⾜し算 MUX レジスタ 番号 計算のしくみ スイッチをクロック発振器 (一定の周期で 0 / 1 を繰り返す, 音叉のようなもの) に変更して,計算が自動で 進むように 66
  23. クロック 発振器 レジスタ 万能演算装置 太⽮印は 64 bit + × −

    ÷ 命令 ⾜し算 MUX レジスタ 番号 計算のしくみ 1. クロック発振器が ON になり, レジスタが後ろの入力を 保存して前に出力 2. 演算装置が計算を実行し,計算 結果がレジスタの後ろまで伝播 3. クロック発振器が OFF になる を 1 クロックとして,これを高速に 繰り返す → 計算機の完成! 67
  24. クロック 発振器 レジスタ 万能演算装置 太⽮印は 64 bit + × −

    ÷ 命令 ⾜し算 MUX レジスタ 番号 計算のしくみ 1. クロック発振器が ON になり, レジスタが後ろの入力を保存して前に 出力 2. 演算装置が計算を実行し,計算結果が レジスタの後ろまで伝播 3. クロック発振器が OFF になる もしクロック発振器が再び ON になるまでに 計算が終わらないと…? 68
  25. クロック 発振器 レジスタ 万能演算装置 太⽮印は 64 bit + × −

    ÷ 命令 ⾜し算 MUX レジスタ 番号 1. クロック発振器が ON になり, レジスタが後ろの入力を保存して前に 出力 2. 演算装置が計算を実行し,計算結果が レジスタの後ろまで伝播 3. クロック発振器が OFF になる もしクロック発振器が再び ON になるまでに 計算が終わらないと…? → 間違った計算結果がレジスタに 保存されてしまう! 69
  26. 参考情報 x86_64 (AVX-512 対応) の (論理) レジスタの数 整数レジスタ : 64

    bit × 16 個 ベクトルレジスタ : 512 bit × 32 個 小数の計算はこっち レジスタで足りない分はメモリに置く必要がある 72
  27. 参考情報 最近の CPU の速度 クロック周波数 : 2 ~ 4 GHz

    3 GHz のとき,1 クロックに光は 10 cm しか進めない! この時間内にすべての計算が終わるように設計され,回路が 詰め込まれている 73
  28. 76

  29. ヒント 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
  30. 計算してみよう レジスタの 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
  31. 色々な演算の速度を知ろう 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
  32. 命令 意味 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
  33. 最近 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
  34. 素数判定するコードが…? 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
  35. 整数除算が遅い場合の対策 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
  36. ところで, 除算が並列にできるってなんだ? 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
  37. 実際のパイプライン処理 参考資料 : Sunny Cove - Microarchitectures - Intel -

    WikiChip Sunny cove block diagram.png by Chipwikia; CC BY-SA 4.0 95
  38. if(rand() % 2) { // ... } else { //

    ... } 1/2 の確率で分岐予測に失敗し, やり直しに何クロックもかかる → 遅い! for(int i = 0; i < 100; i++) { // ... // ... // ... } for 文の条件式はたいてい true → 分岐予測が基本的に成功し,   1 クロック程度で実行できる 98
  39. クイズ 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
  40. ヒント 掛け算と割り算を 回 さっきは並列になったが, 今回は並列にならない 乗算は 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
  41. 計算してみよう 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
  42. コンパイル結果を見る s %= MOD; に対応する命令が 青くなっている movq : 64 bit

    代入 mulq : 64 bit 乗算 shrq : 64 bit 右シフト subq : 64 bit 減算 除算命令が消えている! 112
  43. 除算の最適化 割る数がコンパイル時定数のとき,除算を掛け算 + 右シフトで 行うことがある ざっくりこんな気持ち 整数部 64 bit,小数部 64

    bit の小数で掛け算を行う で割る代わりに, を 整数部 64 bit,小数部 64 bit の 小数で表しておき,これを掛けて整数部を取る 参考資料 : コンパイラによる整数除算最適化の証明 | nu50218 blog x 1/x 113
  44. アセンブラを読もう! 左のコードを -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
  45. アセンブラを読もう! 命令 意味 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
  46. つまり… u64 div(u64 x) { return x / 3; }

    は u64 div(u64 x) { return (u128)x * 0xAAAAAAAAAAAAAAAB >> 65; } に最適化された 117
  47. 豆知識 : 命令の 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
  48. 豆知識 : 乗算命令 C++ で 64 bit 整数の乗算をすると,結果は mod で返ってくる.

    一方,64 bit 乗算命令は… 命令 意味 mulq %rdi %rax * %rdi の上位 64 bit を %rdx に, 下位 64 bit を %rax に 代入 128 bit で計算されている!? 264 122
  49. ご存知でしたか? – https://www.felixcloutier.com/x86/div 64 bit 除算も 128 bit ÷ 64

    bit ができる! (商が 64 bit に収まらないと RE ) 124
  50. 最適化の例 2 冪の掛け算・割り算はビット演算に最適化される 計算 最適化された命令 (例) x / 32 shrq

    $5, %rax x % 32 andl $31, %eax x * 33 salq $5, %rax addq %rdi, %rax 符号なし整数にした方がこの最適化には有利かも 125
  51. 最適化の例 関数呼び出しを削除して,関数の中身を埋め込む (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
  52. クイズ 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
  53. ヒント 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
  54. メモリアクセスの速度 メモリアクセスのレイテンシはだいたいこれくらい 大きさ (目安) レイテンシ (目安) レジスタ 2 KB 1

    クロック メモリ 8 GB ~ ≥100 クロック メモリ制限はたいてい 1 GB だけど… メモリアクセスは遅い! 132
  55. メモリアクセスの速度 メモリアクセスのレイテンシはだいたいこれくらい 大きさ (目安) レイテンシ (目安) レジスタ 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
  56. 48 KB まで L1 キャッシュが効いて,平均 2.3 クロックくらい 1.25 MB まで

    L2 キャッシュが効いて,平均 3.5 クロックくらい 10 MB くらいまで L3 キャッシュが効いて,平均 7 クロックくらい それ以上はどんどん遅くなっていく 136
  57. 0x00000000 0x00000040 0x00000080 0x000000C0 0x00000100 0x00000140 0x00000180 0x000001C0 0x00000200 0x00000240

    0x00000280 0x000002C0 0x00000300 0x00000340 キャッシュメモリのしくみ メモリ・キャッシュメモリは 64 バイト単位 で読み書きする この 64 バイト単位を ライン と呼ぶ 137
  58. キャッシュメモリのしくみ (L1 キャッシュの場合) L1 キャッシュの大きさ : 注 : 大きさ は

    CPU により異なります ラインからなるグループが 64 個ある 64 × 12 × 64 Bytes 12 140
  59. キャッシュメモリのしくみ (L1 キャッシュの場合) 0x12307080 の値が欲しい! → 64 で割って,ライン番号 0x48C1C2 →

    ライン番号 mod 64 は 2 なので,グループ 2 の中を探す 0x12307080 0x12307080 141
  60. キャッシュメモリのしくみ (L2 キャッシュの場合) L2 キャッシュの大きさ : 注 : 大きさ は

    CPU により異なります ラインからなるグループが 個ある 1024 × 20 × 64 Bytes 20 1024 144
  61. クイズ 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
  62. ヒント のワーシャルフロイド法 の入力生成部分は無視できる 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
  63. 実測 #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
  64. 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
  65. どうしてもシーケンシャルアクセスにできないとき 行列の転置 u64 A[3000][3000]; rep(i, 3000) rep(j, i) { swap(A[i][j],

    A[j][i]); } A[i][j] か A[j][i] のどちらかはシーケンシャルアクセスに できない! 156
  66. A[0][0] A[0][7] A[7][0] A[7][7] ブロック化 A[i][j] の の範囲と の範囲を いくつかのブロックに分け,

    ブロック内のアクセスを キャッシュに乗せるテクニック i j 161
  67. 1 2 3 4 5 6 7 8 9 10

    11 12 13 14 15 16 Cache-Oblivious 1. 全体を 4 個のブロックに 分割し,ブロックごとに 転置を行う 2. 各ブロックを再帰的に 4 個のブロックに分割し, これを繰り返す 165
  68. 1 2 3 4 5 6 7 8 9 10

    11 12 13 14 15 16 Cache-Oblivious 1. 全体を 4 個のブロックに分割し, ブロックごとに転置を行う 2. 各ブロックを再帰的に 4 個のブロックに 分割し,これを繰り返す 3. ブロックの大きさが確実にキャッシュに 乗るところで再帰を止めてループで処理 → 「ブロックごとに処理」を様々な 大きさでできるので,どんな大きさの キャッシュも効率的に使える! 166
  69. 1 2 3 4 5 6 7 8 9 10

    11 12 13 14 15 16 Cache-Oblivious ちょっと順序を工夫すると ヒルベルト曲線 ができる → ヒルベルト曲線の順で 処理すると Cache- Oblivious! (ヒルベルト曲線の計算が大変) 167
  70. メモリ確保は遅い! メモリ確保は とされている (?) が,100 クロックくらいは かかってしまう → std::vector のメモリ再確保や平衡二分探索木における

      new など,逐一メモリ確保する操作が遅い! 参考資料 malloc の旅 (glibc 編) – Motohiro KOSAKI 動画版 mallocの動作を追いかける – @kaityo256 O(1) 171
  71. 対策① std::vector の個数を減らす 小さな std::vector をたくさん作ったり,たくさんの std::vector に push_back したりするとメモリ確保が

    たくさん起こって遅い 長さが固定ならスタック領域を使う std::array にしたり, 多次元の std::vector を flatten したりして std::vector の 個数を減らすと良い 173
  72. ベクトル化 並列化可能な繰り返し処理を,専用命令 (SIMD 命令) を用いてまと めて処理する SIMD (Single Instruction Multiple

    Data) 命令 1 回の命令で,複数のデータに対し,並列に同じ処理を行う命令 179
  73. 例 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
  74. 命令 意味 movl $1, %eax %eax に を代入 vpbroadcastd %eax,

    %ymm0 %ymm0 レジスタ (256 bit) を 8 分割, それぞれに %eax を代入 vmovdqa %ymm0, A(%rip) メモリ上で, A のある位置から 始まる 256 bit に %ymm0 を代入 1 命令で 8 要素に同時に代入ができた! 1 181
  75. SIMD 命令の歴史 MMX : 64 bit SSE シリーズ : 128

    bit AVX シリーズ : 256 bit AVX-512 : 512 bit 現在 AtCoder, QOJ (UCup), CodeChef 等で使える 182
  76. 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
  77. 命令 意味 Lat. TP VPADDD 32 bit 整数 16 個を同時に加算

    1 0.5 32 bit 整数 16 個同時に加算が 1 クロックで!? 185
  78. 命令 意味 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
  79. 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
  80. 命令 意味 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
  81. ベクトルレジスタ 浮動小数点数の計算や SIMD 命令で使うのがベクトルレジスタ レジスタ 長さ 意味 XMM レジスタ 128

    bit ZMM レジスタの先頭 1/4 YMM レジスタ 256 bit ZMM レジスタの先頭 1/2 ZMM レジスタ 512 bit 189
  82. ベクトル化されたコードを書く インテル® C++ コンパイラー 17.0 デベロッパー・ガイドおよびリファレンス とかを いろいろ見ながら #include <immintrin.hpp>

    u32 A[1024]; int main() { __m512i T0 = _mm512_loadu_si512((__m512i*)A); // (略) } と頑張って書くのは大変なので…? 191
  83. 最適化オプションを適用 #pragma GCC optimize("Ofast") ↑ これ以下に書かれたコードに -Ofast を適用, -Ofast には

      -ftree-vectorize が含まれる. ( -Ofast は浮動小数の結果が変わって危険とか言われたりするが, そういう最適化は競プロでは大歓迎) 196
  84. ベクトル化の注意点 重い SIMD 命令を使用するとクロック周波数が低下する – Xeon Platinum 8180 の場合 :

    WikiChip 競プロでは 1 コアしか使わないので,クロック低下は 5 〜 10% 程度 198
  85. ベクトル化の注意点 重い SIMD 命令を使用するとクロック周波数が低下する → AtCoder では -march=native がついているが, -march=native

    でコンパイルすると,256 bit 幅の方が 512 bit 幅より速いと判断されて 512 bit 幅を使ってくれない! 1 コア (競プロ用途なら) なら 512 bit 幅の方が速いので, -mperfer-vector-width=512 を指定して 512 bit 幅を 使ってもらう 199
  86. 定数倍高速化の極地 キャッシュとベクトル化をうまく使ってチューニングすると驚異的な 速度になる Intel の x86-simd-sort std::sort より最大 10 倍高速なクイックソート

    Nyaan さんの AVX2 FFT ACL より 2 倍高速な FFT QCFium さんの quick_floyd_warshall 普通に書くより 8 倍高速な Floyd–Warshall 法 200