トヨタ自動車プログラミングコンテスト2024#6（AtCoder Heuristic Contest 034) 解法

解法

最小費用流 (7.9G)

実はこの問題、トラックの移動経路を固定すると、最適解を求めることが出来ます。最適解というのは、$h_{ij} < 0$なのに更に掘り出すケース、あるいはトラックが複数回同じ頂点を通るケース等も対応した上での厳密な最適解です。

うまくグラフを作って最小費用流を流すのですが、きりさんのこのtweetの図がわかりやすいです。

#AHC034

経路を決めたときの最小費用流は、経路の各点に対応する頂点からその位置に対応する頂点への双方向の道を作る感じでやったhttps://t.co/Y3yTpp7uCr pic.twitter.com/UjAUbe2hcD

— きり (@kiri8128) 2024年6月16日

文章で書くと、トラックの移動経路の頂点数を $M$ として、

• 移動経路に対応した頂点が $M$ 個
• グリッドのマス目に対応した頂点が $N^{2}$ 個
• 始点 $S$, 終点 $T$

で計 $M + N^{2} + 2$ 頂点のグラフを作ります。そして、次のように辺を張ります。

• 各グリッドのマス目$(i, j)$ の頂点について、$h_{ij} > 0$ ならば $S$ から (cap, cost) = $(h_{ij}, 0)$ の辺を張る
• 各グリッドのマス目$(i, j)$ の頂点について、$h_{ij} < 0$ ならば $T$ へ $(-h_{ij}, 0)$ の辺を張る
• 各移動経路の頂点について、対応するグリッドのマス目の頂点との間に、$(\infty, 1)$ の辺を双方向に張る
• 各 $i = 1, 2, \cdots, M - 1$について、 $i$ 番目の移動経路の頂点から $i + 1$ 番目の移動経路の頂点へ $(\infty, 1)$ の辺を張る

実際に、このグラフで $S$ から $T$ に最小費用最大流を流すと、

• 移動経路の頂点とグリッドのマス目の頂点の間の辺の流量が 積み込み/積み下ろし の量に対応し、
• 移動経路の頂点間の辺の流量がトラックの運ぶ量に対応する

ことがわかります。

トラックの移動経路としてなにを取るかですが、次の図のように(ジグザグ + 向きを変えたジグザグ)で盤面を二回ずつ通るようなwalkを採用すると7.9G程度のスコアを取ることが出来ます。

実装例: Submission #54667347 - Toyota Programming Contest 2024#6（AtCoder Heuristic Contest 034）

ちなみにこの最小費用流の嬉しいポイントとして、あらゆる解法に対して最後にやって損しない後処理として機能します。つまり実装しても絶対無駄にならないです。

山登り(9.4G)

さらなる改善のためには、よりよい移動経路を探す必要があります。現状だとそもそも完全に経路を決め打ってしまっているので、改善の余地はありそうです。

先述の解法を動かしてみると、運搬コストについてはかなりよさそうで、トラックの移動コスト $100 \times (M - 1)$ の方に改善幅がありそうな気持ちになります。よって、トラックの移動距離を短くするような遷移で山登りをすることを考えます。

この思想を元に、「現在の移動経路について、$i$ 番目の頂点と $i + 3$ 番目の頂点が隣接していたならば、$i + 1$ 番目と $i + 2$ 番目の頂点を削除する」という遷移を考えます。これは上記のジグザグにおいて、曲がる部分をちょっとショートカットして短くする遷移です。

実際にこの遷移を書くと、とんでもなく強力なことがわかります。なんと時間いっぱい山登りするだけで9.4G程度のスコアが取れます。

実装例: Submission #54668251 - Toyota Programming Contest 2024#6（AtCoder Heuristic Contest 034）

なお、更なる改善として自然なのは焼きなましですが、ここで最小費用流(=スコア関数)の速度が問題になります。seed 0で、AC Libraryだと600 ~ 700回程度、後述の強力最小費用流ライブラリを使っても1500回程度しか最小費用流が流せなかったので、焼きなましは厳しかったです。うまくやればこの程度の遷移回数でもなんとかなるんでしょうか？

さらなる改善(9.5G)

コンテスト中に9.48G, 後に9.53Gまで伸ばしましたが、正直もう本質的な改善は行っておらず、泥臭いことをしました。

• より高速な最小費用流を使う
  • AC LibraryはPrimal-Dual法の実装としてそれなりに速いものになっていると思いますが、そもそも最小費用流に対してはNetwork Simplex / Cost Scalingなど高速なアルゴリズムが色々あります。LEMON (C++ library) - Wikipedia がC++の実装として有名で、実際にNetwork Simplexを使うと3倍程度は速くなるようです。
• 多点スタート
  • 高速な最小費用流を使うと実行時間に余裕が出来るので多点スタートをするとスコアが上がります。9.53Gの提出ではジグザグの向きが2通り * 2回ずつ 375ms(=計1500ms)山登り -> 一番良かったやつを更に400ms山登り
• 雰囲気で適当に山登りの遷移を足して、仕上げの山登りで使用

競プロ 乱数 速度調査

疑似乱数の生成速度について軽く調べた。まず結果から紹介する。それぞれ $10^8$ 回乱数を生成してかかった時間をAtCoderのコードテストで計測した。実験コードはここ

Name	Output bit	Time
LCG32	32bit	114ms
mt19937	32bit	325ms
mt19937_64	64bit	389ms
xorshift32	32bit	192ms
xorshift64	64bit	193ms
xoshiro128++	32bit	249ms
xoshiro128**	32bit	238ms
xoshiro256++	64bit	249ms
xoshiro256**	64bit	238ms
pcg32	32bit	134ms
pcg32_fast	32bit	105ms
pcg64	64bit	202ms
pcg64_fast	64bit	157ms
Mwc128XXA32	32bit	92ms
Mwc256XXA64	64bit	108ms

Output bitで分類すると次のようになる。

Name	Output bit	Time
LCG32	32bit	114ms
mt19937	32bit	325ms
xorshift32	32bit	192ms
xoshiro128++	32bit	249ms
xoshiro128**	32bit	238ms
pcg32	32bit	134ms
pcg32_fast	32bit	105ms
Mwc128XXA32	32bit	92ms

Name	Output bit	Time
mt19937_64	64bit	389ms
xorshift64	64bit	193ms
xoshiro256++	64bit	249ms
xoshiro256**	64bit	238ms
pcg64	64bit	202ms
pcg64_fast	64bit	157ms
Mwc256XXA64	64bit	108ms

LCG (線形合同法)

$X = (A \times X + B) \bmod M$ という形で表されるやつ。パラメーターはwikipediaから(あまりよくないパラメーターらしい)

// Reference:
// https://ja.wikipedia.org/wiki/%E7%B7%9A%E5%BD%A2%E5%90%88%E5%90%8C%E6%B3%95
namespace lcg32 {

inline static u32 state = 12345;

u32 next() { return state = 1'103'515'245 * state + 12'345; }

}  // namespace lcg32

mt19937 (メルセンヌ・ツイスター)

C++のSTLのmt19937 / mt19937_64をそのまま

// Reference:
// https://cpprefjp.github.io/reference/random/mt19937.html
namespace mt19937 {

inline static std::mt19937 engine(12345);

u32 next() { return (u32)engine(); }

}  // namespace mt19937

// Reference:
// https://cpprefjp.github.io/reference/random/mt19937_64.html
namespace mt19937_64 {

inline static std::mt19937_64 engine(12345);

u64 next() { return (u64)engine(); }

}  // namespace mt19937_64

xorshift

主観だと競プロで一番使われていそうなやつ。実装はwikipediaから殆どそのまま

// Reference:
// https://ja.wikipedia.org/wiki/Xorshift
namespace xorshift32 {

inline static u32 a = 12345;

u32 next() {
    u32 x = a;
    x ^= x << 13;
    x ^= x >> 17;
    x ^= x << 5;
    return a = x;
}

}  // namespace xorshift32

// Reference:
// https://ja.wikipedia.org/wiki/Xorshift
namespace xorshift64 {

inline static u64 a = 12345;

u64 next() {
    u64 x = a;
    x ^= x << 13;
    x ^= x >> 7;
    x ^= x << 17;
    return a = x;
}

}  // namespace xorshift64

xoshiro++ / xoshiro**

https://prng.di.unimi.it/ に詳しい紹介がある。速くて質もいいらしい。実装はこのサイトから殆どそのまま

// Reference:
// https://prng.di.unimi.it/
namespace xoshiro128plusplus {

inline static u32 s[4] = {123, 234, 345, 567};

u32 next() {
    const u32 result_starstar = std::rotl(s[0] + s[3], 7) + s[0];

    const u32 t = s[1] << 9;

    s[2] ^= s[0];
    s[3] ^= s[1];
    s[1] ^= s[2];
    s[0] ^= s[3];

    s[2] ^= t;

    s[3] = std::rotl(s[3], 11);

    return result_starstar;
}

}  // namespace xoshiro128plusplus

// Reference:
// https://prng.di.unimi.it/
namespace xoshiro128starstar {

inline static u32 s[4] = {123, 234, 345, 567};

u32 next() {
    const u32 result_starstar = std::rotl(s[1] * 5, 7) * 9;

    const u32 t = s[1] << 9;

    s[2] ^= s[0];
    s[3] ^= s[1];
    s[1] ^= s[2];
    s[0] ^= s[3];

    s[2] ^= t;

    s[3] = std::rotl(s[3], 11);

    return result_starstar;
}

}  // namespace xoshiro128starstar

// Reference:
// https://prng.di.unimi.it/
namespace xoshiro256plusplus {

inline static u64 s[4] = {123, 234, 345, 567};

u64 next() {
    const u64 result_starstar = std::rotl(s[0] + s[3], 23) + s[0];

    const u64 t = s[1] << 17;

    s[2] ^= s[0];
    s[3] ^= s[1];
    s[1] ^= s[2];
    s[0] ^= s[3];

    s[2] ^= t;

    s[3] = std::rotl(s[3], 45);

    return result_starstar;
}

}  // namespace xoshiro256plusplus

// Reference:
// https://prng.di.unimi.it/
namespace xoshiro256starstar {

inline static u64 s[4] = {123, 234, 345, 567};

u64 next() {
    const u64 result_starstar = std::rotl(s[1] * 5, 7) * 9;

    const u64 t = s[1] << 17;

    s[2] ^= s[0];
    s[3] ^= s[1];
    s[1] ^= s[2];
    s[0] ^= s[3];

    s[2] ^= t;

    s[3] = std::rotl(s[3], 45);

    return result_starstar;
}

}  // namespace xoshiro256starstar

PGC

https://www.pcg-random.org/index.html に詳しい紹介がある。速くて質もいいらしい。

標準のpcg32 / pcg64と、ストリーム機能(とちょっと短い周期)の代わりに速度に特化したpcg32_fast / pcg64_fastがある。

pcg32 / pgc32_fast は wikipedia の実装を参考に、 pcg64 / pcg64_fast はrustの rand_pcg crateの実装を参考にした。

// Reference:
// https://www.pcg-random.org/download.html
// https://en.wikipedia.org/wiki/Permuted_congruential_generator
namespace pcg32 {

const u64 MULT = 6364136223846793005;
const u64 INC = 1442695040888963407;

inline static u64 state = 123;

u32 next() {
    u64 x = state;
    state = state * MULT + INC;

    u32 count = x >> 59;
    x ^= x >> 18;

    return std::rotr(u32(x >> 27), count);
}

}  // namespace pcg32

// Reference:
// https://www.pcg-random.org/download.html
// https://en.wikipedia.org/wiki/Permuted_congruential_generator
namespace pcg32_fast {

const u64 MULT = 6364136223846793005;

inline static u64 state = 123;

u32 next() {
    u64 x = state;
    state = state * MULT;

    u32 count = x >> 61;
    x ^= x >> 22;

    return (u32)(x >> (22 + count));
}

}  // namespace pcg32_fast

// Reference:
// https://www.pcg-random.org/download.html
// https://crates.io/crates/rand_pcg
namespace pcg64 {

const u128 MULT = u128(2549297995355413924) << 64 | u128(4865540595714422341);
const u128 INC = u128(6364136223846793005) << 64 | u128(1442695040888963407);

inline static u128 state = 123;

u64 next() {
    u128 x = state;
    state = state * MULT + INC;

    u32 rot = x >> 122;
    return std::rotr((u64)(x >> 64) ^ (u64)x, rot);
}

}  // namespace pcg64

// Reference:
// https://www.pcg-random.org/download.html
// https://crates.io/crates/rand_pcg
namespace pcg64_fast {

const u128 MULT = u128(2549297995355413924) << 64 | u128(4865540595714422341);

inline static u128 state = 123;

u64 next() {
    u128 x = state;
    state = state * MULT;

    u32 rot = x >> 122;
    return std::rotr((u64)(x >> 64) ^ (u64)x, rot);
}

}  // namespace pcg64_fast

Mwc128xxa32 / Mwc256xxa64

このブログで紹介されている。質が良くてPCGより速いらしい。

実装は githubから殆どそのまま。変更点はx1とcをまとめている。

// Reference:
// https://tom-kaitchuck.medium.com/designing-a-new-prng-1c4ffd27124d
// https://github.com/tkaitchuck/Mwc256XXA64
namespace mwc128xxa32 {

const u32 MULT = 3487286589;

inline static u32 x2 = 12345;
inline static u32 x3 = 0xcafef00d;
inline static u64 c_x1 = u64(0xd15ea5e5) << 32 | 23456;

u32 next() {
    u64 x = (u64)(x3)*MULT;
    u32 result = (x3 ^ x2) + ((u32)(c_x1) ^ (u32)(x >> 32));
    x3 = x2;
    x2 = (u32)(c_x1);
    c_x1 = x + (c_x1 >> 32);
    return result;
}

}  // namespace mwc128xxa32

// Reference:
// https://tom-kaitchuck.medium.com/designing-a-new-prng-1c4ffd27124d
// https://github.com/tkaitchuck/Mwc256XXA64
namespace mwc256xxa64 {

const u64 MULT = 0xfeb3'4465'7c0a'f413;

inline static u64 x2 = 12345;
inline static u64 x3 = 0xcafef00dd15ea5e5;
inline static u128 c_x1 = u128(0x1405'7B7E'F767'814F) << 64 | 23456;

u64 next() {
    u128 x = (u128)(x3)*MULT;
    u64 result = (x3 ^ x2) + ((u64)(c_x1) ^ (u64)(x >> 64));
    x3 = x2;
    x2 = (u64)(c_x1);
    c_x1 = x + (c_x1 >> 64);
    return result;
}

}  // namespace mwc256xxa64

高速剰余算 div2by1 実装してみた

div2by1というアルゴリズムがある -> https://gmplib.org/~tege/division-paper.pdf

これはBarrett reductionやMontgomery乗算と違い、

• (k, k) -> 2k-bitの乗算器でk-bitの剰余算ができる(Barrett reductionは(2k, 2k) -> 4k-bitの乗算器を必要とする)
• modが偶数でも動作する(Montgomery乗算はmodが奇数の必要がある)

という二つの性質を持つ。今回は次の6種類のプログラムを実装し、$(8 \times 10^{7})! \bmod 998244353$を計算して速度を比較した。

• 32-bit Montgomery
• 32-bit div2by1
• mod < $2^{30}$という制約を仮定し32-bit div2by1をちょっと改造したもの
• これら三種のアルゴリズムをAVX2で高速化したもの

ベンチマークのコード、およびdiv2by1改造版のpythonコードについては、div2by1 改造版 · GitHub

手元(Ryzen 5 5600X)、およびAtCoderのコードテストでの実行時間は次の通り。

	Local	AtCoder
naive	285ms	482ms
naive(const mod)	227ms	298ms
montgomery	202ms	257ms
montgomery AVX	28ms	51ms
div2by1	316ms	458ms
div2by1 AVX	42ms	107ms
my div2by1	292ms	383ms
my div2by1 AVX	52ms	83ms

次のようなことがわかる

• そもそもnaive(=除算命令を使っているはず)が直観よりかなり速い。おそらくIce Lakeから除算が速くなったというやつ(Intel Ice Lakeのプロセッサは整数除算命令がアツい - chroot("/home/hibari"))。でもIce lakeの製造開始は2019年らしい、老人さん？笑
• montgomery + AVX2が速すぎる。4倍速以上になるとは思ってなかった。
• div2by1 あんまり速くない。montgomery速すぎ + AVXで高速化とかするならもうmod 998244353決め打ちで問題なさそう　であることを考えると正直出番がなさそうな

区間mul 区間積 O(log N)

問題

長さNの整数列a_iが与えられます　$Q$個のクエリを処理してください

• given l, r, x: $a_l \cdots a_r$を$x$倍
• given l, r: $a_l \times a_{l+1} \times \cdots \times a_r \bmod 998244353$を出力

$O(\log N)$ per queryで解けるがおそらく$O( \log^{2} N)$と識別不可能。

ちなみに元ネタはこれ(解法は違う): 区間代入/区間積 Θ(logN)/query - noshi91のメモ

$O(\log^{2} N)$ 解法

普通に遅延伝搬segtreeに乗せる。ノードには区間の総積と区間の長さを乗せる。ACL風に書くとS = pair<modint, int>

作用(mapping)の中でpow_modを呼ばないといけないため$O(\log^{2} N)$になる

$O(\log N)$ 解法

作用が可換なので遅延伝搬しない遅延伝搬segtree(何て呼ばれてるんでしょうこれ)が使える。 遅延伝搬しないverは次のようになる。もちろん素直に実装すると$O(\log^{2} N)$になるのだが、よく考えるとどちらも$O(\log N)$になる。

ノードごとに2つのmodint a, bを持つ。初期値はaが区間の総積でbが1

• mul: [l, r]をsegtreeの区間に分割する。分割された区間、およびその区間を子孫に持つすべての区間についてb *= x^([l, r]と自分の区間の共通部分の面積)
• prod: [l, r]をsegtreeの区間に分割する。分割された区間それぞれについて、自分と先祖のbのprodを求め、cとする。そしてa * c^(区間の長さ)を求める。これをすべての区間について掛け合わせる

mulについては次のように高速化する。

• [l, r]を分割した区間: bにかかる係数はx^(2^i)の形になっているので、まとめて$O(\log N)$で前計算できる
• それ以外の区間: 子のbにかかる係数の積を自分のbに掛ければよい

prodについては次のように高速化する。

• cについてはdfsしながらまとめて計算できるので、結局ある数列$d_1, d_2, ..., d_k (k = \log N)$について $d_1 \times d_2^{2} \times d_3^{4} \cdots$ が求められれば良い。これは $d_1 \times square(d_2 \times square(d_3 \times ...)))$という形で計算すれば $O(\log N)$

AHC030 環境構築 振り返り

AHC030に出て、11位でした。

自分は長期マラソンは殆どやったことがなく、2015年に3回topcoder MMに出たのが最初で最後(のはず)でした。 なので今回のAHCで環境整備系も全て一からやることになりました。なのであえてそちらについての感想や振り返りを書きます。

コンテスト中に使ったもの、また時間があったら欲しかったもの、を個人的に重要だと感じた順番で書いています。

ローカルテスター

Psyhoさんも言っています(https://twitter.com/FakePsyho/status/1605639454600806401)が、ローカルテスターが一番重要な環境整備要素でした。自分は適当な100行程度のpythonコードを準備して、改変しながら使っていました。

• 絶対スコア / 相対スコアの一覧を(csv)出力する
• 実行時間の一覧を出力する
• ケースをMでフィルターする

の3つの機能を含んでいて、今回の問題だとこれらは必須だと思います。

また、

• テストケースの並列実行機能

は結局実装しなかったのですが(なんで？)、必須級だと思います。3秒 * 100ケース回すと5分待ちなので、効率がすごい落ちた自覚があります。

個人的には複数の問題に対応できる強力なものを整備するよりは、pythonコードやらなんやらを問題ごとに毎回改変していくのが楽そうだなぁと思っています

Jupyter notebook (IPython notebook)

ちょっとした実験をすぐ書けて、プロットも出来るので非常に便利なツールでした。自分は今回は

• いろんなk, vで正規分布をプロットしてみる
• テストケースのMの分布を調べる
• その他細かい実験/計算

に使いました。インストールが簡単(vscodeならプラグイン入れるだけ)なので使い得だと思います。

　

でかつよCPU

先述の並列テストケース実行とも関係して、CPUは多ければ多いほどいいと思いました。先述のように自分の用意したローカルテスターは並列実行しないのですが、最終盤は結局そのローカルテスターをいろんなMに対して並列に走らせる、とかやっていてCPUが足りない状態になりました。

クラウドにデカいインスタンス借りてsshするのが"正解"なのは間違いないのですが、インスタンス立てたり落としたりするのがどうしても億劫な気持ちになってしまい…　手元のPCが強ければそれが一番良いのは間違いないと思います

でかつよクラウドインスタンス

上と同じ話です。最終盤は強いCPUか強いクラウドインスタンスのどちらかは欲しい

google spread sheet

先述のローカルテスターの出力を張ってスコアの比較に使っていました。 最初のほうは快適だったのですが、画像のシート一覧を見るとわかるように終盤になるにつれて限界になっていきました。改善の余地がありそうです(がspread sheet職人はやりたくない　うーん)

順位表のスナップショット機能 (未実装)

今回の問題は相対スコアなので、順位表のスコアが常時変動します。特に「自分の提出で他の人の点数がどのぐらい変わったか」は、自分の提出がbestを含むかを含んでおり、重要な情報でした。

そのため、提出直前の順位表を保存しておく必要があったのですが、注意力が低く複数回失敗しました。

考えられる対策としては定期的に順位表のスナップショットを取るスクリプトを走らせればいいです。でもどのぐらいの頻度なら怒られないのか、とか、そもそも参加者がみんな個人で定期的にスナップショットを取るというのは変な話なので、公式でスナップショットを提供してくれたら嬉しいなぁと思っています。

マージテクの逆でよく出てくる"2個の木のうち小さいほうを探す"処理ってcoroutineと相性がいいよね

背景

次のような問題を考えます

• 2個の木が与えられます。部分木の頂点数を$n, m$とした時に、$O(\min{(n, m)})$時間で小さいほうの部分木の頂点を列挙してください。

このような問題は「データ構造をマージする一般的なテクの逆」などと呼ばれるテクニックを使う問題で出てきます。具体例としては I - 盆栽 が一番有名だと思います。

冒頭の問題ですが、解法自体は対して難しくなく、「2つの木に並列にBFS/DFSして、どちらかが終わったら打ち切ればいい」というだけの話です。ですがいざ実装をしようとするとなかなか面倒です。しかもウニグラフ等で計算量が壊れがちだったりして厄介です。

実はこの実装はcoroutineと呼ばれる概念と相性が良いです。coroutineはC++だとC++20で入った機能 + 主な用途が並列処理やI/O bottleneckの処理等なので、おそらく競プロでの知名度は低いと思いますが、大体の新しめの言語には実装されている機能です。

実際に冒頭の問題を実装することを考えます。まず、$O(\max{(n, m)})$時間かけていいときの実装例を示します。ただ愚直にdfsをしているだけです。

using Tree = vector<vector<int>>;

void list_vertex(const Tree& tree, int u, int p, vector<int>& result) {
    result.push_back(u);
    for (int v : tree[u]) {
        if (v == p) continue;
        list_vertex(tree, v, p, result);
    }
}

vector<int> small_tree_vertex(const Tree& tree1, const Tree& tree2) {
    vector<int> result1, result2;
    list_vertex(tree1, 0, -1, result1);
    list_vertex(tree2, 0, -1, result2);

    if (result1.size() < result2.size()) {
        return result1;
    } else {
        return result2;
    }
}

これをcoroutineを使って実装すると次のようになります。

using Tree = vector<vector<int>>;

// https://github.com/lewissbaker/cppcoro/blob/master/include/cppcoro/recursive_generator.hpp
cppcoro::recursive_generator<int> list_vertex(const Tree& tree, int u, int p) {
    co_yield u;
    for (int v : tree[u]) {
        if (v == p) continue;
        co_yield list_vertex(tree, v, p);
    }
}

vector<int> small_tree_vertex(const Tree& tree1, const Tree& tree2) {
    vector<int> result1, result2;
    auto co1 = list_vertex(tree1, 0, -1);
    auto co2 = list_vertex(tree2, 0, -1);
    for (auto it1 = co1.begin(), it2 = co2.begin();; it1++, it2++) {
        if (it1 == co1.end()) return result1;
        if (it2 == co2.end()) return result2;
        result1.push_back(*it1);
        result2.push_back(*it2);
    }
}

少しlist_small_tree_vertexがごちゃごちゃしましたが、これで $O(\min{(n, m)})$ 時間で動作します。並列BFSを実装したことがあればなかなか驚きの実装量だと思います。また、C++23ならばstd::views::zipを使えばより簡潔な実装になるはずです。

coroutineというのは、ざっくり言うと「途中で中断と再開」が可能な関数です。実際に、新しいlist_vertex関数は、「頂点を見つけたら(= co_yield uにたどり着いたら)その頂点を返して関数を中断、そしてit++が呼ばれたらdfsをそこから再開」という挙動をします。なので、list_vertexの帰り値を普通のイテレーターのように扱い、どちらかのイテレーターが末尾に到達したらそこまでの結果を返すだけでよいです。

なお、C++だと再帰関数をcoroutineにするにはcppcoro::recursive_generatorのような追加実装が必要なようですが、MITライセンスで公開されているので適切にやれば自分で実装しなくても大丈夫です。

実際に盆栽を解いたコードはこちらです: Submission #49240549 - 東京大学プログラミングコンテスト2014 。冒頭(286行目まで)にこのrecursive_generatorが張り付けられているのでウォっとなりますが、それ以降だけ見ると結構簡潔ではないでしょうか。

-march=native 諸々

概要

-march=nativeについて色々調べた。

話題の発端

こちらのツイートであると思われる。

分からん
これ何かの未定義動作踏んでる？？ pic.twitter.com/ACV8fb8PjQ

— AllDirections (@AllDirections4) 2023年12月21日

シンプルなコードで、しかも-march=nativeを付けた場合のみ壊れる、ということで非常に力がある。自分もこれを機にmarchについて調べてしまった。

そもそもなぜ上記のコードは壊れているのか？

元のコードからC++要素を取り除くと次のようになる。もちろんこのコードも壊れていることがコードテストから確認できる。

#include <cstdio>
#include <cstring>
#include <cassert>

int main() {
    long long a[4] = {1, 1, 1, 0}, b[4];
    memmove(b, a, 4 * sizeof(long long));
    
    for (int i = 0; i < 4; i++) {
      printf("%lld ", b[i]);
    }
    printf("\n");
    return 0;
}

実際にアセンブリを確認すると (godbolt)、vpbroadcastq および vmovdqa で b を {1, 1, 1, 1}で上書きした後 b[3] に対して何もしていないことがわかる。

なお、-march=nativeが実際どう変換されているかはgcc -### -march=native /usr/include/stdlib.hで確認できる。AtCoderだと-march=icelake-server。

これと -march=native を付けた場合のみ壊れるという現象から、GCCのAVX512周りになにかバグがあるのだろうと検討が付く。実際にGCCのissue trackerを眺めると、どうやら今回のバグはこれっぽい 108599 – [12 Regression] Incorrect code generation newer intel architectures 。12.3で修正されているので、AtCoderのGCCがアップデートされればこの問題は解決されそう(いつだろう)。

-march=native -mtune=nativeってそもそも何？

とても雑に言うと

• -march=native: コンパイルしたパソコン(のCPU)専用のa.outを作ってくれという命令。生成されたa.outを他のパソコンにコピーすると、動くかもしれないし動かないかもしれない。
• -mtune=native: コンパイルしたパソコン(のCPU)向けのa.outを作ってくれという命令。生成されたa.outを他のパソコンにコピーすると、動くけどちょっと遅いかもしれない。

という認識。例えば上記のvpbroadcastq命令が動くパソコンは限られるため、-march=nativeを付けないと使用されない。

なお、x86 Options (Using the GNU Compiler Collection (GCC)) にあるように、-march=nativeは-mtune=nativeを含むため、両方指定する必要はない。

Specifying -march=cpu-type implies -mtune=cpu-type, except where noted otherwise.

-march=nativeって効果あるの？

GCC12からは次の理由で格段に効果が上がっている(GCC11までは-O3と組み合わせないと効果が薄い)。

• -march=nativeで解禁される命令の大半は自動ベクトル化 (自動ベクトル化について: gcc での自動ベクトル化 Wiki - yukicoder )向けの命令である。
• GCC11までは-O3を指定しないと自動ベクトル化がonにならなかったが、GCC12からは-O2でonになる。ただし-O3より自動ベクトル化のしきい値が高い(fvect-cost-model=cheap)。

おそらく最も効果があるものの一つはbitsetのand/or/xor/count/any/all等なので、 ABC 329 F で不正を試みる。自明な $O(NQ / w)$ 解をMLE対策に少し工夫して投げると、次のように

• (C++20, -march=native): 1352ms / 4s AC https://atcoder.jp/contests/abc329/submissions/48871083
• (C++17, -march=nativeなし): TLE(時々AC) https://atcoder.jp/contests/abc329/submissions/48871094

約3倍の高速化が確認できる。ただし、C++17でもpragmaをモリモリと付けるとACする

• (C++17, pragmaモリモリ): 1767ms / 4s AC https://atcoder.jp/contests/abc329/submissions/48871873

さらに GCC optimize("Ofast")を付けると、C++20より速くなる。

• (C++17, pragma, Ofast) 971ms https://atcoder.jp/contests/abc329/submissions/48871954

なんか実行時間がめちゃくちゃブレる(インスタンスガチャ？)ので、ブレの範疇な気もする　まったく同じコードを2回投げて 1359ms vs 1907ms とか出た ( https://atcoder.jp/contests/abc329/submissions/48871947, https://atcoder.jp/contests/abc329/submissions/48871969 )

(不正以外で)まともに効果がありそうなのは DP 系だろうか、modintをMontgomery乗算で実装するとSIMD(自動ベクトル化)と相性がいいという小ネタもあり、云々

-march=nativeでpragmaって代替きかないの？

大体聞きそうな気はする、懸念点は

• こだわるとジャッジごとに異なるpragmaを用意する必要がありそうで、ダルい
• -march=...とpragmaが本当に等価なのかわかってない(例えばyukicoderの記事には #pragma GCC optimize ("O3") は -O3 でのコンパイルとは異なるようです とある)
• GCC13.2だとなんかpragma使えないかも？ https://twitter.com/yosupot/status/1730356100363645093

あたりだろうか

結局 -march=native って危険なの？

もちろんまともな根拠はなくただの直観になるが、「わずかに危険だけど、問題になることはほぼない」程度ではないかと思っている

• Safe CFLAGS - Gentoo wiki Gentoo wiki曰く、-march=nativeはおススメである

A recommended default choice for CFLAGS or CXXFLAGS is to use -march=native

• また、march=nativeのように"コンパイルしたパソコン専用に最適化"という概念自体もかなり使われているはず。確かrustだとcargo installでデフォルトで-march=native相当のオプションがonになったはず

• そもそもGCCのバグに出会うのがレアイベント

  • 自分は-march=nativeの有無でどうこうというのは初めて出会った気がする、記憶力がないだけか？
  • AtCoderにGCC12と-march=nativeが導入される以前はノーカンという話もある。

• 一方で、pragmaで大体何とかなりそうだしわざわざ入れる必要がないのではという意見もありそう

皆さんの意見はどうでしょうか(ブン投げ)