【AtCoder】パ研合宿2023 第1日「SpeedRun」D - Bishop

問題

解法

 (x, y) に操作を行うと  (x + K, y + K) または  (x + K, y - K) に移る。

つまり、 (x + y, x - y) に操作を行うと  (x + y + 2 K, x - y) または  (x + y, x - y - 2 K) に移る。

ここで  a = x + y, b = x - y という変数を導入すると

 (a, b) は操作後に  (a + 2 K, b) または  (a, b - 2K) に移るといえる。

これで  a b を独立に考えることができるようになったため、答えは  \left \lceil \frac{|a|}{2K} \right \rceil + \left \lceil \frac{|b|}{2K} \right \rceil になる。

#include "my_template.hpp"
using namespace std;

void solve() {
    I64(K, X1, Y1, X2, Y2);
    i64 x = abs(X1 - X2), y = abs(Y1 - Y2);
    i64 a = abs(x + y), b = abs(x - y);
    i64 ans = ceil(a, 2 * K) + ceil(b, 2 * K);
    print(ans);
    return;
}

int main() {
    solve();
    return 0;
}

感想

無限人通していてかなり厳しい気持ちになった。

45度回転をマンハッタン距離の最大化の文脈でしか使えていないということが判明。 もうちょっと一般化して、 x+y, x-y という変数を導入することで式変形ができるようになる手法、くらいに思っておくと良さそう。

マンハッタン距離の最大化

 (x_1, y_1), (x_2, y_2), ... , (x_n, y_n) が与えられる。 点  (x_q, y_q ) とのマンハッタン距離の最大値について

 (a_i, b_i) = (x_i + y_i, x_i - y_i) とおくと

 \max_{1 \leq i \leq n} ( |x_i - x_q| + |y_i - y_q| )

 = \max_{1 \leq i \leq n} ( \max(x_i - x_q, x_q - x_i ) + \max(y_i - y_q, y_q - y_i) )

 = \max_{1 \leq i \leq n} ( x_i - x_q + y_i - y_q, x_i - x_q + y_q - y_i, x_q - x_i + y_i - y_q, x_q - x_i + y_q - y_i)

 = \max_{1 \leq i \leq n} ( a_i - a_q, b_i - b_q, b_q - b_i, a_q - a_i)

 = \max_{1 \leq i \leq n} ( \max( |a_i - a_q|, |b_i - b_q |))

となり、 a, b それぞれで最大値を求めてから大きい方を選択すれば求められるというテクニックのこと。

またこの  \max( |a_i - a_q|, |b_i - b_q |)  (a_i, b_i)  (a_q, b_q) のチェビシェフ距離と言う。

【AtCoder】ARC173 C - Not Median

問題

解法

以下 0-indexed とする。

 i = 0, N - 1 の場合、区間を一方向に伸ばしながら Segment Tree 上の二分探索を行うなどして求める。

それ以外のとき、各  P_{i} に対し、自身より大きい値か小さい値かしか興味がないので便宜上 +- とおく。

(1) +0+-0- のとき

この区間の中央値は  P_{i} ではないため答えとなり、かつ 3 が最小なのでこれで良い。

(2) -0++0- のとき

+0- を例に取ると -+-+...-+0-+...-+ となる区間は 0 を含みつつ奇数長の区間を選ぶ場合、どこからどこまでを選んでも中央値が 0 になってしまうのでダメ。

逆に初めて --++ が連続するところから 0 まで(場合によっては奇数長になるように反対側の区間を 1 だけ余計に)取ると、その区間は中央値が 0 にならない。 反対側の区間を必要以上に長く取る必要はないので、これを左右それぞれについて計算し、区間長がより小さいものが答えになる。

よって左右それぞれについて愚直に ++-- が出てくる部分まで伸ばしていくことで答えを求めることができる。

これは一見計算量が  O(N ^2) のように見えるが、この愚直がどこからどこまで行われるかについて考えると実は  O(N) である。

 i に対し  P _{i} をプロットした図を考えると (2) となるのは *↘*↘**↗*↗* であり、これらの位置を  i _{0}, i _{1}, ... とおくと 愚直に伸ばしていく操作は隣接する  i で必ず終了するということがわかる。

隣接する  i*↘*↘* の場合を例に取ると 0****0****0****0 であるが、どこでも必ず *↘*↘* が登場していることがわかる。(図がないため説明が難しいが)

よって以上を実装すれば AC できる。(P[i - 1] - P[i]) * (P[i + 1] - P[i]) > 0 などのオーバーフローに気をつける。(一敗)

#include "my_template.hpp"
#include "data_structure/segment_tree.hpp"
#include "algebra/monoid_sum.hpp"
using namespace std;

signed main() {
    INT(N);
    VEC(i64, P, N);
    REP(i, N) P[i]--;
    vector<int> ans(N, -1);
    // calc ans[0], ans[N - 1]
    REP(t, 2) {
        SegmentTree<MonoidSum<int>> seg(N);
        int len = 0;
        REP(i, N) {
            seg.set(P[i], 1);
            len++;
            auto f = [&](int s) -> bool { return s <= len / 2; };
            if (len >= 3 and len % 2 == 1) {
                if (seg.max_right(0, f) != P[0]) {
                    ans[0] = len;
                    break;
                }
            }
        }
        // reverse
        REV(P);
        REV(ans);
    }
    // calc ans[1], ans[2], ... , ans[N - 2]
    REP(i, 1, N - 1) {
        if ((P[i - 1] - P[i]) * (P[i + 1] - P[i]) > 0) {
            // +0+ or -0-
            ans[i] = 3;
        } else {
            // +0- or -0+
            int l = i - 2, r = i + 2;
            while (l >= 0 and (P[l] - P[i]) * (P[l + 1] - P[i]) < 0) l--;
            while (r < N and (P[r] - P[i]) * (P[r - 1] - P[i]) < 0) r++;
            int len = INF<int>;
            if (l >= 0) chmin(len, i - l + 1);
            if (r < N) chmin(len, r - i + 1);
            if (len == INF<int>) continue;
            ans[i] = len / 2 * 2 + 1;
        }
    }
    print(ans);
    return 0;
}

感想

愚直操作が  O(N)目から鱗だった。こういう計算量の見積もり方をしたことがなく面白かった。(解けなかったけど)

ARC173 B(点集合が与えられるので非退化な三角形を可能な限りたくさん作るやつ)もたまたますんなりと解けたが証明はできなかったので、ARC は証明能力が重要だなと思った。

ARC の黄 diff 解けね〜な〜

【Codeforces】Codeforces Round 927 (Div. 3) G. Moving Platforms

問題

解法

現在時刻を  t とすると、ダイクストラ法で各遷移において  L _{i} + S _{i}  x ≡ L _{j} + S _{j}  x \pmod H となる最小の  x > t を求めることができれば良いので線形合同式が解ければ良い。

しかしこれを解くのは意外と難しく、法が素数ではないときに逆元がないが  x が解を持つという場合がある。(例: 9x \equiv 3 \pmod {15} に対し  x \equiv 2 \pmod 5

なので tourist 氏や maspy 氏のコードを参考に ライブラリ を作成した。

結局のところ  ax + by = c が解ければ良いのだが、解  (x, y) が求められたときに  x の周期が  b ではなく b / gcd(a, b) であるという点に注意する。

#include "my_template.hpp"
#include "graph/read_graph.hpp"
#include "math/linear_diophantine.hpp"

using namespace std;

void solve() {
    I64(N, M, H);
    VEC(i64, L, N);
    VEC(i64, S, N);
    auto g = read_graph<int>(N, M);

    vector<i64> dp(N, INF<i64>);
    pqueg<pair<i64, int>> que;
    dp[0] = -1;
    que.emplace(-1, 0);

    while (!que.empty()) {
        auto [t, i] = que.top();
        que.pop();
        if (t != dp[i]) continue;
        FORE(e, g[i]) {
            int j = e.to;
            // L[i] + S[i] * x ≡ L[j] + S[j] * x (mod H) となる x を求める
            // (S[i] - S[j]) * x ≡ L[j] - L[i] (mod H)
            auto [x, mod] = linear_congruence(S[i] - S[j], L[j] - L[i], H);
            if (mod == -1) continue;
            i64 nt = (t + 1);
            if (nt % mod <= x) {
                nt = nt / mod * mod + x;
            } else {
                nt = ceil(nt, mod) * mod + x;
            }
            if (chmin(dp[e.to], nt)) {
                que.emplace(dp[e.to], e.to);
            }
        }
    }
    i64 ans = dp[N - 1] + 1;
    if (ans >= INF<i64> / 2) ans = -1;
    print(ans);
    return;
}

int main() {
    INT(T);
    REP(T) solve();
    return 0;
}

感想

線形合同式を理解し直すことができて良かった。

【AtCoder】ARC169 C - Not So Consecutive

問題

解法?(本題とは関係なし)

まず、dp[i][j][k] = i 個目まで見て末尾が j で k 個続いているものの個数 として各要素の遷移に  O(N) かけて全体  O(N ^ 4) のコードができる。 (下記コードのコメントアウト部分に対応)

簡単に思いつく高速化としてほとんどの adp[i + 1][a][1] に遷移するのでまとめて計算する。

#include "my_template.hpp"
#include "math/static_modint.hpp"
using namespace std;

using mint = mint998;

int main() {
    INT(N);
    VEC(int, A, N);
    REP(i, N) if (A[i] != -1) A[i]--;
    vector dp(N, vector<mint>(N + 2));
    if (A[0] == -1) {
        REP(a, N) dp[a][1] = 1;
    } else {
        dp[A[0]][1] = 1;
    }
    REP(i, 1, N) {
        vector np(N, vector<mint>(N + 2));
        if (A[i] == -1) {
            mint s = 0;
            REP(j, N) {
                REP(k, j + 2) {
                    s += dp[j][k];
                    np[j][1] -= dp[j][k];
                    if (k <= j) np[j][k + 1] += dp[j][k];
                    /*
                    REP(a, N) {
                        if (j == a) {
                            if (k + 1 <= j + 1) np[j][k + 1] += dp[j][k];
                        } else {
                            np[a][1] += dp[j][k];
                        }
                    }
                    */
                }
            }
            REP(a, N) np[a][1] += s;
        } else {
            REP(j, N) {
                REP(k, j + 2) {
                    if (j == A[i]) {
                        if (k <= j) np[j][k + 1] += dp[j][k];
                    } else {
                        np[A[i]][1] += dp[j][k];
                    }
                    /*
                    REP(a, A[i], A[i] + 1) {
                        if (j == a) {
                            if (k + 1 <= j + 1) np[j][k + 1] += dp[j][k];
                        } else {
                            np[a][1] += dp[j][k];
                        }
                    }
                    */
                }
            }
        }
        swap(dp, np);
    }
    mint ans = 0;
    REP(j, N) REP(k, j + 2) ans += dp[j][k];
    print(ans);
    return 0;
}

ここまではバーチャルコンテスト中に考えた解法で、結論としてはこれを  O(N ^ 2) に落とすことはできなかった。

解法

dp[i][j] = i 個目まで見て連続部分列の先頭要素として j から始まっているものの個数 とする。

貰う DP で遷移を考える。

イメージとしては LIS の DP の 2 次元版のように、配列を更新しながら計算していく。 (つまり、dp[i + 1] の計算が終われば dp[i] は不要になるタイプの DP ではないということ)

どういうことかというと、すでに計算が終了している dp[i][j] についても、今見ている dp[i] に遷移できなくなったら個数を 0 に更新するという処理を加える。 つまり、dp[i][j] について、[ ... ][ j ... ということはわかっているが、そこから今見ている要素の 1 つ前まで [ ... ][ j ... j ] で埋めることができなくなったら 0 で更新してその遷移が不可能なことを表現する。

できなくなるというのは、同じ要素が連続しすぎてダメになった場合と、A[i] != -1 なので埋める値が決まってしまった場合である。 同じ要素が連続しすぎてダメになる場合は各 dp[i] について各ループごとに先頭から順番に 1 つずつ増えていくので先頭要素を取り除けば良い。 埋める値が決まってしまった場合は残るのは 1 つなのでその手前の要素を後ろの要素を全て取り除けば良い。

遷移を考慮すると dp 配列の行方向の和と列方向の和をそれぞれ持っておく必要があるので別で用意して、以下のような高速化を行うことができる。

#include "my_template.hpp"
#include "math/static_modint.hpp"
using namespace std;

using mint = mint998;

int main() {
    INT(N);
    VEC(int, A, N);
    REP(i, N) if (A[i] != -1) A[i]--;
    {
        using pr = pair<int, mint>;  // dp[i][j] = x -> (j, x) を格納した deque
        vector<deque<pr>> dp(N + 1);
        vector<mint> dp_sum_i(N + 1);  // dp_sum_i[i] = SUM(dp[i][a] for a)
        vector<mint> dp_sum_a(N + 1);  // dp_sum_a[a] = SUM(dp[i][a] for i)

        dp[0].push_back({N, 1});
        dp_sum_i[0] += 1;
        dp_sum_a[N] += 1;

        REP(i, N) {
            mint dp_sum = 0;
            REP(j, i + 1) dp_sum += dp_sum_i[j];
            vector<mint> np(N + 1, dp_sum);
            REP(a, N) np[a] -= dp_sum_a[a];

            // apply
            REP(a, N) {
                dp[i + 1].push_back({a, np[a]});
                dp_sum_i[i + 1] += np[a];
                dp_sum_a[a] += np[a];
            }

            // shrink
            if (i == 0) {
                dp[0].clear();
                dp_sum_i[0] = 0;
                dp_sum_a[N] = 0;
            }

            REP(j, i + 1) {
                if (LEN(dp[j]) > 0 and dp[j].front().first == i - j) {
                    dp_sum_i[j] -= dp[j].front().second;
                    dp_sum_a[dp[j].front().first] -= dp[j].front().second;
                    dp[j].pop_front();
                }
            }

            if (A[i] != -1) {
                REP(j, i + 2) {
                    while (LEN(dp[j]) > 0 and dp[j].front().first != A[i]) {
                        dp_sum_i[j] -= dp[j].front().second;
                        dp_sum_a[dp[j].front().first] -= dp[j].front().second;
                        dp[j].pop_front();
                    }
                    while (LEN(dp[j]) > 0 and dp[j].back().first != A[i]) {
                        dp_sum_i[j] -= dp[j].back().second;
                        dp_sum_a[dp[j].back().first] -= dp[j].back().second;
                        dp[j].pop_back();
                    }
                }
            }
        }
        mint ans = 0;
        REP(i, N + 1) ans += dp_sum_i[i];
        print(ans);
    }
#if 0
    // 以下の O(N ^ 4) のコードを高速化する
    {
        vector dp(N + 1, vector<mint>(N + 1));
        dp[0][N] = 1;
        REP(i, N) {
            // calc dp[i + 1]
            REP(a, N) REP(j, i + 1) dp[i + 1][a] += SUM(dp[j], mint(0)) - dp[j][a];
            // shrink
            if (i == 0) dp[0][N] = 0;
            REP(j, i + 1) dp[j][i - j] = 0;
            REP(j, i + 2) {
                REP(a, N + 1) {
                    if (A[i] != -1 and a != A[i]) {
                        dp[j][a] = 0;
                    }
                }
            }
        }
        mint ans = 0;
        REP(i, N + 1) ans += SUM(dp[i], mint(0));
        print(ans);
    }
#endif
    return 0;
}

感想

全然高速化ができなくて解説を読んでもかなり理解に苦しんだが、LIS の DP の 2 次元版っぽく考えることでやっと理解できた気がする。

本番中に通せるようになるのはまだまだ先だなと感じる。(そもそもいつか通せるようになるんですかね…?)

おまけ

よく考えると上記のコードで dp_sum_i はいらない。(dp_sum の計算で使っているが、dp_sum_a で代用できる。) 以下のコードと解説の maroon さんのコードを比較してようやく z[j][i] が自分のコードにおける dp[i][j] に対応しているということがわかった。 転置するとだいぶ配列に対する 0 埋め操作が簡単になる。

#include "my_template.hpp"
#include "math/static_modint.hpp"
using namespace std;

using mint = mint998;

int main() {
    INT(N);
    VEC(int, A, N);
    REP(i, N) if (A[i] != -1) A[i]--;
    {
        using pr = pair<int, mint>;  // dp[i][j] = x -> (j, x) を格納した deque
        vector<deque<pr>> dp(N + 1);
        // vector<mint> dp_sum_i(N + 1);  // dp_sum_i[i] = SUM(dp[i][a] for a)
        vector<mint> dp_sum_a(N + 1);  // dp_sum_a[a] = SUM(dp[i][a] for i)

        dp[0].push_back({N, 1});
        // dp_sum_i[0] += 1;
        dp_sum_a[N] += 1;

        REP(i, N) {
            mint dp_sum = 0;
            // REP(j, i + 1) dp_sum += dp_sum_i[j];
            REP(a, N + 1) dp_sum += dp_sum_a[a];
            vector<mint> np(N + 1, dp_sum);
            REP(a, N) np[a] -= dp_sum_a[a];

            // apply
            REP(a, N) {
                dp[i + 1].push_back({a, np[a]});
                // dp_sum_i[i + 1] += np[a];
                dp_sum_a[a] += np[a];
            }

            // shrink
            if (i == 0) {
                dp[0].clear();
                // dp_sum_i[0] = 0;
                dp_sum_a[N] = 0;
            }

            REP(j, i + 1) {
                if (LEN(dp[j]) > 0 and dp[j].front().first == i - j) {
                    // dp_sum_i[j] -= dp[j].front().second;
                    dp_sum_a[dp[j].front().first] -= dp[j].front().second;
                    dp[j].pop_front();
                }
            }

            if (A[i] != -1) {
                REP(j, i + 2) {
                    while (LEN(dp[j]) > 0 and dp[j].front().first != A[i]) {
                        // dp_sum_i[j] -= dp[j].front().second;
                        dp_sum_a[dp[j].front().first] -= dp[j].front().second;
                        dp[j].pop_front();
                    }
                    while (LEN(dp[j]) > 0 and dp[j].back().first != A[i]) {
                        // dp_sum_i[j] -= dp[j].back().second;
                        dp_sum_a[dp[j].back().first] -= dp[j].back().second;
                        dp[j].pop_back();
                    }
                }
            }
        }
        mint ans = 0;
        REP(a, N + 1) ans += dp_sum_a[a];
        // REP(i, N + 1) ans += dp_sum_i[i];
        print(ans);
    }
    return 0;
}

【AtCoder】Introduction to Heuristics Contest A - AtCoder Contest Scheduling

問題

解法

  • 公式解説の内容を(巨大近傍法を除き)実装(124895469点)

  •  T _{0} = 1500, T _{1} = 100 に変更(コンテスト1位の shindannin さんのパラメータを採用)

  •  T を線形に減少するように変更( T = T _{0} ^{1 - t} T _{1} ^{t} から  T = T _{0} (1 - t) + T _{1} t に)

  • 差分計算における prev, next の求め方を std::set を使わずに線形探索に変更(オーダーは悪化しているが、定数倍の差で10倍近く高速化できる)

  •  T の更新を 10000 回に 1 回行う(現在時刻の取得回数を減らして高速化)

  • 戻す操作を update_point, update_swap を使わずに行う(関数の呼び出し回数を減らして高速化)

以上の改善を行った結果、126193376点になった。

// 初期解: ランダム
// 更新: 1点変更 + 2点スワップ (差分計算によってスコア計算を高速化)
// 差分計算における prev, next の求め方を std::set を使わずに線形探索に変更
// T の更新を 10000 回に 1 回行う
// 戻す操作を update_point, update_swap を使わずに行う
// 焼きなまし
// T0 = 1500, T1 = 100 (shindanninさんのパラメータ)
// TL = 1950 [ms]
// t = elapsed() / TL
// T = T0 ^ (1 - t) * T1 ^ t -> T = T0 * (1 - t) + T1 * t
// d = new_score - old_score
// d < 0 (score が減少する場合) でも確率 exp(d / T) で遷移
#include "my_template.hpp"
#include "misc/xor_shift.hpp"
#include "misc/timer.hpp"

using namespace std;

constexpr f64 TL = 1950;
constexpr f64 T0 = 1500;
constexpr f64 T1 = 100;

struct State {
    i64 score;
    vector<int> t;
};

struct Solver {
    int D;
    vector<i64> c;
    vector<vector<i64>> s;
    State st;
    Solver() {
        INT(Dtmp);
        D = Dtmp;
        VEC(i64, ctmp, 26);
        c = ctmp;
        VV(i64, stmp, D, 26);
        s = stmp;
    }
    i64 calc_score(vector<int>& t) {
        vector<int> last(26, -1);
        i64 res = 0;
        REP(i, D) {
            res += s[i][t[i]];
            last[t[i]] = i;
            REP(j, 26) res -= c[j] * (i - last[j]);
        }
        return res;
    }
    void init() {
        st.t.resize(D);
        REP(i, D) st.t[i] = rnd.rand_int(26);
        st.score = calc_score(st.t);
    }
    // 1点変更
    void update_point(int d, int q) {
        // t[d] <- q としたときの score を差分計算
        int old_q = st.t[d];
        st.score -= s[d][old_q];
        {
            int prev = d - 1, next = d + 1;
            while (prev >= 0 and st.t[prev] != st.t[d]) prev--;
            while (next < D and st.t[next] != st.t[d]) next++;

            int d1 = (d - prev), d2 = (next - d);
            st.score += c[old_q] * (d1 * (d1 - 1) / 2 + d2 * (d2 - 1) / 2);
            int d3 = (next - prev);
            st.score -= c[old_q] * (d3 * (d3 - 1) / 2);
        }
        st.t[d] = q;
        {
            int prev = d - 1, next = d + 1;
            while (prev >= 0 and st.t[prev] != st.t[d]) prev--;
            while (next < D and st.t[next] != st.t[d]) next++;

            int d1 = (d - prev), d2 = (next - d);
            st.score -= c[q] * (d1 * (d1 - 1) / 2 + d2 * (d2 - 1) / 2);
            int d3 = (next - prev);
            st.score += c[q] * (d3 * (d3 - 1) / 2);
        }
        st.score += s[d][q];
    }
    // 2点スワップ
    void update_swap(int d1, int d2) {
        int q1 = st.t[d1], q2 = st.t[d2];
        update_point(d1, q2);
        update_point(d2, q1);
    }
    void solve() {
        init();
        int count = 0;
        f64 t, T;
        while (true) {
            if (count % 10000 == 0) {
                t = timer.elapsed() / TL;
                if (t >= 1.0) break;
                T = T0 * (1.0 - t) + T1 * t;
            }
            count++;
            if (rnd.rand_int(2)) {
                // 1点変更
                int d = rnd.rand_int(D);
                int old_q = st.t[d];
                int new_q = rnd.rand_int(26);
                if (old_q == new_q) continue;

                i64 old_score = st.score;
                update_point(d, new_q);
                i64 new_score = st.score;
                if (old_score > new_score and rnd.rand_double() > exp((new_score - old_score) / T)) {
                    // もどす
                    st.t[d] = old_q;
                    st.score = old_score;
                }
            } else {
                // 2点スワップ
                int d1 = rnd.rand_int(D - 1);
                int d2 = rnd.rand_int(d1 + 1, min(D - 1, d1 + 15));
                if (st.t[d1] == st.t[d2]) continue;

                i64 old_score = st.score;
                update_swap(d1, d2);
                i64 new_score = st.score;
                if (old_score > new_score and rnd.rand_double() > exp((new_score - old_score) / T)) {
                    // もどす
                    swap(st.t[d1], st.t[d2]);
                    st.score = old_score;
                }
            }
        }
        show(count);
        show(st.score);
        REP(i, D) print(st.t[i] + 1);
    }
};

int main() {
    timer.reset();
    Solver solver;
    solver.solve();
    return 0;
}

【AtCoder】ARC170 C - Prefix Mex Sequence

問題

解法

「mex を追加する、または mex 以外を追加する」という操作において、種類数という情報だけで遷移元と遷移先を表現できる(情報として十分である)ということに注目する。

dp[i] = 数列の要素が i 種類 とすると

  • mex を追加する場合: 種類数だけで考えたら 1 種類増えるだけ(増えて M + 1 種類を超えることはないが)

  • mex を追加しない場合: 1 種類増えたりしなかったりするが常に mex 以外は M 種類なので遷移先がわかれば係数はわかる

M と N の大小関係に気を使わないと TLE してしまうため、新しく種類数の上限を K として求め、0 から K まで順番に見ていきながら適宜ダメな遷移先をカットするイメージで実装すると良い。

#include "my_template.hpp"
#include "math/static_modint.hpp"
using namespace std;

using mint = mint998;

int main() {
    INT(N, M);
    VEC(int, S, N);
    const int K = min(M + 1, N);  // 種類数の上限
    vector<mint> dp(K + 1, 0);    // dp[i] = 数列の要素が i 種類
    dp[0] = 1;
    REP(i, N) {
        vector<mint> np(K + 1, 0);
        if (S[i] == 1) {
            // mex を入れるので必ず j 種類 -> (j + 1) 種類になる
            REP(j, K + 1) {
                if (j + 1 > M + 1) break;
                if (j + 1 <= K) np[j + 1] += dp[j];
            }
        }
        if (S[i] == 0) {
            // (M + 1) 種類のうち, mex にならない M 種類のどれかを追加
            REP(j, K + 1) {
                // 種類数が増える
                if (j + 1 <= K) np[j + 1] += dp[j] * (M - j);
                // 種類数が増えない
                np[j] += dp[j] * j;
            }
        }
        swap(dp, np);
    }
    mint ans = SUM(dp, mint(0));
    print(ans);
    return 0;
}

感想

種類数に注目すれば十分という発想がなかった。

 S _{i} = 1 である要素数に注目し、適宜  S _{i} = 0 に対応する部分に値を埋めていく方法を考えていたが、埋めた要素の組合せによっては後に破綻することがある。

状態を細かく持てばできるが、状態をまとめることができなかった。

種類数が一種の不変量(?)(その情報だけで閉じているみたいな)になっていると見ることもできて、そういう意味では(広義の)不変量の発見力が足りないとも言える。

【AtCoder】ABC334 G - Christmas Color Grid 2

問題のリンク

LowLinkのアルゴリズムをベースに少し改良を行い実装する。

ところで、LowLinkを改良するアルゴリズムを考えるときには、(実際にはまとめて行っているが)まずDFS Treeをとり、そのDFS Tree上にある辺かそうでない辺(後退辺)であるかで処理を場合明けしてもう一度DFSを行っているとイメージするとわかりやすいと感じた。(DFS Treeを強く意識するのが良い。) また、ordlowdp[i]を後退辺をi回通ったものしてそれぞれdp[0]dp[1]という考え方をすると良いかもしれない。

#include "my_template.hpp"
#include "graph/read_graph.hpp"
#include "graph/low_link.hpp"
#include "math/static_modint.hpp"
using mint = mint998;
using namespace std;

int main() {
    INT(H, W);
    auto [G, id, loc] = read_grid(H, W);
    int N = LEN(G);
    LowLink lnk(G);
    mint ans = 0;
    REP(i, N) ans += lnk.count_components() - 1 + lnk.count_components_remove(i);
    ans /= N;
    print(ans);
    return 0;
}

解法2(分割統治法)

Undo可能UnionFindを使うらしい。(時間のあるときに書く)

感想

LowLinkが橋・関節点列挙のアルゴリズムというのは知っていたが、使い道をよく知らなかったので勉強になった。 今回の場合は無向グラフにおいて頂点を1つ削除したときの連結成分数の変化を前計算  O(V + E) クエリ  O(1) で行っていると言える。