yukicoder No.767 配られたジャパリまん

蟻本に載っている典型問題というのはすぐにわかったが、そこから詳細を詰めるのにかなり時間がかかった。
No.767 配られたジャパリまん - yukicoder

問題

縦H, 横Wの、格子点が(H+1)(W+1)個あるグリッドの上にK個の点がある。K個の点の部分集合2^K個のそれぞれに対し、その点を通行不能にした時に、右もしくは下に移動することを繰り返して(0,0)から(H,W)へ行く方法の総数を求めよ。

解法

通行不能な点が決まっているときは、その点のうちどこを通るかを全通り調べて包除原理をすることにより、点の個数をuとしてO(2^u)時間で計算できる。DPを上手くすることでO(u^2)でもできる。(dp[i] -> dp[j] : 係数はiからjへ行く方法の総数の(-1)倍。ここでは説明しない。参考: E - Ribbons on Tree)
これを2^K通り全てに対して行うことを考える。愚直にやるとO(2^KK^2)時間かかり間に合わないが、一部の計算をまとめて高速化することができる。
DP[S] = (Sの点が禁止される場合における上の包除原理の結果) とすると、実はv ∈ S, v < max SであるvそれぞれについてDP[S /\ {1, ..., v}]からDP[S]へ遷移があるようなDPを考えることができる。これでO(2^KK)時間になる。

解説の方法とは違う。(答えが包除原理をする前の値の交代和みたいな形で書けるので、高速メビウス変換をする。Superset側に減算する。)

実装上の注意点

ソートするので、ソート後の添字と元の添字の変換で間違えないようにする
mod = 10^8 + 7

提出: https://yukicoder.me/submissions/304194 (Rust)

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

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! {
        h: usize,
        w: usize,
        k: usize,
        ab: [(usize, usize); k],
    }
    let mut ab: Vec<(usize, usize, usize)>
        = ab.into_iter().enumerate().map(|(i, (x, y))| (x, y, i)).collect();
    ab.push((0, 0, k));
    ab.push((h, w, k + 1));
    ab.sort_unstable();
    let (fac, invfac) = init_fac();
    let mut trans = vec![vec![ModInt::new(0); k + 2]; k + 2];
    for i in 0 .. k + 2 {
        for j in i + 1 .. k + 2 {
            if ab[i].1 > ab[j].1 { continue; }
            let x = ab[j].1 - ab[i].1;
            let y = ab[j].0 - ab[i].0;
            trans[i][j] = fac[x + y] * invfac[x] * invfac[y];
        }
    }
    let mut dp = vec![ModInt::new(0); 1 << (k + 2)];
    dp[1] = ModInt::new(MOD - 1);
    for bits in 1 .. 1usize << (k + 1) {
        let bits = 1 + 2 * bits;
        let hi = 8 * std::mem::size_of_val(&bits)
            - bits.leading_zeros() as usize - 1;
        for j in 0 .. hi {
            if (bits & 1 << j) == 0 { continue; }
            let prev = bits & (1 << (j + 1)) - 1;
            dp[bits] -= dp[prev] * trans[j][hi];
        }
    }
    let mut ans = vec![ModInt::new(0); 1 << k];
    for bits in 0 .. 1 << k {
        let mut conv = 0;
        for i in 0 .. k {
            if (bits & 1 << i) != 0 { conv |= 1 << ab[i + 1].2; }
        }
        ans[conv] = dp[1 + 2 * bits + (1 << (k + 1))];
    }
    for bits in 0 .. 1 << k {
        puts!("{}\n", ans[bits]);
    }
}

まとめ

この手の包除を頭良くやると指数からO(N^2)になりますよという問題、まだ苦手っぽい。