2019-05-04

いろはちゃんコンテスト Day2: F - 総入れ替え

https://atcoder.jp/contests/iroha2019-day2/tasks/iroha2019_day2_f

問題概要

箱が3つあり、それぞれの箱には100円玉がA1枚、B1枚、C1枚、50円玉がA2枚、B2枚、C2枚入っている。コインが1枚以上入っている箱をひとつ指定して、そこからランダムにコインを1枚取り出すという操作を2人で交互に繰り返すゲームを考える。双方が自分の手に入れる金額の期待値を最大化するように行動するとき、先手が最終的に手に入れる金額の期待値を求めよ。

0 <= A1, A2, B1, B2, C1, C2 <= 10

解法

dp[a1][a2][b1][b2][c1][c2]: 残っているコインの枚数がこれのとき先手が最終的に手に入れる金額の期待値

とすると、先手番の場合たとえば箱Aを選んだとき確率 a / (a+b) で期待値 (dp[a1 -1][a2][b1][b2][c1][c2] + 100) 円となり、確率 b / (a+b) で期待値 (dp[a1][a2 - 1][b1][b2][c1][c2] + 50) 円となる。なのでこのときの期待値は

a / (a+b) * (dp[a1 -1][a2][b1][b2][c1][c2] + 100) + b / (a+b) * (dp[a1][a2 - 1][b1][b2][c1][c2] + 50)

である。他の箱の場合も同様の計算をして最も大きい値のものをdp[a1][a2][b1][b2][c1][c2]として採用すればよい。後手番の場合は期待値の最も低いものをとる（決まった合計金額を2人で取り合うゲームなので、後手の期待値の最大化＝先手の期待値の最小化）。後手のdp計算の場合 +100 とか +50 は入らないことに注意する。あとはdpのインデックスがゼロのときとかをよしなにやる。実装はメモ化再帰のほうがやりやすい（当社比）。

感想

期待値の問題に対してものすごい苦手意識があるが、この問題の場合はそもそもよくあるタイプの2人ゲームなのでそういう感じのメモ化再帰を丁寧にやっていくことを考えるべきだった

コード (C++)

#include <bits/stdc++.h>
using namespace std;
#define REP(i,n) for (int i=0;i<(n);i++)
#define REP2(i,m,n) for (int i=m;i<(n);i++)
typedef long long ll;

int A, B, C, D, E, F, parity;
double mem[11][11][11][11][11][11];

double dfs(int a, int b, int c, int d, int e, int f) {
    if (a == 0 && b == 0 && c == 0 && d == 0 && e == 0 && f == 0) return 0;
    if (a < 0 || b < 0 || c < 0 || d < 0 || e < 0 || f < 0) return 0;
    if (mem[a][b][c][d][e][f] > 0) return mem[a][b][c][d][e][f];
    
    if ((a + b + c + d + e + f) % 2 == parity) {
        double x = (a + b > 0) ? (dfs(a-1, b, c, d, e, f) + 100) * a / (a + b) + (dfs(a, b-1, c, d, e, f) + 50) * b / (a + b) : 0;
        double y = (c + d > 0) ? (dfs(a, b, c-1, d, e, f) + 100) * c / (c + d) + (dfs(a, b, c, d-1, e, f) + 50) * d / (c + d) : 0;
        double z = (e + f > 0) ? (dfs(a, b, c, d, e-1, f) + 100) * e / (e + f) + (dfs(a, b, c, d, e, f-1) + 50) * f / (e + f) : 0;
        return mem[a][b][c][d][e][f] = max({x, y, z});
    } else {
        double x = (a + b > 0) ? dfs(a-1, b, c, d, e, f) * a / (a + b) + dfs(a, b-1, c, d, e, f) * b / (a + b) : -1;
        double y = (c + d > 0) ? dfs(a, b, c-1, d, e, f) * c / (c + d) + dfs(a, b, c, d-1, e, f) * d / (c + d) : -1;
        double z = (e + f > 0) ? dfs(a, b, c, d, e-1, f) * e / (e + f) + dfs(a, b, c, d, e, f-1) * f / (e + f) : -1;
        double ret = 1 << 29;
        if (x != -1) ret = min(ret, x);
        if (y != -1) ret = min(ret, y);
        if (z != -1) ret = min(ret, z);
        return mem[a][b][c][d][e][f] = ret;
    }
}

void solve() {
    cin >> A >> B >> C >> D >> E >> F;
    parity = (A + B + C + D + E + F) % 2;

    cout.precision(20);
    cout << fixed << dfs(A, B, C, D, E, F) << endl;
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    solve();
}

2019-05-03

AtCoder Regular Contest 027: D - ぴょんぴょんトレーニング

問題記録 AtCoder ARC ARC #027 平方分割 DP DP遷移の行列化

https://atcoder.jp/contests/arc027/tasks/arc027_4

問題概要

N個の石が一直線に並んでおり、石x からは 1 ～ Hx 個先のどれかの石にジャンプすることができる。以下の形式のクエリがD個飛んでくるのでそれぞれに解答せよ。

整数s, t (s < t) が与えられる。石sからスタートし、ちょうど石tにたどり着くような移動経路は何通りあるか？

N <= 3 × 10⁵

Hx <= 10

D <= 5000

解法

Hがたかだか10なので、DPの遷移を10×10の行列で表すことができる (dp[11], dp[10], ..., dp[2]) = A ・ (dp[10], dp[9], ..., dp[1]) のような形)。あとはよくあるやつで行列をセグ木とか平方分割とかに乗せればよさそうなのだが、やたら定数倍が重くメモリ制限も小さめなので適当にやるとだめになる (3*10⁵ * 10² * 64 bits で既に240MB)。自分は平方分割でやったが、合成前の個別の行列については陽に持たず必要になるたびいちいち生成するようにした。あとは行列の積を取る方向が微妙にややこしい（大きいインデックスのほうが左に来る）のでがんばる。あと公式解説にあるとおり各クエリに答えるときは行列同士をかけるのではなくベクトルにかけていくようにして時間計算量を抑えた。

感想

虚無みがあったが、DPの遷移を行列でやるやつについての理解は深まったのでよかった

コード (C++)

#include <bits/stdc++.h>
using namespace std;
#define REP(i,n) for (int i=0;i<(n);i++)
#define REP2(i,m,n) for (int i=m;i<(n);i++)
typedef long long ll;

const ll MOD = 1000000007;
const int SEG = 600;
const int M = 10;

int N, Q, s, t;
int H[303030];
ll bucket[SEG][M][M];
ll tmp[M][M];
ll tmp2[M][M];
ll vec[M];
ll vec_tmp[M];

void init_tmp() {
    REP(i, M) tmp[i][i] = 1;
}

void init_vec() {
    vec[0] = 1;
    REP2(i, 1, M) vec[i] = 0;
}

void gen(int i) {
    REP(j, M) REP(k, M) tmp[j][k] = 0;
    REP2(j, 1, M) tmp[j][j-1] = 1;
    REP2(j, 1, M+1) if (i - j >= 0 && H[i-j] >= j) tmp[0][j-1] = 1;
}

void prod1(int i) {
    gen(i);
    REP(j, M) vec_tmp[j] = vec[j];
    REP(j, M) vec[j] = 0;
    REP(j, M) REP(k, M) (vec[j] += tmp[j][k] * vec_tmp[k] % MOD) %= MOD;
}

void prod2(int seg) {
    REP(i, M) vec_tmp[i] = vec[i];
    REP(i, M) vec[i] = 0;
    REP(j, M) REP(k, M) (vec[j] += bucket[seg][j][k] * vec_tmp[k] % MOD) %= MOD;
}

void solve() {
    cin >> N;
    REP(i, N) cin >> H[i];
    cin >> Q;

    REP(i, SEG) REP(j, M) bucket[i][j][j] = 1;
    for (int i = 0; i < N; ++i) {
        gen(i);
        REP(x, M) REP(y, M) tmp2[x][y] = bucket[i/SEG][x][y];
        REP(x, M) REP(y, M) bucket[i/SEG][x][y] = 0;
        REP(x, M) REP(y, M) REP(z, M) (bucket[i/SEG][x][y] += tmp[x][z] * tmp2[z][y]) %= MOD;
    }

    while (Q--) {
        cin >> s >> t;
        --t;
        init_vec();
        for (int i = s / SEG * SEG; i <= t / SEG * SEG; i += SEG) {
            if (i < s || t < i + SEG - 1) {
                for (int j = max(i, s); j <= min(t, i + SEG - 1); ++j) {
                    prod1(j);
                }
            } else {
                prod2(i / SEG);
            }
        }
        cout << vec[0] << "\n";
    }
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    solve();
}

2019-05-03

AtCoder Regular Contest 045: D - みんな仲良し高橋君

AtCoder ARC #045 関節点問題記録 ARC

https://atcoder.jp/contests/arc045/tasks/arc045_d

問題概要

整数Nと、XY平面上の(2N+1)個の点が与えられる。各点はX座標またはY座標が同じ他の1点とペアになることができる。ただしすでにペアをつくっている点が他の点とペアになることはできない。すべての点について「その点だけを除いたとき残された点でN個のペアを作ることができるか」を判定せよ。

N <= 10⁵

解法

まずペアになることが可能な点同士で辺を張ってグラフをつくる。ここで異なる連結成分に属する点同士はペアにできないので、各連結成分は独立に考えてよい。

このグラフの重要な性質として「連結成分の頂点数が偶数であるならば、その連結成分内ではすべての点をペアにすることが可能」というものがある（詳しくは公式解説参照）。よって問題は以下のように言い換えることが可能である。

グラフからある点を除いたとき、すべての連結成分の頂点数が偶数となるか？

自明なパターンから処理していくと、まず頂点数が奇数であるような連結成分が2個以上ある場合、題意を満たすようなペア作りは不可能である（どの点を除いたとしても必ずひとつは奇数の連結成分が残ってしまうため）。

そうでない場合、つまり頂点数が奇数であるような連結成分がただひとつである場合、まず除く点の属する連結成分の頂点数が偶数のときはNGになる（上と同じ理由）。

残るパターンは頂点数が奇数であるような連結成分がただひとつで、かつ除く点の属する連結成分の頂点数が奇数の場合である。このとき、まず除く点が関節点でなければ無条件でOKとなる。なぜならその点を除いても連結成分はバラバラにならないので、結果として残る連結成分はすべて偶数となるからである。逆に関節点の場合、その点を取り除くと連結成分が分裂する。なので分かれた先の各連結成分が偶数であるかを確認する必要がある。この問題でつくったグラフの場合、関節点を取り除いた結果として連結成分が3個以上に分かれることはない（必ず2つになる）ので、関節点検出のときにつくるDFS木上で部分木の個数をかぞえるようなことをすれば分裂後の頂点数の偶奇はわかる（具体的にはある関節点uを取り除いたとき、その連結成分は「DFS木上において関節点uの子であってlow[v]>=ord[u]であるような頂点vを根とする部分木に含まれるすべての点の集合」と「それとu以外のすべての点の集合」の2つの連結成分に分かれる。2つなのでどちらかの偶奇だけがわかればよく、DFS木上では部分木のサイズを数えるほうが簡単）。

感想

なんかややこしくてつらかった　関節点とか橋って問題によって求められる操作が微妙に違っててうまいことライブラリにするのむずない？

コード (D言語)

import std.stdio, std.array, std.string, std.conv, std.algorithm;
import std.typecons, std.range, std.random, std.math, std.container;
import std.numeric, std.bigint, core.bitop;

void main() {
    auto N = readln.chomp.to!int;
    auto P = (2*N+1).iota.map!(i => readln.split.map!(to!int).array).array;

    auto X = new Tuple!(int, int)[][](2*N+1);
    auto Y = new Tuple!(int, int)[][](2*N+1);
    foreach (i; 0..2*N+1) {
        X[P[i][0]-1] ~= tuple(P[i][1]-1, i);
        Y[P[i][1]-1] ~= tuple(P[i][0]-1, i);
    }

    auto G = construct_graph(N, X, Y);
    auto U = construct_uf(N, G);

    if ((2*N+1).iota.map!(i => U.table[i] < 0 && -U.table[i] % 2 == 1).sum != 1) {
        foreach (i; 0..2*N+1) writeln("NG");
        return;
    }

    G.articulation_points;

    debug {
        G.adj.each!writeln;
        G.is_articulation.writeln;
        G.dfs_subtree.writeln;
    }

    foreach (i; 0..2*N+1) {
        if (-U.table[U.find(i)] % 2 == 0) {
            writeln("NG");
        } else if (!G.is_articulation[i]) {
            writeln("OK");
        } else {
            writeln(G.divide_evenly[i] ? "OK" : "NG");
        }
    }
}

class UndirectedGraph {
    int N;
    int[][] adj;
    int[] dfs_subtree;
    bool[] is_articulation;
    bool[] divide_evenly;

    this (int N) {
        this.N = N;
        adj = new int[][](N);
    }

    void add_edge(int u, int v) {
        adj[u] ~= v;
        adj[v] ~= u;
    }

    void articulation_points() {
        auto ord = new int[](N);
        auto low = new int[](N);
        auto dfs_tree = new int[][](N);
        auto used = new bool[](N);
        dfs_subtree = new int[](N);
        is_articulation = new bool[](N);
        divide_evenly = new bool[](N);
        int cnt;

        int dfs(int n, int p) {
            ord[n] = cnt++;
            low[n] = ord[n];
            used[n] = true;
            foreach (m; adj[n]) {
                if (m == p) continue;
                if (used[m]) low[n] = min(low[n], ord[m]);
                else low[n] = min(low[n], dfs(m, n)), dfs_tree[n] ~= m; 
            }
            return low[n];
        }

        int dfs2(int n, int p) {
            dfs_subtree[n] = 1;
            foreach (m; dfs_tree[n]) {
                if (m == p) continue;
                auto s = dfs2(m, n);
                dfs_subtree[n] += s;
            }
            return dfs_subtree[n];
        }

        foreach (i; 0..N) {
            cnt = 0;
            if (!used[i]) {
                dfs(i, -1);
            }
        }

        foreach (i; 0..N) {
            if (ord[i] == 0 && dfs_tree[i].length >= 2) {
                is_articulation[i] = true;
            } else if (ord[i] != 0 && dfs_tree[i].map!(j => (low[j] >= ord[i])).any) {
                is_articulation[i] = true;
            }
        }

        foreach (i; 0..N) {
            if (ord[i] == 0) {
                dfs2(i, -1);
            }
        }

        foreach (i; 0..N) {
            if (!is_articulation[i]) continue;
            if (ord[i] == 0) {
                divide_evenly[i] = dfs_subtree[dfs_tree[i][0]] % 2 == 0;
            } else {
                divide_evenly[i] = dfs_tree[i].filter!(j => low[j] >= ord[i]).map!(j => dfs_subtree[j]).sum % 2 == 0;
            }
        }
    }
}

class UnionFind {
    int N;
    int[] table;
    int[][] group;

    this(int n) {
        N = n;
        table = new int[](N);
        fill(table, -1);
        group = N.iota.map!(i => [i]).array;
    }

    int find(int x) {
        return table[x] < 0 ? x : (table[x] = find(table[x]));
    }

    void unite(int x, int y) {
        x = find(x);
        y = find(y);
        if (x == y) return;
        if (table[x] > table[y]) swap(x, y);
        group[x] ~= group[y];
        group[y] = [];
        table[x] += table[y];
        table[y] = x;
    }
}


UndirectedGraph construct_graph(int N, Tuple!(int, int)[][] X, Tuple!(int, int)[][] Y) {
    auto G = new UndirectedGraph(2*N+1);
    foreach (i; 0..2*N+1) {
        X[i].sort();
        foreach (j; 0..X[i].length.to!int-1) {
            G.add_edge(X[i][j][1], X[i][j+1][1]);
            if (j + 2 < X[i].length) {
                G.add_edge(X[i][j][1], X[i][j+2][1]);
            }
        }
    }
    foreach (i; 0..2*N+1) {
        Y[i].sort();
        foreach (j; 0..Y[i].length.to!int-1) {
            G.add_edge(Y[i][j][1], Y[i][j+1][1]);
            if (j + 2 < Y[i].length) {
                G.add_edge(Y[i][j][1], Y[i][j+2][1]);
            }
        }
    }
    foreach (i; 0..2*N+1) {
        G.adj[i] = G.adj[i].dup.sort().uniq.array;
    }
    return G;
}

UnionFind construct_uf(int N, UndirectedGraph G) {
    auto uf = new UnionFind(2*N+1);
    foreach (i; 0..2*N+1) foreach (j; G.adj[i]) uf.unite(i, j);
    return uf;
}

2019-02-25

Aizu Online Judge 2710: Marked Ancestor

問題記録 AOJ AOJ-ICPC 木 UnionFind 逆順の操作

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2170

問題概要

N頂点の木とQ個のクエリが与えられる。木は常に頂点1が根である。また初めは頂点1だけにマークがされている。クエリにはMとQの2種類があり、Mクエリでは指定された頂点にマークを付ける。Qクエリでは指定された頂点の祖先でマークが付いているもののうちその頂点から最も近い頂点を求める（祖先には自分も含む）。Qクエリで求めた頂点番号の合計を出力せよ。

N, Q <= 10⁵

解法

マークを付ける操作を考えると、これは「マークを付けた頂点を根とする部分木を木から切り離す」操作ともいえる。切り離された部分木に含まれるノードたちにとっては、その時点での最も近いマークされた祖先はその部分木の根である。

このように見るとこれはもともとひとつだった集まりをクエリごとに分断していき、途中途中でその集まりの親を聞かれている問題と捉えることができる。このような操作を要求される問題に対する典型テクとして「操作を逆順で見ていく」というのがあり、この問題にもそれを使うことができる。

なぜ逆順で見るのがいいかというと、「集合を切って分割していく操作」というのは扱いづらいが、これを逆から見た「分割されている集合たちを繋いでいく操作」はかなり扱いやすい性質を持っているからである。具体的にはUnionFindが使える。

というわけでクエリを逆回しに見ていくことにしたいので、まずクエリを先に読んですべてのマークがついたことにしてしまう。そしてこの時点でマークがついていないノードはすべてそのノードの親とuniteする。このようにしてできる集合たちが上述の「部分木を切り離していく操作」を最後まで行ったときのそれぞれの部分木に対応する。あとはここからまたクエリを逆順に見ていく。

逆順に見ていってQクエリが飛んできた場合、指定された頂点が属する集合の根に当たる頂点、すなわち最も深さが小さい頂点を答えればよい。これはUnionFindを適当に改造してそのような情報を持たせることにすれば実装できる。

Mクエリが飛んできた場合、その頂点からマークを剥がすことになる。これはつまり部分木を元の木にくっつける作業で、マークを剥がす頂点とその直接の親をuniteすればよい。ただしその頂点に何度もマークが行われているときは注意が必要で、この場合いちばん最初につけられたマーク以外には意味がないため、他は無視しないといけない。

感想

AOJ-ICPC軽く埋めたろかと思って見たけど普通にむずくて解けなかった

木の頂点をUFで管理するという発想がなかったのは反省

コード(D言語)

import std.stdio, std.array, std.string, std.conv, std.algorithm;
import std.typecons, std.range, std.random, std.math, std.container;
import std.numeric, std.bigint, core.bitop, core.stdc.string;

bool solve() {
    auto s = readln.split.map!(to!int);
    auto N = s[0];
    auto Q = s[1];
    if (N == 0) return false;

    auto G = new int[][](N);
    auto D = new int[](N);
    auto P = new int[](N);
    P[0] = -1;

    foreach (i; 1..N) {
        auto n = readln.chomp.to!int;
        P[i] = n-1;
        G[i] ~= n-1;
        G[n-1] ~= i;
    }

    void dfs(int n, int p, int d) {
        D[n] = d;
        foreach (m; G[n]) if (m != p) dfs(m, n, d+1);
    }

    dfs(0, -1, 0);

    Tuple!(int, int)[] queries;
    auto marked = new int[](N);

    foreach (i; 1..Q+1) {
        auto t = readln.split;
        auto n = t[1].to!int - 1;
        if (t[0] == "M" && !marked[n]) {
            queries ~= tuple(0, n);
            marked[n] = i;
        } else if (t[0] == "Q") {
            queries ~= tuple(1, n);
        }
    }

    auto uf = new UnionFind(N, D);
    foreach (i; 0..N) if (!marked[i]) uf.unite(i, P[i]);
    long ans = 0;

    foreach_reverse (q; queries) {
        if (q[0] == 0) {
            uf.unite(q[1], P[q[1]]);
        } else {
            ans += uf.find(q[1]) + 1;
        }
    }

    ans.writeln;

    return true;
}

void main() {
    while (solve()) {}
}


class UnionFind {
    int N;
    int[] table;
    int[] depth;

    this(int n, int[] d) {
        N = n;
        table = new int[](N);
        fill(table, -1);
        depth = d.dup;
    }

    int find(int x) {
        return table[x] < 0 ? x : (table[x] = find(table[x]));
    }

    void unite(int x, int y) {
        if (x == -1 || y == -1) return;
        x = find(x);
        y = find(y);
        if (x == y) return;
        if (depth[x] > depth[y]) swap(x, y);
        table[x] += table[y];
        table[y] = x;
    }
}

2019-02-14

AtCoder Regular Contest 084: F - XorShift

AtCoder ARC ARC #084 GCD 多項式問題記録

https://atcoder.jp/contests/arc084/tasks/arc084_d

問題概要

N個の非負整数Aiが黒板に書かれている。黒板に書いてある数に対して「数をひとつ選びそれを2倍したものを黒板に書き足す」「数をふたつ選びそれらのbitwise xorをとったものを黒板に書き足す」という操作を好きな順番で何度でも行えるとき、黒板に書かれうる数のうち大きさがX以下であるものが何種類存在するか求めよ。

1 <= N <= 6

1 <= X, Ai <= 2⁴⁰⁰⁰

解法

なぜそうするかはとりあえずおいといて、与えられた数の2進表記を多項式の係数とみなすことにする。 1111 は x³ + x² + x¹ + 1 となり、 101 は x² + 1 となる、いう具合である。このようにすると「ある数を2倍する操作」は「多項式にxをかける操作」とみなすことができるようになる。また「2つの数のxorをとる操作」は「多項式同士の足し算をする操作（ただし係数はmod2で考える）」になる。

ではこのような操作を繰り返していって作ることのできる多項式はどのようなものになるか？実はこれはGCDのような概念を考えることで解決することができる。どういうことかというと、「与えられたN個の多項式に上の操作を繰り返して作れる多項式の集合」を考えたとき、実はそれとまったく同じ集合を「ひとつの多項式だけが含まれる集合」から始めても作ることができる。その「ひとつの多項式」をここではGCDと呼んでいる。（普通の整数のGCDも、上の操作を「整数倍」と「足し算・引き算」に置き換えると同じような性質がある（あるよね？）（なかったらすいません））

GCDの求め方だが、互除法と同じ形でやれる。すなわち次数が大きい方の多項式をA, 次数の小さい方（同じでも良い）をBとして、最大次数が揃うまでBにxをかけてから足す。その結果をCとして再帰的に同じことをし続ける、という形である。これでGCDとなる多項式を求めることができる。

GCDとなる多項式Gが求まったら、あとはそれを使ってどんな数が作れるかを考えればよい。ここでGの最大次数をdとすると、実はd次より上の項はどんな組み合わせでも自由につくることができる。そして(d-1)次より下の項はそれに応じて自動的に一意に決まる。なぜなら x倍の操作でGを左にシフトしていくことはできるが右にシフトさせることはできないからである。もしd次以上の項の係数をフリップさせたいのであればGをそこまで左シフトさせてから足せばいいが、(d-1)次以下に対してはそれができないので自由にいじることはできない。

よって答えとしては、基本的にXのd次以上の項は全部自由に作れるので、そこの部分だけを切り取って2進数と見なした数は全部作ることができる。ただしひとつだけ例外があって、作った結果d次以上の全部の項がXと一致している場合、後ろの桁の結果によってはXをオーバーしている可能性がある。なのでそのケースだけは実際にGから構成してみてXを超えないか確かめる必要がある。もし超えていたらその数は作れないので除く。あとは0も数える対象であるのでそれを忘れず含めるようにする。

感想

わからんが、とにかく書いた

GCDが「集合の要素hoge倍する・または集合の要素同士を足し算（引き算）して集合に加えていく」という操作をするときの最小要素？的な感じになるというのはもっとちゃんと理解しておくべきことっぽい

コード (D言語)

import std.stdio, std.array, std.string, std.conv, std.algorithm;
import std.typecons, std.range, std.random, std.math, std.container;
import std.numeric, std.bigint, core.bitop, core.stdc.string;

immutable long MOD = 998244353;

int[] gcd(int[] a, int[] b) {
    while (true) {
        if (a.length < b.length) swap(a, b);
        if (b.length == 0) return a.dup;
        ulong msb = a.length;

        foreach (i; 0..a.length) {
            if (i < b.length) {
                a[i] ^= b[i];
            }
            if (a[i]) {
                msb = min(i, msb);
            }
        }

        a = a[msb..$];
    }
}

void main() {
    auto s = readln.split;
    auto N = s[0].to!int;
    auto X = s[1].map!(a => to!int(a == '1')).array;
    auto P = N.iota.map!(_ => readln.chomp.map!(a => to!int(a == '1')).array).array;

    auto G = P.front.dup;
    foreach (i; 1..N) G = gcd(G, P[i].dup);

    long ans = 0;

    for (long i = X.length.to!long - G.length.to!long, tmp = 1; i >= 0; i -= 1, tmp = tmp * 2 % MOD) {
        if (X[i]) {
            ans = (ans + tmp) % MOD;
        }
    }


    auto M = new int[](X.length);

    for (int i = 0; i + G.length <= M.length; ++i) {
        if (M[i] != X[i]) {
            for (int j = 0; j < G.length; ++j) {
                M[i+j] ^= G[j];
            }
        }
    }

    if (M <= X) {
        ans = (ans + 1)  % MOD;
    }

    ans.writeln;
}

2019-02-07

DISCO presents ディスカバリーチャンネルコードコンテスト2019 予選: D - チップ・ストーリー　～黄金編～

問題記録 AtCoder 中国剰余定理

https://atcoder.jp/contests/ddcc2019-qual/tasks/ddcc2018_qual_d

問題概要

1以上10¹²以下の秘密の整数Nがある。Nをi進数で表したときの各桁の数字の和Ai (i = 2〜30) がそれぞれ与えられるので、条件を満たすNが存在するか判定せよ。存在する場合にはそのようなNをひとつ示せ。

解法

まず重要な事実として「整数nを(k-1)で割ったあまりと、nをk進数で表したときの桁和を(k-1)で割ったあまりは等しい」というものがある。つまり N = Ai (mod (i-1)) である。

このように見ると、「いろんな進数でのNの桁和」という情報は「いろんな数でNを割ったときのあまり」の情報として使うことができる。そしてこれは中国剰余定理と非常に似た形をしている。

中国剰余定理: あるひとつの整数nを互いに素ないくつかの数 m1, m2, ..., mk で割ったときのあまりがわかっているとき、nをm1×m2× ... ×mkで割ったときのあまりは一意に定まる。

この定理は解があるよ〜という話しかしていないが、実際に復元する手順についてもいろいろやり方が知られている（人のコードを見る限り競プロでは拡張ユークリッドの互除法を用いて実装するのが一般的っぽい（が、自分はそれをよく理解してないのでググって出てきた手計算での復元方法っぽいやつを実装した http://did2.blog64.fc2.com/blog-entry-86.html））。

さてこの定理を使うためにはmodがすべて互いに素でなければならない。というわけでi=2〜30の中からi-1が素数のものだけを抜き出してきて計算する。すると N の mod (2×3×5×...×29) での値を出すことができる(modは1〜29の中のすべての素数の積)。2×3×5×...×29 がだいたい 6.5*10⁹ くらいで、Nは10¹²以下という制約があるので、あとはこのmodの条件を満たす数を総当たりして、逐一条件を満たすかをチェックしていけばよい。

感想

中国剰余定理よく知らなくて解説読んだ後放置してたけどこの問題の範囲ではそんなに難しいことが要求されてるわけではなかった　整数論べんきょうしたいですね

コード(D言語)

import std.stdio, std.array, std.string, std.conv, std.algorithm;
import std.typecons, std.range, std.random, std.math, std.container;
import std.numeric, std.bigint, core.bitop, core.stdc.string;

void main() {
    auto A = 29.iota.map!(_ => readln.chomp.to!long).array;
    A = 0 ~ A;

    const long[] P = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29];
    const long PP = P.reduce!"a * b";
    long ans = 0;

    foreach (p; P) {
        // x * (PP / p) = A[p] (mod p)
        long x = A[p] * powmod((PP / p % p), p - 2, p) % p;
        ans += x * (PP / p) % PP;
        ans %= PP;
    }

    while (ans <= 10L^^12) {
        if (iota(1, 30).map!(i => digit_sum(ans, i+1) == A[i]).all) {
            ans.writeln;
            return;
        }
        ans += PP;
    }

    writeln("invalid");
}

long digit_sum(long x, long base) {
    long ret = 0;
    while (x > 0) {
        ret += x % base;
        x /= base;
    }
    return ret;
}

long powmod(long a, long x, long m) {
    long ret = 1;
    while (x) {
        if (x % 2) ret = ret * a % m;
        a = a * a % m;
        x /= 2;
    }
    return ret;
}

2019-02-02

AtCoder Grand Contest 025: D - Choosing Points

問題記録 AtCoder AGC AGC #025 2部グラフ格子点独立集合

https://atcoder.jp/contests/agc025/tasks/agc025_d

問題概要

正整数N, D1, D2が与えられる。xy平面上で 0 <= x < 2N, 0 <= y < 2N に存在する整数座標の点の集合であって、要素数がN²個でありかつ集合の中のどの2点を取ってもその2点の距離がsqrt(D1)でもsqrt(D2)でもないという条件を満たす集合をひとつ構成せよ。このような集合は必ず存在する。

N <= 300

D1, D2 <= 2*10⁵

解法

2次元平面上で距離が定数Dであるような整数座標の点を繋いでできるグラフは必ず2部グラフになる。sqrt(D1), sqrt(D2)でそれぞれこのような2部グラフを作った上でそれぞれの点を2色で塗り分けると、すべての点は以下の4種類のいずれかひとつに分類されることになる。

D1のグラフで色1, D2のグラフで色1
D1のグラフで色1, D2のグラフで色2
D1のグラフで色2, D2のグラフで色1
D1のグラフで色2, D2のグラフで色2

作ったグラフの性質上、これら4つの集合においては集合内のどの相異なる2点を持ってきてもそれらの距離がsqrt(D1), sqrt(D2)になることはない（その距離になるためには元のグラフで隣接している必要があるが、2部グラフで同色になるならば隣接ではありえない）。

また鳩の巣原理よりこれら4つの分類のうち必ずひとつ以上には (全体の点の数 / 4) 個以上の点が属することになる。この問題で考えるべき格子点は4N²個なので、これらのうち必ずひとつはN²個以上になる。よってそれを答えとすればよい。

なおグラフを作るのを愚直にやるとO(N⁴)になってしまうが、最初に「原点からの距離がsqrt(D1), sqrt(D2)となる整数座標の点」を前計算で列挙しておけば、これは各x座標においてたかだか2点しか存在しないので、各点においてO(N)で距離sqrt(D1), sqrt(D2)の点が列挙できるためO(N³)で済む。

作ったグラフが必ず2部グラフになることの証明は公式解説に詳しい。

感想

わかりません……　格子点をある定数の距離で結ぶと2部グラフになる　は覚えとくとべんりそう

あと参考にした解説記事 (AtCoder Grand Contest 025 - Unhappy Go Lucky!) で印象に残った文章があるので教訓として引用しときます

まずは与えられたグラフの特徴を考察する。独立集合は、二部グラフであれば楽勝で取れる。これが「独立集合」と聞いた瞬間に「二部グラフでは」と行かないと辛い。 AtCoder Grand Contest 025 - Unhappy Go Lucky!

コード (D言語)

import std.stdio, std.array, std.string, std.conv, std.algorithm;
import std.typecons, std.range, std.random, std.math, std.container;
import std.numeric, std.bigint, core.bitop, core.stdc.string;

void main() {
    auto s = readln.split.map!(to!int);
    auto N = s[0];
    auto D = [s[1], s[2]];
    auto G = new int[][][](2, 4*N*N);
    auto C = new int[][](2, 4*N*N);
    C[0][] = C[1][] = -1;

    auto xyd = new int[][](2, 2*N);

    foreach (i; 0..2) {
        xyd[i][] = -1;
        foreach (x; 0..2*N) {
            foreach (y; 0..2*N) {
                if (x*x+y*y == D[i]) {
                    xyd[i][x] = y;
                }
            }
        }
    }


    foreach (i; 0..2) {
        foreach (x1; 0..2*N) {
            foreach (y1; 0..2*N) {
                foreach (x2; x1..2*N) {
                    if (xyd[i][x2-x1] == -1) continue;
                    foreach (k; [-1, 1]) {
                        int y2 = y1 + k * xyd[i][x2-x1];
                        if (y2 < 0 || y2 >= 2*N) continue;
                        G[i][x1*2*N+y1] ~= x2*2*N+y2;
                        G[i][x2*2*N+y2] ~= x1*2*N+y1;
                    }
                }
            }
        }
    }

    foreach (i; 0..2) {
        foreach (x1; 0..2*N) {
            foreach (y1; 0..2*N) {
                if (C[i][x1*2*N+y1] == -1) {
                    dfs(x1*2*N+y1, 0, G[i], C[i]);
                }
            }
        }
    }

    int cnt = 0;

    foreach (i; 0..2) {
        foreach (j; 0..2) {
            if (iota(4*N*N).map!(k => C[i][k] == i && C[j][k] == j).sum < N*N) continue;
            foreach (x; 0..2*N) {
                foreach (y; 0..2*N) {
                    if (C[0][x*2*N+y] == i && C[1][x*2*N+y] == j) {
                        writeln(x, " ", y);
                        ++cnt;
                        if (cnt >= N * N) {
                            return;
                        }
                    }
                }
            }
        }
    }
}

void dfs(int n, int c, ref int[][] G, ref int[] C) {
    if (C[n] != -1) return;
    C[n] = c;
    foreach (m; G[n]) if (C[m] == -1) dfs(m, c^1, G, C);
}