Learning Algorithms

アルゴリズムの勉強メモ

二重辺連結成分分解

二重辺連結成分分解の実装です。

二重頂点連結成分分解についてはこちらの記事をご覧ください。

無向グラフ上で、橋が存在しない部分グラフを一つの頂点につぶすと、橋を辺とする木ができます。以下は、それをそのまま実装しています。

struct LowLink {
        set<pair<int, int>> bridge;
        vector<int> articulation, ord, low;
        vector<bool> used;
        vector<vector<int>> g;
        int n, k = 0;
        LowLink(const vector<vector<int>> &g) : g(g) {
                n = g.size();
                ord.resize(n, 0);
                low.resize(n, 0);
                used.resize(n, false);
        }
        void dfs(int u, int prev) {
                used[u] = true;
                ord[u] = k ++;
                low[u] = ord[u];
                bool is_articulation = false;
                int cnt = 0;
                for (auto v : g[u]) if (v != prev) {
                        if (!used[v]) {
                                cnt ++;
                                dfs(v, u);
                                low[u] = min(low[u], low[v]);
                                if (low[v] > ord[u]) {
                                        bridge.emplace(min(u, v), max(u, v));
                                } 
                                if (prev != -1 && low[v] >= ord[u]) {
                                        is_articulation = true;
                                }
                        } else {
                                low[u] = min(low[u], ord[v]);
                        }
                }
                if (prev == -1 && cnt > 1) is_articulation = true;
                if (is_articulation) articulation.push_back(u);
        }
};
struct TwoEdgeConnectedComponent {
        int n;
        vector<vector<int>> g, tree;
        vector<int> cmp;
        TwoEdgeConnectedComponent(const vector<vector<int>> &g) : g(g) {
                n = (int) g.size();
                cmp.assign(n, -1);
        }
        void build() {
                LowLink lnk(g);
                lnk.dfs(0, -1);
                int k = 0;
                function<void (int, int)> dfs = [&](int u, int prev) {
                        cmp[u] = k;
                        for (auto v : g[u]) if (cmp[v] == -1 && lnk.bridge.count({min(u, v), max(u, v)}) == 0) {
                                dfs(v, u);
                        }
                };
                for (int i = 0; i < n; i ++) if (cmp[i] == -1) {
                        dfs(i, -1);
                        k ++;
                }
                tree.resize(k);
                for (auto e : lnk.bridge) {
                        tree[cmp[e.first]].push_back(cmp[e.second]);
                        tree[cmp[e.second]].push_back(cmp[e.first]);
                }
        }
};

以下のような問題が簡潔に実装できます。
D - 旅行会社高橋君

int main() {
        int n, m;
        scanf("%d%d", &n, &m);
        vector<vector<int>> g(n);
        for (int i = 0; i < m; i ++) {
                int a, b;
                scanf("%d%d", &a, &b);
                a --, b --;
                g[a].push_back(b);
                g[b].push_back(a);
        }
        TwoEdgeConnectedComponent tecc(g);
        tecc.build();
        LCA lca(0, tecc.tree);
        int q;
        scanf("%d", &q);
        while (q --) {
                int a, b, c;
                scanf("%d%d%d", &a, &b, &c);
                a --, b --, c --;
                if (lca.dist(tecc.cmp[a], tecc.cmp[b]) + lca.dist(tecc.cmp[b], tecc.cmp[c]) == lca.dist(tecc.cmp[a], tecc.cmp[c])) {
                        printf("OK\n");
                } else {
                        printf("NG\n");
                }
        }
        return 0;
}

Atcoder Petrozavodsk Contest 001 H. Generalized Insertion Sort

H - Generalized Insertion Sort

解法

まず制約に注目すると、$25000 \geq n \log n$ にしか見えないので、そのような計算回数を実現するアルゴリズムを考えたい気持ちになります。

木に関するアルゴリズムであって、計算量に $\log$ が現れるものは僕は現時点ではそんなに知らなくて、$LCA$ や $HL$ 分解などでした。すると、自然に $HL$ 分解のときに使うパスに近いものを考える発想に至ります。

木の頂点であって、その頂点の部分木がパスであるようなものをパス頂点と呼ぶことにします。パス頂点であって、同一のパスを共有する頂点同士はまとめてソートできないでしょうか。

もしこれができれば、それ以降その木のパス頂点はいじる必要がなくなるので、破壊してしまってよさそうです。パス頂点になる条件を考えるとこの破壊の回数は $\log n$ 回で終了します。

さて、まず明らかにどの操作も根が関わることになるので、根の値に注目します。根の値がいずれかのパス頂点の $index$ と一致するならば、それをそのパスに移動させることができます。ここで、パス全体をまとめてソートするために、最終的な値の列と相対的な順序関係が同じになるように丁寧にソートをしていきます。つまり、パスの下の方から見ていって、挿入しようとする値と、今見ている値を比較して適切な位置に挿入します。もしどの $index$ とも一致しないならば、適当に根にもってきたい値を探してその値の位置に適当に挿入します。

計算量自体は $n$ が小さいことから、すべて適当に書くだけでも通ります。

実装
#include <cstdio>
#include <vector>
#include <algorithm>
#include <functional>
#include <map>
#include <set>
#include <string>
#include <iostream>
#include <cassert>
#include <cmath>
using namespace std;

int main() {
        int n;
        scanf("%d", &n);
        vector<vector<int>> g(n);
        vector<int> par(n);
        for (int i = 1; i < n; i ++) {
                int p;
                scanf("%d", &p);
                g[p].push_back(i);
                par[i] = p;
        }
        vector<int> a(n);
        for (int i = 0; i < n; i ++) {
                scanf("%d", &a[i]);
        }
        vector<int> color(n, 0); //2:need to complete, 1:determined, 0:no considering, -1:dead
        vector<int> ans;
        for (;;) {
                if (color[0] == -1) break;
                vector<bool> path(n);
                function<void (int)> check_path = [&](int u) {
                        path[u] = true;
                        int cnt = 0;
                        for (auto v : g[u]) if (color[v] != -1) {
                                cnt ++;
                                check_path(v);
                                path[u] = path[u] && path[v];
                        }
                        if (cnt > 1) path[u] = false;
                };
                check_path(0);
                set<int> paths;
                for (int i = 0; i < n; i ++) {
                        if (path[i]) {
                                paths.insert(i);
                        }
                }
                for (int i = 0; i < n; i ++) {
                        if (paths.count(a[i])) {
                                color[i] = 2;
                        }
                }
                function<int (int)> find_bottom = [&](int u) {
                        assert(path[u]);
                        for (auto v : g[u]) if (color[v] != -1) {
                                return find_bottom(v);
                        }
                        return u;
                };
                function<void (int, int, int)> insert = [&](int u, int val, int col) {
                        int tmpa = a[u];
                        int tmpc = color[u];
                        a[u] = val;
                        color[u] = col;
                        if (u != 0) insert(par[u], tmpa, tmpc);
                };
                function<bool (int, int, int)> compare = [&](int pos, int in_num, int cur_num) {
                        while (true) {
                                if (pos == in_num) {
                                        return false;
                                } else if (pos == cur_num) {
                                        return true;
                                }
                                if (pos == 0) {
                                        break;
                                }
                                pos = par[pos];
                        }
                        return false;
                };
                vector<int> depth(n);
                function<void (int, int)> get_depth = [&](int u, int d) {
                        depth[u] = d;
                        for (auto v : g[u]) if (color[v] != -1) {
                                get_depth(v, d + 1);
                        }
                };
                while (true) {
                        bool ok = true;
                        for (int i = 0; i < n; i ++) {
                                ok = ok && a[i] == i;
                        }
                        if (ok) { 
                                break;
                        }
                        if (color[0] == 2) {
                                int target = a[0];
                                int cur = find_bottom(target);
                                while (color[cur] == 1 && compare(cur, target, a[cur])) {
                                        cur = par[cur];
                                }
                                insert(cur, target, 1);
                                ans.push_back(cur);
                        } else {
                                get_depth(0, 0);
                                int idx = -1;
                                int ma = -1;
                                for (int i = 0; i < n; i ++) {
                                        if (ma < depth[i] && color[i] == 2) {
                                                ma = depth[i];
                                                idx = i;
                                        }
                                }
                                if (idx == -1) {
                                        break;
                                }
                                insert(idx, a[0], color[0]);
                                ans.push_back(idx);
                        }
                }
                for (int i = 0; i < n; i ++) {
                        if (path[i]) {
                                color[i] = -1;
                        }
                }
        }
        printf("%d\n", (int) ans.size());
        for (int i = 0; i < (int) ans.size(); i ++) {
                printf("%d\n", ans[i]);
        }
        return 0;
}

CS Academy 069 D. Cover the Tree

CS Academy

解法

木を最小個数のパスで被覆する問題。選ぶパスは重複して良いことに注意します(もし重複を許さないなら、それはオイラー閉路を適切に構築してそれを出力する問題です)。

とりあえず、葉は必ずパスの端点として選ぶ必要がありそうです。よって(葉の数 $+ 1$) / $2$ 個のパスは必ず必要です。逆に、ある葉ではない頂点を根として考えると、その根を通過するようなパスばかりを選んでいけるならば、この個数で被覆可能です。というわけで、そのような頂点を求めます。

それは明らかに「葉に関する重心」なので(未証明)それを $C$ とします。これは重心を求めるアルゴリズムを葉の数に書き換えると求められます。

葉の重心を選ぶことで次のようにして葉のマッチングさせることができます。まず $C$ から伸びる部分木の中から葉の数が大きい順に $2$ つ選んでその中から適当に $1$ つずつ葉を選んでマッチングをさせる、というプロセスを繰り返します。これは優先度付きキューで効率よく処理できます。

「葉に関するの重心」が、葉になってしまう場合(頂点数 $2$ の木など)の場合分けが面倒なので、頂点倍加のテクニックを使いました。頂点倍加のテクニックとは、与えられた木に対して、すべての辺の途中に新しい頂点を追加することで、(一般的に言う)重心が必ず $1$ 個になるようにするテクニックです。この問題では、これによって葉以外の重心が必ず存在するようにできます。以上でこの問題が解けました。

実装
#include <cstdio>
#include <vector>
#include <algorithm>
#include <functional>
#include <cassert>
#include <queue>
using namespace std;

int main() {
        int n;
        scanf("%d", &n);
        vector<vector<int>> g(n + n - 1);
        for (int i = 0; i < n - 1; i ++) {
                int a, b;
                scanf("%d %d", &a, &b);
                a --, b --;
                g[a].push_back(n + i);
                g[n + i].push_back(a);
                g[b].push_back(n + i);
                g[n + i].push_back(b);
        }
        n += n - 1; //attention
        vector<int> centroids;
        vector<int> leaves(n);
        int leaves_cnt = 0;
        for (int i = 0; i < n; i ++) {
                if ((int) g[i].size() == 1) {
                        leaves_cnt ++;
                }
        }
        function<void (int, int)> dfs = [&](int u, int prev) {
                leaves[u] = 0;
                bool is_centroid = true;
                for (auto v : g[u]) if (v != prev) {
                        dfs(v, u);
                        leaves[u] += leaves[v];
                        if (leaves[v] > leaves_cnt / 2) is_centroid = false;
                }
                if ((int) g[u].size() == 1) leaves[u] ++;
                if (leaves_cnt - leaves[u] > leaves_cnt / 2) is_centroid = false;
                if ((int) g[u].size() == 1) is_centroid = false;
                if (is_centroid) centroids.push_back(u);
        };
        dfs(0, -1);
        int c = centroids[0];
        vector<vector<int>> leaves_set;
        vector<int> tmp;
        function<void (int, int)> dfs2 = [&](int u, int prev) {
                if ((int) g[u].size() == 1) tmp.push_back(u);
                for (auto v : g[u]) if (v != prev && v != c) {
                        dfs2(v, u);
                }
        };
        for (auto u : g[c]) {
                tmp.clear();
                dfs2(u, -1);
                leaves_set.push_back(tmp);
        }
        swap(leaves_set[0], leaves_set[1]);
        sort(leaves_set.begin(), leaves_set.end(), [&](const vector<int> &a, const vector<int> &b) {
                return (int) a.size() > (int) b.size(); 
        });
        int k = (int) leaves_set.size();
        vector<pair<int, int>> ans;
        priority_queue<pair<int, int>> que; //the number of leaves, index
        for (int i = 0; i < k; i ++) {
                que.push({ (int) leaves_set[i].size(), i });
        }
        while (que.size() >= 2) {
                auto p1 = que.top();
                que.pop();
                auto p2 = que.top();
                que.pop();
                int a = leaves_set[p1.second].back();
                int b = leaves_set[p2.second].back();
                ans.emplace_back(a, b);
                leaves_set[p1.second].pop_back();
                leaves_set[p2.second].pop_back();
                if (p1.first > 1) que.push({ p1.first - 1, p1.second });
                if (p2.first > 1) que.push({ p2.first - 1, p2.second });
        }        
        assert(que.size() <= 1);
        if (que.size() > 0) {
                auto p1 = que.top();
                que.pop();
                ans.emplace_back(leaves_set[p1.second].back(), c);
        }
        printf("%d\n", (int) ans.size());
        for (int i = 0; i < (int) ans.size(); i ++) {
                printf("%d %d\n", ans[i].first + 1, ans[i].second + 1);
        }
}

Atcoder Regular Contest 069 E. Frequency

E - Frequency

他の人の実装を見るとなんだかあっさり実装されているので、より簡単な解法がありそうです。

解法

数列に出現する数 $x$ について、数列全体に $x$ 以上の数が何個あって、その一番左の数の $index$ が何であるのかがわかれば答えがわかりそうです。

それを求めるのは難しいので、構成する数列は必ず広義単調減少になることを利用します。今、位置 $idx$ の数を構成に使っているとすると、次に構成に使う数というのは、区間 $[0, idx)$ の中で最大の値を持つものです。これはセグメント木で求めることができます。ただし、最大の値だけでなく、その $index$ (複数ある場合は一番左の)も返すようにします。

左側にある最大値よりも大きい数が右側にある場合は少し注意が必要で、構成に使う数を今の位置に保持しながら、右の数を処理してから左の数に移動するようにします。これで最適な数列が構成できました。

実装
#include <cstdio>
#include <vector>
#include <algorithm>
#include <functional>
#include <map>
#include <set>
#include <string>
#include <iostream>
#include <cassert>
#include <cmath>
using namespace std;

template<typename Type>
struct SegmentTree {
        vector<Type> data;
        int n;
        SegmentTree(int x) {
                n = pow(2, ceil(log2(x)));
                data.resize(2 * n - 1);
                for (int i = 0; i < 2 * n - 1; i ++) {
                        data[i] = make_pair(0, -1);
                }
        }
        Type merge(Type x, Type y) { //merge function (val, idx) maximum value and minimum idx
                if (x.first > y.first) {
                        return x;
                } else if (x.first < y.first) {
                        return y;
                } else {
                        return make_pair(x.first, min(x.second, y.second));
                }
        }
        void update(int k, Type p) { //update k-th value to p
                k += n - 1;
                data[k] = p;
                while (k > 0) {
                        k = (k - 1) / 2;
                        data[k] = merge(data[k * 2 + 1], data[k * 2 + 2]);
                }
        }
        // [l, r)
        Type query(int a, int b) { return query(a, b, 0, 0, n); }
        Type query(int a, int b, int k, int l, int r) {
                if (r <= a || b <= l) return make_pair(0, -1); //initial value
                if (a <= l && r <= b) return data[k];
                int m = (l + r) / 2;
                return merge(query(a, b, k * 2 + 1, l, m), query(a, b, k * 2 + 2, m, r));
        }
};

int main() {
        int n;
        scanf("%d", &n);
        vector<int> a(n);
        for (int i = 0; i < n; i ++) {
                scanf("%d", &a[i]);
        }
        SegmentTree<pair<int, int>> seg(n);
        for (int i = 0; i < n; i ++) {
                seg.update(i, make_pair(a[i], i));
        }
        map<int, int> cnt;
        for (int i = 0; i < n; i ++) {
                cnt[a[i]] ++;
        }
        sort(a.rbegin(), a.rend());
        a.erase(unique(a.begin(), a.end()), a.end());
        vector<long long> ans(n);
        int idx = n;
        int cur = 0;
        long long total = 0;
        for (;;) {
                pair<int, int> ma = seg.query(0, idx);
                total += cnt[a[cur]];
                if (cur + 1 < a.size()) ans[ma.second] += (long long) (a[cur] - a[cur + 1]) * total;
                else ans[ma.second] += (long long) a[cur] * total;
                cur ++;
                if (cur == a.size()) break;
                if (seg.query(0, ma.second).first == a[cur]) {
                        idx = ma.second;
                }
        }
        for (int i = 0; i < n; i ++) {
                printf("%lld\n", ans[i]);
        }
        return 0;
}

Atcoder Regular Contest 087 E. Prefix-free Game

E - Prefix-free Game

解法

まず文字列は $0$ と $1$ しか含まないので、完全二分木で考えました(この考え方は、一般的な文字列においても使えるらしく $Trie$木と呼ばれているらしいです)。

すると、木の頂点を順に塗りつぶしていく問題に変わります。ある頂点を選んだとき、その頂点を含む部分木全体と、その祖先すべてを塗ることになります。

よって各頂点を選んだときに、どのような部分木が残るのかに注意して $Grundy$ 数を求めればよさそうです。対称性に注意すると、現在の状況から遷移できる状態を考えて、高さ $k - 1$ の完全二分木の $Grundy$ 数 $g(k)$ は次のようになります。

$g(k) = mex\{0,\ g(k - 1),\ g(k - 1) \oplus g(k - 2),\ g(k - 1) \oplus g(k - 2) \oplus g(k - 3),\ ... \}$

このままだと効率よく求めることができませんが、実験をすると実は $k$ を割り切る最大の $2$ の冪であることがわかります(解説動画を見ると、$Grundy$ 数は、まず何も考えずに実験してくださいとおっしゃっていた)。

証明はよくわからないので、誰か教えてください。

$Trie$ 木を知らなくても普通に動的に木を作っていけばよいです。

ちなみによくあるテクニックですが、ある数 $x$ が $2$ で割り切れる回数は、__builtin_ctz(x)または__builtin_ctzll(x)で求められるので、今回の場合 $Grundy$ 数 $g(x)$ は1LL << __builtin_ctzll(x)です。

#include <cstdio>
#include <vector>
#include <algorithm>
#include <functional>
#include <map>
#include <set>
#include <string>
#include <iostream>
#include <cassert>
#include <cmath>
#include <memory>
using namespace std;

struct node_t {
        shared_ptr<node_t> left, right;
        node_t(shared_ptr<node_t> l, shared_ptr<node_t> r) : left(l), right(r) {}
};

using node = shared_ptr<node_t>;

node new_node() {
        return node(new node_t(NULL, NULL));
}

int main() {
        int n;
        long long l;
        scanf("%d %lld", &n, &l);
        vector<string> s(n);
        for (int i = 0; i < n; i ++) {
                cin >> s[i];
        }
        node root = new_node();
        function<void (node, const string&, int)> dfs = [&](node u, const string &str, int depth) {
                if (depth >= (int) str.size()) {
                        return;
                }
                if (str[depth] == '0') {
                        if (u->left == NULL) {
                                u->left = new_node();
                        }
                        dfs(u->left, str, depth + 1);
                } else {
                        if (u->right == NULL) {
                                u->right = new_node();
                        }
                        dfs(u->right, str, depth + 1);
                }
        };
        for (int i = 0; i < n; i ++) {
                dfs(root, s[i], 0);
        }
        function<long long (long long)> get_grundy = [&](long long x) {
                return 1LL << __builtin_ctzll(x);
        };
        long long ans = 0;
        function<void (node, int)> dfs2 = [&](node u, int depth) {
                if ((u->left == NULL && u->right != NULL) || (u->left != NULL && u->right == NULL)) {
                        ans ^= get_grundy(l - depth);
                }
                if (u->left != NULL) {
                        dfs2(u->left, depth + 1);
                }
                if (u->right != NULL) {
                        dfs2(u->right, depth + 1);
                }
        };
        dfs2(root, 0);
        puts(ans ? "Alice" : "Bob");
        return 0;
}

Atcoder Petrozavodsk Contest 001 F. XOR Tree

F - XOR Tree

解法

$XOR$ の性質を使ってかなりきれいに解けます。

まず、パスに関する操作が非常に扱いにくいため、なんとかしたい気持ちになります。

$XOR$ の性質で重要なものとして、任意の数について、それ自身との $XOR$ は$0$ になるというものがあります。これをうまく使えないでしょうか。

ここで、パスの自明な性質として、端の点の次数は $1$ で、途中の点の次数は $2$ であるという性質があります(それはそうで)。このことを考えると、通過するたびに $XOR$ の演算によって打ち消されていくような不変量のようなものを見つけることができそうです。

具体的には、次のようにします。すべての頂点について、その頂点 $v$ から伸びるすべての辺に書かれた数の $XOR$ を計算しておき、それを $a_v$ とおきます。すると、操作を行うパスの途中の点では、$x$ の $XOR$ を $2$ 回とることになり、打ち消しあうことになります。これで操作を、端点のみに任意の数 $x$ を $XOR$ として作用させる操作に置き換えることができました。この値 $a_i$ の値がすべて $0$ になったとき、さらにそのときに限り、題意をみたすことになります。

この置き換えのすごいところは、グラフが連結であるという前提から、もはやグラフの構造は全く気にしなくてよくなったことです。

さて、次にどの $2$ 点を選んでいくかについて考えていきます。この選んだ $2$ 点に辺を張ったグラフを考えてみると、実はこれは適切にやることで必ず森になることが示せます。閉路をつくるメリットがないことを述べればよいわけですが、それは次の図によって示すことができます。

f:id:KokiYamaguchi:20180205171312j:plain

$3$ 頂点に対して、画像左のように $3$ 本の辺を張るとします。この結果、それぞれの頂点に作用させる $XOR$ の値はそれぞれ $xy, yz, zx$ となっています。ところが、$XOR$ の性質より、画像右のような $2$ 本の辺で代用することができるはずです。よって閉路を作って得することはないことがわかり、森になることが示されました。

さて、辺を張って森を作ったとして、そのそれぞれの木に注目すると、(操作の回数はその連結成分のサイズ $-\ 1$) になるはずです。これを条件を満たす辺の張り方にできるかどうかは、実はその頂点に書かれた数のすべての $XOR$ が $0$ であるかどうかを確かめるだけでよいです。

これを示します。まず全体の $XOR$ の値が $0$ でない場合は、いくら操作を行っても $0$ にすることはできません。なぜなら、どんな操作もある $2$ 頂点に同一の値の $XOR$ を作用させることになるため、全体の $XOR$ は不変だからです。逆に、全体の $XOR$ が $0$ であれば、葉の方から順に $0$ にしていくことができるはずです。これで示されました。

これで木の集合として条件を満たすものがわかりました。よってすべての頂点をこの集合に分割することになりますが、もとから $0$ の頂点は無視してよく、$0$ ではない同じ数同士はサイズ $2$ の木として作ってしまった方がよいことがわかります。

これを繰り返すと、取りうる値は $0$ から $15$ であることから、$0$ を無視すると、高々 $k = 15$ 個の頂点のみが残ります。これはいわゆる部分集合を扱う $O(3^k)$ の $bitDP$ で計算できます。

以上でこの問題が解けました。

実装
#include <cstdio>
#include <vector>
#include <algorithm>
using namespace std;

int main() {
        int n;
        scanf("%d", &n);
        vector<int> p(n);
        for (int i = 0; i < n - 1; i ++) {
                int a, b, c;
                scanf("%d %d %d", &a, &b, &c);
                p[a] ^= c;
                p[b] ^= c;
        }
        vector<int> cnt(16, 0);
        for (int i = 0; i < n; i ++) {
                cnt[p[i]] ++;
        }
        int ans = 0;
        int mask = 0;
        for (int i = 1; i < 16; i ++) {
                ans += cnt[i] / 2;
                if (cnt[i] & 1) {
                        mask |= 1 << i;
                }
        }
        vector<bool> ok(1 << 16, 0);
        for (int i = 0; i < (1 << 16); i ++) {
                int check = 0;
                for (int j = 0; j < 16; j ++) {
                        if (i & (1 << j)) {
                                check ^= j;
                        }
                }
                ok[i] = check == 0;
        }
        vector<int> dp(1 << 16, 0);
        for (int i = 0; i < (1 << 16); i ++) {
                if (i & 1) continue; //ignore the nodes whose value is already 0
                int sub = i; //subset of i
                while (sub > 0) {
                        if (ok[sub]) {
                                dp[i] = max(dp[i], dp[i ^ sub] + 1);
                        }
                        sub = (sub - 1) & i;
                }
        }
        printf("%d\n", ans + __builtin_popcount(mask) - dp[mask]);
        return 0;
}

Atcoder Petrozavodsk Contest 001 I. Simple APSP Problem

I - Simple APSP Problem

$2000$ 点の問題ですが、以下のアルゴリズムで使う発想ができれば、解くことができます。

Hirschberg's Algorithmについて - Learning Algorithms

English version is available here.

解法

グリッドを適当に二つに分割することを考えてみます。すると、その $2$ つのグリッドそれぞれに端点を持つ最短経路というのは、$Hirschberg's\ Algorithm$ でも見たように、そのグリッド間の境界を必ず一回だけ通ります。この分割を、白しか存在しない $2$ 行(列)において行うと、その境界を通るものだけまとめて計算することができます。その値は、それぞれのグリッドに含まれる白のマスの個数の積に等しくなります。

以降、この行(列)は考慮しなくて良いので、縮約してしまってよさそうです。実はグリッドのほとんどのマスは白なので、これを繰り返すことで、グリッドのサイズを $O(n) * O(n)$ にまで落とすことができます。

したがって縮約したグリッド上で全点対最短距離を $O(n^4)$ で求めればよさそうです。さて、これでこの問題が解けました。と、言いたいところですが、一つ注意があって、例えば以下のグリッドを縮約するとします。

f:id:KokiYamaguchi:20180205160002j:plain

この場合、画像左に示した部分を縮約することになります。したがって画像右のように、マスに重み(縮約されているマスの個数)がついていなければいけないはずです。

f:id:KokiYamaguchi:20180205160017j:plain

これは、行と列それぞれ独立に、白マスだけが続く長さと黒マスの出現を保存していき、新しくマッピングし直した座標で、それらの長さの積をとったものを重みとすることで実装できます。

縮約したグリッド上でのマス $A$ (重み $W_A$)とマス $B$ (重み $W_B$)の最短経路 $d$ は、縮約前の $A$ に含まれるすべてのマスから、縮約前の $B$ に含まれるすべてのマスへの最短距離になっているので、$d * W_A * W_B$ を全点対について計算して総和をとります。最後に重複を除くために $2$ の逆元をかけて、上で求めていた結果と足し合わせるとそれが答えです。

実装
#include <cstdio>
#include <vector>
#include <algorithm>
#include <functional>
#include <map>
#include <set>
#include <string>
#include <iostream>
#include <cassert>
#include <cmath>
#include <queue>
using namespace std;

const int MOD = 1e9 + 7;

struct state { int y, x, step; };
static const int dx[] = {1, 0, -1, 0}, dy[] = {0, 1, 0, -1};

int main() {
        long long h, w;
        scanf("%lld %lld", &h, &w);
        int n;
        scanf("%d", &n);
        vector<int> w_cnt(w, 0), h_cnt(h, 0);
        vector<bool> w_exist(w, false), h_exist(h, false);
        vector<pair<int, int>> black;
        for (int i = 0; i < n; i ++) {
                int y, x;
                scanf("%d %d", &y, &x);
                w_cnt[x] ++;
                h_cnt[y] ++;
                w_exist[x] = true;
                h_exist[y] = true;
                black.emplace_back(x, y);
        }
        for (int i = 1; i < w; i ++) w_cnt[i] += w_cnt[i - 1];
        for (int i = 1; i < h; i ++) h_cnt[i] += h_cnt[i - 1];
        long long ans = 0;
        //precalc
        for (int i = 0; i < w - 1; i ++) {
                if (!w_exist[i] && !w_exist[i + 1]) {
                        long long left = (long long) (i + 1) * h % MOD - w_cnt[i];
                        long long right = (long long) (w - (i + 1)) * h % MOD - (n - w_cnt[i]);
                        ans += left * right;
                        ans %= MOD;
                }
        }
        for (int i = 0; i < h - 1; i ++) {
                if (!h_exist[i] && !h_exist[i + 1]) {
                        long long left = (long long) (i + 1) * w % MOD - h_cnt[i];
                        long long right = (long long) (h - (i + 1)) * w % MOD - (n - h_cnt[i]);
                        ans += left * right;
                        ans %= MOD;
                }
        }
        //compress
        map<int, int> newx, newy;
        vector<pair<long long, bool>> widths, heights; //(length, is_white)
        {
                int cnt = 0;
                for (int i = 0; i < w; i ++) {
                        if (!w_exist[i]) {
                                cnt ++;
                        } else {
                                if (cnt) {
                                        widths.emplace_back(cnt, true);
                                        cnt = 0;
                                }
                                newx[i] = (int) widths.size();
                                widths.emplace_back(1, false);
                        }
                }
                if (cnt) widths.emplace_back(cnt, true);
        }
        {
                int cnt = 0;
                for (int i = 0; i < h; i ++) {
                        if (!h_exist[i]) {
                                cnt ++;
                        } else {
                                if (cnt) {
                                        heights.emplace_back(cnt, true);
                                        cnt = 0;
                                }
                                newy[i] = (int) heights.size();
                                heights.emplace_back(1, false);
                        }
                }
                if (cnt) heights.emplace_back(cnt, true);
        }
        //re-write the grid
        int neww = (int) widths.size();
        int newh = (int) heights.size();
        vector<vector<long long>> s(newh, vector<long long> (neww, 1)); //-1 when it's black, weight when it's white
        for (auto b : black) {
                s[newy[b.second]][newx[b.first]] = -1;
        }
        for (int i = 0; i < newh; i ++) {
                for (int j = 0; j < neww; j ++) {
                        if (heights[i].second || widths[j].second) {
                                s[i][j] = heights[i].first * widths[j].first % MOD;
                        }
                }
        }
        //BFS
        long long sum = 0;
        for (int sy = 0; sy < newh; sy ++) {
                for (int sx = 0; sx < neww; sx ++) {
                        if (s[sy][sx] == -1) continue;
                        long long res = s[sy][sx];
                        vector<vector<bool>> used(newh, vector<bool>(neww, false));
                        queue<state> q;
                        q.push({sy, sx, 0});
                        used[sy][sx] = true;
                        while (!q.empty()) {
                                state p = q.front(); q.pop();
                                if (p.y != sy || p.x != sx) {
                                        assert(s[p.y][p.x] != -1);
                                        sum += (long long) p.step * res % MOD * s[p.y][p.x] % MOD;
                                        sum %= MOD;
                                }
                                for (int d = 0; d < 4; d ++) {
                                        int xx = p.x + dx[d], yy = p.y + dy[d];
                                        if (xx < 0 || xx >= neww || yy < 0 || yy >= newh) continue;
                                        if (used[yy][xx] || s[yy][xx] == -1) continue;
                                        used[yy][xx] = true;
                                        q.push({yy, xx, p.step + 1});
                                }
                        }
                }
        }
        sum *= (MOD + 1) / 2;
        sum %= MOD;
        ans += sum;
        ans %= MOD;
        printf("%lld\n", ans);
        return 0;
}