Learning Algorithms

アルゴリズムの勉強メモ

AOJ Prime Caves

AOJ Prime Caves

感想

実装がめんどくさい

解法

どう見てもDPするだけ。$\ m \leq 10^6\ $なので、$\ 1000 * 1000\ $の正方形のマスを用意して考えればよい。すごく面倒だが、各数を投げるとその配列上での位置を返す関数を作る。これは各正方形の左上の数が、偶数の二乗になっていることを考慮して実装した。ここに$\ 10^6\ $以下の素数をすべて投げてそこにフラグを立てておく。

あとは、$\ n\ $が表す位置から下に向かって三角形を広げるようにDPをする。下限は計算するのが面倒なので、一番下まで計算してしまって良い。これで前半の答えがわかり、さらにこの最大値をとるような位置の中で、その数が素数であるようなものの最大値をとれば、それが後半の答えである。このとき、位置からもとの数がなんだったのかを知る必要があるので、さっきの関数の中で、ついでに位置から数がわかるように配列に書き込んでおくとよい。$\ O(m)\ $。

実装
#include "bits/stdc++.h"
using namespace std;

int s[1000][1000];
int dp[1000][1000];
int p[1000][1000];

pair<int, int> get_pos(int n) {
        if (n == 1) return make_pair(500, 499);
        int top_left = 1;
        int k = 1;
        while (n > top_left) {
                top_left = (2 * k) * (2 * k);
                k ++;
        }
        k --;
        int top_right = top_left - 2 * k + 1;
        int bottom_right = top_right - 2 * k + 1;
        int bottom_left = bottom_right - 2 * k + 1;
        int y, x;
        if (n <= bottom_left) {
                y = 499 - k + 1 + (2 * k - 1 - (bottom_left - n));
                x = 499 - k + 1;
        } else if (n <= bottom_right) {
                y = 499 + k;
                x = 499 - k + 1 + (n - bottom_left);
        } else if (n <= top_right) {
                y = 499 - k + 1 + (top_right - n);
                x = 499 + k;
        } else { 
                y = 499 - k + 1;
                x = 499 - k + 1 + (top_left - n);
        }
        p[y][x] = n;
        return make_pair(y, x);
}

int N = 1000000;
vector<int> primes;
vector<bool> is_prime(N + 1, true);
void init() {
        is_prime[0] = is_prime[1] = false;
        for (int i = 2; i <= N; i ++) {
                if (is_prime[i]) {
                        primes.push_back(i);
                        for (int j = i + i; j <= N; j += i) is_prime[j] = false;
                }
        }
}

int main() {
        init();
        int m, n;
        while (scanf("%d%d", &m, &n) && n) {
                for (int i = 0; i < 1000; i ++) for (int j = 0; j < 1000; j ++) s[i][j] = 0;
                for (auto p : primes) {
                        if (p > m) break;
                        int y, x;
                        tie(y, x) = get_pos(p);
                        s[y][x] = 1;
                }
                int dp[1000][1000] = { 0 };
                int sy, sx;
                tie(sy, sx) = get_pos(n);
                dp[sy][sx] = s[sy][sx];
                for (int dy = 0; sy + dy + 1 < 1000; dy ++) {
                        for (int dx = -dy; dx <= dy && sx + dx < 1000; dx ++) {
                                dp[sy + dy + 1][sx + dx - 1] = max(dp[sy + dy + 1][sx + dx - 1], s[sy + dy + 1][sx + dx - 1] + dp[sy + dy][sx + dx]);
                                dp[sy + dy + 1][sx + dx]     = max(dp[sy + dy + 1][sx + dx]    , s[sy + dy + 1][sx + dx]     + dp[sy + dy][sx + dx]);
                                dp[sy + dy + 1][sx + dx + 1] = max(dp[sy + dy + 1][sx + dx + 1], s[sy + dy + 1][sx + dx + 1] + dp[sy + dy][sx + dx]);
                        }
                }
                int ans = 0;
                for (int i = 0; i < 1000; i ++) {
                        for (int j = 0; j < 1000; j ++) {
                                ans = max(ans, dp[i][j]);
                        }
                }
                if (ans == 0) { 
                        cout << 0 << ' ' << 0 << endl;
                        continue;
                }
                int res = -1;
                for (int i = 0; i < 1000; i ++) {
                        for (int j = 0; j < 1000; j ++) {
                                if (dp[i][j] == ans && s[i][j]) {
                                        res = max(res, p[i][j]);
                                }
                        }
                }
                cout << ans << ' ' << res << endl;
        }
        return 0;
}

Atcoder Grand Contest 004 F. Namori

Atcoder Grand Contest 004 F. Namori

F: Namori - AtCoder Grand Contest 004 | AtCoder

感想

解にたどり着くまでに考えることがたくさんある。ケースに応じて適切なアルゴリズムに落とし込んでいくのが、おもしろい。

解法

まず連続する同色の頂点を反転させるという操作は、(白、白)->(黒、黒)と(黒、黒)->(白、白)の2パターンがあるため、考えるのが難しい。そこで、どんな操作も(白、黒)->(黒、白)という操作に変えることができれば、扱うのが容易になるだろう。そこで、二部グラフの要領で白黒に塗ってしまったグラフを考え、そのグラフとの$\ XOR\ $をとったような色の状態を考えることにする。すると、ある連続する2点を選んだ時に、その一方のみが$\ XOR\ $をとったようなものになるはずなので、(白、黒)->(黒、白)という操作のみに落とし込むことができる。

さらにこれを「反転の操作」ではなく、「黒の位置を白の位置に動かす操作」に見ることができることに気付く。すると結局、すべての黒い部分を白い部分に移動させるために必要な操作回数を求める問題になる。

この発想を思いつくのは難しいのだが、下の類題を見れば、割と典型であることに気付くだろう。下の問題は反転の操作を駒を動かす操作に見立てるところまで同じである(反転と$\ XOR\ $のテク)。

F: Prime Flip - AtCoder Regular Contest 080 | AtCoder

ここまでくれば、$\ m = n - 1\ $の場合は易しい問題に変わる。まずグラフが木なので、二部グラフの要領で黒を塗っていく。そのあと、この黒を白に移動させていくのだが、まずこの個数が一致していなければ明らかに不可能であることがわかる。逆に一致していれば有限回の操作で可能である。

さて、グラフが木であることを利用すると、ある辺に注目したときに、どちらかの部分木を見て、そこに黒の頂点が$\ B\ $個、白の頂点が$\ W\ $個あったとすると、無駄なく移動させるためには、この辺を経由する回数というのは明らかに$\ |B - W|\ $回であることが言える。よって各辺についてこの値を計算し、その総和をとれば答えである。これは簡単な$\ DFS\ $によって$\ O(n)\ $で計算可能である。

次に、$\ m = n\ $の場合を考える。これはNamoriグラフと呼ばれ、一つだけ閉路を含む木のような構造である。この閉路が偶閉路なのか奇閉路なのかで場合分けをする。

まずグラフが偶閉路をもつ場合を考える。この場合は木の場合と同様に(グラフが二部グラフなので)白黒に塗っておくことができる。これらの白黒の個数が一致していない場合はやはり不可能である。

さて、閉路上にない辺は、上のようにしてどちらかの部分木の白黒の数を数えて$\ |B - W|\ $なる数を加算していけば良い。閉路上の辺については、どこか一つの辺の移動量を固定しないと計算できない。そこで、ある閉路上の辺$\ (s, t)\ $について、ここを移動する黒の数を$\ 0\ $と仮定し、この辺を切ってしまうことにする。

このあとに、上と同じように移動量を考えていくのだが、閉路上の辺については、$\ B - W\ $の値を一旦vector flowなどに保管しておくことにする。そして、このflowから閉路上の移動量がどうなっているのかを求めよう。

$\ s\ $と$\ t\ $の間の移動量を$\ 0\ $に仮定した場合の移動量がflowに格納されている。したがって、この場合の答えは、flowのそれぞれの値の絶対値の総和をとればよい。仮に$\ s, t\ $間の移動量を$\ 1\ $にした場合は、このflowのそれぞれの数から$\ 1\ $を引いた数の絶対値をとったものの総和をとり、それに$\ 1\ $を足すことで答えになる。一般に$\ s, t\ $の移動量が$\ x\ $だったとすると、答えは$\ \sum |flow_i - x| + |x|\ $となる。この式をぐっと睨むとこれは下に凸だと気付くので、三分探索して最小値をとればそれが答えになる。この計算量は$\ O(n \log n)\ $である。

ある辺が閉路上の辺であるかどうかの判定は、ある頂点を根にして$\ DFS\ $しているときに、ある頂点以下に$\ s, t\ $がそれぞれ含まれているかどうかをみる配列を用意しておき、どちらか一方のみ含まれているならば、その辺は閉路上の辺であることがわかる。

三分探索は、終わらせ方がよくわからなかったので、幅が$\ 3\ $とかになったら終了して、あとは幅$\ 20\ $くらいで適当に最小値を見つけるようにした。

最後に、グラフが奇閉路をもつ場合を考える。この場合は、そもそも二部グラフではないので白黒に塗り分けることができないのだが、閉路は一つなので、塗り分け不可能な部分は必ず一つである。その辺が結ぶ2点上では、(白、白)->(黒、黒)あるいは(黒、黒)->(白、白)の操作をしてよいものとして考えることにする。ここで白黒の個数調整ができると考えると、初期状態で白と黒の偶奇が一致していれば可能であることがいえる。

この個数調整を回数を$\ cnt\ $としておくと、$\ s, t\ $においては黒$\ cnt\ $分多くみて考え、この辺を切ってしまえば、あとは木の場合と全く同様に$\ O(n)\ $で処理することができる。

解法がめちゃくちゃ長くなった。

実装
#include &quot;bits/stdc++.h&quot;
using namespace std;

vector<int> g[101010];
int color[101010];
int black, white, type, cnt;
int s, t;
long long ans;
bool in_s[101010], in_t[101010];
vector<long long> flow;

void check(int v, int prev, int d) {
        color[v] = d;
        if (d & 1) black ++;
        else white ++;
        for (auto u : g[v]) if (u != prev) { 
                if (color[u] == -1) {
                        check(u, v, 1 - d);
                } else {
                        s = v, t = u;
                        if (color[u] == 1 - d) type = 1;
                        else type = 2;
                }
        }
}

pair<int, int> dfs(int v, int prev, int d) {
        int b = 0, w = 0;
        for (auto u : g[v]) if (u != prev) {
                pair<int, int> get = dfs(u, v, d + 1);
                b += get.first;
                w += get.second;
        }
        if (d & 1) b ++;
        else w ++;
        if (type == 2 && (v == s || v == t)) b -= cnt;
        ans += abs(b - w);
        return make_pair(b, w);
}

pair<int, int> dfs2(int v, int prev, int d) {
        int b = 0, w = 0;
        in_s[v] = false, in_t[v] = false;
        in_s[v] |= v == s;
        in_t[v] |= v == t;
        for (auto u : g[v]) if (u != prev) {
                pair<int, int> get = dfs2(u, v, d + 1);
                in_s[v] |= in_s[u];
                in_t[v] |= in_t[u];
                b += get.first;
                w += get.second;
        }
        if (d & 1) b ++;
        else w ++;
        if ((in_s[v] && in_t[v]) || (!in_s[v] && !in_t[v])) ans += abs(b - w);
        else flow.push_back(b - w);
        return make_pair(b, w);
}

int main() {
        int n, m;
        scanf(&quot;%d%d&quot;, &n, &m);
        for (int i = 0; i < m; i ++) {
                int a, b;
                scanf(&quot;%d%d&quot;, &a, &b);
                a --, b --;
                g[a].push_back(b);
                g[b].push_back(a);
        }
        type = 0;
        memset(color, -1, sizeof color);
        check(0, -1, 0);
        if ((type == 0 && black != white) ||
            (type == 1 && black != white) ||
            (type == 2 && black % 2 != white % 2)) {
                puts(&quot;-1&quot;);
                return 0;
        }
        if (type == 0) dfs(0, -1, 0);
        g[s].erase(remove(g[s].begin(), g[s].end(), t), g[s].end());
        g[t].erase(remove(g[t].begin(), g[t].end(), s), g[t].end());
        if (type == 1) {
                dfs2(0, -1, 0);
                int lb = -101010, ub = 101010;
                while (ub - lb > 3) {
                        int left = (lb * 2 + ub) / 3;
                        int right = (lb + ub * 2) / 3;
                        long long left_ans = 0, right_ans = 0;
                        for (auto f : flow) {
                                left_ans += abs(left + f);
                                right_ans += abs(right + f);
                        }
                        left_ans += abs(left); //self
                        right_ans += abs(right); //self
                        if (left_ans < right_ans) ub = right;
                        else lb = left;
                }
                long long res = 1LL << 60;
                for (int i = lb - 10; i < lb + 10; i ++) {
                        long long sum = 0;
                        for (auto f : flow) sum += abs(i + f);
                        sum += abs(i); //self
                        res = min(res, sum);
                }
                ans += res;
        }
        if (type == 2) {
                cnt = (black - white) / 2;
                dfs(0, -1, 0);
                ans += abs(cnt);
        }
        printf(&quot;%lld\n&quot;, ans);
        return 0;
}

Atcoder Grand Contest 014 E. Blue and Red Tree

Atcoder Grand Contest 014 E. Blue and Red Tree

E: Blue and Red Tree - AtCoder Grand Contest 014 | AtCoder

感想

難しいがおもしろい。解説を読んで、少し人のコードを参考にして書いた。

解法

辺の色が異なる2つの木を重ね合せて考える。そのグラフ上で、ある2頂点$\ a, b\ $に対して青と赤の辺がどちらも張られているとする。この2頂点に関する操作はどのタイミングでも可能なので、あとの方でやった方がよい。

さらにこれらの2点間は、他の点での操作をするときに、いつでも移動可能であるため、この2点に辺を張る操作を「最後にする」という意味をもたせて、縮約して考えてみる。この新しくできたグラフの上でまた同様のことを考えていく。 もしこの縮約を繰り返していって一つの頂点にまとめてしまえるならば、その縮約をしていった逆の順番で、辺を構成していくことが可能であるといえる。

逆に頂点が2個以上残っているある段階でそれ以上縮約できない状態になったとしたら、小さい図を書いて実験してみると、なんとなくNOになりそうだとわかる。木である性質を使って色々場合分けしたら証明できると思う。ちゃんとした証明はeditorialで。

この縮約はある部分でやったからといって他の部分でできなくなるということはないので貪欲にやってよい。そのためにはいわゆるマージテクを使う。

簡単に辺を削除できるように、setで辺の情報をもつようにした。さらに、辺の色は、重複しているときのみ考えればよいため、色は実質無視して良い。

実装の方針としては、辺の重複が現れるたびに、queueなどのデータ構造にその頂点のペアを保管しておく。そこから取り出した頂点のペアが$\ (a, b)\ $だとすると、その頂点につながる辺の数が少ない方を、多い方にマージするようにする。多い方を$\ a\ $とすると、まず$\ a\ $と$\ b\ $の間の辺を取り除き、$\ b\ $とつながる各頂点を$\ a\ $とつなぎ、最後に$\ b\ $からの情報をすべてクリアする。

ただし、注意すべきことは、ある時点でpushされた頂点のペアが、popされるまでの間に縮約によって別の頂点に変化していることがある。それに対処するために、UnionFindのようなものを使う。すなわち、根が自身がマージされたノードであるような構造を持っておけば、矛盾なくマージされた情報がわかる。

細かいことだが、もう一つ着目するところがあって、途中にあるif (g[a].size() < g[b].size()) swap(a, b);は必要に思われるが、実はこれを書かなくても高速に動く(というかあまりswapが実行されない)。実はemplace(a, c)で追加されるときにまだ十分にマージされていないcの方が小さい確率が十分高く、その結果bothの中の(a, b)についてもbの方が小さい確率が十分に高いからである。

実装
#include "bits/stdc++.h"
using namespace std;

#define mp      make_pair

int par[101010];

int find(int x) { return par[x] == x ? x : par[x] = find(par[x]); }

int main() {
        int n;
        scanf("%d", &n);
        set<int> g[101010];
        queue<pair<int, int>> both;
        for (int i = 0; i < n; i ++) par[i] = i;
        for (int i = 0; i < 2 * n - 2; i ++) {
                int a, b;
                scanf("%d%d", &a, &b);
                a --, b --;
                if (g[a].count(b)) { 
                        both.emplace(a, b);
                } else { 
                        g[a].insert(b);
                        g[b].insert(a);
                }
        }
        int ans = 0;
        while (!both.empty()) {
                ans ++;
                int a, b;
                tie(a, b) = both.front();
                both.pop();
                a = find(a), b = find(b);
                //if (g[a].size() < g[b].size()) swap(a, b);
                a = par[a], b = par[b];
                if (a == b) continue;
                par[b] = a;
                g[a].erase(b);
                g[b].erase(a);
                for (auto c : g[b]) {
                        g[c].erase(b);
                        if (g[a].count(c)) { 
                                both.emplace(a, c);
                        } else { 
                                g[a].insert(c);
                                g[c].insert(a);
                        }
                }
                g[b].clear();
        }
        puts(ans == n - 1 ? "YES" : "NO");
        return 0;
}

行列ライブラリ

行列ライブラリ

競プロ用に行列ライブラリを作った。

特別に綺麗だとか高速だとかそういうことはないけれど、もし利用したい人がいれば、勝手に利用していただいても構いません。多分動きます。

Typeは、値が小数になる場合と、逆行列が必要なときにdoubleにしておく。その他の場合はintでよい。
int GetRank(Matrix a)は任意の行列のランクを返す。
bool Inv(Matrix a, Matrix &inv)は第二引数に逆行列として受け取りたい、aと同じ大きさの空の行列を参照で渡して使う。もし逆行列が存在するならばtrueを返し、第二引数の中が逆行列に書き換えられる。存在しないならばfalseを返す。
AddSubは常に引数を2つ渡すようにすればよい。
いずれの関数も、定義できない演算をしようとすると、assertで怒られる。

また追記、修正するかもしれない。

typedef double Type;
typedef vector<vector<Type>> Matrix;
int GetRank(Matrix a) {
        int h = a.size(), w = a[0].size();
        int res = 0, now = 0;
        for (int i = 0; i < h; i ++) {
                Type ma = 0.0;
                int pivot;
                for (int j = i; j < h; j ++) {
                        if (a[j][now] > ma) {
                                ma = a[j][now];
                                pivot = j;
                        }
                }
                if (ma == 0.0) {
                        now ++;
                        if (now == w) break;
                        i --;
                        continue;
                }
                if (pivot != i) {
                        for (int j = 0; j < w; j ++) { 
                                swap(a[i][j], a[pivot][j]);
                        }

                }
                Type tmp = 1.0 / a[i][now];
                for (int j = 0; j < w; j ++) a[i][j] *= tmp;
                for (int j = 0; j < h; j ++) {
                        if (i != j) {
                                Type tmp2 = a[j][now];
                                for (int k = 0; k < w; k ++) {
                                        a[j][k] -= a[i][k] * tmp2;
                                }
                        }
                }
                res ++;
        }
        return res;
}
bool Inv(Matrix a, Matrix &inv) {
        assert(a.size() == a[0].size() && inv.size() == inv[0].size());
        int n = a.size();
        for (int i = 0; i < n; i ++) {
                for (int j = 0; j < n; j ++) {
                        inv[i][j] = (i == j ? 1.0 : 0.0);
                }
        }
        for (int i = 0; i < n; i ++) {
                Type ma = 0.0;
                int pivot;
                for (int j = i; j < n; j ++) { 
                        if (a[j][i] > ma) {
                                ma = a[j][i];
                                pivot = j;
                        }
                }
                if (ma == 0.0) return false;
                if (pivot != i) {
                        for (int j = 0; j < n; j ++) { 
                                swap(a[i][j], a[pivot][j]);
                                swap(inv[i][j], inv[pivot][j]);
                        }

                }
                Type tmp = 1.0 / a[i][i];
                for (int j = 0; j < n; j ++) {
                        a[i][j] *= tmp;
                        inv[i][j] *= tmp;
                }
                for (int j = 0; j < n; j ++) {
                        if (i != j) {
                                Type tmp2 = a[j][i];
                                for (int k = 0; k < n; k ++) {
                                        a[j][k] -= a[i][k] * tmp2;
                                        inv[j][k] -= inv[i][k] * tmp2;
                                }
                        }
                }
        }
        return true;
}
Matrix Add(const Matrix &a, const Matrix &b, bool minus = false) {
        assert(a.size() == b.size() && a[0].size() == b[0].size());
        int h = a.size(), w = a[0].size();
        Matrix c(h, vector<Type> (w));
        for (int i = 0; i < h; i ++) {
                for (int j = 0; j < w; j ++) {
                        c[i][j] = a[i][j] + (minus ? -1 : 1) * b[i][j];
                }
        }
        return c;
}
Matrix Sub(const Matrix &a, const Matrix &b) {
        return Add(a, b, true);
}
Matrix Mul(const Matrix &a, const Matrix &b) {
        assert(a[0].size() == b.size());
        Matrix c(a.size(), vector<Type> (b[0].size()));
        for (int i = 0; i < a.size(); i ++) {
                for (int k = 0; k < b.size(); k ++) {
                        for (int j = 0; j < b[0].size(); j ++) {
                                c[i][j] = (c[i][j] + a[i][k] * b[k][j]);
                        }
                }
        }
        return c;
}
Matrix Pow(Matrix a, long long n) {
        assert(a.size() == a[0].size());
        Matrix b(a.size(), vector<Type> (a.size()));
        for (int i = 0; i < a.size(); i ++) {
                b[i][i] = 1;
        }
        while (n > 0) {
                if (n & 1) b = Mul(b, a);
                a = Mul(a, a);
                n >>= 1;
        }
        return b;
}
void PrintMatrix(const Matrix &a) {
        int h = a.size(), w = a[0].size();
        for (int i = 0; i < h; i ++) {
                for (int j = 0; j < w; j ++) {
                        cout << a[i][j] << ' ';
                }
                cout << endl;
        }
}

AOJ Graph Automata Player

AOJ Graph Automata Player

Graph Automata Player | Aizu Online Judge

感想

ICPCのチーム練習で解いた。好きな感じの問題なので、じーーっと考察していると解法がわかった。

解法

ある頂点に注目して考える。ある状態から次の状態に遷移するとき、その頂点から移動可能な点に書かれた数の総和が奇数であるとき、その頂点は$\ 1\ $になる、という風に読み替えるとわかりやすい。よく考えると、この頂点から(自分を含む)(他の各頂点に対して接続しているかどうか)の情報列と、(各頂点の値)の情報列の内積のようなもの(?)をとれば、その総和が計算できる。このことから、隣接行列を$\ G\ $、時刻$\ t\ $での状態を$\ S_{t}\ $、時刻$\ t - 1\ $での状態を$\ S_{t - 1}\ $とすると、$\ S_{t} = G * S_{t - 1}\ $と表すことができる。

よって、$\ G\ $の逆行列が存在するならば、明らかに$\ S_{-T} = (G^{-1})^T * S_{0}\ $が成り立つので、これを計算するだけでよい。

存在しない場合は、解の存在性を見ればよい。すなわち、$\ G^T * S_{-T} = S_{0}\ $という方程式を満たす解$\ S_{-T}\ $が存在するかどうかを判定できればよい。線形代数の教科書を取り出して復習してみると、$\ Ax = b\ $の解が存在するための必要十分条件が$\ rank([A : b]) = rank(A)\ $が成り立つことであるとわかるため、そのように判定する。解が存在すればambiguous、存在しないならばnoneである。

ただし、行列の中身は常に$\ mod\ 2\ $で見ておき、小数が現れないようにすることができる。

実装
#include "bits/stdc++.h"
using namespace std;

typedef int Type;
typedef vector<vector<Type>> Matrix;

int GetRank(Matrix a) {
        int h = a.size(), w = a[0].size();
        int res = 0, now = 0;
        for (int i = 0; i < h; i ++) {
                Type ma = 0.0;
                int pivot;
                for (int j = i; j < h; j ++) {
                        if (a[j][now] > ma) {
                                ma = a[j][now];
                                pivot = j;
                        }
                }
                if (ma == 0.0) {
                        now ++;
                        if (now == w) break;
                        i --;
                        continue;
                }
                if (pivot != i) {
                        for (int j = 0; j < w; j ++) { 
                                swap(a[i][j], a[pivot][j]);
                        }

                }
                Type tmp = 1.0 / a[i][now];
                for (int j = 0; j < w; j ++) a[i][j] *= tmp;
                for (int j = 0; j < h; j ++) {
                        if (i != j) {
                                Type tmp2 = a[j][now];
                                for (int k = 0; k < w; k ++) {
                                        a[j][k] = ((a[j][k] - a[i][k] * tmp2) + 2) % 2;
                                }
                        } }
                res ++;
        }
        return res;
}
bool Inv(Matrix a, Matrix &inv) { //square
        int n = a.size();
        for (int i = 0; i < n; i ++) {
                for (int j = 0; j < n; j ++) {
                        inv[i][j] = (i == j ? 1.0 : 0.0);
                }
        }
        for (int i = 0; i < n; i ++) {
                Type ma = 0.0;
                int pivot;
                for (int j = i; j < n; j ++) { 
                        if (a[j][i] > ma) {
                                ma = a[j][i];
                                pivot = j;
                        }
                }
                if (ma == 0.0) return false; //no inverse matrix
                if (pivot != i) {
                        for (int j = 0; j < n; j ++) { 
                                swap(a[i][j], a[pivot][j]);
                                swap(inv[i][j], inv[pivot][j]);
                        }

                }
                Type tmp = 1.0 / a[i][i];
                for (int j = 0; j < n; j ++) {
                        a[i][j] *= tmp;
                        inv[i][j] *= tmp;
                }
                for (int j = 0; j < n; j ++) {
                        if (i != j) {
                                Type tmp2 = a[j][i];
                                for (int k = 0; k < n; k ++) {
                                        a[j][k] = ((a[j][k] - a[i][k] * tmp2) + 2 ) % 2;
                                        inv[j][k] = ((inv[j][k] - inv[i][k] * tmp2) + 2) % 2;
                                }
                        }
                }
        }
        return true;
}

Matrix Mul(const Matrix &a, const Matrix &b) {
        Matrix c(a.size(), vector<Type> (b[0].size()));
        for (int i = 0; i < a.size(); i ++) {
                for (int k = 0; k < b.size(); k ++) {
                        for (int j = 0; j < b[0].size(); j ++) {
                                c[i][j] = ((c[i][j] + a[i][k] * b[k][j])) % 2;
                        }
                }
        }
        return c;
}
Matrix Pow(Matrix a, long long n) { //a^n, square
        Matrix b(a.size(), vector<Type> (a.size()));
        for (int i = 0; i < a.size(); i ++) {
                b[i][i] = 1;
        }
        while (n > 0) {
                if (n & 1) b = Mul(b, a);
                a = Mul(a, a);
                n >>= 1;
        }
        return b;
}
void PrintMatrix(const Matrix &a) {
        int h = a.size(), w = a[0].size();
        for (int i = 0; i < h; i ++) {
                for (int j = 0; j < w; j ++) {
                        cout << a[i][j] << ' ';
                }
                cout << endl;
        }
}
int main() {
        int n;
        scanf("%d", &n);
        Matrix g(n, vector<int> (n));
        for (int i = 0; i < n; i ++) {
                for (int j = 0; j < n; j ++) {
                        scanf("%d", &g[i][j]);
                }
        }
        Matrix v(n, vector<int> (1));
        for (int i = 0; i < n; i ++) scanf("%d", &v[i][0]);
        long long t;
        scanf("%lld", &t);
        Matrix inv(n, vector<int> (n));
        if (Inv(g, inv)) {
                inv = Pow(inv, t);
                auto ans = Mul(inv, v);
                for (int i = 0; i < n; i ++) cout << ans[i][0] << (i == n - 1 ? '\n' : ' ');
        } else {
                auto g_t = Pow(g, t);
                Matrix g_tv(n, vector<int> (n + 1));
                for (int i = 0; i < n; i ++) {
                        for (int j = 0; j < n; j ++) {
                                g_tv[i][j] = g_t[i][j];
                        }
                }
                for (int i = 0; i < n; i ++) g_tv[i][n] = v[i][0];
                if (GetRank(g_tv) == GetRank(g_t)) {
                        cout << "ambiguous" << endl;
                } else {
                        cout << "none" << endl;
                }
        }
        return 0;
}

CODE FESTIVAL 2017 qual B C. 3 Steps

CODE FESTIVAL 2017 qual B C. 3 Steps

C: 3 Steps - CODE FESTIVAL 2017 qual B | AtCoder

感想

qual Aには出れなかったので、qual Bから参戦していた。40分で3完し、そのあとはひたすらDのバグをつぶそうとしていた。

総合243位、日本人98位なので本戦には行けないが、得意な感じの問題がCで出てくれたのでパフォーマンスは悪くなかった。そのCについて解法をまとめておく。

解法

移動をスタートする頂点を$\ S\ $に固定して考える。まず頂点$\ S\ $には$\ 0\ $と書き、$\ S\ $から進んだ距離を各頂点に書いていくことにする。ただし、$\ 2\ $と書いた次の頂点には、$\ 1\ $と書くことにする(つまり交互に書く)。なぜなら、本来$\ 3\ $と書かれるべき頂点には、$\ S\ $から辺を張れるので、それは$\ 1\ $とみなせるからである。すると、$\ S\ $を除くすべての頂点が$\ 1\ $または$\ 2\ $で塗られることになる。ただし、これは少し嘘で、塗り方に矛盾が生まれる場合があって、それは奇閉路が存在するときに限る。$\ 1\ $と書かれた頂点には辺を張ることが可能であり、$\ 2\ $と書かれた頂点には辺を張ることが不可能である。先に述べたように、奇閉路があるときは、すべての頂点を$\ 1\ $で塗れるので、すべての頂点に辺を張ることが可能である。
グダグダ書いたが結局二部グラフかどうかを見て、適当に完全グラフまたは二部グラフを作って辺を数えればよい。$\ O(n)\ $。答えは最大$\ 10^{10}\ $程度にまで大きくなるのでオーバーフローに気をつけること。

実装
#include "bits/stdc++.h"
using namespace std;

void dfs(int v, int prev, int color, vector<int> &used, vector<vector<int>> &g, bool &can) {
        used[v] = color;
        for (auto u : g[v]) if (u != prev) {
                if (used[u] == -1) dfs(u, v, 1 - color, used, g, can);
                else if (used[u] != 1 - color) can = true;
        }
}

int main() {
        int n, m;
        cin >> n >> m;
        vector<vector<int>> g(n);
        for (int i = 0; i < m; i ++) {
                int a, b;
                cin >> a >> b;
                a --, b --;
                g[a].push_back(b);
                g[b].push_back(a);
        }
        vector<int> used(n, -1);
        bool can = false;
        dfs(0, -1, 0, used, g, can);
        long long ans;
        if (can) ans = (long long) n * (n - 1) / 2;
        else {
                int black = 0;
                for (int i = 0; i < n; i ++) black += used[i];
                ans = (long long) black * (n - black);
        }
        cout << ans - m << endl;
        return 0;
}

CODE FESTIVAL 2017 qual B D. 101 to 010

CODE FESTIVAL 2017 qual B D. 101 to 010

D: 101 to 010 - CODE FESTIVAL 2017 qual B | AtCoder

感想

コンテスト中はあやふやな実装でバグらせた。
DEGwerさんの提出コードを見ると、驚くほどわかりやすく実装されていたのでさすがだと思った。左右端の処理もきれい。
考えていることをしっかり整理してコードに落とす力が足りていない。

解法

左から普通にDPをする。DPの更新は101となっている部分に到達した時に、左側と右側の$\ 1\ $の個数をみていくようなDPで構わない。一見すると$\ O(n^2)\ $に見えるが、よく考えると各$\ 1\ $は高々$\ 2\ $回しか見ないので、$\ O(n)\ $で高速に計算される。

実装
#include "bits/stdc++.h"
using namespace std;

int dp[505050];

int main() {
        int n; 
        scanf("%d", &n);
        string s; 
        cin >> s;
        s = "0" + s + "0";
        for (int i = 1; i < s.size() - 1; i ++) {
                dp[i + 1] = max(dp[i + 1], dp[i]);
                if (s[i - 1] == '1' && s[i] == '0' && s[i + 1] == '1') {
                        int prev = i - 1, next = i + 1;
                        while (s[prev] != '0') {
                                dp[i + 1] = max(dp[i + 1], dp[prev - 1] + (i - prev));
                                prev --;
                        }
                        while (s[next] != '0') {
                                dp[next] = max(dp[next], dp[i - 2] + (next - i));
                                next ++;
                        }
                }
        }
        printf("%d\n", dp[s.size() - 1]);
        return 0;
}