diverta 2019 Programming Contest F - Edge Ordering

5時間45分かけてようやく解けた。本番で解きたかった。
F - Edge Ordering

問題

N頂点M辺の無向グラフがあり、辺には0からM-1までの番号が付いている。0番からN-2番までの辺は全域木をなすことが保証されている。
これらに1からMまでの整数の重みを付けた重み付き無向グラフを考える。ある重み付けについて、それの最小全域木が0番からN-2番までの辺でできた木であるような場合、その重み付けを妥当な重み付けと呼ぶ。妥当な重み付けすべてについて、そのときの最小全域木の重みの総和を求めよ。

解法

Tを、0番からN-2番までの辺からなる全域木とする。
まず、妥当な重み付けがどのような性質を満たすべきかを考える。N-1番からM-1番までの各辺eについて、以下が成り立つことが必要十分であることがわかる:

• e = (u, v)としたとき、Tの上のu-v間のパスをe_1, ..., e_kとしたとき、e_1, ..., e_kの重みは全てeの重み未満である。

ここで、本来の問題を簡略化した問題を考える。「妥当な重み付けは何通りあるだろうか?」
これには F - Leftmost Ball と同じテクニックが使える。各部分集合に対して、以下の値がわかればよい:

• Sのどれかの辺よりも大きい重みがついていなければならない(どれかの辺に依存する)辺の集合をとする。S + Uを「トポロジカルソート」(上の大小関係を満たすように1から|S|+|U|までの整数で重み付けすること)する方法の個数の、全ての並び替えの個数(|S|+|U|)!に対する割合はいくらか?

これは以下のようなDPでできる:

• DP1[{}] = 1
• DP1[ ] += DP1[S] / ints[ ] (ここで、 ints[S] := (Sの辺とSのどれかの辺に依存する辺の個数) = (上の|S| + |U|))

DP1[{}] = 1は自明である(空を並べる方法は1通りであり、これはトポロジカルソート)。Sに新しい辺eを追加するときを考える。Sの辺に依存せずeに依存するような辺はints[ ] - 1 - ints[S]本ある。それらをもともとあるints[S]本の辺の中に織り交ぜ、それらの先頭にeをつける方法の割合は1 / (全体の長さ)である。

以上で妥当な重み付けの総数のM!に対する割合が、したがって妥当な重み付けの総数がわかった。本来の問題に戻る。
妥当な重み付けについて、0番からN-2番までの辺についた重みの合計の総和が知りたいのだった。これを求めるために、以下の部分問題を解く必要がある。

• 自分を含めたx人の人が横一列に並んでおり、自分の左から数えた位置の期待値はkである(1 <= k <= x)。この列にy人がランダムに入る。ランダムに入った後の自分の左から数えた位置の期待値を求めよ。

これを考える。まずkが固定値である場合(自分が必ずk番目にいる場合)を考える。y人それぞれについて、自分の左に入る確率はk / (x + 1)である。よって、自分の位置の増分の期待値はy * k / (x + 1)であり、よって位置の期待値はk * (x + y + 1) / (x + 1)である。
これはkの一次式なので、kが固定値でない場合にも成り立つ。

元の問題に戻る。簡略化した問題のDPを少し変更することで元の問題を解きたい。DP2[S] := (上のS + Uを並べた際の、Sの辺の順位の和 (妥当でない場合は0) の期待値)と置くと、以下のようなDPでできる:

• DP2[{}] = 0
• DP2[ ] += (DP2[S] * ints[ ] / (ints[S] + 1) + DP1[S] * (|S| + 1)) / ints[ ]

DP2[{}] = 0は自明である(空列において、順位の和は0)。Sに辺eを追加するにあたり、もともとあるints[S]本の前に必ず追加されるのはeの1本で、残りのints[ ] - 1 - ints[S]本はランダムに織り交ぜられる。よって、S内の辺fに対して、もともとの順位の期待値がuだったとすると、新しい順位の期待値は1 + f * ints[ ] / (ints[S] + 1)である。S内に辺は|S|本あるので、期待値は|S| + (もともとの期待値) * ints[ ] / (ints[S] + 1)になる。

残りはintsの計算である。各辺eについて以下の操作ができればよい:

• e = (u, v)とする。Tの上のu-v間のパスをe_1, ..., e_kとしたとき、S = {e_1, ..., e_k}としてなるすべてのTに対してints[T] += 1をする。

このような操作は高速ゼータ変換で行うことができる。

以上でこの問題を解くことができた。計算量はO(N * 2^N)である。

登場する典型

• ゼータ変換
  • intsの計算に必要
  • ある部分集合について、「なるすべてのTに対してdp[T] += 1; 」という操作を実行したくなったら考えるべき
• 確率
  • それっぽい方法はたくさんあるが、大体は誤り
  • 確率を扱うときには、どこが独立か、どこの独立性を無視して良いかというのを正確に理解する
    • 例えば今回は、妥当な重み付けが0番〜N-2番の辺を軽くする方に偏っているので、妥当な重み付けに限ったとき、例えば(0番の重み) < (N番の重み)となる条件付き確率は1/2よりも大きい
• 期待値と確率が入り混じったDPの遷移
  • 期待値の変換が線形なら、確率と「確率 * 期待値」をもっておけばうまくいく

実装上の注意点

• 解説では場合の数を求めているが、この解法では確率と期待値を求めている。どちらも等価だが、確率の方はうまく考えないと混乱するので注意。

提出: #5391914 (Rust)

// ModInt省略

fn dfs(r: usize, v: usize, par: usize,
       u: usize, used: &mut [Vec<usize>], g: &[Vec<(usize, usize)>]) {
    used[r][v] = u;
    for &(w, idx) in &g[v] {
        if par == w { continue; }
        dfs(r, w, v, u | 1 << idx, used, g);
    }
}
 
// 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 construct_ints(n: usize, ab: &[(usize, usize)]) -> Vec<usize> {
    let m = ab.len();
    let mut g = vec![vec![]; n];
    for i in 0..n - 1 {
        let (a, b) = ab[i];
        g[a].push((b, i));
        g[b].push((a, i));
    }
    let mut used = vec![vec![0; n]; n];
    for i in 0..n {
        dfs(i, i, n, 0, &mut used, &g);
    }
    // cnt, ints
    let mut cnt = vec![0usize; 1 << (n - 1)];
    let mut ints = vec![0usize; 1 << (n - 1)];
    for i in 0..m {
        let (a, b) = ab[i];
        cnt[used[a][b]] += 1;
    }
    for i in 0..n - 1 {
        for bits in 0..1 << (n - 1) {
            if (bits & 1 << i) == 0 {
                cnt[bits | 1 << i] += cnt[bits];
            }
        }
    }
    for bits in 0..1 << (n - 1) {
        ints[bits] = m - cnt[(1 << (n - 1)) - 1 - bits];
    }
    ints
}
 
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 (fac, invfac) = fact_init(m + 1);
    let mut inv = vec![ModInt::new(0); m + 1];
    for i in 1..m + 1 {
        inv[i] = fac[i - 1] * invfac[i];
    }
    let ints = construct_ints(n, &ab);
    let mut dp1 = vec![ModInt::new(0); 1 << (n - 1)];
    let mut dp2 = vec![ModInt::new(0); 1 << (n - 1)];
    dp1[0] = ModInt::new(1);
    for bits in 0..1 << (n - 1) {
        for j in 0..n - 1 {
            if (bits & 1 << j) == 0 {
                let x = ints[bits] as i64;
                let y = ints[bits ^ 1 << j] as i64 - x - 1;
 
                let intercalate = inv[(x + y + 1) as usize];
                dp1[bits ^ 1 << j] += dp1[bits] * intercalate;
                let ratio = ModInt::new(x + y + 1) * inv[x as usize + 1];
                dp2[bits ^ 1 << j] += intercalate *
                    (ratio * dp2[bits] + dp1[bits] * (bits.count_ones() as i64 + 1));
            }
        }
    }
    puts!("{}\n", fac[m] * dp2[(1 << (n - 1)) - 1]);
}

まとめ

最後のDP2を計算するパートを思いつくのに一番時間がかかった。解けた後意気揚々と解説を見たらその一番時間がかかるパートが「ほぼ先ほどと同様のDP」で済まされていてキレかけた(は?)。多分諦めて解説を見ていたらキレていたと思う。