競プロ典型 90 問 006 - Smallest Subsequence（★5）

AtCoder 競プロ典型90問競プロ典型90問とその類題競プロ典型90問難易度5 「次の要素」へのポインタを求める Greedy 辞書順文字列問題部分列各アルファベットごとに考える前処理

辞書順最小なものを求めるとき、しばしば貪欲法が有効ですね！

問題概要

英小文字のみからなる長さ $N$ の文字列 $S$ が与えられます。

$S$ の長さ $K$ の部分文字列であって、辞書順最小のものを求めてください。

制約

$1 \le K \le N \le 10^{5}$

辞書順最小 → 貪欲法！

「辞書順最小のものを求めよ」と言われたら、とにかく貪欲法！！！

まずは求めたい文字列の先頭の文字について、次のように考えます。

もし先頭の文字が a であって、長さの部分文字列が存在するならば、先頭の文字は a であると考えてよい
- 文字列 $S$ にそもそも文字 a がないとダメ
- 文字 a があったとしても、その後ろの文字数が $K-1$ 文字未満の場合はダメ
そのような部分文字列が存在しないならば、先頭の文字を b にできるかを考える
それもダメならば、先頭の文字を c にできるかを考える
それもダメならば、先頭の文字を d にできるかを考える
...

2 文字目以降の都合など、お構いなしで OK！！！！！

とにかく 2 文字目以降が全部 z になったとしても、1 文字目が小さい方が勝ちです！！！

たとえば $S$ = "bbbbazzzz" で $K = 5$ のとき、答えは "azzzz" です。

ちなみに、同じ文字が複数個あるならば、より先頭側にある文字を選ぶようにします。たとえば $S$ = "zzzazzzazzz" で $K = 5$ のとき、答えは "aazzz" です。先頭の a を選ぶとき、 $2$ 個あるうちの前の方の a を選びます。

具体的なアルゴリズム

以上の着想を具体的にアルゴリズムとしてまとめると、次のようになります。

文字列 $S$ の中で、最後にとった文字の添字を $j$ とする (j = -1 と初期化しておく)

各に対して
- 文字 a, b, ... の順に見て行って、次の条件を満たす最小のを求める
  - S[k] がその文字である
  - $k \ge j$ である
  - S[k:] は $N - i$ 文字以上ある
- そのような $k$ が存在するならば、答えに S[k] を加える、存在しなければ次の文字を見る

計算量を見積もってみましょう。 $k$ を探索する部分には $O(N)$ の計算量を要することから、 $M$ をアルファベットの種類数として $O(MNK)$ となります。

このままでは間に合いません。

前処理

前処理として、次の配列 nex を求めましょう。この配列は極めて汎用性が高いので、実にさまざまな問題で活用できます！！！

いっそライブラリ化しても良いように思います。

github.com

nex[i][c] ← 文字列 $S$ の $i$ 文字目以降で、文字 $c$ が出現する最小の添字 (存在しない場合は $N$ )

たとえば、 $S$ = "abca" のとき

nex[0][a] = 0
nex[0][b] = 1
nex[0][c] = 2
nex[0][d] = -1
nex[1][a] = 3
nex[1][b] = 1
nex[1][c] = 2
nex[2][a] = 3
nex[2][b] = -1
nex[2][c] = 2
nex[3][a] = 3
nex[3][b] = -1
nex[3][c] = -1

というようになります。

この配列は、文字列 $S$ の後ろから眺めていくことにより、次のコードのように $O(MS)$ の計算量で求められます。

またこの配列 nex を使うと、先ほどの貪欲法において $k$ を探索する部分の計算量は $O(1)$ にできますので、全体の計算量は $O(MS)$ へと改善されます。

コード (C++ と Python3)

C++

#include <iostream>
#include <vector>
#include <string>
using namespace std;

// res[i][c] := i 文字目以降で最初に文字 c が登場する index
// 存在しないときは N
vector<vector<int> > calc_next(const string &S) {
    // 文字列 S の長さ
    int N = (int)S.size();

    // 答え
    vector<vector<int> > res(N + 1, vector<int>(26, N));

    // 後ろから見る
    for (int i = N - 1; i >= 0; --i) {
        // i + 1 文字目以降の結果を一旦 i 文字にコピー
        res[i] = res[i + 1];

        // i 文字目の情報を反映させる
        res[i][S[i] - 'a'] = i;
    }
    return res;
}

int main() {
    // 入力
    int N, K;
    string S;
    cin >> N >> K >> S;
    
    // 前処理
    string res = "";
    auto nex = calc_next(S);

    // 貪欲法
    int j = -1;
    for (int i = 0; i < K; ++i) {
        for (char c = 'a'; c <= 'z'; ++c) {
            int k = nex[j + 1][c - 'a'];

            // ちゃんと K 文字が作れれば OK
            if (N - k >= K - i) {
                res += c;
                j = k;
                break;
            }
        }
    }
    cout << res << endl;
}

Python3

# res[i][c] := i 文字目以降で最初に文字 c が登場する index
# 存在しないときは N
def calc_next(S):
    # 文字列 S の長さ
    N = len(S)

    # 答え
    res = [[N] * 26 for _ in range(N + 1)]

    # 後ろから見る
    for i in range(N - 1, -1, -1):
        # i + 1 文字目以降の結果を一旦 i 文字にコピー
        for j in range(26):
            res[i][j] = res[i + 1][j]

        # i 文字目の情報を反映させる
        res[i][ord(S[i]) - ord('a')] = i

    # 答えを返す
    return res

# 入力
N, K = map(int, input().split())
S = input()

# 前処理
res = ''
nex = calc_next(S)

# 貪欲法
j = -1
for i in range(K):
    for ordc in range(26):
        k = nex[j + 1][ordc]

        # ちゃんと K 文字が作れれば OK
        if N - k >= K - i:
            res += chr(ord('a') + ordc)
            j = k
            break
print(res)

2021-10-08

競プロ典型 90 問 005 - Restricted Digits（★7）

競プロ典型90問 AtCoder 競プロ典型90問とその類題競プロ典型90問難易度7 漸化式 DP 行列行列累乗ダブリング桁DP 整数を「10倍してaを足す」で捉える DP状態:あまり sparseな問題 DP高速化:FFT FFT きたまさ法数え上げ問題 DP高速化:ダブリング

きたまさ法によく似たタイプの DP ダブリング高速化！

ただかなり難しい問題だと思うので、004 まで順調に解いていた方が、この問題で挫折しないように注意！！！

無理に解こうとせずに飛ばすのも一案だと思います......

問題概要

$K$ 種類の数字 $c_{1}, c_{2}, \dots, c_{K}$ のみを使うことで作れる $N$ 桁の正の整数のうち、 $B$ の倍数は何個あるか、 $1000000007$ で割ったあまりを答えてください。

制約

$1 \le K \le 9$
$1 \le N \le 10^{18}$
$2 \le B \le 1000$

まずは DP

今回の問題のように、「〜という数字を使って整数を作り、それが $B$ の倍数になるようにする」というタイプの問題では、DP が有効な印象があります。

drken1215.hatenablog.com

そのような DP では、整数というものを、

$10$ 倍する
$a$ を足す

という $2$ つの手続きを繰り返すことで作り上げるものと考えます。たとえば " $4649$ " という整数は

整数 $0$ からスタートして
$10$ 倍して $4$ を足すと、 $4$ になる
$10$ 倍して $6$ を足すと、 $46$ になる
$10$ 倍して $4$ を足すと、 $464$ になる
$10$ 倍して $9$ を足すと、 $4649$ になる

というように作れます。これを踏まえて今回の問題は、次の配列 dp を考えることで (計算量を無視すれば) 解けます。

dp[i][r] ← 所定の $K$ 種類の数字のみを用いた $i$ 桁の整数であって、 $B$ で割ったあまりが $r$ であるようなものの個数

このように、「 $B$ の倍数の個数」を聞かれているときに、「 $B$ で割ったあまり」を定義に含めた DP を考えると有効なことは多々あります。

この DP の遷移式は次のように考えられます。

dp[i + 1][(r * 10 + c[k]) % B] += dp[i][r]

この時点で計算量は $O(NBK)$ になります。このままでは当然 TLE。

行列累乗することで $O(B^{3} \log N)$ にはなりますが、 $B \le 1000$ なのでそれでも間に合いません。このようなとき、ダブリング解法を考えることが有効な場合があるようです。

ダブリング解法を適用すると、 $O(B^{2} \log N)$ になります。

ダブリングの考え方

今回は配列 dp[N] (dp[N][0], dp[N][1], ..., dp[N][B-1] という $B$ 個の要素からなる配列) を求めたいということになります。

上の「dp[i + 1][(r * 10 + c[k]) % B] += dp[i][r]」という漸化式は、配列 dp[0] からスタートして、dp[0] → dp[1] → dp[2] → ... → dp[N] という順に計算していくものと解釈できます。しかしこのままでは $O(N)$ 回のステップを要することになります。

しかし今回は、ダブリングすることで、 $O(\log N)$ ステップの更新で dp[N] が求められるのです！

そこで用いられるアイデアは、繰り返し二乗法と同じものです。繰り返し二乗法とは、たとえば $3^{64}$ を計算したいときに「 $3$ を $64$ 回かける」とやるのではなく、

$3^{1}$ からスタートして
$3^{2}$ (= $3^{1} \times 3^{1}$ )
$3^{4}$ (= $3^{2} \times 3^{2}$ )
$3^{8}$ (= $3^{4} \times 3^{4}$ )
$3^{16}$ (= $3^{8} \times 3^{8}$ )
$3^{32}$ (= $3^{16} \times 3^{16}$ )
$3^{64}$ (= $3^{32} \times 3^{32}$ )

という順に計算することで、わずか $6$ 回のかけざんで求められます！

この方法だと一見 $3^{N}$ を計算するのに $N$ が $2$ の冪乗である必要があるように思われるかもしれません。しかしたとえば $3^{100}$ を計算するときにも、

$100 = 64 + 32 + 4$

であることに着目して、

$3^{100} = 3^{64} \times 3^{32} \times 3^{4}$

というように計算できます。念のために、一般に通用する方法を整理しましょう。 $100$ は二進法で表すと $1100100$ となります。このことから

$100 = 2^{2} + 2^{5} + 2^{6}$

であることが導かれます。一般に $3^{N}$ を計算する場合には、 $N$ を二進法で表すことで求められます。このような方法をダブリングと呼ぶことがあります。

今回の DP の場合

以上の $3^{100}$ を計算したような方法を使えるためには、

配列 dp[i] と配列 dp[j] とから、配列 dp[i + j] を計算する

ことが高速にできればよいことになります。これさえできれば、 $3^{N}$ を計算するときのようなダブリング手法によって、配列 dp[N] を $O(\log N)$ 回のステップで計算できます。

ここで、dp[i + j] とは $i+j$ 桁の数を考えていることになります。これを「前半の $i$ 桁分」と「後半の $j$ 桁分」とに分けて考えてみましょう。

前半 $i$ 桁分を $B$ で割ったあまりが $p$ (dp[i][p] 通りあります)
後半 $j$ 桁分を $B$ で割ったあまりが $q$ (dp[j][q] 通りあります)

であるとすると、そのような $i+j$ 桁の整数を $B$ で割ったあまりは

$(p \times 10^{j} + q)$ % $B$

となります。つまり、

dp[i + j][(p * tj + q) % B] += dp[i][p] * dp[j][q]

という遷移で表せます。ここで $10^{j}$ を $B$ で割ったあまりを tj と書いています。これはまさに、「配列 dp[i] と配列 dp[j] とから配列 dp[i + j] を計算する遷移式」です。そしてこの遷移は $O(B^{2})$ の計算量でできます。

よって全体としては、ダブリングによって $O(\log N)$ 回の遷移を行うことになりますので、 $O(B^{2} \log N)$ の計算量で解けます。

コード (C++ と PyPy3)

C++

#include <iostream>
#include <vector>
using namespace std;

// MOD
constexpr int MOD = 1000000007;

// N の対数
constexpr int LOG = 62;

int main() {
    // 入力
    long long N, B, K;
    cin >> N >> B >> K;
    vector<int> C(K);
    for (auto& ci : C) cin >> ci;

    // dp[i] と dp[j] を掛け合わせて dp[i+j] を得る処理
    // tj: 10^j を B で割ったあまり
    auto mul = [&](const vector<long long>& dpi,
                   const vector<long long>& dpj,
                   long long tj) -> vector<long long> {
        vector<long long> res(B, 0);
        for (long long p = 0; p < B; ++p) {
            for (long long q = 0; q < B; ++q) {
                res[(p * tj + q) % B] += dpi[p] * dpj[q];
                res[(p * tj + q) % B] %= MOD;
            }
        }
        return res;
    };

    // ten[i]: 10^(2^i) を B で割ったあまり
    vector<long long> ten(LOG, 10);
    for (int i = 1; i < LOG; ++i) {
        ten[i] = (ten[i - 1] * ten[i - 1]) % B;
    }

    // dp[2^i][r] を doubling[i][r] で書くことにする
    vector<vector<long long>> doubling(LOG, vector<long long>(B, 0));

    // 初期化 (doubleing[0] := dp[1])
    for (int k = 0; k < K; ++k) {
        doubling[0][C[k] % B] += 1;
    }

    // ダブリング
    for (int i = 1; i < LOG; ++i) {
        doubling[i] = mul(doubling[i - 1], doubling[i - 1], ten[i - 1]);
    }

    // ダブリングした結果をもとに答えを求める
    vector<long long> res(B, 0);
    res[0] = 1;
    for (int i = 0; i < LOG; ++i) {
        // N を二の冪乗の積で表すときに、2^i を含むかどうか
        if (N & (1LL << i)) {
            res = mul(res, doubling[i], ten[i]);
        }
    }
    cout << res[0] << endl;
}

Python3 (PyPy3)

# MOD
MOD = 1000000007

# N の対数
LOG = 62

# 入力
N, B, K = map(int, input().split())
C = list(map(int, input().split()))

# dp[i] と dp[j] を掛け合わせて dp[i+j] を得る関数
# tj: 10^j を B で割ったあまり
def mul(dpi, dpj, tj):
    res = [0] * B
    for p in range(B):
        for q in range(B):
            res[(p * tj + q) % B] += dpi[p] * dpj[q]
            res[(p * tj + q) % B] %= MOD
    return res

# ten[i]: 10^(2^i) を B で割ったあまり
ten = [10] * LOG
for i in range(1, LOG):
    ten[i] = (ten[i - 1] * ten[i - 1]) % B

# dp[2^i][r] を doubling[i][r] と書くことにする
doubling = [[0] * B for _ in range(LOG)]

# 初期化 (doubling[0] = dp[1])
for k in range(K):
    doubling[0][C[k] % B] += 1

# ダブリング
for i in range(1, LOG):
    doubling[i] = mul(doubling[i - 1], doubling[i - 1], ten[i - 1])

# ダブリングした結果をもとに答えを求める
res = [0] * B
res[0] = 1
for i in range(LOG):
    if N & (1 << i):
        res = mul(res, doubling[i], ten[i])
print(res[0])

2021-08-14

AtCoder ABC 213 H - Stroll (赤色, 600 点)

AtCoder AtCoder600点 ABC-H 赤色diff FFT 分割統治法分割統治FFT DP高速化:FFT DP高速化グラフ問題数え上げ問題「総和=K」を扱う

分割統治 FFT 面白いね！！
公式解説がとても丁寧なので、備忘録程度に

問題へのリンク

問題概要

頂点数 $N$ のが与えられます。 $M$ 組の頂点対に対して、次のように無向辺を張っていきます。 $i$ 組めの頂点対に対しては

長さ 1 の辺が $p_{i, 1}$ 本
長さ 2 の辺が $p_{i, 2}$ 本
...
長さ $T$ の辺が $p_{i, T}$ 本

というように辺を張っていきます。このように作られたグラフにおいて、頂点 0 から出発して頂点 0 へと戻ってくる長さ $T$ のウォークの本数を 998244353 で割ったあまりを求めてください。

制約

$2 \le N \le 10$
$1 \le T \le 4 \times 10^{4}$

まずは DP

いかにも計算量的に間に合わないことはわかりつつも、まずは愚直な DP を立てることが大事な気がする。

dp[v][t] ← 距離 $t$ だけ進んで、頂点 $v$ に到達するような場合の数

このとき、 $v$ を終点にもつような各辺 $e = (u, v)$ に対して、

dp[v][t] += $\displaystyle \sum_{i = 0}^{t-1}$ dp[u][t - i] × p[e][i]

という遷移が立てられる。この時点で $O(MT^{2})$ にはなっているけど、とても間に合わないということで悩んでいた。そのジレンマは公式解説にまさに書かれていた。

僕もこの式が畳み込みっぽいとは思っていて、でも dp[t] の計算に dp[0], dp[1], ..., dp[t-1] も必要で難しいな...となった。それと類似の状況は下に貼る問題でも生じていた。畳み込みの添字に i < j という順序関係があって難しいという問題だった。そしてそのときに FFT が役に立ったことを思い出してはいた。なので今回も FFT 的にできないかな...という気持ちはあった。

drken1215.hatenablog.com

ただそのあとを詰められなかった。コンテスト後に解説を見て、分割統治 FFT (オンライン FFT) という考え方を知った。上に貼った問題もオンライン FFT とは違うけど、似た感じなのかなと思った。分割統治 FFT の解説については公式解説がとても丁寧なのでそちらに。

備忘録として、統治部分の式を。

t = mid, mid + 1, ..., right - 1 に対して

dp[v][t] += $\displaystyle \sum_{i = {\rm left}}^{{\rm mid}}$ dp[u][i] × p[e][t - i]

計算量は $O(M T (\log T)^{2})$ となる。なんか順序っぽい構造が入った FFT は分割統治法を使うとうまくいくことがあるよ、ということで頭に留めたいと思う。

コード

#include <bits/stdc++.h>
using namespace std;

#include "atcoder/convolution.hpp"
#include "atcoder/modint.hpp"
using namespace atcoder;
using mint = modint998244353;

int main() {
    // 入力
    int N, M, T;
    cin >> N >> M >> T;
    vector<int> A(M * 2), B(M * 2);
    vector<vector<mint>> P(M * 2, vector<mint>(T + 1));
    for (int i = 0; i < M; ++i) {
        cin >> A[i] >> B[i];
        --A[i], --B[i];
        B[i + M] = A[i], A[i + M] = B[i];
        for (int t = 1; t <= T; ++t) {
            long long p;
            cin >> p;
            P[i + M][t] = P[i][t] = p;
        }
    }

    // 分割統治 FFT
    vector<vector<mint>> dp(N, vector<mint>(T + 1, 0));
    dp[0][0] = 1;
    auto rec = [&](auto self, int left, int right) -> void {
        if (right - left <= 1) return;

        int mid = (left + right) / 2;

        // まず左半分を更新
        self(self, left, mid);

        // 左半分から右半分への遷移を更新
        for (int e = 0; e < M * 2; ++e) {
            int u = A[e], v = B[e];

            vector<mint> L(mid - left, 0), R(right - left, 0);
            for (int t = left; t < mid; ++t) L[t - left] = dp[u][t];
            for (int t = 0; t < right - left; ++t) R[t] = P[e][t]; 
            
            auto seki = convolution(L, R);
            for (int t = mid; t < right; ++t) {
                dp[v][t] += seki[t - left];
            }
        }
        
        // 最後に右半分を更新
        self(self, mid, right);
    };
    rec(rec, 0, T + 1);
    cout << dp[0][T].val() << endl;
}

2021-08-12

AOJ 2345 Network Reliability (JAG 冬コン 2011 F) (700 点)

AOJ JAG JAG春冬コングラフ問題指数探索系問題グラフの頂点集合の部分集合を走査する確率 bitDP O(3^N) 連続量問題グラフのconnectivity

ABC 213 G の類題に挙がっていたのでやってみた！

問題へのリンク

問題概要

頂点数 $N$ 、辺数 $M$ の単純無向グラフが与えられます。また $0$ 以上 $100$ 以下の整数 $p$ が与えられます。

グラフの各辺は $p$ % の確率でそれぞれ独立に消失します。グラフ全体が連結である確率を求めてください。

制約

$1 \le N \le 14$
$0 \le M \le 100$
$0 \le P \le 100$

考えたこと

$M \le 100$ なので、辺の消失パターン ( $2^{M}$ 通りある) すべてを調べる方法では間に合わない。

しかしながら、このようにグラフ上の指数探索系の問題では、「ある頂点 $v$ を一つ固定して、その周りの様子で場合分けする」のが定石だと思われる。ここでは、頂点 $v$ を含む連結成分 $S$ がどのようになるかで場合分けする考え方で問題を解いていこう。ここで、

$f(S)$ = もとのグラフのうち、頂点集合を $S$ のみに限ったグラフを考えたとき、それが連結になる確率

とする。これを求める漸化式を立てていく。なお答えは $f(V)$ と表される。

$f(S)$ を求める漸化式

$f(S)$ をいきなり求めることが難しいので、余事象を考えて全体から引くことを考えよう。つまり、 $v \in T \subsetneq S$ であるような $T$ をすべて考えて引いていく方針で考えてみよう。

$f(S)$ を求めるために、 $v \in S$ となる頂点 $v$ を選び、 $S$ 内部で $v$ を含む連結成分 $T$ で場合分けして考えることにすると、

$f(S) = 1.0 - \sum_{v \in T \subsetneq S} f(T) \times p^{|C(T, S-T)|}$

となる。ここで $C(X, Y)$ は、頂点集合 $X$ と頂点集合 $Y$ との間にある辺集合を表します。

これは「部分集合の部分集合」を考えていくような遷移になっているので、 $O(3^{N})$ な bit DP となります。なお、 $C$ については、頂点集合 $S$ に含まれる辺の本数を $C(S)$ として、

$C(X, Y) = e(X \cup Y) - e(X) - e(Y)$

と求められます。 $e$ については $O(N 2^{N})$ で計算できます。

drken1215.hatenablog.com

コード

#include <bits/stdc++.h>
using namespace std;

int main() {
    // 入力
    int N, M;
    double P;
    cin >> N >> M >> P;
    P /= 100;
    vector<vector<bool>> G(N, vector<bool>(N, false));
    for (int i = 0; i < M; ++i) {
        int a, b;
        cin >> a >> b;
        --a, --b;
        G[a][b] = G[b][a] = true;
    }

    // P の冪乗を確保しておく
    vector<double> power(M + 1, 1);
    for (int i = 0; i < M; ++i) power[i+1] = power[i] * P;

    // e(S) を求める
    vector<long long> e(1<<N, 0);
    for (int u = 0; u < N; ++u)  {
        for (int v = 0; v < N; ++v) {
            if (G[u][v]) {
                int S = (1<<u) | (1<<v);
                e[S] = 1;
            }
        }
    }
    for (int i = 0; i < N; ++i) {
        for (int S = 0; S < (1<<N); ++S) {
            if (S & (1<<i)) {
                e[S] += e[S ^ (1<<i)];
            }
        }
    }

    // O(3^N) DP
    vector<double> f(1<<N, 1.0);
    for (int S = 1; S < (1<<N); ++S) {
        // S に含まれる頂点 v を 1 つ選ぶ
        int v = -1;
        for (v = 0; v < N; ++v) if (S & (1<<v)) break;

        // S の部分集合を走査
        for (int T = S - 1; T >= 0; --T) {
            T &= S;

            if (T & (1 << v)) {
                int c = e[S] - e[T] - e[S - T];
                f[S] -= f[T] * power[c];
            }
        }
    }
    cout << fixed << setprecision(10) << f[(1<<N)-1] << endl;
}

2021-08-12

AtCoder ABC 213 G - Connectivity 2 (橙色, 600 点)

AtCoder AtCoder600点 ABC-G 橙色diff 部分グラフ列挙問題数え上げ問題グラフ問題指数探索系問題 bitDP 高速ゼータ変換グラフのconnectivity グラフ・盤面・数列の個数の数え上げ一般化した問題を解く包除原理除原理グラフの頂点集合の部分集合を走査する各kに対して場合分け SubsetConvolution 高速畳み込み計算累積和 O(3^N) 補集合を考える

面白かった！！　より一般化した問題 (グラフの頂点集合の任意の部分集合に対して、それらを連結にする辺の選び方の数え上げ) を考えた方が考えやすいね。

問題へのリンク

問題概要

頂点数 $N$ 、辺数 $M$ の単純無向グラフ $G$ が与えられます。 $G$ の辺集合の部分集合 ( $2^{M}$ 通りある) たちを考える。

各頂点 $k = 1, 2, \dots, N-1$ に対して、頂点 $0$ と頂点 $k$ とが連結であるような部分集合が何個あるかを 998244353 で割ったあまりを求めよ。

制約

$1 \le N \le 17$

連結グラフの数え上げ問題へ

頂点数の小さい無向グラフが与えられたときに、全頂点が連結であるような部分グラフ (辺集合の部分集合) の個数を数え上げる問題は割と典型な気がする。そしてそれを求める過程では、頂点集合のすべての部分集合 $S$ に対しても、 $S$ が連結であるような部分グラフ ( $S$ 内部の辺のみを考える) の個数をすべて求めることになるのだ。

つまり、次の関数 $f(S)$ が求められることになる (後述)。

$f(S)$ = 頂点集合 $S$ ( $\subset V$ ) について、その内部にある辺たちの部分集合のうち、 $S$ 全体が連結であるようなものの個数

そして今回の問題は $S$ として、とくに $S = (0, k)$ をとれば、よさそうだ。ただし注意点として、 $f(S)$ はあくまで頂点集合 $S$ 内部の辺のみを考えている。今回求めたいのは、もとのグラフ全体の辺たちの部分集合のうち、 $S$ 全体を連結にするようなものの個数なのだ。

そこで、次のような場合分けをしよう。 $S$ を含む連結成分で場合分けするのだ。

$S$ を含む連結成分をなる頂点の部分集合が $T$ であるとき、
$V - T$ に含まれる辺の本数を $e(V - T)$ とすると
そのような場合の数は通り
- $T$ と $V - T$ との間に辺を引くことは考えなくてよい

というふうに考えられる。これを各 $T$ に対して総和をとればよいということになる。以上より、もとの問題は $f(T)$ を求める問題へと帰着された。

なお、上のことを実現するためには、各部分集合 $S$ に対して $e(S)$ を計算する必要もある。愚直にやると $O(M 2^{N})$ になる。今回はこれでも十分間に合うが、高速ゼータ変換を用いると $O(N 2^{N})$ になる。

高速ゼータ変換

具体的には、頂点集合の部分集合 $S$ に対して

$h(S)$ = $S$ のサイズが 2 であり、それらがグラフ $G$ のある辺の両端点になっているとき 1、それ以外のとき 0

という関数 $h$ を定義すると

$e(S) = \sum_{T \subset S} h(T)$

と表せることに注意しよう。これはまさに高速ゼータ変換を使える形になっている。下のコードのような感じで in-place にできる。ここで、初期状態の $e$ は $h$ を表しているものとする。これによって、 $e(S)$ の計算は $O(N 2^{N})$ となる。

// 初期状態の e は h を表すものとする
// e に対して in-place に累積和をとるような更新をする
for (int i = 0; i < N; ++i) {
    for (int S = 0; S < (1<<N); ++S) {
        if (S & (1 << i)) {
             e[S] += e[S ^ (1 << i)];
        }
    }
}

$f(S)$ を求める

それでは、次の $f(S)$ を求めていきます。

$f(S)$ = 頂点集合 $S$ ( $\subset V$ ) について、その内部にある辺たちの部分集合のうち、 $S$ 全体が連結であるようなものの個数

このように、グラフの頂点集合の部分集合に対する何かを計算するときには、1 つ頂点 $v$ を選んでその頂点 $v$ に対する挙動で場合分するという方針がしばしば有効なイメージがある (例：彩色数)。

ここでも $v \in S$ を一つとる。また、 $f(S)$ を求めるためには、逆に条件を満たさないものを数え上げて引くことにする。条件を満たさないとき、 $v$ を含む連結成分は $S$ の真部分集合のいずれかになる。それを $T$ としたときは

$f(T) \times 2^{e(S - T)}$ 通り

となる。よって

$\displaystyle f(S) = 2^{e(S)} - \sum_{T \subsetneq S, v \in T} f(T) \times 2^{e(S - T)}$

と求められることがわかった (除原理)。これは、 $f(S)$ が $O(3^{N})$ な bitDP で計算できることを意味している！！！　あとは丁寧に実装していけば OK。

なお、Subset Convolution という超高度な技を使えば $O(N^{2} 2^{N})$ で計算できるみたい (あまりよく分かってない...)。

コード

計算量は

$e(S)$ を求める： $O(N 2^{N})$ (高速ゼータ変換)
$f(S)$ を求める： $O(3^{N})$ (Subset Convolution を用いれば $O(N^{2} 2^{N})$

となる。

#include <bits/stdc++.h>
using namespace std;

// modint
template<int MOD> struct Fp {
    long long val;
    constexpr Fp(long long v = 0) noexcept : val(v % MOD) {
        if (val < 0) val += MOD;
    }
    constexpr int getmod() const { return MOD; }
    constexpr Fp operator - () const noexcept {
        return val ? MOD - val : 0;
    }
    constexpr Fp operator + (const Fp& r) const noexcept { return Fp(*this) += r; }
    constexpr Fp operator - (const Fp& r) const noexcept { return Fp(*this) -= r; }
    constexpr Fp operator * (const Fp& r) const noexcept { return Fp(*this) *= r; }
    constexpr Fp operator / (const Fp& r) const noexcept { return Fp(*this) /= r; }
    constexpr Fp& operator += (const Fp& r) noexcept {
        val += r.val;
        if (val >= MOD) val -= MOD;
        return *this;
    }
    constexpr Fp& operator -= (const Fp& r) noexcept {
        val -= r.val;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr Fp& operator *= (const Fp& r) noexcept {
        val = val * r.val % MOD;
        return *this;
    }
    constexpr Fp& operator /= (const Fp& r) noexcept {
        long long a = r.val, b = MOD, u = 1, v = 0;
        while (b) {
            long long t = a / b;
            a -= t * b, swap(a, b);
            u -= t * v, swap(u, v);
        }
        val = val * u % MOD;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr bool operator == (const Fp& r) const noexcept {
        return this->val == r.val;
    }
    constexpr bool operator != (const Fp& r) const noexcept {
        return this->val != r.val;
    }
    friend constexpr istream& operator >> (istream& is, Fp<MOD>& x) noexcept {
        is >> x.val;
        x.val %= MOD;
        if (x.val < 0) x.val += MOD;
        return is;
    }
    friend constexpr ostream& operator << (ostream& os, const Fp<MOD>& x) noexcept {
        return os << x.val;
    }
    friend constexpr Fp<MOD> modpow(const Fp<MOD>& r, long long n) noexcept {
        if (n == 0) return 1;
        if (n < 0) return modpow(modinv(r), -n);
        auto t = modpow(r, n / 2);
        t = t * t;
        if (n & 1) t = t * r;
        return t;
    }
    friend constexpr Fp<MOD> modinv(const Fp<MOD>& r) noexcept {
        long long a = r.val, b = MOD, u = 1, v = 0;
        while (b) {
            long long t = a / b;
            a -= t * b, swap(a, b);
            u -= t * v, swap(u, v);
        }
        return Fp<MOD>(u);
    }
};

const int MOD = 998244353;
using mint = Fp<MOD>;





int main() {
    // 入力
    int N, M;
    cin >> N >> M;
    vector<vector<bool>> G(N, vector<bool>(N, false));
    for (int i = 0; i < M; ++i) {
        int a, b;
        cin >> a >> b;
        --a, --b;
        G[a][b] = G[b][a] = true;
    }

    // 2 の冪乗を確保しておく
    vector<mint> two(N*N + 1, 1);
    for (int i = 0; i < N*N; ++i) two[i+1] = two[i] * 2;

    // e(S) を求める
    vector<long long> e(1<<N, 0);
    for (int u = 0; u < N; ++u)  {
        for (int v = 0; v < N; ++v) {
            if (G[u][v]) {
                int S = (1<<u) | (1<<v);
                e[S] = 1;
            }
        }
    }
    for (int i = 0; i < N; ++i) {
        for (int S = 0; S < (1<<N); ++S) {
            if (S & (1<<i)) {
                e[S] += e[S ^ (1<<i)];
            }
        }
    }

    // f(S) を求める
    vector<mint> f(1<<N, 0);
    for (int S = 0; S < (1<<N); ++S) {
        if (S == 0) continue;
        
        // S に含まれる頂点 v を 1 つ選ぶ
        int v = -1;
        for (v = 0; v < N; ++v) if (S & (1<<v)) break;

        // 全体
        f[S] = two[e[S]];
        
        // S の部分集合を走査していく
        for (int T = S - 1; T >= 0; --T) {
            T &= S;

            // T が v を含む場合のみ
            if (T & (1 << v)) {
                f[S] -= f[T] * two[e[S ^ T]];
            }
        }
    }

    // 集計する
    int V = (1<<N) - 1;
    for (int k = 1; k < N; ++k) {
        mint res = 0;
        for (int S = 0; S < (1<<N); ++S) {
            if (!(S & (1<<0))) continue;
            if (!(S & (1<<k))) continue;
            res += f[S] * two[e[V ^ S]];
        }
        cout << res << endl;
    }   
}

けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

競プロ典型 90 問 006 - Smallest Subsequence（★5）

問題概要

制約

辞書順最小 → 貪欲法！

具体的なアルゴリズム

前処理

コード (C++ と Python3)

C++

Python3

競プロ典型 90 問 005 - Restricted Digits（★7）

問題概要

制約

まずは DP

ダブリングの考え方

今回の DP の場合

コード (C++ と PyPy3)

C++

Python3 (PyPy3)

AtCoder ABC 213 H - Stroll (赤色, 600 点)

問題概要

制約

まずは DP

コード

AOJ 2345 Network Reliability (JAG 冬コン 2011 F) (700 点)

問題概要

制約

考えたこと

$f(S)$ を求める漸化式

コード

AtCoder ABC 213 G - Connectivity 2 (橙色, 600 点)

問題概要

制約

連結グラフの数え上げ問題へ

高速ゼータ変換

$f(S)$ を求める

コード