Learning Algorithms

アルゴリズムの勉強メモ

Atcoder Regular Contest 013 C. 笑いをとれるかな?

Atcoder Regular Contest 013 C. 笑いをとれるかな?

C: 笑いをとれるかな? - AtCoder Regular Contest 013 | AtCoder

感想

友人が問題文の内容がおもしろいと言っていたのでやったみた。確かにおもしろかった。

解法

各座標軸を基準にみて端点同士の点より内側にある部分は決して食べることができない。よってこれは大きなハバネロが一つ埋め込まれているものとして考えてよく、それぞれの座標のハバネロの位置の最小値と最大値がわかれば、あとはその外側の幅の部分を食べていけることになる。

この幅が高々$\ 6 * N\ $個あって、それらを山とみれば、単純な$\ Nim\ $になっているので、$\ XOR\ $をとると答えがわかる。

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

static const int INF = 0x3f3f3f3f;

int main() {
        int n;
        scanf("%d", &n);
        int ans = 0;
        for (int i = 0; i < n; i ++) {
                int X, Y, Z;
                int m;
                scanf("%d%d%d%d", &X, &Y, &Z, &m);
                int min_x = INF, max_x = 0, min_y = INF, max_y = 0, min_z = INF, max_z = 0;
                for (int i = 0; i < m; i ++) {
                        int x, y, z;
                        scanf("%d%d%d", &x, &y, &z);
                        min_x = min(min_x, x);
                        max_x = max(max_x, x);
                        min_y = min(min_y, y);
                        max_y = max(max_y, y);
                        min_z = min(min_z, z);
                        max_z = max(max_z, z);
                }
                max_x = X - max_x - 1;
                max_y = Y - max_y - 1;
                max_z = Z - max_z - 1;
                int res = 0;
                res = min_x ^ max_x ^ min_y ^ max_y ^ min_z ^ max_z;
                ans ^= res;
        }
        puts(ans ? "WIN" : "LOSE");
        return 0;
}

Atcoder Grand Contest 011 E. Increasing Numbers

Atcoder Grand Contest 011 E. Increasing Numbers

E: Increasing Numbers - AtCoder Grand Contest 011 | AtCoder

感想

数論の問題などでたまに見かけるレピュニットの性質をうまく利用した問題で少しおもしろかった。

解法

増加的な数を9個以下のレピュニットの和で表してやると解ける。

レピュニットを$\ R\ $とすると、$\ 9R + 1\ $がきれいになることを利用するのは下の問題でも同じである。

Problem 132 - Project Euler

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

//----------------------------------------Library---------------------------------------------

class Bigint {
public:
        static const int BASE = 100000000, LEN = 8;
        bool negative;
        std::vector<int> a;
        Bigint& normalize();
public:
        Bigint(int x = 0);
        Bigint(const std::string& s);
        Bigint& operator = (const Bigint& x);
        Bigint& operator = (int x);
        Bigint& operator = (const std::string& s);
        const bool operator < (const Bigint& x) const;
        const bool operator > (const Bigint& x) const;
        const bool operator <= (const Bigint& x) const;
        const bool operator >= (const Bigint& x) const;
        const bool operator != (const Bigint& x) const;
        const bool operator == (const Bigint& x) const;
        Bigint operator -() const;
        Bigint& operator += (const Bigint& x);
        Bigint& operator -= (const Bigint& x);
        Bigint& operator *= (const Bigint& x);
        Bigint& operator /= (const Bigint& x);
        Bigint& operator %= (const Bigint& x);
        const Bigint operator + (const Bigint& x) const;
        const Bigint operator - (const Bigint& x) const;
        const Bigint operator * (const Bigint& x) const;
        const Bigint operator / (const Bigint& x) const;
        const Bigint operator % (const Bigint& x) const;
        friend std::pair<Bigint, Bigint> divmod(const Bigint& lhs, const Bigint& rhs);
        friend std::istream& operator >> (std::ostream& is, Bigint& x);
        friend std::ostream& operator << (std::ostream& os, const Bigint& x);
        friend const Bigint abs(Bigint x);
};
Bigint& Bigint::normalize() {
        int i = a.size()-1;
        while (i >= 0 && a[i] == 0) i --;
        a.resize(i + 1);
        if (a.size() == 0) negative = false;
        return *this;
}
Bigint::Bigint(int x) : negative(x < 0) {
        x = abs(x);
        for (; x > 0; x /= BASE) a.push_back(x % BASE);
}
Bigint::Bigint(const std::string& s): negative(false) {
        int p = 0;
        if (s[p] == '-') { p ++; negative = true; }
        else if (s[p] == '+') p ++;
        for (int i = s.size() - 1, v = BASE; i >= p; i --, v *= 10) {
                int x = s[i] - '0';
                if (x < 0 || 9 < x) {
                        std::cerr << "error: parse error:" << s << std::endl;
                        exit(1);
                } 
                if (v == BASE) {
                        v = 1;
                        a.push_back(x);
                } else {
                        a.back() += x * v;
                }
        }
        normalize();
}
Bigint& Bigint::operator = (const Bigint& x) {
        negative = x.negative;
        a = x.a;
        return *this;
}
Bigint& Bigint::operator = (int x) { return *this = Bigint(x); }
Bigint& Bigint::operator = (const std::string& s) { return *this = Bigint(s); }
const bool Bigint::operator < (const Bigint& x) const {
        if (negative != x.negative) return negative < x.negative;
        if (a.size() != x.a.size()) return (a.size() < x.a.size()) ^ negative;
        for(int i = a.size()-1; i >= 0; i --)
                if (a[i] != x.a[i]) return (a[i] < x.a[i]) ^ negative;
        return false;
}
const bool Bigint::operator > (const Bigint& x) const { return x < (*this); }
const bool Bigint::operator <= (const Bigint& x) const { return !(x < (*this)); }
const bool Bigint::operator >= (const Bigint& x) const { return !((*this) < x); }
const bool Bigint::operator != (const Bigint& x) const { return (*this) < x || x < (*this); }
const bool Bigint::operator == (const Bigint& x) const { return !((*this) < x || x < (*this)); }
Bigint Bigint::operator -() const {
        Bigint ret(*this);
        if (a.size()) ret.negative = !ret.negative;
        return ret;
}
Bigint& Bigint::operator += (const Bigint& x) {
        if (negative != x.negative) return *this -= -x;
        if (a.size() < x.a.size()) a.resize(x.a.size());
        int up = 0;
        for (int i = 0; i < a.size(); i ++) {
                a[i] += (i < x.a.size() ? x.a[i] : 0) + up;
                up = a[i] / BASE;
                a[i] %= BASE;
        }
        if (up) a.push_back(1);
        return *this;
}
Bigint& Bigint::operator -= (const Bigint& x) {
        if (negative != x.negative) return *this += -x;
        std::vector<int> b(x.a);
        if ((*this < x) ^ negative) {
                a.swap(b);
                negative = !negative;
        }
        for (int i = 0, up = 0; i < a.size(); i ++) {
                a[i] += BASE - (i < b.size() ? b[i] : 0) + up;
                up = a[i] / BASE - 1;
                a[i] %= BASE;
        }
        return this->normalize();
}
Bigint& Bigint::operator *= (const Bigint& x) {
        negative ^= x.negative;
        std::vector<int> c(a.size() * x.a.size() + 1);
        for (int i = 0; i < a.size(); i ++) {
                long long tmp = 0;
                for (int j = 0; j < x.a.size(); j ++) {
                        long long v = (long long)a[i] * x.a[j] + c[i+j] + tmp;
                        tmp = v / BASE;
                        c[i + j] = (int)(v % BASE);
                }
                if (tmp) c[i + x.a.size()] += (int)tmp;
        }
        a.swap(c);
        return this->normalize();
}
Bigint& Bigint::operator /= (const Bigint& x) {
        return *this = divmod(*this, x).first;
}
Bigint& Bigint::operator %= (const Bigint& x) {
        return *this = divmod(*this, x).second;
}
const Bigint Bigint::operator + (const Bigint& x) const {
        Bigint res(*this); return res += x;
}
const Bigint Bigint::operator - (const Bigint& x) const {
        Bigint res(*this); return res -= x;
}
const Bigint Bigint::operator * (const Bigint& x) const {
        Bigint res(*this); return res *= x;
}
const Bigint Bigint::operator / (const Bigint& x) const {
        Bigint res(*this); return res /= x;
}
const Bigint Bigint::operator % (const Bigint& x) const {
        Bigint res(*this); return res %= x;
}
std::pair<Bigint, Bigint> divmod(const Bigint& lhs, const Bigint& rhs) {
        if (!rhs.a.size()) {
                std::cerr<<"error: division by zero"<<std::endl;
                exit(1);
        }
        Bigint x(abs(rhs)), q, r;
        for (int i = lhs.a.size() - 1; i >= 0; i --) {
                r = r * Bigint::BASE + lhs.a[i];
                int head = 0, tail = Bigint::BASE;
                if (r >= x) {
                        while (head + 1 < tail) {
                                int mid = (head + tail) / 2;
                                if (x * Bigint(mid) > r) tail = mid;
                                else head = mid;
                        }
                        r -= x * head;
                }
                q.a.push_back(head);
        }
        reverse(q.a.begin(), q.a.end());
        bool neg = lhs.negative ^ lhs.negative;
        q.negative = neg;
        r.negative = neg;
        return std::make_pair(q.normalize(), r.normalize());
}
std::istream& operator >> (std::istream& is, Bigint& x) {
        std::string tmp;
        is >> tmp;
        x = Bigint(tmp);
        return is;
}
std::ostream& operator << (std::ostream& os, const Bigint& x) {
        if (x.negative) os << '-';
        if (!x.a.size()) os << 0;
        else os << x.a.back();
        for (int i = x.a.size()-2; i >= 0; i --) {
                os.width(Bigint::LEN);
                os.fill('0');
                os << x.a[i];
        }
        return os;
}
const Bigint abs(Bigint x) {
        x.negative = false;
        return x;
}
std::string str(Bigint x) {
        stringstream st;
        st << x;
        return st.str();
}

//------------------------------------------------------------------------------------------



bool check(Bigint n, int k) {
        string s = str((Bigint)9 * n + (Bigint)(9 * k));
        int sum = 0;
        for (int i = 0; i < s.size(); i ++) sum += s[i] - '0';
        return sum <= 9 * k;
}

int main() {
        Bigint n;
        cin >> n;
        int lb = 0, ub = 500000;
        while (ub - lb > 1) {
                int mid = (lb + ub) / 2;
                if (check(n, mid)) ub = mid;
                else lb = mid;
        }         
        cout << ub << endl;
        return 0;
}


下の実装はboostの多倍長整数ライブラリによるすっきりした実装。boostはの多倍長演算は遅いっぽく(要検証)、ほとんど$\ TLE\ $してしまった。

実装2
#include "bits/stdc++.h"
 
#include <boost/multiprecision/cpp_int.hpp>
#define Bigint boost::multiprecision::cpp_int
 
bool check(Bigint n, int k) {
        Bigint t = 9 * (n + k);
        std::string s = t.str();
        int sum = 0;
        for (int i = 0; i < s.size(); i ++) sum += s[i] - '0';
        return sum <= 9 * k;
}
 
int main() {
        Bigint n;
        std::cin >> n;
        int lb = 0, ub = 500000;
        while (ub - lb > 1) {
                int mid = (lb + ub) / 2;
                if (check(n, mid)) ub = mid;
                else lb = mid;
        }         
        std::cout << ub << std::endl;
        return 0;
}

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;
}