DISCO presents ディスカバリーチャンネル コードコンテスト2019 本戦 D - DISCO!

まあ実家だよね。
D - DISCO!

問題

'D', 'I', 'S', 'C', 'O'のみからなる文字列Sが与えられる。このとき、以下のQ個の質問に答えよ:

  • SのL_i文字目からR_i文字目までからなる部分文字列を考える。この文字列に部分列として含まれる"DISCO"は何通りあるか? つまり、単調増加な添字の5つ組であって、それぞれの位置の文字を並べると"DISCO"になるようなものは何通りあるか? 答えを2^32で割った余りを求めよ。

解法

まずは「Sの"DISCO"という部分文字列を数え上げる問題」を考えよう。これは以下のような6 * |S|エントリのDPでできる:

  • DP[i][j]: i文字目まで見たときに、"DISCO"の何文字目までできているか (0 <= i <= |S|, 0 <= j <= 5)
  • 遷移: i文字目を取るか取らないかの2通り。ただしi文字目を取れるのはS[i] == "DISCO"[j]のときだけ。

こうすると、例えばS[i] = 'D'の時の遷移は

(DP[i][0] ... DP[i][5]) =
(DP[i - 1][0] ... DP[i - 1][5])
(1 0 0 0 0 0
1 1 0 0 0 0
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1)

と、行ベクトルに右から行列を掛ける形で表せる。最終的に欲しい値は、(DP[0][0] ... DP[0][5]) = (1 0 0 0 0 0)としたときのDP[|S|][5]である。
よって、文字S[i]に対応する行列をA[i]と置くと、(1 0 0 0 0 0)A[1]A[2] ... A[|S|]の第5成分が欲しい値である。
今は1文字目から|S|文字目までの遷移を考えたが、欲しいのは行列積なので、任意の区間[L, R]に対しA[L] A[L + 1] ... A[R] (の0行5列成分)が分かればよい。これは行列をセグメント木に載せればできる。計算量は愚直にやるとO(K^3 (|S| + Q) log |S|)であり、制限時間が長いのでこれで十分である。

登場する典型

  • DPの遷移が線形変換であるときに、高速化をする定跡
    • 行列累乗、複数の変換を行列積に置き換えることでできる高速化(セグメント木に載せるなど)

実装上の注意点

  • 制限時間が13秒と長いので、間に合う中で一番簡単な実装を選択する。今回の場合、登場する行列が全て可逆なので、行列積バージョンの累積和を取るO(K^2(|S| + Q))の解法を書くこともできるが、セグメント木で殴るO(K^3 (|S| + Q) log |S|)の解法でもよい。
    • なお、遷移を両側から行うために、"DISCO"の部分文字列10通りについて遷移を手書きする方法があるが、より簡単な実装を考えるべきである。
  • Rustでオーバーフローする演算を行いたいときは、core::num::Wrappingを使うか基本型のwrapping_***メソッドを使う。
    • 間違えてi32で演算しないようにする(1敗)

提出: Submission #4039621 - DISCO presents ディスカバリーチャンネル コードコンテスト2019 本戦 (Rust)

// SegTree省略

extern crate core;
use core::num::Wrapping;
const SZ: usize = 6;
type Mat = [[Wrapping<u32>; 6]; 6];
 
fn mul(a: Mat, b: Mat) -> Mat {
    let mut ans = [[Wrapping(0); 6]; 6];
    for i in 0 .. SZ {
        for j in 0 .. SZ {
            for k in 0 .. SZ {
                ans[i][k] += a[i][j] * b[j][k];
            }
        }
    }
    ans
}
 
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 mut unit: Mat = [[Wrapping(0); 6]; 6];
    for i in 0 .. SZ { unit[i][i] = Wrapping(1); }
    let mut st = SegTree::new(n, |x, y| mul(x, y), unit);
    for i in 0 .. s.len() {
        let idx = match s[i] {
            'D' => 0,
            'I' => 1,
            'S' => 2,
            'C' => 3,
            'O' => 4,
            _ => unreachable!(),
        };
        let mut m = unit;
        m[idx][idx + 1] = Wrapping(1);
        st.update(i, m);
    }
    for (l, r) in lr {
        let ans = st.query(l, r);
        puts!("{}\n", ans[0][5]);
    }
}

まとめ

この手の解ける問題で方針をミスらなかったり実装を単純化したりするのは極めて大事。

ARC058 F - 文字列大好きいろはちゃん / Iroha Loves Strings (解説を見た)

実装が大変すぎる。
F - 文字列大好きいろはちゃん / Iroha Loves Strings

問題

文字列の列s_1, ..., s_Nが与えられる。これの部分列を順序を保ったまま連結して、長さKの文字列を作ることを考える。ありえる文字列の辞書順最小値を求めよ。
制約

  • 1 <= N <= 2000
  • 1 <= K <= 10000
  • 1 <= |s_i| <= K
  • \sum |s_i| <= 10^6
  • アルファベットの大きさは26

解法

以下の解法では、辞書順を以下で定義する。

  • どちらかがどちらかのprefixである場合、長い方が小さい。
  • そうでない場合は普通の辞書順。

O(NK^2)のDPは色々考えられるが、どれもTLEする。そのため、うまいDPを考える必要がある。
以下のようにすると都合の良い性質が得られる。

  • DP[i][j]: i番目の文字列までを見て、j文字の文字列を作る時の最小値。ただし、(i+1)番目以降の文字列をどう使っても最終的に長さKの文字列を作れない場合は無効な値とする。

こうすることで、DP[i][j]のうち辞書順最小(ただし長い方が小さい)のものをCS[i]とおけば、CS[i]のprefix以外から遷移するのは無駄であるため、各jについてDP[i][j]はCS[i]のprefixか無効な値かのいずれかである。空間計算量はO(NK)である。
この上の遷移CS[i] -> CS[i + 1], DP[i][..] -> DP[i + 1][..]を高速化する方法を考える。ここで必要になるのは以下のタイプの比較である。

  • 添字u, vおよびブール値f1, f2が与えられた時、CS[i][..u] + (f1 ? S[i] : "") < CS[i][..v] + (f2 ? S[i] : "")か?

これはS[i] + CS[i]にz_algorithmを適用しておけば高速に計算できる。比較は以下のように考えるとわかりやすい多少はマシになる。

  • S[i]が先頭に来ているので、S[i]とS[i]のsuffixおよびCS[i]のsuffixのLCP(最長共通prefix)の長さは高速に計算できる。
  • それを利用して、S[i]とS[i]のスライスおよびCS[i]のスライスの辞書順比較も高速に計算できる(下記コードのcmp_slice)
  • 「CS[i][..u] < CS[i][..v] + S[i]か?」や「CS[i][..u] + S[i] < CS[i][..v] + S[i]か?」という問題を解く時には、先頭から見ていってどちらかのS[i]との境目で区切って比較する
    • 例: u > vのときCS[i].[..u] < CS[i][..v] + S[i]か? という問いに答えることを考える。先頭v文字は同一。なのでCS[i][v..u]とS[i]の比較結果が分かれば良い。これは上のデータ構造を使えば高速にできる。
    • 例2: u < v, u + |S[i]| > vのときCS[i][..u] + S[i] < CS[i][..v] + S[i]か?という問いに答えることを考える。先頭u文字は同一。u..v文字目を見て(<=> S[i][..v-u]とCS[i][u..v]を比較して)等しくなければそれが答えで、等しい場合はさらにv文字目以降を見る(<=> S[i][v-u..]とS[i]を比較する)。全ての比較がS[i]のprefixとの比較になっていることに注意。

上を頑張って実装するとACできる。計算量はO(NK)である。

登場する典型

実装上の注意点

  • 力任せに実装するとバグらせたり頭が爆発したりするので、極力共通化・単純化する
  • DPの更新を2ステップに分ける都合上、それぞれで必要な比較が微妙に違うので注意!
    • 1回目はCS[i + 1]を求めるためのものなので、辞書順で小さくなるかどうかを判定する (cmpsの返り値の0番目)
    • 2回目はDP[i][j]を求めるためのものなので、CS[i + 1]のprefixであるかどうかを判定する (cmpsの返り値の1番目)

提出: Submission #4032098 - AtCoder Regular Contest 058 (Rust)

/*
 * Z algorithm. Calculates an array a[i] = |lcp(s, s[i...|s|])|,
 * where s is the given string.
 * If n = s.length(), the returned array has length n + 1.
 * E.g. z_algorithm("ababa") = {5, 0, 3, 0, 1, 0}
 * Reference: http://snuke.hatenablog.com/entry/2014/12/03/214243
 * Verified by: AtCoder ARC055-C (http://arc055.contest.atcoder.jp/submissions/1061771)
 */
fn z_algorithm<T: PartialEq>(s: &[T]) -> Vec<usize> {
    let n = s.len();
    let mut ret = vec![0; n + 1];
    ret[0] = n;
    let mut i = 1; let mut j = 0;
    while i < n {
        while i + j < n && s[j] == s[i + j] { j += 1; }
        ret[i] = j;
        if j == 0 { i += 1; continue; }
        let mut k = 1;
        while i + k < n && k + ret[k] < j {
            ret[i + k] = ret[k];
            k += 1;
        }
        i += k; j -= k;
    }
    ret
}

fn calc_poss(s: &[Vec<char>], k: usize) -> Vec<Vec<bool>> {
    let n = s.len();
    let mut poss = vec![vec![false; k + 1]; n + 1];
    poss[n][0] = true;
    for i in (0 .. n).rev() {
        let len = s[i].len();
        for j in 0 .. k - len + 1 {
            poss[i][j + len] |= poss[i + 1][j];
        }
        for j in 0 .. k + 1 {
            poss[i][j] |= poss[i + 1][j];
        }
    }
    poss
}

fn solve() {
    use std::cmp::Ordering::*;
    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,
        s: [chars; n],
    }
    let poss = calc_poss(&s, k);
    let mut cs = vec![Vec::new(); n + 1];
    let mut dp = vec![vec![false; k + 1]; n + 1];
    dp[0][0] = true;
    for i in 0 .. n {
        let len = s[i].len();
        let mut css = s[i].clone();
        css.extend_from_slice(&cs[i]);
        let zarr = z_algorithm(&css);
        // cs[i + 1] is implicitly represented by
        // cs[i][..cs_len] + if s_appended { s[i] } else { "" },
        // where (cs_len, s_appended) = cs_next.
        let mut cs_next = (0, false);
        // compares css[a .. a + blen] and css[b .. b + blen].
        // Returns (comparison, is_prefix).
        // If one is a prefix of the other, the longer one is smaller.
        // Either a or b must be 0.
        let cmp_slice = |a: usize, alen: usize, b: usize, blen: usize| {
            let common_len = match (a, b) {
                (0, _) => zarr[b],
                (_, 0) => zarr[a],
                _ => unreachable!(),
            };
            if min(alen, blen) <= common_len {
                (blen.cmp(&alen), true) // longer is smaller
            } else {
                (css[a + common_len].cmp(&css[b + common_len]), false)
            }
        };
        // Compares cs[i][..j] + (if fappend { s[i] } else { "" }) and the current cs[i + 1].
        // If one is a prefix of the other, the longer one wins (is smaller).
        let cmps = |(j, fappend): (usize, bool),
        (cs_len, s_appended): (usize, bool)| {
            match (fappend, s_appended) {
                (true, true) if j == cs_len => (Equal, true),
                (true, true) if j > cs_len + len =>
                    cmp_slice(len + cs_len, j - cs_len, 0, len),
                (true, true) if cs_len > j + len =>
                    cmp_slice(0, len, len + j, cs_len - j),
                (true, true) if j > cs_len => {
                    let res = cmp_slice(len + cs_len, j - cs_len, 0, j - cs_len);
                    if res.0 != Equal { return res }
                    cmp_slice(0, len, j - cs_len, len - (j - cs_len))
                },
                (true, true) if cs_len > j => {
                    let res = cmp_slice(0, cs_len - j, len + j, cs_len - j);
                    if res.0 != Equal { return res }
                    cmp_slice(cs_len - j, len - (cs_len - j), 0, len)
                },
                (true, true) => unreachable!(),

                (true, false) if cs_len <= j => (Less, true),
                (true, false) => cmp_slice(0, len, len + j, cs_len - j),

                (false, true) if j <= cs_len => (Greater, true),
                (false, true) => cmp_slice(len + cs_len, j - cs_len, 0, len),

                (false, false) => (cs_len.cmp(&j), true),
            }
        };
        // get min
        for j in 0 .. k + 1 {
            if !poss[i + 1][k - j] || !dp[i][j] { continue; }
            if cmps((j, false), cs_next).0 == Less {
                cs_next = (j, false);
            }
        }
        for j in 0 .. k - len + 1 {
            if !poss[i + 1][k - j - len] || !dp[i][j] { continue; }
            if cmps((j, true), cs_next).0 == Less {
                cs_next = (j, true);
            }
        }
        // update dp
        for j in 0 .. k + 1 {
            if !poss[i + 1][k - j] || !dp[i][j] { continue; }
            dp[i + 1][j] |= cmps((j, false), cs_next).1;
        }
        for j in 0 .. k - len + 1 {
            if !poss[i + 1][k - j - len] || !dp[i][j] { continue; }
            dp[i + 1][j + len] |= cmps((j, true), cs_next).1;
        }
        cs[i + 1] = cs[i][..cs_next.0].to_vec();
        if cs_next.1 {
            cs[i + 1].extend_from_slice(&s[i]);
        }
    }

    puts!("{}\n", cs[n][..k].iter().cloned().collect::<String>());
}

まとめ

実装の簡略化を頑張ったが、それでもかなり複雑になってしまった。実装が上手く行かず鬱病になったが、某氏にこれは2度と実装したくない問題worst5に入ると言われたので気分が楽になった(は?)

KEYENCE Programming Contest 2019 E - Connecting Cities

色々な解き方があって本当に面白い問題だと思う。
E - Connecting Cities

問題

整数N, Dと長さNの整数列A_iが与えられる。以下の完全グラフの最小全域木の重みを求めよ。

  • 頂点数はN
  • 頂点iと頂点jを結ぶコストはD|i - j| + A_i + A_j

解法

Code Festival 2017 J (Tree MST) と 距離空間上のボロノイ図 - とこはるのまとめの方針に従い、2N頂点のグラフを作り、そのグラフのヴォロノイ図を計算することで解を求める。
頂点1, ..., 2Nを持つグラフに以下のように辺を張る:

  • (i, i + N)間に重みA_i (1 <= i <= N)
  • (i + N, i + N + 1)間に重みD (1 <= i <= N - 1)

このグラフのヴォロノイ図を求めればよい。以下のようなことを考えながら実装する。

  • 頂点1, ..., Nから色1, ..., Nの液体が同じ速さで流れ始める
  • 色の違う液体は触れると触れた場所がその場で固まり、互いに混じり合わない。触れた場所以外はそのまま流れ続ける。
  • 最終的に全ての液体の動きが止まった時の液体の配置がヴォロノイ図である

これを実装するとダイクストラ法に似た実装になる。以下のデータを管理すれば良い。

  • 「色iの液体がt秒後に頂点vに到達した」というイベント。時系列順に処理するので優先度付きキューがよい。
  • 「頂点vは色iで塗られた、あるいはどの色でも塗られていない」という状態。色に順番があることにすれば、一つの頂点が2色以上で同時に塗られることはないと思って良い。

最終的なヴォロノイ図において、色と色の境界はちょうどN-1個である。それぞれの境界に対応する辺を集めたものがMSTである。
計算量は優先度付きキューを使うところが一番重く、O(N log N)である。

登場する典型

実装上の注意点

  • 今回は境界がちょうどN-1個なので、Union-Find木などを使う必要はない
  • 2色が合流するイベントを適切に処理する。色の順番を見て同じイベントを2回数えないようにするのもよいし、あえて2回数えて後で2で割るのもよい。

提出: Submission #4014722 - KEYENCE Programming Contest 2019 / キーエンス プログラミング コンテスト 2019 (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,
        d: i64,
        a: [i64; n],
    }
    let mut g = vec![Vec::new(); 2 * n];
    for i in 0 .. n {
        g[i].push((n + i, a[i]));
        g[n + i].push((i, a[i]));
    }
    for i in 0 .. n - 1 {
        g[n + i].push((n + i + 1, d));
        g[n + i + 1].push((n + i, d));
    }
    let mut que = BinaryHeap::new();
    const INF: i64 = 1 << 50;
    let mut dp = vec![(INF, 2 * n); 2 * n];
    for i in 0 .. n {
        que.push((0, i, i));
    }
    let mut tot = 0;
    while let Some((dist, v, w)) = que.pop() {
        let dist = -dist;
        if dp[w].0 == INF {
            dp[w] = (dist, v);
            for &(u, cost) in g[w].iter() { 
                que.push((-dist - cost, v, u));
            }
        } else {
            // dp[w].1 and v meets
            if v != dp[w].1 {
                tot += dp[w].0 + dist;
            }
        }
    }
    puts!("{}\n", tot / 2);
}

まとめ

最小全域木は構造が綺麗なので、解法がたくさんあるのが面白い。(この他にも想定解法の枝刈りKruskal、Boruvka、分割統治Kruskalなどがある。)

KEYENCE Programming Contest 2019 F - Paper Cutting

母関数おじさんにとっては多分実家。
F - Paper Cutting

問題

長方形の中に縦線がH個、横線がW個ある。これからK本の辺を選んで、順番にその線に従って長方形をカットする。
各カットについて、そのスコアはカット後の分割された長方形の個数である。P(H + W, K)通りあるカットの方法について、スコアの合計の総和を求めよ。

解法

全てのカットが一様ランダムに選ばれるとした時の、スコアの合計の期待値を求めれば良い。
和の期待値は期待値の和なので、i回目のカットで縦線/横線を選んだときのスコアの期待値が分かれば良い。(iを0-indexedと思うことにして、0<=i<=K-1とする。)
i回目のカットが縦の場合を考える。i回目のカットの前に縦にj回カット、横に(i-j)回カットされていたとする。そのときの確率はC(H-1,j)*C(W,i-j)/C(H+W-1,i)で、スコアは(j+2)(i-j+1)である。これの積を全てのj (0 <= j <= i)について足し合わせれば良い。
ここで、係数がjの多項式なので、これを母関数を使って表すことにすると、C(H-1,j)(j+2)は(1/x) * (d/dx((1+x)^{h-1}x^2) = (1+x)^{h-2}((h-1)x+2(1+x))におけるx^jの係数、C(W,u)*(u+1)は(d/dx)((1+x)^w x) = (1+x)^{w-1}(wx + 1 + x)におけるx^uの係数である。求めたい和は畳み込みなので、結局これらの積(1+x)^{h+w-3}((h-1)x+2(1+x))wx + 1 + x)におけるx^iの係数である。これはあらかじめ階乗や階乗の逆数のテーブルを用意しておけば定数時間で計算できる。
i回目のカットが横の場合も同様である。i回目のカットが縦/横になる確率はそれぞれH/(H+W), W/(H+W)なので、これを掛けた値を合計すれば良い。

登場する典型

  • 組み合わせの問題を母関数を使って表す

実装上の注意点

  • H, Wが10^7とかなり大きいので、ループの中にO(log mod)の処理を乗せてはいけない!
    • mod逆元など
    • 定数時間でもmodは重いので、極力減らす!

提出:
Submission #4009374 - KEYENCE Programming Contest 2019 / キーエンス プログラミング コンテスト 2019
(Rust)

// ModInt省略

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! {
        h: usize,
        w: usize,
        k: usize,
    }
    let (fac, invfac) = fact_init(h + w + 1);
    let comb = |x, y| {
        if x < y { ModInt::new(0) }
        else { fac[x] * invfac[x - y] * invfac[y] }
    };
    let invcomb = |x, y| {
        if x < y { ModInt::new(0) }
        else { invfac[x] * fac[x - y] * fac[y] }
    };
    let mut tot = ModInt::new(0);
    let hp = ModInt::new(h as i64) * ModInt::new((h + w) as i64).inv();
    let wp = ModInt::new(w as i64) * ModInt::new((h + w) as i64).inv();
    for i in 0 .. k {
        let c0 = comb(h + w - 2, i) * 2;
        let c1 = if i > 0 { comb(h + w - 2, i - 1) } else { 0.into() };
        let c2 = if i > 0 && h + w >= 3 { comb(h + w - 3, i - 1) } else { 0.into() };
        let c3 = if i > 1 && h + w >= 3 { comb(h + w - 3, i - 2) } else { 0.into() };
        let fac = invcomb(h + w - 1, i);
        let mut coef = c0;
        if i > 0 {
            coef += c1 * (w as i64 + 1) * 2;
            if h > 0 {
                coef += c2 * (h as i64 - 1);
                if i > 1 {
                    coef += c3 * (h as i64 - 1) * (w as i64 + 1);
                }
            }
        }
        coef *= fac;
        tot += coef * hp;
        let mut coef = c0;
        if i > 0 {
            coef += c1 * (h as i64 + 1) * 2;
            if w > 0 {
                coef += c2 * (w as i64 - 1);
                if i > 1 {
                    coef += c3 * (h as i64 + 1) * (w as i64 - 1);
                }
            }
        }
        coef *= fac;
        tot += coef * wp;
    }
    puts!("{}\n", tot * fac[h + w] * invfac[h + w - k]);
}

まとめ

Eで1時間くらい使ってしまったせいで、これを40分で解けたのにコンテスト終了に間に合わせられなかった。

ARC072 F - Dam

1.7年前(2017/4/22)の本番では解法の見当もつかなかったが、改めて考えてみたら分かった。
F - Dam

問題

数Lと数列v_1, ..., v_N, t_1, ..., t_Nが与えられる。このとき、各1 <= i <= Nに対して以下の問いに答えよ。

  • 許容量Lリットルのダムに水を流入させることを考える。1日目にt_1度の水v_1リットルが流入し、その後水を任意量捨て、2日目にt_2度の水v_2リットルが流入し、…ということをi日目まで行う。各流入後に水を捨てる量をうまく調節して、i日目の流入のあとダムに水がLリットルある状態での温度を最大化したい。最大値を求めよ。

1 <= N <= 5 * 10^5
1 <= L <= 10^9
1 <= v_i <= 10^9
0 <= t_i <= 10^9

解法

i日目に水を捨てたあとの状態で、水がvリットル残っている状態での温度の最大値をf(i, v)と置くことにする。f(i, v)はvについて単調減少である。
g(i, v)をi日目に水を捨てる前の温度の最大値とすると、g(i, v) = (f(i - 1, v - v_i) * (v - v_i) + t_i * v_i) / v, f(i, v) = max_{w >= v} g(i, w)である。
ここで、f, gの代わりにp(i, v) = vf(i, v), q(i, v) = vg(i, v)を考えることにすると、q(i, v) = p(i - 1, v - v_i) + t_i * v_i なのでq(i, v)のグラフはp(i - 1, v)のグラフに(0, 0)と(v_i, t_i * v_i)を端点とする線分を連結させたものになる。またg(i, v)からf(i, v)を作るときに行われるchmaxのような処理は、p, qを使って表すとp(i, v) = v * max_{w >= v} (q(i, w) / w)である。qからpを作る操作は原点から見ていって傾きが増加したらグラフを書き換える、という操作を貪欲にやることで行える。図にすると以下のようになる。(Desmos | Graphing Calculatorで作成)

f:id:koba-e964:20190112160034p:plain
chmax前
f:id:koba-e964:20190112160010p:plain
chmax後
この操作を効率的にやるためには、dequeもしくはstackがあればよい。
計算量はO(N)である。

登場する典型

  • 最適値をパラメータを使って表す

実装上の注意点

  • stackを使って実装することもできるが、解説のようにdequeを使う方がよい(今回はstackを使った)

提出: Submission #3980185 - AtCoder Regular Contest 072 (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,
        l: f64,
        tv: [(f64, f64); n],
    }
    // Let x be the sum of v for all entry (v, _, _).
    // An entry (v, p, acc) means
    // (x - v, acc - p) and (x, acc) are connected by a segment.
    // We are to find max_y{(f(x) - f(x - y)) / y}.
    let mut pool: Vec<(f64, f64, f64)> = Vec::new();
    // Sentinel. l liter of water with temperature -1.0
    pool.push((l, -1.0, -l));
    // The pointer is at pool[cur].0 + pos
    // This pointer should always be at x - l, where x is described above.
    let mut cur = 0;
    let mut pos = 0.0;
    for (t, v) in tv {
        let mut last = pool[pool.len() - 1].2;
        pool.push((v, t * v, last + t * v));
        // Advances pointer by v
        let mut rem = v;
        while rem > 0.0 {
            if pool[cur].0 - pos < rem {
                rem -= pool[cur].0 - pos;
                cur += 1;
                pos = 0.0;
            } else {
                pos += rem;
                break;
            }
        }
        assert!(pos <= pool[cur].0 + 1e-9);
        // Greedy chmax
        while pool.len() >= 2 {
            let (v1, p1, acc1) = pool[pool.len() - 1];
            if v1 >= l { break; }
            let (v2, p2, _acc2) = pool[pool.len() - 2];
            if p2 / v2 > p1 / v1 {
                pool.pop();
                pool.pop();
                if v1 + v2 >= l {
                    pool.clear();
                    let newp = p1 + p2 / v2 * (l - v1);
                    pool.push((l, newp, newp));
                    // x = l. Pointer's position is reset to 0.0.
                    cur = 0;
                    pos = 0.0;
                } else {
                    pool.push((v1 + v2, p1 + p2, acc1));
                    if cur == pool.len() {
                        cur -= 1;
                        pos += pool[cur].0;
                    }
                }
            } else { break; }
        }
        let (v, p, acc) = pool[cur];
        let acclast = pool.last().unwrap().2;
        puts!("{:.15}\n", (acclast - acc + (p / v) * (v - pos)) / l);
    }
}

まとめ

実装の単純化が少し不足していた。

ARC086 F - Shift and Decrement

最初面倒 + 間違っている方針を考えていたが、某にもっと簡単な方法があると言われたので簡単な方法に気づくことができた。
F - Shift and Decrement

問題

N個の非負整数A_1, ..., A_Nが黒板に書かれている。以下の操作をK回まで行う。

  • 以下の2個から一つ選ぶ。
    • 操作A: 黒板の全ての整数xについて、それをfloor(x/2)で置き換える。
    • 操作B: 黒板の全ての整数xについて、それをx - 1で置き換える。これは全ての整数が1以上であるときにのみ行える。

最終的に黒板に書かれている非負整数の組み合わせとしてあり得るものをmod (10^9 + 7)で出力せよ。

解法

まず、どちらの操作も広義単調増加な変換しかしないので、操作によって整数の大きさの順番が変わることはない。よって、Aは多重集合ではなく数列として扱うことができる。Aは昇順にソートされていると仮定してよい。
どんな状態へも、[(操作B1回以下) -> (操作A)]の有限回の繰り返し -> 操作Bの有限回の繰り返しで到達することができる。最後の操作Aの前に操作Bを2回行う場合、それを直後の操作Aの後の操作B1回に置き換えられるためである。最適(操作回数最小)なものがこの形に限ることもほぼ明らか。
1を長さNのすべて1からなる数列とし、減算や乗算を項ごとに定めることにする。[(操作B1回以下) -> (操作A)]の有限回の繰り返しでAからvに到達したとすると、そこからはv, v - 1, v - 21, ..., v -k1に行くことができる。
逆に、到達可能な状態は全てこのように表せるので、上の(v, k)を全列挙すれば良い。これにはメモ化再帰が使える。
さて、このように(v, k)が列挙されたとき、それをまとめることを考えよう。vの各要素からvの最小値を引いてできる配列をd(v)としたとき、v1 - s1 = v2 - t1が成立するのはd(v1) = d(v2), v1[1] - s = v2[1] - tのとき、及びそのときに限る。(仮定からv1もv2もソートされていることに注意。)
以上から、d(v)が同じ状態をまとめると、それぞれについてあり得る最終状態はd(v) + u1 (uは然るべき区間のunionの上を動く)であることが分かる。各d(v)について、uが動く集合(下記コードのrangesのunion)の大きさを求めればよい。これは区間の出現・消失イベントの起きる時刻でソートすればできる。
計算量はどうなるだろうか? 最初のメモ化再帰で訪れる頂点の個数を見積もる。操作Aを使う個数を固定し、最後に操作Aを使った状態だけを考えると、現れる数列はA/2^b - 1からA/2^bの間に収まるので、高々N通りである。探索順を工夫すれば、各状態を2回以上探索しないようにできるため、探索すべき状態の個数はO(N * log max A)程度である。1状態ごとにかかる時間はO(N)なので、全体の計算量はO(N^2 * log max A)である。

登場する典型

  • 最適戦略についての考察
  • メモ化再帰
  • 区間のunionの大きさ

実装上の注意点

  • traverseの順番に気をつける
    • kが大きい順にtraverseすることで、ある状態に2回目以降訪れた場合にそこから先の探索が起きないことが保証できる(ダイクストラ法に近い)

提出: Submission #3974898 - AtCoder Regular Contest 086 (Rust)

fn chmax<T: Ord>(a: &mut T, b: T) -> bool {
    if *a >= b { return false; }
    *a = b;
    true
}

fn dfs(a: Vec<i64>, k: i64, memo: &mut HashMap<Vec<i64>, i64>) {
    {
        let entry = memo.entry(a.clone()).or_insert(-1);
        if !chmax(entry, k) { return }
    }
    if k == 0 { return };
    let half = a.iter().map(|&x| x / 2).collect();
    dfs(half, k - 1, memo);
    if k > 1 && a[0] > 0 {
        let half = a.iter().map(|&x| (x - 1) / 2).collect();
        dfs(half, k - 2, memo);
    }
}

const MOD: i64 = 1_000_000_007;

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: i64,
        a: [i64; n],
    }

    let mut a = a;
    a.sort();

    let mut memo = HashMap::new();
    dfs(a, k, &mut memo);
    // compress
    let mut comp = HashMap::new();
    for (a, k) in memo {
        let mi = a[0];
        let mut diff = vec![0; n - 1];
        for i in 0 .. n - 1 { diff[i] = a[i + 1] - mi; }
        let range = (max(0, mi - k), mi);
        let entry = comp.entry(diff).or_insert(Vec::new());
        entry.push(range);
    }
    let mut tot = 0;
    for (_, ranges) in comp {
        let mut events = Vec::new();
        for (x, y) in ranges {
            events.push((x, 1));
            events.push((y + 1, -1));
        }
        events.sort();
        let mut cur = 0;
        let mut last = 0;
        for (t, delta) in events {
            if cur > 0 {
                tot += t - last;
                tot %= MOD;
            }
            last = t;
            cur += delta;
        }
    }
    puts!("{}\n", tot);
}

まとめ

本番誰も解けなかったのが意外である。ちなみに私は時間ギリギリにEを通して狂喜してた。

ARC080 F - Prime Flip

最初見た時こんなの無理だろと思ったが、期間をおいて考え直したら分かった。
F - Prime Flip

問題

無限枚のカードがある。カードには1番から番号が振られている。最初x[1], ..., x[N]番目のカードは表向きで、それ以外は裏向きである。
以下の操作を何回か行うことができる。

  • 3以上の素数pを選ぶ。連続するp枚のカードをひっくり返す。

最終的に全てのカードを裏向きにする時、操作回数の最小値を求めよ。

  • 1 <= N <= 100
  • xはソートされている

解法

まず、範囲クエリなのでimos法を使い、Z/2Zの元を要素として持つ配列aに以下の操作を繰り返して全て0にする問題に変換する:

  • 3以上の素数pと正の整数uを選ぶ。a[u] += 1, a[u + p] += 1を同時に行う。

ここで、a[u]とa[v]を同時に1から0にするのに必要なコストを考えよう。v - uが3以上の素数なら1である。v - uが偶数の時、素数がたくさんあることを考えるとp1 - p2 = v - uなる素数p1, p2があるので、コストは2である。v - uが素数でない奇数の時、同様にコストは3である。
よってa[u] = 1であるような奇数と偶数をマッチングし、マッチングできたものをコスト1でマッチング、できなかったものをコスト2または3でマッチングするという実行可能解が存在することが分かる。以降最適解の中には必ずこの形のものがあることを証明する。
最適解でflipする整数全体を頂点集合とし、同時にflipする整数の間に辺を設けた無向グラフGを考える。このグラフは単純であるとしてよい。次数3以上の頂点は、上の議論を用いて次数のうち2個を別の整数にすることができるので、Gの頂点は全て次数2以下としてよい。
次数2以下の頂点からなるグラフはパスかサイクルかその直和である。サイクルは無駄なのでないとしてよい。Gの各パスについて、上の議論およびGの最適性からそれらは長さ3以下である。

登場する典型

  • Z/2Z上の範囲加算クエリをたくさん飛ばせるので、imosして2点更新クエリにする
  • 素数はたくさんあるので、2, 3個組み合わせれば全ての整数を表現できる
  • 最適解が特定の形を持つことを示し、探索すべき範囲を減らす

実装上の注意点

  • bipartite_matchingの内部で、n = 0の時にg[0]を参照しないようにする(2敗(は?))
  • Dinic法で殴る際は、辺の重みを適切に設定する(この場合全て1) (1敗)

提出: Submission #3957344 - AtCoder Regular Contest 080 (Rust)

fn prime_sieve(n: usize) -> Vec<bool> {
    let mut ans = vec![true; n];
    ans[0] = false;
    ans[1] = false;
    for i in 2 .. n {
        if !ans[i] { continue; }
        for j in 2 .. (n - 1) / i + 1 {
            ans[i * j] = false;
        }
    }
    ans
}

fn bipartite_matching(g: &[Vec<bool>]) -> usize {
    let n = g.len();
    if n == 0 { return 0; }
    let m = g[0].len();
    let mut to = vec![None; m];
    let mut visited = vec![false; n];
    let mut ans = 0;
    fn augment(v: usize, g: &[Vec<bool>],
               visited: &mut [bool], to: &mut [Option<usize>])
               -> bool {
        if visited[v] { return false; }
        visited[v] = true;
        for i in 0 .. g[v].len() {
            if !g[v][i] { continue; }
            if let Some(w) = to[i] {
                if augment(w, g, visited, to) {
                    to[i] = Some(v); return true;
                }
            } else {
                to[i] = Some(v); return true;
            }
        }
        false
    }
    for i in 0 .. n {
        for v in visited.iter_mut() { *v = false; }
        if augment(i, &g, &mut visited, &mut to) { ans += 1; }
    }
    ans
}

fn solve() {
    let out = std::io::stdout();
    let mut out = BufWriter::new(out.lock());
    macro_rules! puts {
        ($($format:tt)*) => (write!(out,$($format)*).unwrap());
    }
    input! { x: [i64] }
    const W: usize = 1 << 24;
    let is_prime = prime_sieve(W);
    let mut hs = HashSet::new();
    for x in x {
        if hs.contains(&x) {
            hs.remove(&x);
        } else {
            hs.insert(x);
        }
        let x = x + 1;
        if hs.contains(&x) {
            hs.remove(&x);
        } else {
            hs.insert(x);
        }
    }
    let mut even = Vec::new();
    let mut odd = Vec::new();
    for x in hs {
        if x % 2 == 0 { even.push(x); }
        else { odd.push(x); }
    }
    let n = even.len();
    let m = odd.len();
    let mut g = vec![vec![false; m]; n];
    for i in 0 .. n {
        for j in 0 .. m {
            if is_prime[(even[i] - odd[j]).abs() as usize] {
                g[i][j] = true;
            }
        }
    }
    let size = bipartite_matching(&g);
    let mut tot = 3 * ((n - size) % 2);
    tot += 2 * ((n - size) / 2);
    tot += 2 * ((m - size) / 2);
    tot += size;
    puts!("{}\n", tot);
}

まとめ

それっぽい方針を思いついた後の正当性の証明に手こずった。