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について色々調べた。

話題の発端

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

シンプルなコードで、しかも-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 および vmovdqab{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対策に少し工夫して投げると、次のように

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

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

なんか実行時間がめちゃくちゃブレる(インスタンスガチャ?)ので、ブレの範疇な気もする まったく同じコードを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 って危険なの?

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

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で大体何とかなりそうだしわざわざ入れる必要がないのではという意見もありそう

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

FHC 2023 Final 反省会会場 / 並列化研究

概要 / 言い訳タイム

FHC 2023 Finalの順位表を見ると、私がBをダウンロードだけして提出していないことがわかると思います。 そもそもカクタス(F)をシバけないとどうしようもないセットではあったのですが、それでもBを出していれば一応5位で入賞であり、このムーブは奇妙です。

これ

実際何が起きたのかというと、普通にTLEしました。AC率からもわかるように最悪ケースを作るのが難しい問題ではないのですが、$N \times M \le 1{,}000{,}000$に対して $N = M = 300$がほぼ最悪ケースだと勘違いしました。このコンテストは世界top25のコンテストです。

反省

そもそもこのミスはリカバリー可能なはずでした。FHCは手元実行なので、適当にケースごとに並列化すれば容易に高速化できるはずです。この準備はしようしようと思っていたのですが、面倒でやらなかったらやられてしまいました。

というわけで、この記事はFHC用に並列化環境を整備する話になります。

成果物

成果物をFHC 2023 Qual A1問題に適応したのがこちらになります。

qual-A1.cpp · GitHub

非ライブラリ部分は次のようになります。

using namespace std;

void solve(auto input_end, auto output) {
    int s, d, k;
    cin >> s >> d >> k;
    input_end();

    int buns = 2 * (s + d);
    int patties = s + 2 * d;

    output([&] {
        if (buns < k + 1 || patties < k) {
            cout << "NO" << endl;
        } else {
            cout << "YES" << endl;
        }
        cout << flush;
    });
};

int main() {
    int t;
    cin >> t;
    fhc_solve([&](auto i, auto o){ solve(i, o); }, t);
    return 0;
}

可能な限り、「いつものように普通にsolve関数を書いたら、勝手に並列化される」に近いものを目指しました。

  • input_end()を入力が終わった後に呼ばないといけない
  • output()にラムダを渡して、そこですべての出力を行わないといけない
  • outputの最後で必ずflushしないといけない
  • multi-threadを使っているので、グローバル変数を書き換えたりすると、何が起きるかわからない
  • 実行のたびに出力をファイルに全部保存して消してないので、出力が大きいとやばそう

などが残った制約です。代償としてライブラリ側はfreopenなどを乱用したコードになりました。

実行すると次のようになります

$ ./A1/main < A1/big.in > A1/big.out                    
Start FHC solver: tmp = "/tmp/output-14960984700273349145", parallel = 12
[#0] Start case: 1 / 79
[#0] End case: 1 (0 ms)
[#1] Start case: 2 / 79
[#1] End case: 2 (0 ms)
[#0] Start case: 3 / 79
[#0] End case: 3 (0 ms)
[#2] Start case: 4 / 79
[#2] End case: 4 (0 ms)
[#3] Start case: 5 / 79
[#10] Start case: 6 / 79
[#3] End case: 5 (0 ms)
[#1] Start case: 7 / 79
[#10] End case: 6 (0 ms)
[#5] Start case: 8 / 79
[#1] End case: 7 (0 ms)
:

私のPCは(論理)12コアなので、12並列でケースが実行されます。

Bに再挑戦

実際にBに再挑戦してみました、80s -> 28sなので、おおよそ3倍弱の高速化のようです。80sってそもそも間に合ってね?については、コンテスト中の6分間であわてて定数倍高速化した後のコードだからです(本当は一番最初のコードで試したかったのですが、上書きしていたため入手できませんでした…)。

旧
./B/main_single < B/test.in > B/test3.out  80.11s user 0.03s system 99% cpu 1:20.14 total

新
./B/main < B/test.in > B/test3.out  148.96s user 0.15s system 527% cpu 28.245 total

また、ログは次のような感じです。最大ケースのCase 31(N = M = 1000)に引っ張られていることがわかります。

:
[#4] End case: 98 (13 ms)
[#4] Start case: 99 / 100
[#4] End case: 99 (5 ms)
[#4] Start case: 100 / 100
[#4] End case: 100 (14 ms)
[#8] End case: 45 (1874 ms)
[#7] End case: 37 (11847 ms)
[#3] End case: 38 (11849 ms)
[#2] End case: 36 (11859 ms)
[#0] End case: 39 (11978 ms)
[#10] End case: 35 (12436 ms)
[#6] End case: 41 (12013 ms)
[#9] End case: 40 (12544 ms)
[#11] End case: 42 (12237 ms)
[#5] End case: 43 (11407 ms)
[#1] End case: 31 (28236 ms)
./B/main < B/test.in > B/test3.out  148.96s user 0.15s system 527% cpu 28.245 total

Eにも挑戦

また、結構実行時間がやばかったEに対してもやってみたところ… MLEしました。

並列なしでメモリを4GBぐらい使うとんでもプログラムなので、12並列ならばさもありなんです。

3並列ならば、高速化が確認できました(155s -> 70s)。しかし、いざ本番でMLEで突然死、は困るので、この弱点の対応法は悩みどころです。48GB余裕なぐらいメモリ増設すればいいだけでは? お安い対策としては、そもそもクラウドに巨大インタンスを借りてそこでやるなどが考えられます。ルール的にOKなのかは微妙ですが…

./E/main_single < ./E/big.in > ./E/big2.out  139.51s user 15.66s system 99% cpu 2:35.24 total
./E/main < ./E/big.in > ./E/big2.out  183.57s user 5.30s system 269% cpu 1:10.12 total

ログを確認したところ、一番時間のかかるケースは10s程度だったので、クラウド等にガチ強力インスタンス借りれば10s切れそうではある

[#0] End case: 15 (10019 ms)

遅延Segtree3

大嘘昔話

実は「segtreeというのはモノイドを載せられて、lazysegtreeはそれにいい感じの作用が行えて…」のようにsegtreeが抽象化されたのは割と最近です。 昔はなんかsegtreeって大体実装一緒だな…と思いながら、みな自分のstarryskytree.cppを毎回コピペして適当に書き直して使っていました。

もちろん適切にクラス等を使って最強のsegtreeを作れば勝ちまくりモテまくりであることには皆薄々気づいており、私もそのような夢を追い求める若者の一人でした。 その名残がこちらです。

2023年

時は2023年、昔はC++11の機能は新しいといわれていましたが、今ではC++20がどこのジャッジでも使えます(正確には使えないジャッジは引退しました)。C++20と言えばconcept、今回は昔を思いながら、conceptの勉強がてら抽象化segtreeに挑戦してみました。

上の遅延SegTree / 遅延Segtree2にあるように、大体実装方法は

  • structを自分で定義してそれをsegtreeに渡す
  • lambda/関数 を演算の個数だけ用意して一気にsegtreeに渡す

の2種類だと思うんですが、今回は両方できるようにしてみました。2018年ならいざ知らず同様の実装がそこら中にあると思う。

コード

こちらです。

#include <vector>
#include <iostream>
#include <numeric>

template <class T>
concept monoid = requires (T& x, typename T::S s) {
    { x.op(s, s) } -> std::same_as<typename T::S>;
    { x.e() } -> std::same_as<typename T::S>;
};

template <monoid M>
struct SegTree {
    using S = M::S;

    M m;
    std::vector<S> v;

    SegTree(M _m, std::vector<S> _v) : m(_m), v(_v) {
    }

    S all_prod() {
        // TODO optimize :)
        S val = m.e();
        for (auto x : v) {
            val = m.op(val, x);
        }
        return val;
    }
};

template <class T, class OP, class E>
struct LambdaMonoid {
    using S = T;
    S op(S a, S b) { return _op(a, b); }
    S e() { return _e(); }

    LambdaMonoid(OP op, E e) : _op(op), _e(e) {}
  private:
    OP _op;
    E _e;
};
template <class OP, class E>
LambdaMonoid(OP op2, E e2)->LambdaMonoid<decltype(e2()), OP, E>;

// --- ここまでライブラリ ---

struct AddInt {
    using S = int;
    S op(S a, S b) { return a + b; }
    S e() { return S(0); }
};

int main() {
    std::vector<int> v(10);
    std::iota(v.begin(), v.end(), 0);
    SegTree seg0(AddInt(), v);

    SegTree seg1(
        LambdaMonoid([&](int a, int b) { return a + b; }, [&]() { return 0; }),
        v);

    std::cout << "sum(1..10) = " << seg0.all_prod() << std::endl;
    std::cout << "sum(1..10) = " << seg1.all_prod() << std::endl;

    return 0;
}

タイトルに遅延segtreeとありますが、シンプルに詐欺です。普通のsegtreeしか試してみていません。

解説

まず、最初の数行がいきなり重要です

template <class T>
concept monoid = requires (T& x, typename T::S s) {
    { x.op(s, s) } -> std::same_as<typename T::S>;
    { x.e() } -> std::same_as<typename T::S>;
};

template <monoid M>
struct SegTree {
    using S = M::S;

    M m;
    :

これは、monoid conceptを定義し、struct SegTreeがこのmonoid conceptを満たすMしか受け取れないようにしています。Mmonoid conceptを満たすとは、

  • using S = hoge として値の型が定義されている
  • op(S, S) -> Sをメンバとして持つ
  • e() -> Sをメンバとして持つ

という、大体atcoder libraryと同じ定義です。例えば下のほうにあるstruct AddIntmonoid conceptを満たします。

struct AddInt {
    using S = int;
    S op(S a, S b) { return a + b; }
    S e() { return S(0); }
};

なので、こういうAddInt構造体を用意して、main関数の最初で行われているように

    std::vector<int> v(10);
    std::iota(v.begin(), v.end(), 0);
    SegTree seg0(AddInt(), v);

と書けば[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]の格納されたsegtreeが出来ます。

また、わざわざstruct AddIntのようにstructを定義しなくても、ラムダ式を書き並べるだけでsegtreeを使うこともできます。

    SegTree seg1(
        LambdaMonoid([&](int a, int b) { return a + b; }, [&]() { return 0; }),
        v);

これがどういう仕組みかというと、ラムダ式を受け取ってmonoid structのようにふるまうLambdaMonoidもライブラリ側に用意しておきました。

template <class T, class OP, class E>
struct LambdaMonoid {
    using S = T;
    S op(S a, S b) { return _op(a, b); }
    S e() { return _e(); }

    LambdaMonoid(OP op, E e) : _op(op), _e(e) {}
  private:
    OP _op;
    E _e;
};
template <class OP, class E>
LambdaMonoid(OP op2, E e2)->LambdaMonoid<decltype(e2()), OP, E>;

応用編として、「pair<S, S>を使うことで一般のモノイドをReverse可能なモノイドに変換するやつ」とかが実現できると思っています。なんのこっちゃという話ですが、平衡二分木でこういうのが欲しくなります。

template <monoid M> struct AttachReverse {
    using S = std::pair<M::S, M::S>;
    S op(S a, S b) { return {m.op(a, b), m.op(b, a)}; }
    S e() { return {m.e(), m.e()}; }

    S rev(S a) { return {a.second, a.first}; }

    AttachReverse(M _m) : m(_m) {}

  private:
    M m;
};

結論

全てが懐古

Weighted balanced binary tree の平衡条件をWolfram Engineで

注: ガチ誰得記事

Weighted balanced binary treeという平衡二分探索木があります。詳細はwikipedia(Weight-balanced tree - Wikipedia)が詳しいのですが、ざっくり言うと

  • 「子のサイズが親のサイズの少なくとも $\alpha$ 倍」
  • もしくは同値な表現として「左右の子のサイズが高々 $\frac{1}{\alpha} - 1$ 倍しか違わない」

という平衡条件の木です。ここで、サイズというのは葉にのみ値を載せる木ならば葉の個数、全部の頂点に値を載せる木ならば(頂点の個数) + 1とします。

$\alpha$ を大きくすればするほど回転が増える代わりに強く平衡(=木の高さが低くなる)し、$\alpha$ を小さくすればするほど回転が減る代わりに平衡が弱くなります。

このサイズ情報はまず間違いなく競技プログラミングだと管理する(K番目のアクセスに必要なため)ので、赤黒木等と違い追加メモリが実質0と言う利点があります。カロリー0!

この木は強そうなのですが、うまく動く証明が難しいです。葉にのみ値を載せる木の場合のmerge関数はたとえば次のようになります。

Node* merge(Node* l, Node* r) {
    if (is_balanced(l->size, r->size)) return make_node(l, r);
    if (l->size > r->size) {
        auto x = l->l;
        auto y = merge(l->r, r);
        if (is_balanced(x->size, y->size)) return make_node(x, y);
        if (is_balanced(x->size, y->l->size) && is_balanced(x->size + y->l->size, y->r->size)) return make_node(make_node(x, y->l), y->r);
        return make_node(make_node(x, y->l->l), make_node(y->l->r, y->r));
    } else {
        ...
    }
}

このコードはどういうことかというと、lとrをマージしたい時、

  • lとrをそのままマージして問題ない: マージする
  • lがrに対してデカすぎてマージできない場合は、l->rとrを再帰的にmergeする。そしてl->lをx、merge(l->r, r)をyとする。
  • a = x, b = y->l->l, c = y->l->r, d = y->rとして、下図の3つの木のうち、少なくともどれか一つはうまくバランスしている。

という、大変魔法のようなマージ関数になります。

もちろん、任意の $\alpha$ についてこのmerge関数が動くわけではなく、wikipediaにあるように $0 \lt \alpha \lt 1 - \frac{1}{\sqrt{2}}$ ならば動きます。

今日はこの事実をWolfram Engine (Wolfram Engine) を利用して確かめていきましょう。

  • 変数 t を $\frac{1}{\alpha} - 1$ とします、つまり左右の子は高々 t 倍しかサイズが違わないという意味です。
  • 変数 a, b, c, d を上図の a, b, c, d のサイズとします
  • 変数 er のサイズとします。つまり l のサイズは a + b + c + d - e です。

まず、次のことが言えます

  • lr より t 倍以上大きい: (a+b+c+d-e) > t e
  • 元々 l->ll->rはバランスしていた: a <= t(b+c+d-e) && t a >= (b+c+d-e)
  • bcはバランスしている: b <= t c && t b >= c
  • b + cdはバランスしている: (b+c) <= t d && t(b+c) >= d

そしてこの条件下で

  • 図の一番上のマージが失敗する = (ab+c+dがバランスしていない): (t a < b+c+d || a > t(b+c+d))
  • 二番目のマージが失敗する = (ab+cがバランスしていない OR a+b+cdがバランスしていない): (t a < b+c || a > t(b+c) || t(a+b+c) < d || a+b+c > t d)
  • 三番目のマージが失敗する = (ab / cd / a+bc+d のどれかがバランスしていない): (t a < b || a > t b || t(a+b) < c+d || a+b > t(c+d) || t c < d || c > t d )

これを全てまとめてwolfram engineに投げると次のようになります

In[1]:= Resolve[Exists[{a, b, c, d, e}, 1 <= t && 0 <= a && 0 <= b && 0 <= c && 0 <= d && 0 <= e && (a+b+c+d-e) > t e && a <= t(b+c+d-e) && t a >= (b+c+d-e
) && b <= t c && t b >= c && (b+c) <= t d && t(b+c) >= d && (t a < b+c+d || a > t(b+c+d)) && (t a < b+c || a > t(b+c) || t(a+b+c) < d || a+b+c > t d) && (t
 a < b || a > t b || t(a+b) < c+d || a+b > t(c+d) || t c < d || c > t d )], t]

Out[1]= 1 <= t < 1 + Sqrt[2]

というわけで、$1 + \sqrt{2} \le t$、つまり$\alpha \le \frac{1}{2 + \sqrt{2}} = 1 - \frac{1}{\sqrt{2}}$ならば正しく動くことがわかりました。

結論: wolfram engineってすごい

データ構造をマージする一般的なテクを高速化するやや一般的なテク

知る人ぞ知る謎の知識みたいになってる気がしたのでメモ。

サイズ $A$ とサイズ $B$ のデータ構造 ($A < B$) が $O(A \log B)$ 時間でマージできる場合、合計 $O(N \log ^2 N)$ 時間で $N$ 個の要素がマージできます。競プロだと std::setstd::priority_queue などです。

ここで、もしもマージが $O(A (\log B - \log A + 1))$ に高速化できたならば、合計の計算量が $O(N \log N)$ になります。証明は普通のデータ構造をマージする一般的なテクとほぼ同じで、ある要素を含むデータ構造のサイズが $1 \to x \to y \to N$ と成長したなら、この要素にかかる計算量は $(\log x - \log 1 + 1) + (\log y - \log x + 1) + (\log N - \log z + 1) = O(\log N)$ 、といったノリです。

でも$O(A (\log B - \log A + 1))$ なんて謎の計算量なんて…と思うかもしれませんが、実はかなり多くのデータ構造がこの計算量を達成できます。すごいですね。

平衡二分木の merge / split が登場してややこしいので、データ構造をマージする操作を meld と呼ぶことにします

二分ヒープ

std::priority_queue の中身です。(自分で std::priority_queue 相当のものを実装する必要がありますが)実はpriority_queue + マージテクは $\log$ 一個です。でも $O(N \log ^2 N)$ のマージがそもそも爆速だからうまみが少なそう。

まず、二分ヒープを線形時間で構築するアルゴリズムを理解する必要があります。たとえば下のコードのようになります。

void down_heap(vector<int>& d, int u) {
    int n = int(d.size());
    while (2 * u + 1 < n) {
        int v = (2 * u + 1);
        if (v + 1 < n && h.d[v] < h.d[v + 1]) v++;
        if (h.d[u] >= h.d[v]) break;
        swap(h.d[u], h.d[v]);
        u = v;
    }
}

void build_heap(vector<int>& d) {
    int n = int(d.size());
    for (int i = n - 1; i >= 0; i--) {
        down_heap(d, i);
    }
    return h;
}

「ある頂点について、葉の方向に自分より大きい要素とswapしていく」down_heap関数を、葉から順に全ての頂点に呼ぶと二分ヒープが構築できます。計算量は各頂点の「葉までの距離」の総和になり、式をコネコネすると線形時間であることがわかります。

meld関数のコードは次のようになります。

void meld(vector<int>& h, vector<int>& other) {
    if (h.len() < other.len()) swap(h, other);
    if (other.empty()) return;

    int l = int(h.size()), r = l + int(other.size()) - 1;

    h.insert(h.end(), other.begin(), other.end());

    while (l) {
        l = (l - 1) / 2;
        r = (r - 1) / 2;
        for (int i = r; i >= l; i--) {
            down_heap(h, i);
        }
    }
}

これは何をしているのかというと、小さいほうのヒープの要素を適当に大きいヒープの末尾に追加した後、「新しく追加した要素を子孫に持つ要素」全てに対して先述のdown_heap関数を呼んでいるだけです。 計算量ですが、各頂点についてmin(other.size(), 葉までの距離)であること、そしてこれの総和が上記の $O(A (\log B - \log A + 1))$ になることが示せます。

std::set (by merge-split baseの平衡二分木)

std::setも、自作の平衡二分木で実装することで $O(A (\log B - \log A + 1))$が達成できます。まず、一般に merge-split base の平衡二分木ならなんでも可能な(ハズの)方法を紹介しますが、おそらく定数倍はとても悪いです。

$A \le \sqrt{B}$ の場合は愚直に $A$ 回 insertすればよいので、$\sqrt{B} < A$ としてよい。

  • サイズ $B$ の木を サイズ $B / A$ (ぐらい)の木に $A$ 分割する。
    • 愚直にsplitを使うと $O(A \log B)$ だが、例えば8分割したいなら (2分割 => それぞれを2分割 => それぞれを2分割)、のように、なるべくsplitする木のサイズが小さくなるように切っていくと $O(A (\log B - \log A + 1))$ になる
  • 分割した木それぞれに対して、サイズ $A$ の木の(対応する値の範囲の)要素を追加していく
    • 追加する要素が $B / A$ 個以下なら愚直にinsert
    • 追加する要素が $B / A$ 個以上なら(両方の木をvectorにする) => (std::merge) => (vectorを長さ $B / A$ のブロックに分割する) => (それぞれのブロックから、線形時間で新しい木を構築する)で、線形時間でマージする
  • 分割した木たちを(splitと同様に、いい感じに)マージして新しい木を作る。
    • サイズが $(B / A)$以上$2(B / A)$以下の木たちのマージになり、splitと同様の計算量解析が可能

で$O(A (\log B - \log A + 1))$ になっている…ハズ

std::set (by merge-split baseの多くの平衡二分木)

葉にだけ値を持たせる赤黒木など、サイズ $L, R (L < R)$の木のmerge (not meld)が $O(\log R - \log L + 1)$であることが保証可能な平衡二分木は割と簡単に実装できます。(参考: https://www2.ioi-jp.org/camp/2012/2012-sp-tasks/2012-sp-day4-copypaste-slides.pdf) このような平衡二分木なら、次のような割と素朴なmeld関数が上手くいくはずです。

Node* meld(Node* n, deque<int>& q) {
    while(q.size() && q.front() == n->val) q.pop_front(); // multisetにしたい場合はこの行を消す
    if (n->val < q.front()) return n;
    if (is_leaf(n)) {
        vector<int> v;
        while (q.size() && q.front() < n->val) {
            v.push_back(q.front()); q.pop_front();
        }
        return merge(build_tree(v), n); // build_tree(v): |v|時間で木を構築する関数
    }
    Node* l = meld(n->l, q);
    Node* r = meld(n->r, q);
    return merge(l, r); // (1)
}

Node* meld(Node* n, Node* m) {
    if (n->size < m->size) swap(n, m);
    deque<int> q = to_deque(m); // to_deque(m): mの要素を舐めて(sortedな)dequeを線形時間で生成する関数
    n = meld(n, q);
    return merge(n, build_tree(q));
}

まず、meld関数が呼ばれる回数は、n->sizeが $B / A$以上かどうかで場合分けすれば示せます。

  • n->size が $B / A$ 以上 : そもそもこういう頂点は $O(A)$ 個しかない
  • n->size が $B / A$ 未満: 高さが $O(\log (B / A)) = O(\log B - \log A + 1)$ なので合計 $O(A (\log B - \log A + 1))$

コード中の(1)の部分の計算量が本質です。ここで、merge関数の計算量により、meld前後で木のサイズが $k$ 倍に増えていたら (1) の部分の計算量が $O(\log k)$ であることを利用します。

  • 木のサイズが $2$ 倍以上に増えるような頂点数は?: 高々 $O(A)$
  • 木のサイズが $4$ 倍以上に増えるような頂点数は?: 高々 $O(A / 2)$
  • :

より、(*)の部分の計算量の合計は(meldが呼ばれた回数に加えて)高々 $O(A)$ です。

std::set (by splay-tree)

は愚直に小さいほうの集合を昇順(or 降順)にinsertしていくだけで $O(A (\log B - \log A + 1)$ って噂を聞いたんですが、本当かわかりません。

std::set (by treap)

は割と容易に $O(A (\log B - \log A + 1)$ のmeldが書けるって聞きました。いかがでしたか