PAKEN Programming Camp 2019 Day 1 - Speedrun O: りんごだいぼうけん

関係ないところを高速化しようとしてハマった。
りんごだいぼうけん

問題

H * W のグリッド上にリンゴが M 個、スタートとゴールがそれぞれ 1 個、壁が何個かある。リンゴを K 個以上回収してスタートからゴールまで最短距離で行け。

解法

M 個のリンゴ間、およびスタートからリンゴ、リンゴからゴールまでの距離を計算し、bitDP を行う。
距離の計算は、単純にダイクストラ法を使うと計算量が O(HWM log (HW)) ~= 4 * 10^8 になり落ちてしまうので、01BFS で線形 (O(HWM) ~= 2 * 10^7) にする。
bitDP は O(2^M M^2) ~= 4 * 10^8 かかるが、これでも通る。(おそらくキャッシュに乗りやすいため)

登場する典型

  • 01BFS
    • O(HWM log(HW)) だと落ちてしまうので
  • bitDP
    • 3s あれば N <= 20 くらいでも O(2^N N^2) で ok

実装上の注意点

頻繁に動く添字はなるべく多次元配列の最後の添字にするなどして、キャッシュに乗りやすいように書く

提出: https://onlinejudge.u-aizu.ac.jp/beta/review.html#PakenCamp19Day1/4077654 (Rust)

/*
 * Dijkstra's algorithm.
 * Verified by: AtCoder ABC035 (http://abc035.contest.atcoder.jp/submissions/676539)
 */

struct Dijkstra {
    edges: Vec<Vec<(usize, i64)>>, // adjacent list representation
}

/*
 * Code from https://doc.rust-lang.org/std/collections/binary_heap/
 */
#[derive(Copy, Clone, Eq, PartialEq)]
struct State {
    cost: i64,
    position: usize,
}

// The priority queue depends on `Ord`.
// Explicitly implement the trait so the queue becomes a min-heap
// instead of a max-heap.
impl Ord for State {
    fn cmp(&self, other: &State) -> Ordering {
        // Notice that the we flip the ordering here
        match other.cost.cmp(&self.cost) {
            std::cmp::Ordering::Equal => other.position.cmp(&self.position),
            x => x,
        }
    }
}

// `PartialOrd` needs to be implemented as well.
impl PartialOrd for State {
    fn partial_cmp(&self, other: &State) -> Option<Ordering> {
        Some(self.cmp(other))
    }
}


impl Dijkstra {
    fn new(n: usize) -> Self {
        Dijkstra { edges: vec![Vec::new(); n] }
    }
    fn add_edge(&mut self, from: usize, to: usize, cost: i64) {
        self.edges[from].push((to, cost));
    }
    /*
     * This function returns a Vec consisting of the distances from vertex source.
     */
    fn solve(&self, source: usize, inf: i64) -> Vec<i64> {
        let n = self.edges.len();
        let mut d = vec![inf; n];
        let mut que = std::collections::VecDeque::new();
        que.push_back(State {cost: 0, position: source});
        while let Some(State {cost, position: pos}) = que.pop_front() {
            if d[pos] <= cost {
                continue;
            }
            d[pos] = cost;
            for adj in &self.edges[pos] {
                if adj.1 == 0 {
                    que.push_front(State {cost: cost + adj.1, position: adj.0});
                } else {
                    que.push_back(State {cost: cost + adj.1, position: adj.0});
                }
            }
        }
        return d;
    }
}

fn main() {
    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,
        a: [chars; h],
    }
    let mut dijk = Dijkstra::new(h * w);
    let dxy = [(1, 0), (0, 1), (-1, 0), (0, -1)];
    let mut st = 0;
    let mut en = 0;
    let mut apples = vec![];
    for i in 0..h {
        for j in 0..w {
            let v = i * w + j;
            match a[i][j] {
                'a' => apples.push(v),
                's' => st = v,
                'e' => en = v,
                _ => (),
            };
            if a[i][j] == '#' { continue; }
            for d in 0..4 {
                let (dx, dy) = dxy[d];
                let nx = i as i32 + dx;
                let ny = j as i32 + dy;
                if nx < 0 || nx >= h as i32 || ny < 0 || ny >= w as i32 {
                    continue;
                }
                let nx = nx as usize;
                let ny = ny as usize;
                if a[nx][ny] == '#' { continue; }
                dijk.add_edge(v, nx * w + ny, 1);
            }
        }
    }
    let m = apples.len();
    const INF: i64 = 1 << 20;
    let mut sdist = vec![0; m];
    {
        let tmp = dijk.solve(st, INF);
        for i in 0..m {
            sdist[i] = tmp[apples[i]];
        }
    }
    let mut tdist = vec![0; m];
    {
        let tmp = dijk.solve(en, INF);
        for i in 0..m {
            tdist[i] = tmp[apples[i]];
        }
    }
    let mut mat = [[0; 20]; 20];
    for i in 0..m {
        let tmp = dijk.solve(apples[i], INF);
        for j in 0..m {
            mat[i][j] = tmp[apples[j]] as i32;
        }
    }
    let mut dp = vec![[INF as i32; 20]; 1 << m];
    for i in 0..m {
        dp[1 << i][i] = sdist[i] as i32;
    }
    for bits in 0usize..1 << m {
        for j in 0..m {
            if (bits & 1 << j) == 0 { continue; }
            for i in 0..m {
                if (bits & 1 << i) != 0 { continue; }
                let to = bits ^ 1 << i;
                dp[to][i] = min(dp[to][i], dp[bits][j] + mat[j][i]);
            }
        }
    }
    let mut mi = INF as i32;
    for bits in 0usize..1 << m {
        if bits.count_ones() < k as u32 { continue; }
        for i in 0..m {
            if (bits & 1 << i) == 0 { continue; }
            mi = min(mi, dp[bits][i] + tdist[i] as i32);
        }
    }
    puts!("{}\n", if mi >= INF as i32 { -1 } else { mi });
}

まとめ

これでダイクストラではなく bitDP を高速化しようとずっと頑張っていたのは、勘が鈍っているとしか。

PAKEN Programming Camp 2019 Day 1 - Speedrun M: カードはおやつに入りますか?(Are Cards Snacks?)

ちょっとハマったのでメモ。
カードはおやつに入りますか?(Are Cards Snacks?)

問題

N 枚のカードがあり、i 枚目には整数 A[i] が書かれている。square1001 君はカードのうち何枚かを選んで合計 K にしたいので、最小の枚数のカードを除去して邪魔せよ。

解法

愚直にやったら (自分が残すカードの集合, square1001 君がとるカードの集合) の全探索で O(3^N) になる。これでは無理なので高速化を行う。
ある集合 S に対して、「S の部分集合であって、和をとると K になるものの個数」が高速に求めたい気分になるので、dp[S] = (S の和が K なら 1, そうでなければ 0) とした配列 dp に対して高速ゼータ変換を行えばよい。

登場する典型

高速ゼータ変換

実装上の注意点

なし

提出:
https://onlinejudge.u-aizu.ac.jp/beta/review.html#PakenCamp19Day1/4077478
(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: i64,
        a: [i64; n],
    }
    let mut dp = vec![0; 1 << n];
    for bits in 0..1 << n {
        let mut sum = 0;
        for i in 0..n {
            if (bits & 1 << i) != 0 {
                sum += a[i];
            }
        }
        dp[bits] = if sum == k { 1 } else { 0 };
    }
    for i in 0..n {
        for bits in 0..1 << n {
            if (bits & 1 << i) != 0 {
                continue;
            }
            dp[bits | 1 << i] += dp[bits];
        }
    }
    let mut mi = n;
    for bits in 0usize..1 << n {
        let bc = n - bits.count_ones() as usize;
        if dp[bits] == 0 {
            mi = min(mi, bc);
        }
    }
    puts!("{}\n", mi);
}

まとめ

O(3^N) で解きたくなったら高速ゼータ変換を疑うのが良いかもとは思った。(過適合か…?)

第一回日本最強プログラマー学生選手権-予選- E - Card Collector

いい加減マトロイドの基本をスラスラと思い出せるようにならないと非常にまずいと感じた。
E - Card Collector

問題

H 行 W 列のマス目の上に N 枚のカードが置かれている。i番目のカードは R[i] 行 C[i] 列に置かれており、整数 A[i] が書かれている。(同じマス目に2枚以上のカードが置かれる可能性もある。)
以下のようにカードを取るとき、取ったカードに書かれている整数の合計値としてありうるものの最大値を求めよ。

  • 各行から最大1枚取る。
  • その後、各列から最大1枚取る。その際、直前のステップで取ったカードは取れない。

解法

「マトロイドの上の貪欲法」の一言で片付けられるが、片付けずマトロイドについての定理をこの問題での場合に特殊化して記述する。

M をカード全体の集合とする。カードの集合 A に含まれるカードを全部取れるとき、集合 A を独立集合と呼ぶことにし、集合族 F を F := {A \subseteq M | 集合 A は独立集合} で定める。
以下の補題を示す。

(補題: 増加公理) 独立集合 A, B に対して、|A| > |B| が成り立つ場合、ある x \in A \setminus B が存在して  B \cup \{x\} は独立集合。

(証明)*1
 A = A0 \coprod A1, B = B0 \coprod B1 とし、A0, B0 は行ごとに最大1個、A1, B1 は列ごとに最大1個のカードが含まれるような集合とする。 |A| > |B| であるため、|A0| > |B0| か |A1| > |B1| のどちらかは成立する。一般性を失わず、|A0| > |B0| として良い。
A0 でカバーされていて B0 でカバーされていない行が存在するので、ある x \in A0 \setminus B0 が存在して、 B0 \cup \{x\} は行ごとに最大1個のカードしか含まない。

  1. そのような x であって、B1 に含まれないものが存在する場合
    x は B に含まれず  (B0 \cup \{x\}) \coprod B1 は独立集合である。
  2. そのような x がすべて B1 に含まれる場合
    そのような要素を B1 側から B0 側に持ってくることで、 B = A0 \coprod (B1 \setminus (A0 \setminus B0)) とできる。 C := B1 \setminus (A0 \setminus B0) と置くと |A1| > |C| であるため、列に対して同じ議論を行うことで、ある y \in A1 \setminus C が存在して  B \cup \{y\} は独立集合であることがわかる。(終)

増加公理を利用して、以下の貪欲を行う。

  1. 集合 S を  S := \emptyset で初期化する。
  2. 重み A[c] の大きい順にカードを取ることを試みる。カード c を取ろうとする時、 S \cup \{c\} が独立集合であれば、S に c を追加する。そうでなければ何もしない。
  3. 全てのカードをチェックした後の S が最適解である。

これが最適解を与えることの証明を与える。
(証明)*2
仮に S と異なる集合 T が最適解であったとする。両方とも極大なので、増加公理から |S| = |T| である。S の要素を A の値が大きい順に s_0, ..., s_{k - 1}, T の要素を A の値が大きい順に t_0, ..., t_{k - 1} とする。
初めて A[s_i] < A[t_i] が成り立つ i を j とする。C := {t_0, ..., t_j}, D := {s_0, ..., s_{j - 1}} と置くと、増加公理から C のある元 t_u があって  D \cup \{t_u\} は独立。
ここで、A[t_u] >= A[t_j] > A[s_j] であるため、貪欲で s_j が選ばれたという事実と矛盾する。(終)

最後に、ある集合 S が独立集合であるかどうかの判定を行うアルゴリズムを考える。*3 左側 H 頂点、右側 W 頂点の二部グラフを考え、カード i を、左側の R[i] と右側の C[i] を結ぶ重み付きの辺とする。取ったカードたちによる部分グラフを考えた時、各連結成分に対して (辺の本数) <= (頂点の個数) であることが、独立であることと同値なので、各連結成分に対して辺の本数と頂点の個数を管理すればよい。これは UnionFind などを用いてできる。

登場する典型

実装上の注意点

とくになし。

提出: #7329715 (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, h: usize, w: usize,
        rca: [(usize1, usize1, i64); n],
    }
    let mut ans = vec![];
    for (r, c, a) in rca {
        ans.push((a, r, c));
    }
    ans.sort();
    ans.reverse();
    let mut uf = UnionFind::new(h + w);
    let mut tol = vec![1; h + w];
    let mut tot = 0;
    for &(val, r, c) in &ans {
        let v = c + h;
        if !uf.is_same_set(r, v) {
            let a = uf.root(r);
            let b = uf.root(v);
            let sum = tol[a] + tol[b];
            if sum >= 1 {
                uf.unite(a, b);
                tot += val;
                let root = uf.root(a);
                tol[root] = sum - 1;
            }
        } else {
            let a = uf.root(r);
            if tol[a] >= 1 {
                tot += val;
                tol[a] -= 1;
            }
        }
    }
    puts!("{}\n", tot);
}

まとめ

各行最大1個、各列最大1個のマトロイドは、交差が二部マッチングからなる独立性システムになるという事実がよく知られているが、それが交差ではなく合併になっても効率的に解けるという点で面白かった。

*1:行ごと最大1個、列ごと最大1個が明らかにマトロイドであり、(M, F) はこれらの合併なのでマトロイドである、という証明をインライン展開した。

*2:この証明は任意のマトロイドに適用できる。

*3:独立性オラクルとしてこのような単純なアルゴリズムが使えるのはおそらくこの問題特有で、一般にはマトロイドの合併の独立性オラクルの実装は難しい気がする。識者のツッコミ待ちです

AtCoderで赤になるまでにやったこと

AtCoderで赤タッチしたのももう1.5年も前ですが書きます。覚えていないこと・抜け・漏れなどのオンパレードですがご容赦ください。

お前は誰

kobae964 - AtCoderです。2017年11月に初めてAtCoderでredになりました。


なんで今更これ書いたの

当時は面倒だったので書きませんでした。その後人々が変色時に記事を書いていくようになり、見かけた記事は大体全部読んでいました。
微笑ましいと思ったもののAtCoderの赤でそれをやった人を見たことがなかったので、せっかくなので書こうと思いました。すでに誰か書いていたらごめんなさい。

注意事項

筆者はAGCよりもARCの方が得意なため、AGCでの立ち回りは強い人に聞いてください。

赤になるまでにやったこと

問題を解いた

2000問くらい解きました。


使った知識・メタ典型

赤タッチするなら 蟻本 + 赤の人ならほとんどが知っている知識の半分 くらいが必要になると思います。
以下では蟻本上級編以上のレベルの内容を列挙していきます。列挙の漏れは100%確実にあります。

使えた知識 (これが問われたら30分以内にわかる)
知っていたが使いこなせなかった知識 (これが問われたら1時間以上かかる)
  • 包除原理 (2^n系包除, 約数系包除, Zeta/Hadamard変換全般, ポリアの数え上げ定理)
  • 行列累乗 (遷移を行列にエンコードするのが苦手だった)
  • スライド最小値
  • 区間スケジューリング
  • Convex Hull Trick
  • 分割統治
  • 複雑なDP
    • 実家DP (配列をセグメント木で管理して、一点更新/区間和などに高速に対応する)
    • 累積和を取りながらDP (今のAtCoderでは結構見る)
    • 挿入DP
    • 桁DP
    • K種類の物を順番に埋めていく時にindexが最後に埋めた物以外の(K-1)個できるやつ (問題例)
    • …列挙していったらキリがないなこれ
  • 最大長方形 (stackを使うやつ)
  • マンハッタン距離は回す (問題例)
  • 小さい場合に手許でbrute forceして分析
  • Prim法 (は?)
  • 誤読して時間を溶かさない
  • 変な方針を思いついた場合に、変であることに素早く気付く
    • 変な貪欲に突っ走らない
知らなかったけど後で必要になった知識
  • 集合をクラスタに分けるO(3^n)系DP (O(3^n)時間、O(2^n)空間で、dp[S] <- dp[T] * sub[S - T]みたいにするやつ)
  • 構築の手筋
    • 2ベキ
    • 必要条件を列挙して可能性を絞るテクニック
  • 必要条件を集めたらいつのまにか十分条件になってました系の考察
  • 計算量周りのテクニック
  • 木についてのテク (重心分解, HL分解)
  • MSTについての性質(マトロイド性、貪欲が効くことなど、問題例)
  • MSTのアルゴリズム (Boruvka法、Voronoi図を作る方法)
使わなかった知識
  • 枝刈り全探索
  • 幾何
  • 難しいデータ構造 (AtCoder以外では使った)
    • 遅延セグメント木
    • Sparse table
    • 平衡二分木
    • 平方分割

まとめ

赤を目指すには

  • 知識を増やす (汎用性の高い方法で覚えておくとよい)
  • なんとなくの考察 (第一感) の正確性を上げる (一発で正しい答えにたどり着ける確率を上げる)
  • 嘘解法についての嗅覚を鋭くする
  • 式変形を高速に、正確に行う
  • 知らないことでも実験して考える

とかが必要になる気がします。*1 いずれまた赤から落ちると思うのでそうなったらこれの通りにまた頑張ります。

ここまで読んでくださりありがとうございました。

赤になって得られたもの (追記 2019/7/1)

強い人と知り合いになれる確率がグッと高まる

赤になることに限定せず、強くなることによるメリットとして、いわゆる「いつもの勢」と呼ばれる、国内オンサイトコンテスト常連の人たちについての話をしない訳にはいかないでしょう。
赤になった年の翌年、入社した会社で強い人と知り合い、それがきっかけとなりいつもの勢の人々と知り合いになりました。いつもの勢の輪に入れているとは思いませんが、少なくとも観測できるようにはなったと言っていいです。彼らの競プロに対する考え方は、非常に参考になります。ストイックに精進したい人にはいい環境だと思います。

OpenCupに出場できる確率が高まる

これもコネが増えたことによる帰結です。
OpenCupについて、chokudaiさんなどが言及しているのをたまに目にする人がいるかもしれません。OpenCupというのは、ある期間日曜午後に (不定期に) 行われる、5時間チーム戦内輪コンテストです。*2 これも社にいた強い人に教えてもらい、幸運にも組んでいただける人を見つけることができました。精進したい人にはオススメです。

ARCがunratedになる

ARCみたいな面白いコンテストにunratedで出て新戦法を試すことができるのはメリットだと思います。多分tourist出しの練習をした人絶対いると思いますよ

Q. ユーザによっては一度に限り変更できる場合があるようですが、なぜ私は出来ないのでしょうか?

先日ユーザ名を一度だけ変更できる機能が実装されましたが、最大レートが2800以上のユーザはユーザ名の変更ができません。草

f:id:koba-e964:20190701021801p:plain
max rating >= 2800になると現れる文言

モチベーション維持の方法 (追記 2019/7/1)

ブランク期間の長さにペナルティをつける

自分の心の中のレート := 本来のレート - 100 / 1ヶ月 * ブランク期間 くらいにすると、ratedを避けるモチベがなくなっていいっすよ
これはAtCoder固有の事情ですが、ある程度の参加回数になると1回の大失敗で下げられるレートの上限が121 (= -800 * log_2(0.9)) になります。つまり36日ratedに出ないともう何やっても悪くなりません。無敵だな!

ポエム (ここから下は読む必要なし) (追記 2019/7/1)

私は高校時代、数学オリンピックに出ていましたが、ろくに何も精進しなくて本戦 (予選の次にある2番目のやつ) 止まりで、春合宿出場者やIMO出場者に対して羨望に近い憧れを抱いていました。そんな私が、高校時代から知ってはいたものの手を出せずにいた競技プログラミングに手を出したのが大学1年の時です。当時はTopCoder *3 がまだ盛んで、参加して緑落ちしては嘆き、黄色に昇格しては狂喜乱舞し、というのをやっていました。おそらく当時はUnionFindやbitDPすらまともに書けなかったと思います。レッドコーダーたちについては、尊敬の念をもって見ており、自分がなれるとは夢にも思っていませんでした。
転機と呼べるものがいつあったか、正確には思い出せませんが、大学3年くらいの時だと思います。演習で競プロみたいな内容のものがあり、競プロの面白さに目覚めました。そのときはyukicoderあたりを埋めまくって典型を覚えていったと思います。
本格的に取り組み始めたのは院生になってからだと思います。研究そっちのけで (は?) 競プロに打ち込み、M1の11月にオレンジになりました。そこからはより一層のめり込み、友達との食事の時間も考察に明け暮れる始末で (は?)、ratedのたびに2800以上のパフォーマンスを取っては赤パフォだと喜び、2400以下のパフォーマンスを取っては落ち込み、という繰り返しをしました。
途中からARCでパフォーマンス3200を取れるようになってきたので、このまま続けていればいずれは赤になれる、と徐々に確信できるようになりました。が、実際になれたときの喜びはひとしおでした。高校時代憧れの的でしかなかった国際科学オリンピック出場者に一歩だけ近づけたのだという思いが湧きました。

この先 (追記 2019/7/1)

誰か強くなる方法を教えてください m(_ _)m

*1:誰も「赤になる方法」みたいな記事を書きたがらないように見えるのは、こういう情報量ゼロの内容しかないからでは、と思っています。

*2:内輪なので参加するのはちょっと面倒です

*3:今は表記がTopcoderになっていますが、2012年頃はTopCoderだったと思います

階乗 mod pの高速な実装 (Rust)

階乗 mod 素数 - memo のO(p^{1/2} log p)解法の、Rustでの実装を与えた。説明は元記事の方で十分にされているので、この記事での説明は期待しないでほしい。また、実装にはところどころ雑な部分があるので、今後refineするかもしれない。

問題

n! mod pを計算せよ。
No.502 階乗を計算するだけ - yukicoder

コードの仕様

MOD: i64: mod。今回は10^9 + 7

登場する典型

元記事にない補足事項

元記事における整数vは、この実装では2^15に決め打ちしている。これは、v * v >= p = 10^9 + 7を満たす最小の2ベキである。

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

// mod_int::ModInt 省略

macro_rules! define_mod {
    ($struct_name: ident, $modulo: expr) => {
        #[derive(Copy, Clone, PartialEq, Eq, PartialOrd, Ord, Hash)]
        struct $struct_name {}
        impl mod_int::Mod for $struct_name { fn m() -> i64 { $modulo } }
    }
}
const MOD: i64 = 1_000_000_007;
define_mod!(P, MOD);
type ModInt = mod_int::ModInt<P>;

/// FFT (in-place, verified as NTT only)
/// R: Ring + Copy
/// Verified by: https://codeforces.com/contest/1096/submission/47672373
mod fft {
    use std::ops::*;
    /// n should be a power of 2. zeta is a primitive n-th root of unity.
    /// one is unity
    /// Note that the result should be multiplied by 1/sqrt(n).
    pub fn transform<R>(f: &mut [R], zeta: R, one: R)
        where R: Copy +
        Add<Output = R> +
        Sub<Output = R> +
        Mul<Output = R> {
        let n = f.len();
        assert!(n.is_power_of_two());
        {
            let mut i = 0;
            for j in 1 .. n - 1 {
                let mut k = n >> 1;
                loop {
                    i ^= k;
                    if k <= i { break; }
                    k >>= 1;
                }
                if j < i { f.swap(i, j); }
            }
        }
        let mut zetapow = Vec::new();
        {
            let mut m = 1;
            let mut cur = zeta;
            while m < n {
                zetapow.push(cur);
                cur = cur * cur;
                m *= 2;
            }
        }
        let mut m = 1;
        while m < n {
            let base = zetapow.pop().unwrap();
            let mut r = 0;
            while r < n {
                let mut w = one;
                for s in r .. r + m {
                    let u = f[s];
                    let d = f[s + m] * w;
                    f[s] = u + d;
                    f[s + m] = u - d;
                    w = w * base;
                }
                r += 2 * m;
            }
            m *= 2;
        }
    }
}

mod arbitrary_mod {
    use mod_int;
    use fft;
    const MOD1: i64 = 1012924417;
    const MOD2: i64 = 1224736769;
    const MOD3: i64 = 1007681537;
    const G1: i64 = 5;
    const G2: i64 = 3;
    const G3: i64 = 3;
    define_mod!(P1, MOD1);
    define_mod!(P2, MOD2);
    define_mod!(P3, MOD3);

    fn zmod(mut a: i64, b: i64) -> i64 {
        a %= b;
        if a < 0 {
            a += b;
        }
        a
    }
    fn ext_gcd(mut a: i64, mut b: i64) -> (i64, i64, i64) {
        let mut x = 0;
        let mut y = 1;
        let mut u = 1;
        let mut v = 0;
        while a != 0 {
            let q = b / a;
            x -= q * u;
            std::mem::swap(&mut x, &mut u);
            y -= q * v;
            std::mem::swap(&mut y, &mut v);
            b -= q * a;
            std::mem::swap(&mut b, &mut a);
        }
        (b, x, y)
    }
    fn invmod(a: i64, b: i64) -> i64 {
        let x = ext_gcd(a, b).1;
        zmod(x, b)
    }

    // This function is ported from http://math314.hateblo.jp/entry/2015/05/07/014908
    fn garner(mut mr: Vec<(i64, i64)>, mo: i64) -> i64 {
        mr.push((mo, 0));

        let mut coffs = vec![1; mr.len()];
        let mut constants = vec![0; mr.len()];
        for i in 0..mr.len() - 1 {
            let v = zmod(mr[i].1 - constants[i], mr[i].0) * invmod(coffs[i], mr[i].0) % mr[i].0;
            assert!(v >= 0);
            for j in i + 1..mr.len() {
                constants[j] += coffs[j] * v % mr[j].0;
                constants[j] %= mr[j].0;
                coffs[j] = coffs[j] * mr[i].0 % mr[j].0;
            }
        }
        constants[mr.len() - 1]
    }

    // f *= g, g is destroyed
    fn convolution_friendly<P: mod_int::Mod>(a: &[i64], b: &[i64], gen: i64) -> Vec<i64> {
        use mod_int::ModInt;
        let d = a.len();
        let mut f = vec![ModInt::<P>::new(0); d];
        let mut g = vec![ModInt::<P>::new(0); d];
        for i in 0..d {
            f[i] = a[i].into();
            g[i] = b[i].into();
        }
        let zeta = ModInt::new(gen).pow((P::m() - 1) / d as i64);
        fft::transform(&mut f, zeta, ModInt::new(1));
        fft::transform(&mut g, zeta, ModInt::new(1));
        for i in 0..d {
            f[i] *= g[i];
        }
        fft::transform(&mut f, zeta.inv(), ModInt::new(1));
        let inv = ModInt::new(d as i64).inv();
        let mut ans = vec![0; d];
        for i in 0..d {
            ans[i] = (f[i] * inv).x;
        }
        ans
    }


    pub fn arbmod_convolution(a: &mut [i64], b: &mut [i64], mo: i64)
                          -> Vec<i64> {
        use ::mod_int::Mod;
        let d = a.len();
        assert!(d.is_power_of_two());
        assert_eq!(d, b.len());
        for x in a.iter_mut() {
            *x = zmod(*x, mo);
        }
        for x in b.iter_mut() {
            *x = zmod(*x, mo);
        }
        let x = convolution_friendly::<P1>(&a, &b, G1);
        let y = convolution_friendly::<P2>(&a, &b, G2);
        let z = convolution_friendly::<P3>(&a, &b, G3);

        let mut ret = vec![0; d];
        let mut mr = [(0, 0); 3];
        for i in 0..d {
            mr[0] = (P1::m(), x[i]);
            mr[1] = (P2::m(), y[i]);
            mr[2] = (P3::m(), z[i]);
            ret[i] = garner(mr.to_vec(), mo);
        }
        ret
    }
}

// f *= g, g is not destroyed
fn convolution(f: &mut [i64], g: &mut [i64]) {
    let ans = arbitrary_mod::arbmod_convolution(f, g, MOD);
    for i in 0..f.len() {
        f[i] = ans[i];
    }
}

fn grow(d: i64, v: i64, mut h: Vec<i64>,
        invfac: &[ModInt]) -> Vec<i64> {
    assert_eq!(h.len() as i64, d + 1);
    let dd = d as usize;
    let dm = ModInt::new(d);
    let vm = ModInt::new(v);

    let mut aux = vec![1; dd];

    let mut f = vec![0; 4 * dd];
    let mut g = vec![0; 4 * dd];
    for i in 0..dd + 1 {
        f[i] = (invfac[i] * invfac[dd - i] * h[i]).x;
        if (dd + i) % 2 != 0 {
            f[i] = if f[i] == 0 { 0 } else { MOD - f[i] };
        }
    }
    let oldf = f.clone();
    for (idx, &a) in [dm + 1, dm * vm.inv(), dm * vm.inv() + dm + 1].iter().enumerate() {
        for i in 0..4 * dd { f[i] = oldf[i]; }
        for i in 0..4 * dd { g[i] = 0; }
        for i in 1..2 * dd + 2 {
            g[i] = (a - d + i as i64 - 1).inv().x;
        }
        convolution(&mut f, &mut g);
        let mut prod = 1;
        for i in 0..dd + 1 {
            prod = prod * (a - i as i64).x % MOD;
            assert_ne!(prod, 0);
        }
        for i in 0..dd + 1 {
            f[dd + i + 1] = f[dd + i + 1] * prod % MOD;
            prod = prod * (a + i as i64 + 1).x % MOD;
            prod = prod * (a - d + i as i64).inv().x % MOD;
        }
        match idx {
            1 => {
                for i in 0..dd + 1 {
                    h[i] = h[i] * f[dd + 1 + i] % MOD;
                }
            }
            0 => {
                for i in 0..dd {
                    aux[i] = f[dd + 1 + i];
                }
            }
            2 => {
                for i in 0..dd {
                    aux[i] = aux[i] * f[dd + 1 + i] % MOD;
                }
            }
            _ => unreachable!(),
        }
    }
    h.extend_from_slice(&aux);
    h
}

fn gen_seq(d: i64, v: i64) -> Vec<i64> {
    assert!(d > 0 && (d as u64).is_power_of_two());
    let dd = d as usize;

    // precompute factorial and its inv
    let mut fac = vec![ModInt::new(0); 2 * dd + 1];
    let mut invfac = vec![ModInt::new(0); 2 * dd + 1];
    fac[0] = ModInt::new(1);
    for i in 1..2 * dd + 1 {
        fac[i] = fac[i - 1] * (i as i64);
    }
    invfac[2 * dd] = fac[2 * dd].inv();
    for i in (0..2 * dd).rev() {
        invfac[i] = invfac[i + 1] * (i as i64 + 1);
    }
    let mut size = 1;
    // Initialized with [g_1(0), g_1(v)].
    let mut seq = vec![1.into(), (v + 1).into()];
    while size < d {
        seq = grow(size, v, seq, &invfac);
        size *= 2;
    }
    assert_eq!(size, d);
    seq
}

fn fact(n: i64) -> ModInt {
    let d = 1 << 15;
    let aux = gen_seq(d, d);
    // eprintln!("{:?}", aux);
    let mut ans = ModInt::new(1);
    let lim = min(d, (n + 1) / d);
    for i in 0..lim {
        ans *= aux[i as usize];
    }
    for i in lim * d..n {
        ans *= i + 1;
    }
    ans
}

// Uses techniques described in https://min-25.hatenablog.com/entry/2017/04/10/215046.
// Bostan, A., Gaudry, P., & Schost, É. (2007). Linear Recurrences with Polynomial Coefficients and Application to Integer Factorization and Cartier–Manin Operator. SIAM Journal on Computing, 36(6), 1777–1806. https://doi.org/10.1137/s0097539704443793
fn main() {
    input!(n: i64);
    if n >= MOD {
        println!("0");
    } else {
        println!("{}", fact(n));
    }
}

まとめ

  • 400行も書いたので流石に疲れた
  • だれか可変modのジャッジのリンクをください

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 \subseteq \{0, \ldots, N - 2\}に対して、以下の値がわかればよい:

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

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

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

DP1[{}] = 1は自明である(空を並べる方法は1通りであり、これはトポロジカルソート)。Sに新しい辺eを追加するときを考える。Sの辺に依存せずeに依存するような辺はints[  S \cup \{e\} ] - 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[  S \cup \{e\}] += (DP2[S] * ints[  S \cup \{e\} ] / (ints[S] + 1) + DP1[S] * (|S| + 1)) / ints[  S \cup \{e\} ]

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

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

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

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

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

登場する典型

  • ゼータ変換
    • intsの計算に必要
    • ある部分集合 S \subseteq \{0, \ldots, k - 1\}について、「 T \supseteq Sなるすべての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」で済まされていてキレかけた(は?)。多分諦めて解説を見ていたらキレていたと思う。

Codeforces Global Round 2 H. Triple

本番Eまで解いてF, G, Hのどれを解くか、という状態になったときに一番可能性の高かったこの問題に粘着した。結局時間内には解けなかったが、悔しかったので全力で頑張ったら解けた。
https://codeforces.com/contest/1119/problem/H

問題

整数kが与えられる。数列がn個ある。どの数列も長さ(x + y + z)で、要素は0から2^k-1までの整数である。i番目の数列s[i]はa[i]をx個、b[i]をy個、c[i]をz個含む。
t = 0, 1, ..., 2^k - 1に対して、以下の個数を数え上げて998244353で割った余りを求めよ:

  • 長さnの整数列(u[0], ..., u[n - 1])であって、s[0][u[0]] xor s[1][u[1]] xor ... xor s[n - 1][u[n - 1]] = tを満たすもの。

制約

  • 1 <= n <= 10^5
  • 1 <= k <= 17
  • 0 <= a[i], b[i], c[i] <= 2^k - 1
  • 0 <= x, y, z <= 10^9

解法

各列s[i]について長さ2^kの頻度表を作り、それをv[i]と呼ぶ。v[i]のアダマール変換をh[i]とすれば、求めるものはh[i]の要素ごとの積の逆アダマール変換である。
これを愚直にやるとO(n * k * 2^k)時間かかってしまうので高速化する。一般性を失わず、a[i] = 0としてよい。(b[i]とc[i]にa[i]をxorしておき、最終的な結果にa[i]たちのxorをxorする。) こうすると、アダマール変換後の数列h[i]に登場する数はx ± y ± zという形の4通りのみとなる。それぞれがいつ起こるかは以下の通り:

  • h[i][j] = x - y - z: popcount(b[i] & j)が奇数、popcount(c[i]&j)が奇数
  • h[i][j] = x + y - z: popcount(b[i] & j)が偶数、popcount(c[i]&j)が奇数
  • h[i][j] = x - y + z: popcount(b[i] & j)が奇数、popcount(c[i]&j)が偶数
  • h[i][j] = x + y + z: popcount(b[i] & j)が偶数、popcount(c[i]&j)が偶数

各jに対して、これらが現れる回数を数え上げればよい。
実はこれもアダマール変換を用いることでできる。例えば長さ2^kの配列pを、p[j] = (popcount(b[i] & j)もpopcount(c[i]&j)も奇数 ? 1 : 0)で定めると、配列pは以下によってつくられる長さ2^kの配列qのアダマール変換として得られる:

  • qを0で初期化, q[0] += 1/4, q[b[i] ^ c[i]] += 1/4, q[b[i]] -= 1/4, q[c[i]] -= 1/4

これを利用すれば、毎回配列のたかだか16箇所に足し算をし、最後にまとめてアダマール変換を行うことで、x ± y ± zに乗せるべき指数がわかる。

計算量はO(2^k * k + n)である。

登場する典型

実装上の注意点

  • 手数が多い上に、アダマール変換を使う都合上所々で係数の調整が必要なので、注意して実装する。
  • (x-y-z) % MODみたいな演算をバグらせない
    • MOD < 10^9なので2 * MODを足すだけでは不十分、3*MODを足す (1敗)

提出: #52511941 (Rust)

// ModInt省略

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,
        x: i64, y: i64, z: i64,
        abc: [(usize, usize, usize); n],
    }
    let inv2 = ModInt::new(2).inv();
    let mut base = 0;
    // bias: 4 *
    let mut tbl = vec![vec![0i64; 1 << k]; 4];
    for &(a, b, c) in &abc {
        base ^= a;
        let b = a ^ b;
        let c = a ^ c;
        tbl[0][0] += 4;
        tbl[1][0] += 2;
        tbl[1][b] -= 2;
        tbl[2][0] += 2;
        tbl[2][c] -= 2;
        tbl[3][0] += 1;
        tbl[3][b ^ c] += 1;
        tbl[3][b] -= 1;
        tbl[3][c] -= 1;
    }
    for bits in 0..1 << k {
        for i in 0..2 {
            for j in 0..4 {
                if (j & 1 << i) != 0 {
                    tbl[j ^ 1 << i][bits] -= tbl[j][bits];
                }
            }
        }
    }
    // Hadamard
    for c in 0..4 {
        for i in 0..k {
            for bits in 0..1 << k {
                if (bits & 1 << i) == 0 {
                    let x = tbl[c][bits];
                    let y = tbl[c][bits | 1 << i];
                    tbl[c][bits] = x + y;
                    tbl[c][bits | 1 << i] = x - y;
                }
            }
        }
        for bits in 0..1 << k {
            tbl[c][bits] >>= 2;
        }
    }
    let mut g = vec![ModInt::new(0); 4];
    for kind in 0..4 {
        let mut t = x;
        t += if (kind & 1) == 0 { y } else { -y };
        t += if (kind & 2) == 0 { z } else { -z };
        g[kind] = ModInt::new(t + 3 * MOD);
    }
    let mut prod = vec![ModInt::new(1); 1 << k];
    for bits in 0..1 << k {
        for c in 0..4 {
            prod[bits] *= g[c].pow(tbl[c][bits]);
        }
    }
    // Hadamard
    for i in 0..k {
        for bits in 0..1 << k {
            if (bits & 1 << i) == 0 {
                let x = prod[bits];
                let y = prod[bits | 1 << i];
                prod[bits] = x + y;
                prod[bits | 1 << i] = x - y;
            }
        }
    }
    let fac = inv2.pow(k as i64);
    for bits in 0..1 << k {
        prod[bits] *= fac;
    }
    for bits in 0..1 << k {
        puts!("{}{}", prod[bits ^ base], if bits + 1 == (1 << k) { "\n" } else { " " });
    }
}

まとめ

p[j] = (popcount(b[i] & j)もpopcount(c[i]&j)も奇数 ? 1 : 0)の(逆)アダマール変換が4箇所以下がhotな配列という性質が見えにくかった。
解けたときはドーパミンが多量に分泌されて気持ちよくなった。しかし時間内に解けないのはやはり問題なので、次は時間内に解きたい。