yukicoder No.803 Very Limited Xor Subset

この手の問題は最近割と見る印象がある。
No.803 Very Limited Xor Subset - yukicoder

問題

長さNの整数列Aがある。以下の条件を満たす部分列Sの個数をmod 10^9 + 7で求めよ:

  • Sのxor和はX
  • Sは以下のM個の条件(i = 1, ..., M)を満たす
    • L[i]番目からR[i]番目の要素の中でSに含まれるものの個数を2で割った余りはtype[i]である

制約

  • 1 <= N <= 300
  • 0 <= M <= min(300, N(N - 1) / 2)
  • 0 <= X, A[i] <= 10^9
  • 1 <= L[i] <= R[i] <= N
  • 入力は全て整数

解法

A[i]やXを長さ30の0/1ベクトルとみなす。選ぶ個数についての制約は次のようにエンコードできる:

  • j個目の制約用に、A[i]やXに列を追加する(30+j列目になる)。A[i]については、iがL[j]からR[j]に含まれていたらその列は1にする。Xについては、type[j]にする。

こうした場合、元の問題で全ての制約を満たすことと、この問題でAの部分列Sの要素ごとのxorがXになっていることは同値。そのようなSの個数を数え上げればよい。

線形代数を使うと、求めるSの個数はA[i]たちで生成できる空間にXが要素として含まれていたら|ker A| = 2^{dim (ker A)}、含まれていなければ0であることがわかる。この手の計算はGauss-Jordanの消去法ででき、計算量はO(N * M * min(N, M))である。

登場する典型

  • 使用するかしないかをGF(2)の要素と思って線形代数をする
    • Gauss-Jordanの消去法
  • rank A + dim (ker A) = dim cod A
    • 特定の0/1ベクトルをxor和として実現する方法を数えたい場合にかなりよく見るテク
  • 制約をエンコードするために列を追加

実装上の注意点

  • Gauss-Jordanの消去法はそこそこ間違えやすいので気をつける
  • 数え上げだが実質的にはdim kerを求める問題なので、サンプルがあっているからといって過信しない

提出: #325170 (Rust)

// ModInt省略

// returns (rank(a), dim ker(a))
fn gauss_elim_gf2(a: &mut [Vec<i32>], b: &mut [i32]) -> Option<(usize, usize)> {
    let n = a.len();
    let m = b.len();
    let mut row = 0;
    for col in 0..m {
        if row >= n { break; }
        let mut row2 = None;
        for i in row..n {
            if a[i][col] == 1 {
                row2 = Some(i);
                break;
            }
        }
        if row2.is_none() { continue; }
        let row2 = row2.unwrap();
        a.swap(row, row2);
        for k in row + 1..n {
            if a[k][col] == 1 {
                for j in 0..m {
                    a[k][j] ^= a[row][j];
                }
            }
        }
        if b[col] == 1 {
            for j in 0..m {
                b[j] ^= a[row][j];
            }
        }
        row += 1;
    }
    for i in 0..m {
        if b[i] != 0 { return None; }
    }
    Some((row, n - row))
}

fn solve() {
    let out = std::io::stdout();
    let mut out = BufWriter::new(out.lock());
    macro_rules! puts {
        ($($format:tt)*) => (write!(out,$($format)*).unwrap());
    }
    input! {
        n: usize,
        m: usize,
        x: i64,
        a: [i64; n],
        tlr: [(i32, usize1, usize); m],
    }
    let mut mat = vec![vec![0; 30 + m]; n];
    for i in 0..n {
        for j in 0..30 {
            mat[i][j] = if (a[i] & 1 << j) != 0 { 1 } else { 0 };
        }
    }
    let mut b = vec![0; 30 + m];
    for i in 0..30 {
        b[i] = if (x & 1 << i) != 0 { 1 } else { 0 };
    }
    for i in 0..m {
        let (t, l, r) = tlr[i];
        for j in l..r {
            mat[j][30 + i] = 1;
        }
        b[30 + i] = t;
    }
    match gauss_elim_gf2(&mut mat, &mut b) {
        Some((_, ker)) => {
            puts!("{}\n", ModInt::new(2).pow(ker as i64));
        },
        None => {
            puts!("0\n");
        },
    }
}

まとめ

この手の線形代数をやる問題は結構すき (解けるので)

yukicoder No.802 だいたい等差数列

こういうのはかなり得意だと思っている。
No.802 だいたい等差数列 - yukicoder

問題

以下の条件を満たす長さNの整数列Aの個数を数え、10^9+7で割った余りを求めよ。

  • 1 <= A[1] <= ... <= A[N] <= M
  • i = 1, ..., N - 1に対して D1 <= A[i + 1] - A[i] <= D2

制約

  • 2 <= N <= 3 * 10^5
  • 1 <= M <= 10^6
  • 0 <= D1 <= D2 <= M

解法

まず愚直DPを考える。DP[i][j] := (Aの第i項まで見たとき、A[i]=jとなるような数列の個数)としたとき、以下のような遷移を持つDPを行えばよい。

  • init: DP[1][j] = 1
  • step: DP[i + 1][j + k] += DP[i][j] (D1 <= k <= D2)

これのDP[N][1] + ... DP[N][M]が答えである。これを実行することは計算量が多すぎて不可能なので、高速化を行う。
DP配列を多項式と考えることにして、f(i) = \sum DP[i][j] * x^jと置く。(最終的に興味があるのはx^M以下の項なので、x^{M+1}以上の項は無視する。また多項式級数(項が無限個続くもの)を断りなく同一視する。)
最初はf(1) = x + x^2 + ... + x^M = x/(1 - x) + O(x^{M + 1})である。
またstepを多項式の言葉に直すとf(i + 1) = f(i) * (x^{D_1} + ... + x^{D_2}) = f(i) * (x^{D_1} - x^{D_2 + 1}) / (1 - x)である。
また、多項式A + Bx + ... + Cx^Mに1 / (1 - x)を掛けたもの(のマクローリン展開)のx^Mの係数はA + B + ... + Cである。
よってこれらを組み合わせると、x * (x^{D_1} - x^{D_2 + 1})^{N - 1} / (1 - x)^{N + 1} (をマクローリン展開したもの)におけるx^Mの係数が答えであることがわかる。
これはどう計算すべきだろうか? まず、(x^{D_1} - x^{D_2 + 1})^{N - 1} = \sum_{i=0}^{N - 1} C(N - 1, i) * (-1)^i * x^{D_1 * (N - 1 - i) + (D_2 + 1) * i)}なので、「ある整数kに対してx^k / (1 - x)^{N + 1}におけるx^Mの係数」を高速に計算できればよいことがわかる(このようなクエリはN回実行される)。
これは1 / (1 - x)^{N + 1}におけるx^{M - k}の係数に他ならないので、M - k < 0であれば0で、M - k >= 0であればC(M - k + N, N)である。

登場する典型

  • 母関数
    • 愚直DPの高速化

実装上の注意点

  • 二項係数の扱いに気をつける
  • 母関数を使う時は式変形が複雑になりがちなので、式変形が間違っていないかどうかサンプルなどを使って確かめる

提出: #324980 (Rust)

fn fact_init(w: usize) -> (Vec<ModInt>, Vec<ModInt>) {
    let mut fac = vec![ModInt::new(1); w];
    let mut invfac = vec![0.into(); w];
    for i in 1 .. w {
        fac[i] = fac[i - 1] * i as i64;
    }
    invfac[w - 1] = fac[w - 1].inv();
    for i in (0 .. w - 1).rev() {
        invfac[i] = invfac[i + 1] * (i as i64 + 1);
    }
    (fac, invfac)
}

fn solve() {
    let out = std::io::stdout();
    let mut out = BufWriter::new(out.lock());
    macro_rules! puts {
        ($($format:tt)*) => (write!(out,$($format)*).unwrap());
    }
    input! {
        n: usize,
        m: usize,
        d1: usize,
        d2: usize,
    }
    let n = n - 1;
    let (fac, invfac) = fact_init(n + m + 2);
    let mut tot = ModInt::new(0);
    for i in 0..n + 1 {
        let deg = d1 * (n - i) + (d2 + 1) * i + 1;
        if deg > m { continue; }
        let rem = m - deg;
        let tmp = fac[rem + n + 1] * invfac[n + 1] * invfac[rem]
            * fac[n] * invfac[n - i] * invfac[i];
        if i % 2 == 0 {
            tot += tmp;
        } else {
            tot -= tmp;
        }
    }
    puts!("{}\n", tot);
}

まとめ

母関数周りも今後定跡になりそうな気がする。

早稲田大学プログラミングコンテスト2019 D - Choose Your Characters

さすがにこれが解けないのは反省。
D - Choose Your Characters

問題

N種類のキャラクターがいるゲームで対戦を行う。N種類のキャラクターには1からNまでの番号がついている。またそれらの間には相性という概念が定まっており、「キャラクターaはキャラクターbに対して有利(キャラクターbはキャラクターaに対して不利)」という関係がM個存在する。
このゲームでA君とB君が対戦する。A君にどうしても勝ちたいB君は、以下のような条件で特訓を行うことにした。

  • 特訓するキャラクターは[L, R] (閉区間)の範囲の番号のキャラクターに集中して行う。少なければ少ないほどよい。
  • A君がどのキャラクターaを選んでも、[L, R]内のキャラクターbであって、a != bでbはaに対して不利ではないようなbが存在する
    • 少なくとも五分五分以上にはなるようなbが存在する

条件を満たす閉区間[L, R]は存在するか? 存在する場合は長さがもっとも短いものを、複数存在する場合はLが最小のものを出力せよ。存在しない場合は-1を出力せよ。

制約

  • 2 <= N <= 100000
  • 1 <= M <= 200000
  • 同じ条件は1回まで出現する
  • 正反対の条件は出現しない

解法

二分探索を使うと、区間長R - L + 1を固定した時に、条件を満たすようなLが存在するかどうかをO(N)時間程度で判定できればよいことがわかる。
これは可能で、以下のようにすれば良い。

まず、各キャラクターbについて、B君がbを選んだ時にA君に選んでほしくないキャラクター(bが不利になるキャラクターとb自身)を管理する。
以降B君が選ぶ区間[L, R]を伸び縮みさせることを考える。[L, R]の中のキャラクターに対するA君に選んで欲しくないキャラクターが、何回数えられるかを管理する。明らかに、選んで欲しくないキャラクターとしてR - L + 1回登場するキャラクターが存在する場合、今見ている範囲[L, R]ではそのキャラクターには対応できない。
長さlenで条件を満たす区間が存在するかどうか確かめる場合には、[0, len - 1] -> [0, len] -> [1, len] -> [1, len + 1] -> ... のように、右端と左端を交互に右にずらしていけばよい。あるステップで区間に入ったり区間から除外される頂点をvとすれば、そのステップの計算量はO(vを選んだ時にA君に選んでほしくないキャラクターの個数)である。その合計はO(M)なので、長さlenでOKかどうかの判定もO(M)時間である。

2分探索なので全体の計算量はO(M log N)。

登場する典型

  • 尺取り法
    • セグメント木を使いたいが計算速度が間に合わないときは、尺取り法やスライド最小値でなんとかできないかを考える
  • 頻度表
    • setを使いたくなったら頻度表でなんとかできないかを考える

実装上の注意点

  • 尺取りなので気をつける
    • 2分探索時とLを求める時の2回行うので、関数にした方がよい
  • 頻度表なのでデータの整合性に気をつける
    • 変換用の関数を作るとよい

提出: #4543465 (Rust)

fn op(len: usize, a: &mut [usize], count: &mut usize, idx: usize, newval: usize) {
    if a[idx] == len {
        *count -= 1;
    }
    a[idx] = newval;
    if newval == len {
        *count += 1;
    }
}

fn sweep(g: &[Vec<usize>], len: usize) -> Option<usize> {
    let n = g.len();
    // sweep
    let mut freq = vec![0; n];
    let mut cnt = 0;
    for i in 0..len {
        for &v in &g[i] { 
            let val = freq[v] + 1;
            op(len, &mut freq, &mut cnt, v, val);
        }
    }
    for i in 0..n - len + 1 {
        if cnt == 0 {
            return Some(i);
        }
        // progress
        if i < n - len {
            for &v in &g[i + len] {
                let val = freq[v] + 1;
                    op(len, &mut freq, &mut cnt, v, val);
                }
            for &v in &g[i] {
                let val = freq[v] - 1;
                op(len, &mut freq, &mut cnt, v, val);
            }
        }
    }
    None
}

fn solve() {
    let out = std::io::stdout();
    let mut out = BufWriter::new(out.lock());
    macro_rules! puts {
        ($($format:tt)*) => (write!(out,$($format)*).unwrap());
    }
    input! {
        n: usize,
        m: usize,
        ab: [(usize1, usize1); m],
    }
    let mut g = vec![vec![]; n];
    for &(a, b) in &ab {
        g[b].push(a);
    }
    for i in 0..n {
        g[i].push(i);
    }
    let mut pass = n + 1;
    let mut fail = 1;
    while pass - fail > 1 {
        let mid = (pass + fail) / 2;
        if sweep(&g, mid).is_some() {
            pass = mid;
        } else {
            fail = mid;
        }
    }
    // len = pass
    if pass > n {
        puts!("-1\n");
        return;
    }
    let l = sweep(&g, pass).unwrap();
    puts!("{} {}\n", l + 1, l + pass);
}

まとめ

本番は「セグ木に平衡2分探索木を乗せて、毎回O(N log^2 N)時間かけて判定すればいいから全体ではO(N log^3 N)!w」みたいなことを考えて途方にくれていたが、セグ木も平衡2分探索木も自然により良い方法が存在したので、考えるべきだった。

Mujin Programming Challenge 2017 C - Robot and String (解説を見た)

自分で考察したときは手掛かりが何もなかった。
C - Robot and String

問題

英小文字からなる文字列Sが与えられる。以下の形のクエリにQ回答えよ。

  • Sの部分文字列S[L[i]...R[i] + 1] (半開区間、L[i]文字目からR[i]文字目までを抜き出した文字列)について、以下の規則で簡約していった結果が空文字列になるか判定せよ:
    • 文字列の中に隣接する同じ文字がある場合に、それらの中で一番左のものについて、その同じ2文字を辞書順で次の文字1文字に変換する。ただしzzの変換先は空文字列とする。

制約

  • 1 <= |S| <= 5 * 10^5
  • 1 <= Q <= 10^5

解法

まず、各iについて、S[i...x]が空文字列に簡約されるような最小のx > iを求めたい気分になる。文字列の簡約は左から優先して行われるので、空文字列から始めて、i文字目から文字を追加して簡約して…、という処理をj文字目まで行った結果は、S[i..j]を簡約した結果と同じになる。
i文字目から文字の追加・簡約を繰り返すとき、文字列が最初に空文字列に簡約されるのはいつだろうか? そのような場合はもともと"z"だった文字列に'z'を追加した場合、およびその場合に限られる。よって、各iについてS[i..x]が"z"に簡約されるような最小のx > iを求めたい気分になる。この考察を繰り返すと、各iについてS[i..x]が各文字1文字に簡約されるような最小のx > iを求めたい気分になる。

実は以上の情報がDPをする上で十分である。DP[i] := (S[i...x]が空文字列に簡約されるような最小のx > i), DPS[c][i] := (S[i..x]がc 1文字に簡約されるような最小のx > i)とすると、以下の更新を行えばいいことがわかる(更新はすべてminで、式の中に一つでもinvalidな値がある場合は実行されない):

  • DPS[c][i] <- DPS[c][DP[i]]
  • DPS[c + 1][i] <- DPS[c][DPS[c][i]] (c < 'z', ダブリングに似ている)
  • DP[i] <- DPS['z'][DPS['z'][i]]

これらの更新はどこから始めればよいか? まず、i = n - 1のときは更新なしで正しい答えが手に入る。i = n - 2のときは、DP[n - 1]と各文字cについてDPS[c][n - 1]がわかっていれば正しい答えが手に入る。この考察を帰納的に行うことで、iの降順に更新を行えば良いことがわかる。

以上でDP[i]がわかった。あとはこれをダブリングすれば、各クエリに対してO(log |S|)時間で答えを求めることができる。このダブリングの手法は例えばLCAの計算などで有名なので、ここでは省略する。

登場する典型

  • ダブリング
  • 再帰的に考察を進め、自己完結したら止める

実装上の注意点

  • 添字の扱い(上限n - 1かnか)に気をつける

提出: #4397603 (Rust)

// This solution was written after the author read the editorial.
fn solve() {
    let out = std::io::stdout();
    let mut out = BufWriter::new(out.lock());
    macro_rules! puts {
        ($($format:tt)*) => (write!(out,$($format)*).unwrap());
    }
    input! {
        s: chars,
        lr: [(usize1, usize)],
    }
    let n = s.len();
    let s: Vec<_> = s.into_iter().map(|x| (x as u8 - b'a') as usize).collect();
    let mut dp = vec![n + 1; n + 1];
    let mut dps = vec![vec![n + 1; 26]; n + 1];
    for i in (0..n).rev() {
        let mut que = BinaryHeap::new();
        que.push((n - i - 1, s[i]));
        while let Some((pidx, ch)) = que.pop() {
            let idx = n - pidx;
            if ch <= 25 {
                if dps[i][ch] <= idx { continue; }
                dps[i][ch] = idx;
                if dps[idx][ch] <= n {
                    que.push((n - dps[idx][ch], ch + 1));
                }
            } else {
                assert_eq!(ch, 26);
                if dp[i] <= idx { continue; }
                dp[i] = idx;
                for j in 0..26 {
                    if dps[idx][j] <= n {
                        que.push((n - dps[idx][j], j));
                    }
                }
            }
        }
    }
    const B: usize = 20;
    let mut dbl = vec![vec![n + 1; n]; B];
    for i in 0..n {
        dbl[0][i] = dp[i];
    }
    for j in 1..B {
        for i in 0..n {
            let stage = dbl[j - 1][i];
            if stage < n {
                dbl[j][i] = dbl[j - 1][stage];
            }
        }
    }
    for (mut l, r) in lr {
        for i in (0..B).rev() {
            if dbl[i][l] <= r {
                l = dbl[i][l];
            }
            if l == n { break; }
        }
        puts!("{}\n", if l == r { "Yes" } else { "No" });
    }
}

まとめ

「1文字に簡約できる範囲を持つDP」という解説の方法が天才すぎて、それっぽい考察の道筋をこじつけるのに苦労した。

CODE FESTIVAL 2017 Elimination Tournament Round 2 A - Colorful MST

あさぷろで解いたときは解けなかった。面白い問題。
A - Colorful MST

問題

N頂点M辺のグラフがある。各頂点には色1から色Kまでのどれかの色が塗られているか、あるいはまだ色が塗られていない。
以下の操作を行ったときにできるグラフの最小全域木の重みの最小値を求めよ。どのようにしても最小全域木を作れない場合は-1を出力せよ。

  1. まだ色が塗られていない頂点に色1から色Kまでのどれかの色を塗る。別の頂点には別の色を塗ってよい。
  2. 同じ色の頂点を1点に潰して、K頂点のグラフを作る。このとき、対応する色の頂点が存在しないような色については、その頂点は0本の辺が接続することになる。

制約

  • 1 <= K <= N <= 10^5
  • 1 <= M <= 10^5
  • 1 <= (辺の重み) <= 10^9
  • 与えられるグラフは、単純とも連結とも限らない

解法

適当に(K - 1)辺(実行可能解)を取った時、その辺集合が全域木となるような色の塗り方は存在するだろうか? 存在するためには以下の条件が必要十分:

  • すでに塗られている頂点を色ごとに1点に潰した時、潰した後の点集合が森になっている。

(必要性): どのような塗り方をしても、塗られていない状態よりは潰れるので、森になりにくい。よって自明
(十分性): 塗り方を作ればよい。まだ塗られていない頂点に接続する辺1つにつき、まだ使われていない色を1つずつ使っていけばよい。
以上より、すでに塗られている頂点を同一視しながらKruskal法を実行するだけで解ける。

登場する典型

  • 実行可能解の条件を考える

実装上の注意点

  • オーバーフローの危険があるのでlong longを使う
  • ちょうど(K - 1)個の辺を選ぶ
    • 1敗 (は?)

提出: Submission #4343453 - CODE FESTIVAL 2017 Elimination Tournament Round 2 (Parallel) (C++)

const int N = 112222;
int a[N], b[N], w[N];

typedef pair<ll, PI> PLPI;

int main(void) {
  ios::sync_with_stdio(false);
  cin.tie(0);
  int n, m, k;
  cin >> n >> m >> k;
  VI c(n);
  REP(i, 0, n) {
    cin >> c[i];
  }
  vector<PLPI> pool;
  REP(i, 0, m) {
    cin >> a[i] >> b[i] >> w[i];
    a[i]--, b[i]--;
    pool.push_back(PLPI(w[i], PI(a[i], b[i])));
  }
  sort(pool.begin(), pool.end());
  UnionFind uf(k + n);
  int conn = k;
  ll tot = 0;
  REP(i, 0, pool.size()) {
    if (conn <= 1) break;
    ll ww = pool[i].first;
    int x = pool[i].second.first;
    int y = pool[i].second.second;
    x = c[x] == 0 ? k + x : c[x] - 1;
    y = c[y] == 0 ? k + y : c[y] - 1;
    assert (x >= 0);
    assert (y >= 0);
    if (not uf.is_same_set(x, y)) {
      tot += ww;
      uf.unite(x, y);
      conn--;
    }
  }
  cout << (conn > 1 ? -1 : tot) << "\n";
}

まとめ

実行可能解の条件を考えるのも割とよく見る典型という気がする。

CODE FESTIVAL 2017 Exhibition A - Awkward

当時Code Festivalに参加していた時は、何をすればいいのか見当もつかず、出題陣のライブ解説も理解不能だったが、今解き直したらかなり簡単に思えた。A - Awkward

問題

N頂点の木が与えられる。木の頂点の並べ方はN!通りあるが、その中で以下を満たすものの個数を10^9 + 7で割った余りを求めよ。

  • どの2個の隣接する頂点も、隣同士にならない。

解法

包除原理を使う。以下の値をx = 0, ..., N - 1について数え上げることができればよい。

  • N - 1個の辺のうちx個を選ぶ。それらの辺については制約が破られている(隣同士になっている)ような並べ方の総数の合計。
    • 例えば入力例3 (N=5, b = [1, 1, 3, 3])では120, 192, 96, 16, 0となる。問題の答えは120 - 192 + 96 - 16 + 0 = 8。

これを求めるにはどうすれば良いだろうか? まず、辺の集合Sをfixしたときに、その集合Sについては制約が破られているような並べ方の総数を求めよう。辺集合に枝分かれがある(同じ頂点にSの3本以上の辺が接続している)ときは0通りである。それ以外の時を考える。Sを辺集合とした木の部分グラフは、1点と長さ1以上のパスの寄せ集めになっている。パスは列の中で一続きの部分列になること、およびどちらを列の先頭に持ってくるかが2通りあることを考えると、場合の数は(ものの個数)! * 2^(長さ1以上のパスの個数)であることがわかる。
次に、集合Sではなく集合の大きさx = |S|をfixしたときの和を計算したい。「(ものの個数)!」についてはあとでまとめて掛けることにして、「2^(長さ1以上のパスの個数)」の合計をもとめることにすると、これは以下のようなDPで計算することができる。

  • DP[v][x][k]: 頂点vを根とする部分木において、連結成分の個数がx個あり、vからは辺がk本伸びているときの、2^(長さ1以上のパスの個数)の合計

これは二乗の木 DP - (iwi) { 反省します - TopCoder部のテクニックを使えば、O(N^2)で計算することができる。

登場する典型

  • 包除原理
  • 2乗の木DP

実装上の注意点

  • 遷移が複雑なのでバグらせないようにする
  • 余計なループを回すとO(N^3)になるので回さないようにする
    • RustだとDP配列をそのまま返してしまうのが賢明な気がする

提出: Submission #4342066 - CODE FESTIVAL 2017 Exhibition (Parallel) (Rust)

// ModInt省略
// Depends on ModInt.rs
fn fact_init(w: usize) -> (Vec<ModInt>, Vec<ModInt>) {
    let mut fac = vec![ModInt::new(1); w];
    let mut invfac = vec![0.into(); w];
    for i in 1 .. w {
        fac[i] = fac[i - 1] * i as i64;
    }
    invfac[w - 1] = fac[w - 1].inv();
    for i in (0 .. w - 1).rev() {
        invfac[i] = invfac[i + 1] * (i as i64 + 1);
    }
    (fac, invfac)
}

fn dfs(ch: &[Vec<usize>], v: usize) -> Vec<[ModInt; 3]> {
    let mut dp = vec![[ModInt::new(0); 3]; 2];
    dp[1][0] = ModInt::new(1);
    for &w in &ch[v] {
        let sub = dfs(ch, w);
        let sub_sz = sub.len() - 1;
        let cur_sz = dp.len() - 1;
        let mut next_dp = vec![[ModInt::new(0); 3]; cur_sz + sub_sz + 1];
        for i in 1..sub_sz + 1 {
            for j in 1..cur_sz + 1 {
                // add one edge
                for k in 0..2 {
                    next_dp[i + j - 1][1] += sub[i][k] * dp[j][0];
                    next_dp[i + j - 1][2] += sub[i][k] * dp[j][1];
                }
                // don't add an edge
                for k in 0..3 {
                    for l in 0..3 {
                        next_dp[i + j][l] += sub[i][k] * dp[j][l]
                            * if k >= 1 { 2 } else { 1 };
                    }
                }
            }
        }
        dp = next_dp;
    }
    dp
}

fn solve() {
    let out = std::io::stdout();
    let mut out = BufWriter::new(out.lock());
    macro_rules! puts {
        ($($format:tt)*) => (write!(out,$($format)*).unwrap());
    }
    input! {
        n: usize,
        b: [usize1; n - 1],
    }
    let mut ch = vec![vec![]; n];
    for i in 1..n {
        ch[b[i - 1]].push(i);
    }
    const W: usize = 3000;
    let mut invtbl = vec![ModInt::new(0); W];
    for i in 1..W {
        invtbl[i] = ModInt::new(i as i64).inv();
    }
    let (fac, _invfac) = fact_init(W);
    let dp = dfs(&ch, 0);
    let mut tot = ModInt::new(0);
    let mut factor = ModInt::new(1);
    for i in (0..n + 1).rev() {
        let mut t = ModInt::new(0);
        for k in 0..3 {
            t += dp[i][k] * fac[i]
                * if k >= 1 { 2 } else { 1 };
        }
        tot += t * factor;
        factor = -factor;
    }
    puts!("{}\n", tot);
}

まとめ

当時に比べこの手の問題が多く出題されているなど、界隈が進歩しているので、今出題するとしたら800点くらいになりそうな気がする。

CODE FESTIVAL 2018 Final I - Homework

F問題は解説を見ないとわからなかったがI問題は自力で解けた。
I - Homework

問題

荷物がN個ある。i番目の荷物はコスト2^A[i]、価値B[i]である。合計価値をK以上にするために必要なコストの合計を求めよ。

制約

  • 入力は全て整数
  • 1 <= N <= 10^5
  • 0 <= A[i] <= 30
  • 1 <= B[i] <= 10^9
  • 1 <= K <= \sum B[i]

解法

二分探索をすることにして、問題を「合計コストx以下で合計価値をK以上にできるか?」という形に変形する。

ここで、xの偶奇とコスト1の荷物をどうするかに着目しよう。xが奇数であれば、少なくとも1個はコスト1の荷物を取るべきである。またxが偶数であっても奇数であっても、最初に取る荷物以外の荷物は2個ずつ取っていくべきである。これを使うと、xの各ビットを下から見て、立っていたら一番価値の高い荷物を取り、そのあとそれ以外の荷物を2個まとめてコストが2倍の荷物として扱うとよいことがわかる。
これの計算量は、ほとんどの荷物の個数が毎ターン半分になること、最大1個の荷物がそれ以外に次のターンに持ち越されることを考えると、O(N + (ターン数)) = O(N + max A[i])である。

登場する典型

  • 二分探索
  • 1, 2, 2^2, ...という等比数列の典型パターン
    • 今回は下位ビットに着目するといけた

実装上の注意点

  • 二分探索の範囲
    • 上限はN * 2^30以上にすべき
    • 上限に応じてxの大きさも変わるので、ターン数も30ではなくもっと大きくしないといけない

提出: Submission #4325140 - CODE FESTIVAL 2018 Final (Parallel) (Rust)

const B: usize = 31;
const C: usize = 50;
 
// x: max time
fn can_pack(tbl: &[Vec<i64>], k: i64, x: i64) -> bool {
    let mut rem = k;
    let mut cur_bag = tbl[0].clone();
    for i in 0..C {
        cur_bag.sort();
        if (x >> i & 1) == 1 && cur_bag.len() >= 1 {
            rem -= cur_bag.pop().unwrap();
        }
        let mut nxt = vec![];
        while let Some(mut fst) = cur_bag.pop() {
            if let Some(snd) = cur_bag.pop() {
                fst += snd;
            }
            nxt.push(fst);
        }
        if i + 1 < B {
            nxt.extend_from_slice(&&tbl[i + 1]);
        }
        cur_bag = nxt;
    }
    rem <= 0
}
 
fn solve() {
    let out = std::io::stdout();
    let mut out = BufWriter::new(out.lock());
    macro_rules! puts {
        ($format:expr) => (write!(out,$format).unwrap());
        ($format:expr, $($args:expr),+) => (write!(out,$format,$($args),*).unwrap())
    }
    input! {
        n: usize,
        k: i64,
        ab: [(usize, i64); n],
    }
    let mut tbl = vec![vec![]; B];
    for (a, b) in ab {
        tbl[a].push(b);
    }
    let mut pass = 1 << C;
    let mut fail = 0;
    while pass - fail > 1 {
        let mid = (pass + fail) / 2;
        if can_pack(&tbl, k, mid) {
            pass = mid;
        } else {
            fail = mid;
        }
    }
    puts!("{}\n", pass);
}

まとめ

2の冪乗系コストの問題は、たまに見ると面白い。