Codeforces Round #669 (Div. 2) D. Discrete Centrifugal Jumps

考察・実装合計で 2 時間かかった上に、無駄な実装をしてしまった。
Problem - D - Codeforces

問題

長さ n の数列 h が与えられる。1 番目の要素から n 番目の要素まで離散的ジャンプだけを使って移動したい。i < j に対して i から j へのジャンプが離散的であるとは、max(h[i+1], ..., h[j-1]) < min(h[i], h[j]) または min(h[i+1], ... h[j-1]) > max(h[i], h[j]) が成立することをいう。離散的ジャンプの最小回数を求めよ。

制約

  • 2 <= n <= 3 * 10^5
  • 1 <= h[i] <= 10^9

解法

愚直な O(N^2) DP (dp[i] := i から n までの最短手数) を考える。i から j へジャンプできるのは j が i+1 から貪欲に取った単調増加列か単調減少列に含まれる場合に限られる。i の降順に見て、貪欲な単調増加列・単調減少列を stack で管理すれば遷移すべき頂点の集合はわかる。遷移の個数の合計が問題だが、一般には O(N^2) であるため、RMinQ に対応できるセグメント木などで stack の中身に対応する要素を持てば、計算量は O(N log N) 時間となる。(実は、editorial によれば遷移の個数は 4N 以下であるため、このような仕掛けは不要である。)

登場する典型

  • stack を使って、現在の地点から見た貪欲な単調増加列・単調減少列を持つ
    • stack の上を二分探索して、どこまでが遷移としてありえるかを調べる
  • うまい数え方で、遷移の個数が合計 O(N) であることを示す

実装上の注意点

  • 二分探索で off-by-one error が多発しやすいので、stack がどの順番でデータを持つか、二分探索でどこを見るべきかを正確に詰める
    • 今回の場合、降順にデータを持っている stack s に対しては s[j] >= h[i] なる最大の j が欲しいものであり、j が 0 以上なら j 以降の全ての要素を、j が -1 なら 0 以降の全ての要素が遷移先としてありえるので、セグメント木へのクエリは [max(1, j) - 1, |st|) となる。

提出: #93371509 (Rust)

// SegTree, upper_bound omitted

fn solve() {
    let out = std::io::stdout();
    let mut out = BufWriter::new(out.lock());
    macro_rules! puts {
        ($($format:tt)*) => (let _ = write!(out,$($format)*););
    }
    input! {
        n: usize,
        h: [i64; n],
    }
    let mut coord = h.clone();
    coord.sort(); coord.dedup();
    let mut dp = vec![0; n];
    let mut up: Vec<Reverse<(i64, usize)>> = vec![];
    const INF: i64 = 1 << 50;
    let mut ups = SegTree::new(n + 1, min, INF);
    let mut down: Vec<(i64, usize)> = vec![];
    let mut downs = SegTree::new(n + 1, min, INF);
    up.push(Reverse((h[n - 1], n - 1)));
    ups.update(0, 0);
    down.push((h[n - 1], n - 1));
    downs.update(0, 0);
    for i in (0..n - 1).rev() {
        let idx = up.upper_bound(&Reverse((h[i], 0)));
        let idx = max(idx, 1) - 1;
        let val = ups.query(idx, up.len());
        let mut mi = val;
        let idx = down.upper_bound(&(h[i], 1 << 30));
        let idx = max(idx, 1) - 1;
        let val = downs.query(idx, down.len());
        mi = min(mi, val);
        dp[i] = mi + 1;
        // upd
        while let Some(Reverse((x, idx))) = up.pop() {
            if x <= h[i] {
                continue;
            }
            up.push(Reverse((x, idx)));
            break;
        }
        up.push(Reverse((h[i], i)));
        ups.update(up.len() - 1, dp[i]);
        while let Some((x, idx)) = down.pop() {
            if x >= h[i] {
                continue;
            }
            down.push((x, idx));
            break;
        }
        down.push((h[i], i));
        downs.update(down.len() - 1, dp[i]);
    }
    puts!("{}\n", dp[0]);
}

まとめ

遷移の個数が O(N) であることを示せず、O(N^2) を O(N log N) に高速化する典型テクニックの適用でかなり詰まってしまった。反省。

yukicoder No.1031 いたずら好きなお姉ちゃん

解説より速い方法で AC できたので記念に。
No.1031 いたずら好きなお姉ちゃん - yukicoder

問題

長さ N の順列 p が与えられる。以下をちょうど1回実行して得られる順列として、あり得るものの個数を計算せよ:
1. non-empty な区間 [l, r] を選ぶ。
2. [l, r] の範囲内の max と min について、それら 2 要素を入れ替える。

制約

  • 1 <= N <= 100000
  • p は 1, 2, ..., N の順列

解法

max と min の組み合わせとしてそれぞれ p[l], p[r] が選ばれることと、以下の条件は同値:

  • p の l 番目と r 番目の間に、p[r] より小さい要素や p[l] より大きい要素が存在しない。

r を全探索して、l としてあり得る値の個数を高速に求める方針をとる。
まず、各 r について、r を含む p[r] 以上の要素の範囲を求める。これは SparseTable と二分探索を用いて O(N log N) 時間でできる。この範囲を R(r) と呼ぶことにする。
次に、各 r について、l としてあり得る要素の個数を求める。これは l in R(r) かつ p の l 番目と r 番目の間に p[l] より大きい要素が存在しないことと同値。
まず l < r となるもののみを数える。r を左からスキャンして、p[r] から左側に貪欲に取っていった単調増加部分列 (の添字列) をスタックを使って管理する (要素を追加・削除) ことにすると、このスタックに載っていてかつ R(r) の要素でもあるようなものが全て l としてあり得る。スタックの操作がならし O(1) 時間、スタックに載っているある値以上の要素の個数を調べるのが二分探索で O(log N) 時間でできるため、合計 O(N log N) 時間である。
l > r となるものを数えるには、右からスキャンすればよく、そのときも同様の解析ができる。

以上から、合計 O(N log N) 時間でできた。

登場する典型

  • SparseTable
    • RMQ + 二分探索のような、セグメント木で愚直に実行すると O(N log^2 N) 時間かかる処理を O(N log N) 時間に抑えるときに重宝する
  • スタック
    • 「現在の要素から貪欲に取っていったときの増加列」などを作りたいとき、スタックで合計 O(N) 時間で作れる

実装上の注意点

  • off-by-one error に注意する

提出: #471138 (Rust) No.1031 いたずら好きなお姉ちゃん - yukicoder (Rust)

fn solve() {
    let out = std::io::stdout();
    let mut out = BufWriter::new(out.lock());
    macro_rules! puts {
        ($($format:tt)*) => (let _ = write!(out,$($format)*););
    }
    input! {
        n: usize,
        p: [usize; n],
    }
    let mut tot = 0i64;
    let spt = SparseTable::new(&p, min);
    let mut lohi = vec![(0, 0); n];
    for i in 0..n {
        let mut fail = n;
        let mut pass = i;
        while fail - pass > 1 {
            let mid = (pass + fail) / 2;
            if spt.query(i, mid) < p[i] {
                fail = mid;
            } else {
                pass = mid;
            }
        }
        let hi = pass;
        // bias: -1
        fail = 0;
        pass = i + 1;
        while pass - fail > 1 {
            let mid = (pass + fail) / 2;
            if spt.query(mid - 1, i) < p[i] {
                fail = mid;
            } else {
                pass = mid;
            }
        }
        let lo = pass - 1;
        lohi[i] = (lo, hi);
    }
    let mut st: Vec<(usize, usize)> = vec![];
    // left -> right
    for i in 0..n {
        let lo = lohi[i].0;
        while let Some(x) = st.pop() {
            if x.0 > p[i] {
                st.push(x);
                break;
            }
        }
        // binsect bias = -1
        let mut pass = st.len() + 1;
        let mut fail = 0;
        while pass - fail > 1 {
            let mid = (pass + fail) / 2;
            if st[mid - 1].1 >= lo {
                pass = mid;
            } else {
                fail = mid;
            }
        }
        tot += (st.len() + 1 - pass) as i64;
        st.push((p[i], i));
    }
    // right -> left
    st.clear();
    for i in (0..n).rev() {
        let hi = lohi[i].1;
        while let Some(x) = st.pop() {
            if x.0 > p[i] {
                st.push(x);
                break;
            }
        }
        // binsect bias = -1
        let mut pass = st.len() + 1;
        let mut fail = 0;
        while pass - fail > 1 {
            let mid = (pass + fail) / 2;
            if st[mid - 1].1 <= hi {
                pass = mid;
            } else {
                fail = mid;
            }
        }
        tot += (st.len() + 1 - pass) as i64;
        st.push((p[i], i));
    }
    
    puts!("{}\n", tot);
}

まとめ

区間系が苦手すぎる。

yukicoder No.93 ペガサス

挿入 DP の問題を集中特訓していたのに、どうしても挿入 DP 解法が思いつかず別解で通してしまった。
No.93 ペガサス - yukicoder

問題

N 点の順列 p であって以下を満たすものの個数を mod (10^9 + 7) で求めよ:

  • 1 <= i <= N - 2 なるすべての i に対し、|p[i] - p[i + 2]| != 1

制約

  • 1 <= N <= 1000

解法

包除原理を使う。
まず、ある部分集合 V \subseteq \{1, 2, \ldots, N - 2\} について、任意の  i \in V に対して |p[i] - p[i + 2]| = 1 であるようなの順列の個数を求める。このとき、i が 2 個間隔で繋がっている場所には連続した整数が入る。(例: |p[3] - p[5]| = 1 かつ |p[5] - p[7]| = 1 かつ |p[7] - p[9]| = 1 の場合、p[3], p[5], p[7], p[9] には連続した 4 個の整数が昇順または降順に入る。)このような連続した制約のことを streak と呼び、制約で言及される要素の個数を streak の長さと呼ぶ。V を決めたとき、streak の長さが a_0, a_1, ..., a_{k-1} とすると、場合の数は k! 2^{a_i の中で 2 以上のものの個数} である。
包除原理の係数は (-1)^|V| = (-1)^{n+k} であるため、(-1)^{n+k} k! 2^{a_i の中で 2 以上のものの個数} の和が求められれば良い。V の要素を 2, 4, 6, 8, ..., 1, 3, 5, 7, ... と調べて行くことで、自然に (i) streak を伸ばす (ii) 今の streak を打ち切って新しい長さ 1 の streak を始める の2遷移を持つ DP ができる。(2 * floor(N / 2) から 1 への streak を繋ぐことはできないので、その場合は (i) の遷移をなくす)
最終的に必要な状態は (場所, streak の個数, 現在の streak の長さが 2 以上か) である。計算量は O(N^2) 時間、O(N^2) 空間である。

登場する典型

  • 包除原理
  • 持つべき状態数を削減する
    • 数式を見て、現在の streak の長さの区別が 2 通りでよいこと、streak の個数が必要で局所的に計算できないこと、などを判断する (たとえば個数の部分が 2^k であれば、1 個増える遷移を行うたびに 2 倍すればよく、個数を持つ必要はない)

実装上の注意点

  • 正しい式を詰める
    • 最初 N!/(a_0!a_1!...a_{k-1}!) * 2^{a_i の中で 2 以上のものの個数} だと思ったので (なんで?)
    • 手拍子で立式するのをやめる
    • 実装する前に小さい例で試す

提出: #477257 (Rust) No.93 ペガサス - yukicoder (Rust)

// Tags: insertion-dp, inclusion-exclusion-principle
fn solve() {
    let out = std::io::stdout();
    let mut out = BufWriter::new(out.lock());
    macro_rules! puts {
        ($($format:tt)*) => (let _ = write!(out,$($format)*););
    }
    input!(n: usize);
    let (fac, _invfac) = fact_init(n + 1);
    // dp[pos][#cluster][is >= 2]
    let mut dp = vec![vec![[ModInt::new(0); 2]; n + 1]; n + 1];
    if n % 2 == 0 {
        dp[1][1][0] -= 1;
    } else {
        dp[1][1][0] += 1;
    }
    for i in 1..n {
        let cut = i == n / 2 || i == n;
        for j in 0..n {
            let val = dp[i][j][1] * 2 + dp[i][j][0];
            dp[i + 1][j + 1][0] -= val;
        }
        if !cut {
            for j in 0..n + 1 {
                let val = dp[i][j][0] + dp[i][j][1];
                dp[i + 1][j][1] += val;
            }
        }
    }
    let mut tot = ModInt::new(0);
    for i in 0..n + 1 {
        let tmp = (dp[n][i][0] + dp[n][i][1] * 2) * fac[i];
        tot += tmp;
    }
    puts!("{}\n", tot);
}

まとめ

挿入 DP をマスターしたい。

AGC043 E - Topology

本番は ABD しか解けず無事死亡した模様。物理学科の人たちが E を通していて本気でびっくりした。
E - Topology

問題

整数 N と、それぞれの部分集合  S \subseteq \{0, \ldots, N-1\} に対して 0 または 1 の整数 a[S] が与えられる。
以下の条件を満たす xy 平面内の 閉曲線 C を構成せよ。

  • 任意の  S \subseteq \{0, \ldots, N-1\} に対し、以下が成立する:
    • {(i + 0.5, 0.5) | i in S} の点を通過せずに、閉曲線 C を全て y < 0 の領域に移動できる <=> a[S] = 1

答えは以下の形式で与えること。
長さ 250000 以下の点列 (x_i, y_i) (0 <= i <= L) であって、以下を満たすもの:

  • 全て整数座標で、0 <= x_i <= N かつ 0 <= y_i <= 1
  • 隣り合う2点は距離1
  • 閉曲線である。つまり、(x_0, y_0) = (x_L, y_L)

解法

自明な条件として、a[0] = 0であれば非合理。また  S \subseteq T なる S,T に対して、 a[S] = 0 かつ a[T] = 1 の場合も非合理。
これ以外の場合に条件を満たす閉曲線 C を構成する。

与えられた集合 S に対して、{(i + 0.5, 0.5) | i in S} の点を通らずに閉曲線 C を y < 0 の領域に移動できない場合に、C を S によってブロックされるということにする。
S にはブロックされるが、 S の真の部分集合についてはブロックされないような (0,0) を基点とする閉曲線 f(S) を構成する。

  • S が 1 点集合 S = {v} の場合は、(v + 0,5, 0.5) の周りを回る閉曲線 (0, 0) -> ... -> (v, 0) -> (v, 1) -> (v + 1, 1) -> (v + 1, 0) -> ... -> (0, 0) が条件を満たす。
  • S が 2 点以上の集合の場合、適当に一点取って v とし、C1 := f(v), C2 := f(S - {v}) とおく。このとき、交換子 C1 C2 C1^{-1} C2^{-1} は条件を満たす。R^2 - {(i + 0.5, 0.5) | i in S - {v}} において C1 は 1 点とホモトピックであり、w in S - {v} としたとき R^2 - {(i + 0.5, 0.5) | i in S - {v}} において C2 = f(S - {v}) は 1 点とホモトピックである。(帰納法で示せる)

a[S] = 0 である全ての S に対して f(S) を計算し、それらを繋げたものを C とすれば、C は与えられた条件を満たす。

C の長さを見積もりたい。S が 1 点集合のとき |f(S)| <= 2N + 2 である。S が k 点集合のとき、|f(S)| <= 2(2N + 2 + f(S - {v})) であることから帰納法で |f(S)| <= (2N + 2) * (2^k - 1) が言える。これらの和の最悪ケースは、 \sum_{k = 0}^N (2N + 2) 2^k C(N, k) = (2N + 2) 3^N を上回らない。N = 8 のとき (2N + 2) 3^N = 118098 である。

登場する典型

実装上の注意点

  • バグを埋め込みやすく、出力された閉曲線が正しいかどうかもチェックするのが難しいため、小さい例 (N = 2, a = 1110 など) で慎重に意図通りか確かめる。
  • あまり場合分けをしないようにする。この手の AtCoder の構成ゲーは出力のリミットがかなりゆるく設定されていることが多いので、なあなあでも通る (は?)

提出: #11076006 (Rust)

fn conn<T: Eq + Copy + std::fmt::Debug>(ans: &mut Vec<T>, ops: &[T]) {
    assert_eq!(ans.last(), ops.first());
    for i in 1..ops.len() {
        ans.push(ops[i]);
    }
}

fn calc(n: usize, bits: usize) -> Vec<(i64, i64)> {
    if bits == 0 {
        return vec![(0, 0)];
    }
    let mut ma = 0;
    for i in 0..n {
        if (bits & 1 << i) != 0 {
            ma = i;
        }
    }
    let mut sub = calc(n, bits ^ 1 << ma);
    let mut t = vec![];
    for i in 0..ma + 1 {
        t.push((i as i64, 0));
    }
    t.push((ma as i64, 1));
    t.push((ma as i64 + 1, 1));
    for i in (0..ma + 2).rev() {
        t.push((i as i64, 0));
    }
    if bits == 1 << ma {
        return t;
    }
    let mut ans = vec![(0, 0)];
    // commutator
    conn(&mut ans, &sub);
    conn(&mut ans, &t);
    sub.reverse();
    t.reverse();
    conn(&mut ans, &sub);
    conn(&mut ans, &t);
    ans
}

fn solve() {
    let out = std::io::stdout();
    let mut out = BufWriter::new(out.lock());
    macro_rules! puts {
        ($($format:tt)*) => (let _ = write!(out,$($format)*););
    }
    input! {
        n: usize,
        a: chars,
    }
    for i in 0..1 << n {
        for j in 0..1 << n {
            if (i & j) == i && a[i] == '0' && a[j] == '1' {
                puts!("Impossible\n");
                return;
            }
        }
    }
    if a[0] == '0' {
        puts!("Impossible\n");
        return;
    }
    let mut ans = vec![(0, 0)];
    for i in 1..1 << n {
        if a[i] == '1' {
            continue;
        }
        conn(&mut ans, &calc(n, i));
    }
    assert!(ans.len() - 1 <= 250_000);
    puts!("Possible\n");
    puts!("{}\n", ans.len() - 1);
    for (x, y) in ans {
        puts!("{} {}\n", x, y);
    }
}

まとめ

何で群論知っていながら詳細詰めるのに1時間半もかかったの?

第6回 ドワンゴからの挑戦状 予選 C - Cookie Distribution

かなり好きなタイプの問題。
C - Cookie Distribution

問題

N 人の人にクッキーを配りたい。i = 1, ..., K に対して、以下の操作を順番に行う。

  • N 人の中から a[i] 人を等確率にランダムに選び、その人たちにクッキーを 1 個ずつ渡す。

最終的に 人 i がもらったクッキーの枚数を c[i] とし、 c[1] * c[2] * ... c[N] を 嬉しさと呼ぶことにする。
嬉しさの期待値に C(N, a[1]) * ... * C(N, a[K]) を掛けたものは整数になるが、これを mod (10^9 + 7) で求めよ。

制約

  • 1 <= K <= 20
  • 1 <= N <= 1000
  • 1 <= a[i] <= N

解法

dp[i][j] := (i 回目まで処理を行った後の、c[1] * ... c[j] の期待値 (j = 0 の場合は 1)) と置く。
dp[i + 1][j] を dp[i][...] を使って表したい。ここで、c[1] * ... * c[j] を数える代わりに、人 1 から 人 j までがもらったクッキーの j 個組を数えることにする。 (両者は同じものである)
a[i + 1] 個のクッキーを配った後のこの量は、0 <= s <= j となる s ごとに以下の値を求めて合計すればわかる:

  • j 個のうち s 個はもともと配られていたもの、残りの (j - s) 個は今配られた a[i + 1] 個のうちのどれかとしたときの、j 個組の個数

j 個のうち s 個の場所の選び方が C(j, s) 通り、個数が dp[i][s] (対称性より)、a[i + 1] 個のクッキーのうち (j - s) 個が残りの場所に配られる確率は C(n - (j - s), n - a[i + 1]) / C(n, a[i + 1]) である。これらを全て掛けて dp[i][s] * C(j, s) * C(n - (j - s), n - a[i + 1]) / C(n, a[i + 1]) が求める量である。
これを全ての s について足し合わせれば良い。計算量は、DP テーブルのエントリ数が O(KN) 個、各エントリの計算に O(N) 時間かかるので、合計 O(KN^2) である。

なお、dp[i][s] -> dp[i + 1][j] の遷移に登場する変数が s, j - s の形でしか登場しないので、任意 mod の畳み込みを用いることで O(KN log N) にすることもできる。(参考: AtCoder AGC #005 : F - Many Easy Problems - kmjp's blog)

登場する典型

  • 積の期待値を求めたい場合に、組合せ論的に扱いやすいものに置き換えて考える
  • 積の期待値 = 期待値の積 ではないので注意
    • 今回の場合、c[i] 同士は全然独立ではないので、この式は成り立たない

実装上の注意点

とくになし

提出: #9422487 (Rust)

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, k: usize,
        a: [usize; k],
    }
    let (fac, invfac) = fact_init(n + 1);
    let comb = |x, y| {
        if x < y {
            ModInt::new(0)
        } else {
            fac[x] * invfac[y] * invfac[x - y]
        }
    };
    let mut dp = vec![vec![ModInt::new(0); n + 1]; k + 1];
    dp[0][0] += 1;
    for i in 0..k {
        for j in 0..n + 1 {
            let val = dp[i][j];
            dp[i + 1][j] += val;
            for l in 0..j {
                if a[i] + l >= j {
                    let val = dp[i][l] * comb(j, l) * comb(n + l - j, a[i] + l - j);
                    dp[i + 1][j] += val * invfac[n] * fac[a[i]] * fac[n - a[i]];
                }
            }
        }
    }
    let mut ans = dp[k][n];
    for i in 0..k {
        ans *= comb(n, a[i]);
    }
    puts!("{}\n", ans);
}

まとめ

個人的には D が比較的難しかった (solved 数を見る限り C > D だったらしいが)。おかげで E にあまり時間が掛けられなかった。ゆるせね〜 (責任転嫁)

yukicoder No.963 門松列列(2)

最近流行りの例のアレ。
No.963 門松列列(2) - yukicoder

問題

長さ N の交代順列の個数を、mod 1012924417 (= 483 * 2^21 + 1) で割ったあまりを求めよ。
交代順列: 1 から N の並び替えであって、任意の i (2 <= i <= N - 1) に対して a[i - 1] < a[i] と a[i] > a[i + 1] が同値であるもの。

  • 2 <= N <= 202020

解法

N=i のときの問題の答えを dp[i] と置き、そのうち左端で a[1] < a[2] が成り立っているものだけの個数を dp2[i] と置く。n >= 2 のとき dp2[n] = dp[n] / 2 であり、0 <= n <= 1 のときは dp2[n] = dp[n] = 1 とする。
そうすると、挿入 DP の考え方で、n >= 2 のとき dp[n] = \sum_{j = 0}^{n - 1} dp2[j] * dp2[n - 1 - j] * C(n - 1, j) が言える。(n を j 番目に入れるとき、n の左側に置く数の選び方は C(n - 1, j) 通り、n の左右の並べ方はそれぞれ dp2[j], dp2[n - 1 - j] 通り。)
ep[i] := dp2[i] / i! と置くことにすると、n >= 2 のとき ep[n] = \sum_{j = 0}^{n - 1} ep[j] * ep[n - 1 - j] / n が言える。これは、左側からだんだん値が決まっていく畳み込みなので、分割統治しながら FFT をしていけば解ける。(O(N log^2 N) 時間)
さらに考察を進めることもできる。E := \sum_{n = 0}^{\infty} ep[i] x^i と置く。このとき、漸化式を E の式で表しなおすと、E = 1 + x/2 + \int_0^x E(t)^2/2 dt となる。両辺を微分すると  E' = (1 + E^2) / 2 となり、これを初期条件 E(0) = 1 の下で解くと E(x) = tan(x/2 + pi/4) = (1 + sin x) / cos x となる。sin, cos のテイラー展開は簡単にわかり、形式的冪級数の商は O(N log N) 程度で計算する方法があるので、全体で O(N log N) 時間で解ける。

登場する典型

実装上の注意点

  • FFT で元にする範囲、計算した結果を書き込む範囲に気をつける。今回の場合は以下の2つ:
    • [0, x) * [0, x) を [x, 2x) に書き込む
    • 2x <= lo のときに [0, 2x) * [lo, lo + x) の2倍を [lo + x, lo + 2x) に書き込む。2倍するのは [0, 2x) 側が左右両方に現れうるためで、2x まで見ないといけないのは lo -> lo + 2x - 1 みたいな遷移があり得るため。なお、FFT の幅は 2x でよい。長さ 2x のところで1周するが、1周した後の値は x 番目の左側にしか来ないため。
  • 初期値に気をつける
    • 今回は ep[0] = ep[1] = 1

提出: #414769 (Rust)
(動的 FFT)

fn rec(lo: usize, hi: usize, dp: &mut [ModInt], ep: &mut [ModInt],
       fac: &[ModInt], invfac: &[ModInt]) {
    let n = hi - lo;
    debug_assert!(n.is_power_of_two());
    if (lo, hi) == (0, 2) {
        dp[0] += 1;
        dp[1] += 1;
        ep[0] += 1;
        ep[1] += 1;
        return;
    }
    if n == 1 {
        return;
    }
    let mid = (lo + hi) / 2;
    rec(lo, mid, dp, ep, fac, invfac);
    // FFT
    let zeta = ModInt::new(5).pow((MOD - 1) / n as i64);
    let mut tmp = vec![ModInt::new(0); n];
    let mut tmp2 = vec![ModInt::new(0); n];
    for i in lo..mid {
        tmp[i - lo] = ep[i];
    }
    // Difference can be anything in [1, n - 1].
    for i in 0..n {
        tmp2[i] = ep[i];
    }
    fft::transform(&mut tmp, zeta, 1.into());
    fft::transform(&mut tmp2, zeta, 1.into());
    let mut invn = ModInt::new(n as i64).inv();
    // If not overlapping, multiply by two.
    if lo != 0 {
        invn *= 2;
    }
    for i in 0..n {
        tmp[i] = tmp[i] * tmp2[i] * invn;
    }
    fft::transform(&mut tmp, zeta.inv(), 1.into());
    for i in mid..hi {
        dp[i] += tmp[i - lo - 1] * fac[i - 1];
        ep[i] += tmp[i - lo - 1] * fac[i - 1] * invfac[i] * invfac[2];
    }
    rec(mid, hi, dp, ep, 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);
    const W: usize = 1 << 18;
    let (fac, invfac) = fact_init(W);
    let mut dp = vec![ModInt::new(0); W];
    let mut ep = vec![ModInt::new(0); W];
    rec(0, W, &mut dp, &mut ep, &fac, &invfac);
    //debugln!("{:?}", dp);
    //debugln!("{:?}", ep);
    puts!("{}\n", dp[n]);
}

まとめ

solved 数少なすぎ && そろそろ形式的冪級数ライブラリを整備しなければ…。

yukicoder No.951 【本日限定】1枚頼むともう1枚無料!

典型のはずなのにかなりハマった。
No.951 【本日限定】1枚頼むともう1枚無料! - yukicoder

問題

N 種類のピザがあり、i 番目のピザのコストは P[i] 円、美味しさは D[i] である。各ピザはちょうど1枚だけ売られている。以下の行動を好きな回数繰り返すことができる。

  • ピザを買う。そのピザの値段を P[i] 円としたとき、P[i] 円以下の残っているピザを一つ選び、0 円で買うことができる。

買ったピザの美味しさの和を最大にしたい。支払いを K 円以内にするとき、実現できる美味しさの和の最大値を求めよ。

解法

最適戦略としてどのようなものがありえるかを考察する。買ったピザの値段を高い順に並べたものを X_0, X_1, ..., X_{a - 1} とすると、この中で 0 円にすべきなのは X_1, X_3, ... である。値段の低い順番に見た時、以下のようなルールでピザを選んでいくのと同値である。

  • 「正規の金額でピザを買う」「0 円でピザを買う」「ピザを無視する」の3種類の行動を取ることができる。「0 円でピザを買う」は「0 円にできる権利」を使わないと選べない。最初「0 円にできる権利」を 1 回使うことができる。「0 円にできる権利」を使うと 1 個減り、正規の金額でピザを買うと「0 円にできる権利」は 1 回使える状態に戻る。最終的に「0 円にできる権利」は 1 回使える状態になっていないといけない。(最後の購入は正規の金額でなければならない)

0 円にできる権利はどの地点でも最大 1 回しか残っていない。よって、状態は [現在位置][使った金額][権利の残り回数] の O(NK) 状態である。

登場する典型

  • 最適戦略の形を考える
  • ナップサック DP などの、更新の方向が決まっている DP について、次元を一つ落とす
    • 遷移が DP[i][j] -> DP[i + 1][k] (常に j < k) の形をしている場合、逆からループすることで1次元で済ませることができる。
    • 今回の場合は、(使った金額, 権利の残り回数) の辞書順で考えるとよい

実装上の注意点

  • (P[i], D[i]) のソートを忘れない

提出:
#413247 (Rust, 高速化前)
#413250 (Rust, 高速化後)

(高速化後)

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, k: usize,
        pd: [(usize, i64); n],
    }
    let mut pd = pd;
    pd.sort();
    const INF: i64 = 1 << 50;
    let mut dp = vec![[-INF; 2]; k + 1];
    dp[0] = [0, 0];
    // (j, 0) -> (j, 1), (j, 1) -> (j + p, 0) の遷移を一挙にやりたい。
    // (j, l) が大きくなる方向にしか遷移しないため、
    // (j, l) の辞書順の逆に遷移させればよい.
    for i in 0..n {
        let (p, d) = pd[i];
        for j in (0..k + 1).rev() {
            if j + p <= k {
                dp[j + p][0] = max(dp[j + p][0], dp[j][1] + d);
            }
            dp[j][1] = max(dp[j][1], dp[j][0] + d);
        }
    }
    let mut ma = 0;
    for i in 0..k + 1 {
        ma = max(ma, dp[i][0]);
    }
    puts!("{}\n", ma);
}

まとめ

こういう DP 配列を潰す系の高速化は、常識な気もする一方でどこにも書かれていない気がする。誰かまとめてほしい