ABC116 D - Various Sushi

こういう、適切な貪欲をやる問題が苦手すぎる…。
D - Various Sushi

問題

N個の寿司があり、i個目の寿司の種類はt_i、美味しさはd_iである。この中からK個選び、美味しさd_iの総和とt_iの種類数の2乗の和を最大化したい。最大値を求めよ。

解法

寿司を美味しさの順に並べたとき、各種類について一番美味しさが高いものをinitiator、それ以外のものをfollowerと呼ぶことにする。
最適な選び方を一つ固定する。これはどういう形をしているだろうか? まず、各種類について、その種類の寿司を取ったのであれば必ずinitiatorは取っている。また、選んだfollowerよりも価値の高い選んでいないinitiatorがある場合、followerをinitiatorと取り替えることによりスコアが高くなるので、そのようなことはないとして良い。
よって、最適戦略は上からinitiatorを何個か取って、そのあとにfollowerを何個か取った形をしている。上から取るinitiatorの個数を全探索すればよい。
initatorを新しく取るたびに一番美味しさの低いfollowerを捨てる操作が高速にできればよい。これは優先度キュー、vectorなどでできる。

登場する典型

  • 貪欲
  • 最適解の形を考える

実装上の注意点

  • 新しいinitiatorを追加するたびに、2 * (今まで取ったinitiatorの個数) + 1だけスコアに加算するのを忘れない!

提出: Submission #4059339 - AtCoder Beginner Contest 116 (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,
        td: [(usize1, i64); n],
    }
    let mut dt = vec![];
    for (t, d) in td {
        dt.push((d, t));
    }
    dt.sort(); dt.reverse();
    let mut que = BinaryHeap::new();
    let mut seen = vec![false; n];
    let mut tot = 0;
    let mut kind = 0;
    for i in 0 .. k {
        let (d, t) = dt[i];
        tot += d;
        if seen[t] {
            que.push((-d, t));
        } else {
            seen[t] = true;
            tot += 2 * kind + 1;
            kind += 1;
        }
    }
    let mut ma = tot;
    for i in k .. n {
        let (d, t) = dt[i];
        // No point of picking it. Skip.
        if seen[t] { continue; }
        let profit = d + 2 * kind + 1;
        if let Some((tap, _ris)) = que.pop() {
            let tap = -tap;
            tot += profit - tap;
            seen[t] = true;
            kind += 1;
        }
        ma = max(ma, tot);
    }
    puts!("{}\n", ma);
}

まとめ

最近のABCは結構定跡の把握漏れに気付かせてくれる気がする。

DISCO presents ディスカバリーチャンネル コードコンテスト2019 本戦 参加記

楽しいコンテストだった

来社

寝坊して9:24頃到着した(受付は9:00-9:50)。運良くコンセントのある席が残っていてよかった。

本戦問題コード部門

DISCO presents ディスカバリーチャンネル コードコンテスト2019 本戦 - AtCoder
A (00m00s-07m56s): 「300にしては難しくないか?」と思った。氷以外が必ず2マス以上連続で現れるのはわかっていたが、それをうまく実装の簡略化に結びつけられなかった。
B (07m56s-18m27s): idsigmaさんのツイートを以前見ていたおかげでその時の考察を再利用できた。結局大きい方から順にpushしていって転倒数を消費することを考えればよい。大きい方をpushするとオーバーしてしまう場合は小さい方をpushする。


D (18m27s-33m45s): 実家。CとDを見てどちらに行くか迷ったが、yukicoderでhamkoしゃんが似た問題(No.510 二次漸化式 - yukicoder)を出していたので解法を秒で思いついたため、最初にDを解くことにした。これが通せて暫定3位になった時点で作戦勝ちだと思った。
C (33m45s-105m23s): 幾何。当然円の周囲に触れるところで反射させるのが一番効率的だし、ある地点まで光を飛ばせるならそこより近いところにも飛ばせるので、反射で届く(柱にブロックされない)一番遠いところを計算するだけ。ただ、幾何ライブラリを整備しておらずその場で書いたため、単純なバグに長時間苦しんだ。これは負けですね…
E (105m23s-120m00s): 構築。15分しかないのでダメみたいですね…。終了です。

本戦終了時点で凍結含め最悪18位だったので、最終的に10位になれてよかった。Cに時間をかけたせいでEみたいな面白いけど時間がかかる問題に手を出せなかったのは悲しい。

本戦問題装置実装部門

無理でした おわり

コンテスト後

いつもの勢の方々とご飯に行った。kyuriさんと一緒にyosupoさんにEの解を教えてもらったりした

まとめ

反省して幾何ライブラリを整備します (少なくとも、コピペできるようなオンラインの幾何ライブラリを見繕っておきます)

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);
    }
}

まとめ

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