第3回 ドワンゴからの挑戦状 本選: B - ニワンゴくんの約数

AtCoder ドワコン ドワコン3 本戦 約数 Mo's Algorithm クエリ平方分割

https://atcoder.jp/contests/dwacon2017-honsen/tasks/dwango2017final_b

問題概要

N要素の正整数の数列Xが与えられる。Q個の区間[L, R]が与えられるので、それぞれについて「その区間に含まれるすべてのXの要素の積をとったとき、その数の約数は何個か」という問いに答えよ。

N, Q <= 105

max(Xi) <= 105

解法

約数の個数は各素因数について（出現回数＋１）の積をとることで求められる。なので問題で要求されている操作そのものはシンプルで、愚直にやるならばXの要素を素因数分解したうえでそれぞれの区間において各素数の出現回数の和を取ってかけ合わせればよい。

そう考えると各素数の出現回数テーブルを持って累積和をすればよいのではと一瞬思えるが、105までの素数は104個くらいあるのでそれでは計算量的に厳しい。

ということで別の道を探してmax(Xi)の小ささに注目してみる。max(Xi)がたかだか105なので、あるひとつのXiにおいてsqrt(105)より大きい素数はたかだか1回しか現れない（もし2回以上現れたら明らかに105を越えてしまう）。また105以下の素数は104個くらいあるがsqrt(105)以下の素数は70個くらいしかない。まとめると以下のような性質の違いがsqrt(105)を境に現れるといえる。

  1. sqrt(105)以下の素因数 → 種類は少ないがひとつのXiにおいて複数種・複数回現れる可能性がある
  2. sqrt(105)より大きい素因数 → 種類は多いがひとつのXiにおいてたかだか一種類・一回までしか現れない

以上の性質の違いからこれらを場合分けして独立に解くことにするとうまくいく。

sqrt(105)以下の素因数に対して

上述のとおりこのような素因数は70個くらいしかないので単純な累積和（全部の素数の出現回数をカウントするテーブルを作って累積和）を適用できる。

sqrt(105)より大きい素因数に対して

ひとつのXiにおいてたかだか一種類・一回までしか現れないという性質から、もしいまある区間の答えがわかっているならば、その区間の端をひとつだけずらしたときの答えというのが高速に計算できる（たとえば[L, R]の右端をひとつ伸ばして[L, R+1]を見たとき、X_(R+1)の「大きい素因数」の出現回数が1増える。もしもともとの出現回数が3回だったとするとこれによって答えは「元の答え / 3 × 4」に更新できる）。

このように区間をひとつずらしたときの答えの差分更新が高速にできるとき、Mo's algorithm が使える可能性がある。このアルゴリズムの詳細については以下の記事がとてもわかりやすかった。

具体的な話はこの記事を読んでもらうとして自分の理解をざっくりまとめておくと、Mo's algorithmは「区間クエリのいい感じの並び替え方」についてのアルゴリズムである。具体的には「区間の右端か左端をひとつずつずらして差分更新していって最終的にクエリで飛んでくる全部の区間の答えを求めたい、というとき、Mo's algorithm にしたがって区間を並べ替えるとずらしの回数が(N+Q)sqrt(N)に抑えられる」というもの（Nは数列の長さ、Qはクエリ数）。

このケースはまさにMoに適した形の構造をしている。各素因数が何回出てきたかカウントする用の105のテーブルひとつと現在の答えを保持する変数ひとつを持ったうえでMoが返してくるとおりの順番で更新を行っていけばよい。

感想

うしさんありがとう

コード (C++)

Mo の部分は上の記事のものをそのまま貼っただけ（うしさんありがとう）

#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 SQ = 350;
int N, Q;
int A[101010];
int B[101010][100];
int L[101010];
int R[101010];
ll ans[101010];
vector<int> P;
ll tmp = 1;
ll cnt[101010];
ll inv[101010];

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

// from https://ei1333.hateblo.jp/entry/2017/09/11/211011
struct Mo
{
    vector< int > left, right, order;
    vector< bool > v;
    int width;
    int nl, nr, ptr;

    Mo(int n) : width((int) sqrt(n)), nl(0), nr(0), ptr(0), v(n) {}

    void insert(int l, int r) /* [l, r) */
    {
        left.push_back(l);
        right.push_back(r);
    }

    /* ソート */
    void build()
    {
        order.resize(left.size());
        iota(begin(order), end(order), 0);
        sort(begin(order), end(order), [&](int a, int b)
        {
            if(left[a] / width != left[b] / width) return left[a] < left[b];
            return right[a] < right[b];
        });
    }

    /* クエリを 1 つぶんすすめて, クエリのidを返す */
    int process()
    {
        if(ptr == order.size()) return (-1);
        const auto id = order[ptr];
        while(nl > left[id]) distribute(--nl);
        while(nr < right[id]) distribute(nr++);
        while(nl < left[id]) distribute(nl++);
        while(nr > right[id]) distribute(--nr);
        return (order[ptr++]);
    }

    inline void distribute(int idx)
    {
        v[idx].flip();
        if(v[idx]) add(idx);
        else del(idx);
    }

    void add(int idx) {
        if (A[idx] == 1) return;
        tmp = tmp * inv[cnt[A[idx]] + 1] % MOD;
        cnt[A[idx]] += 1;
        tmp = tmp * (cnt[A[idx]] + 1) % MOD;
    }

    void del(int idx) {
        if (A[idx] == 1) return;
        tmp = tmp * inv[cnt[A[idx]] + 1] % MOD;
        cnt[A[idx]] -= 1;
        tmp = tmp * (cnt[A[idx]] + 1) % MOD;
    }
};

void enum_small_primes() {
    for (int i = 2; i < SQ; ++i) {
        bool ok = true;
        for (auto p: P) if (i % p == 0) {
            ok = false;
            break;
        }
        if (ok) {
            P.push_back(i);
        }
    }
}

void solve_small() {
    REP(i, N) REP(j, (int)P.size()) while (A[i] % P[j] == 0) {
        A[i] /= P[j];
        B[i+1][j] += 1;
    }
    REP(i, N) REP(j, (int)P.size()) {
        B[i+1][j] += B[i][j];
    }
    REP(i, Q) REP(j, (int)P.size()) {
        ans[i] = ans[i] * (B[R[i]+1][j] - B[L[i]][j] + 1) % MOD;
    }
}

void solve_large() {
    Mo mo = Mo(N);
    REP(i, Q) {
        mo.insert(L[i], R[i] + 1);
    }
    mo.build();
    REP(i, Q) {
        int idx = mo.process();
        ans[idx] = ans[idx] * tmp % MOD;
    }
}

void solve() {
    cin >> N >> Q;
    REP(i, N) cin >> A[i];
    REP(i, Q) cin >> L[i] >> R[i], --L[i], --R[i];
    REP(i, Q) ans[i] = 1;
    REP(i, 101010) inv[i] = powmod(i, MOD-2);

    enum_small_primes();
    solve_small();
    solve_large();

    REP(i, Q) cout << ans[i] << "\n";
}

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

早稲田大学プログラミングコンテスト2019: F - RPG

問題記録 フロー WUPC2019

https://atcoder.jp/contests/wupc2019/tasks/wupc2019_f

問題概要

N頂点M辺の有向無閉路グラフが与えられる。グラフの各頂点は必ず「戦闘」もしくは「空き地」のいずれかに分類される。また空き地の頂点iには整数のコストCiを支払うことで回復所を設置することができる。頂点1から頂点Nまで移動するとき、どのような経路を通ったとしても2つの戦闘頂点間を移動する間には必ずひとつ以上の回復所を通過するように回復所を設置したい。設置に必要なコスト合計の最小値を求めよ。

N <= 500

M <= 1000

解法

フローをやる。辺の張り方はいろいろあると思うが、自分は次のようにやった。

  • 頂点を倍加して各頂点に「入頂点」「出頂点」をつくる
  • 戦闘頂点ではsourceから出頂点へ、入頂点からsinkへそれぞれ容量INFで辺を張る
  • 空き地頂点では自分の入頂点から自分の出頂点に対して容量Ciで辺を張る
  • あとは与えられたグラフどおりに容量INFで出頂点から入頂点への辺を張っていく

あとはこのグラフにフローを流して最大流を求めれば答えが出る。この問題で求められているような回復所の設置の仕方はすべての戦闘頂点の出頂点がsource側、入頂点がsink側になるようグラフをカットするものと捉えることができて、INFの辺は実質カットできないものと見なすとあとは空き地頂点の容量Ciの辺を切る（＝そこに回復所を設置する）しかないので、結果最小カット（＝最大流）が答えとなる、という感じ（たぶん）

感想

カットで頂点をsource側sink側の2種類に分けるという考え方重要っぽい　フローよくわかってないのでなんか変なこと言ってたら教えてください

こいつ最近F問題の記事ばっか書いてるな

コード (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;
 
immutable long MOD = 10^^9 + 7;
immutable long INF = 10L^^15;

void main() {
    auto s = readln.split.map!(to!int);
    auto N = s[0];
    auto M = s[1];
    auto C = 0 ~ readln.split.map!(to!(long)).array ~ 0;
    auto G = new int[][](N);
    foreach (i; 0..M) {
        s = readln.split.map!(to!int);
        G[s[0]-1] ~= s[1]-1;
    }
    
    int source = 0;
    int sink = N - 1;
    auto ff = new FordFulkerson!long(N*2, source, sink);

    foreach (i; 1..N-1) {
        if (C[i] == -1) {
            ff.add_edge(source, i, INF);
            ff.add_edge(i+N, sink, INF);
        } else {
            ff.add_edge(i+N, i, C[i]);
        }

        foreach (j; G[i]) {
            ff.add_edge(i, j+N, INF);
        }
    }

    long ans = ff.run;
    writeln( ans < INF ? ans : -1 );
}

class FordFulkerson(T) {
    int N, source, sink;
    int[][] adj;
    T[][] flow;
    bool[] used;

    this(int n, int s, int t) {
        N = n;
        source = s;
        sink = t;
        assert (s >= 0 && s < N && t >= 0 && t < N);
        adj = new int[][](N);
        flow = new T[][](N, N);
        used = new bool[](N);
    }

    void add_edge(int from, int to, T cap) {
        adj[from] ~= to;
        adj[to] ~= from;
        flow[from][to] = cap;
    }

    T dfs(int v, T min_cap) {
        if (v == sink)
            return min_cap;
        if (used[v])
            return 0;
        used[v] = true;
        foreach (to; adj[v]) {
            if (!used[to] && flow[v][to] > 0) {
                auto bottleneck = dfs(to, min(min_cap, flow[v][to]));
                if (bottleneck == 0) continue;
                flow[v][to] -= bottleneck;
                flow[to][v] += bottleneck;
                return bottleneck;
            }
        }
        return 0;
    }

    T run() {
        T ret = 0;
        while (true) {
            foreach (i; 0..N) used[i] = false;
            T f = dfs(source, T.max);
            if (f > 0)
                ret += f;
            else
                return ret;
        }
    }
}

CPSCO2019 Session3: F - Flexible Permutation

問題記録 CPSCO2019 DP 順列のDP

https://atcoder.jp/contests/cpsco2019-s3/tasks/cpsco2019_s3_f

問題概要

１～Nの順列Pのうち

  • Pi > i であるような i がA個
  • Pi < i であるような i がB個
  • Pi = i であるような i がN - A - B個

という条件をすべて満たすものが何通りあるか求めよ。

N <= 300

解法

O(N3)解

想定解はO(N3)で、これはwriterによる解説がわかりやすいのでそちらを見るのが早い。ざっくり言うと (1～iを1～Piのどこかに置いた、>がa個、<がb個) を状態にとるDPで、(i+1)への遷移ではまず(i+1)をP(i+1)に置いてから既に置いてある数のうちどれかひとつと位置を入れ替える（もしくは入れ替えを行わない）ことを考える。数を小さい順に見ていっていることからP(i+1)に置かれる数は必ず < になるし、手前に置かれることになる(i+1)は必ず > になる。なので考慮すべき遷移は「> と入れ替える」「< と入れ替える」「入れ替えない」の3通りのみとなり、O(1)で計算できる。よって状態O(N3), 遷移O(1)の全体O(N3)で答えが出る。

O(N2)解

先程の解説でも言及されている通りこの問題にはO(N2)解が存在する。O(N3)解では = の数も考慮するためにO(N3)となっているが、これを捨てることで高速化できる。まず = にする(N-A-B)個の場所を先に決め打ってしまうと、あとは(A+B)要素の順列において「>がA個」「<がB個」「=が0個」の数をかぞえる問題になる（最初に決め打った分は最後にnC(n-a-b)をかければよい）。

DP自体も基本的にはO(N3)版と同じようなことをする。ただし=を捨てたので状態は (1～iを1～Piのどこかに置いた、>がa個) のN2に落ちる。また=が0個なので遷移の際に「入れ替えない」という選択肢はない。ただしどうしても=の存在する状態を経由しないと生み出せない順列があり（例えば2143は213からでないと作れない）、その分を考慮する必要がある。そしてこのようなケースは新たな数を一度に2つ追加してやることでうまく拾うことができる。つまり=がひとつもない状態の中のどこかにひとつ=を生やして、それを(i+2)と入れ替えることにする、という方法である。これは i -> (i+2) への遷移になる。またどの位置を=にするかで (i+1) の自由度があるのでその分を掛けるようにする。

感想

O(N2)でだいぶ悩んだ　説明できている気がしない

コード (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;
 
immutable long MOD = 10^^9 + 7;

void main() {
    auto s = readln.split.map!(to!int);
    auto N = s[0];
    auto A = s[1];
    auto B = s[2];

    if (A == 0 && B == 0) {
        writeln(1);
        return;
    }
    
    if (A == 0 || B == 0) {
        writeln(0);
        return;
    }

    auto dp = new long[][](A+B+1, A+B+1);
    dp[2][1] = 1;

    foreach (i; 2..A+B) {
        foreach (j; 0..i+1) {
            (dp[i+1][j] += dp[i][j] * j % MOD) %= MOD;
            (dp[i+1][j+1] += dp[i][j] * (i - j) % MOD) %= MOD;
            if (i + 2 <= A + B) {
                (dp[i+2][j+1] += dp[i][j] * (i + 1) % MOD) %= MOD;
            }
        }
    }
    
    writeln( comb(N, A+B) * dp[A+B][A] % MOD );
}

long f(int n) {
    long ret = 1;
    foreach (i; 1..n+1) ret = ret * i % MOD;
    return ret;
}

long comb(int n, int r) {
    return f(n) * powmod(f(n-r), MOD-2, MOD) % MOD * powmod(f(r), MOD-2, MOD) % MOD;
}

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

CPSCO2019 Session1: F - Fruits in Season

問題記録 CPSCO2019 二分探索

https://atcoder.jp/contests/cpsco2019-s1/tasks/cpsco2019_s1_f

問題概要

N個の果物があり、これを1日目からN日目まで毎日ひとつずつ食べていく。i番目の果物はj日目に食べるとAi-abs(j-ti)×Biの満足度が得られる。またある順番で果物を食べたとき、全体の満足度は各日に得られた満足度の最小値となる。自由に食べる順番を決められるとき、全体の満足度の最大値を求めよ。

N <= 2×104

解法

二分探索の見た目をしてるのでそうすると、二分探索の中身では「ある満足度xが与えられる。すべての果物の満足度がx以上になるように果物を各日に割り当てることは可能か？」という問題を解くことになる。そして満足度の式から各果物がx以上となるような日は連続した区間になっているので、この問題はさらに次のように言い換えられる。

  • 区間[l, r]がN個与えられる。それぞれの区間に対してl<=i<=rなる整数iをひとつだけ割り当てるとき、整数1～Nを重複なくすべての区間に割り当てることは可能か？

この問題は貪欲に解くことが可能で、割り当てる数を前から順番に見ていって、「その数を含む区間のうち右端が最小のもの」をその数に対して割り当てていけばよい。priority queueとかで実装すると楽。途中でできなくなったらNG。

感想

にぶたんの中身でだいぶ迷ってしまった　区間やりくり系の問題なんかいろいろバリエーションがある気がしてあまり整理できていない

あと本番では式の形を見た瞬間にアッCHTだわからない！！となって一瞬で布団に入って普通に寝てしまった（最悪）

コード (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 A = N.iota.map!(_ => readln.split.map!(to!long).array).array;

    bool check(long x) {
        auto S = new Tuple!(long, long)[](N);
        foreach (i; 0..N) {
            long t = A[i][0];
            long a = A[i][1];
            long b = A[i][2];
            if (a < x) return false;
            long l = (-a + t*b + x) / b + ((-a + t*b + x) % b != 0);
            long r = (a + t*b - x) / b;
            S[i] = tuple(l, r);
        }
        S.sort!();
        auto pq = new BinaryHeap!(Array!(long), "a > b");
        int p = 0;
        foreach (i; 0..N) {
            while (p < N && S[p][0] <= i + 1) pq.insert(S[p][1]), p++;
            if (pq.empty || pq.front < i + 1) return false;
            pq.removeFront;
        }
        return true;
    }

    long hi = 10L^^18;
    long lo = -10L^^18;

    while (hi - lo > 1) {
        long mid = (hi + lo) / 2;
        (check(mid) ? lo : hi) = mid;
    }

    lo.writeln;
}

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

問題記録 いろはコンDay2 期待値 2人ゲーム

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

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 × 105

Hx <= 10

D <= 5000

解法

Hがたかだか10なので、DPの遷移を10×10の行列で表すことができる (dp[11], dp[10], ..., dp[2]) = A ・ (dp[10], dp[9], ..., dp[1]) のような形)。あとはよくあるやつで行列をセグ木とか平方分割とかに乗せればよさそうなのだが、やたら定数倍が重くメモリ制限も小さめなので適当にやるとだめになる (3*105 * 102 * 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();
}

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 <= 105

解法

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

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

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

自明なパターンから処理していくと、まず頂点数が奇数であるような連結成分が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;
}