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

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

AtCoder ABC 180 F - Unbranched (橙色, 600 点)

結構いろんな考え方のできる問題!

問題概要

頂点数  N、辺数  M の無向グラフであって、次の条件を満たすものの個数を 1000000007 で割ったあまりを求めよ。

  • 自己ループを持たない
  • すべての頂点の次数が 2 以下である
  • 各連結成分のサイズを並べたとき、その最大値がちょうど L である

制約

  •  M \le N \le 300

考えたこと

まず、自己ループを持たず、すべての頂点の次数が 2 以下であるとは、以下と同値になる。

  • 連結成分がすべて「パス」か「サイクル」で構成される
    • 多重辺はアリなので、サイズ 2 のサイクルはアリなことに注意

ざっくり次のような DP でできそう

  • dp[ n ][ m ] :=  N 頂点のうち  n 個を用いて、 m 本の辺が張られた状態にするような方法の個数

あとは、このような DP の具体的な方法と、「連結成分の最大サイズが  L」という条件の扱い方について考えていく。

 

解法 (1):「サイズ L 以下」から「サイズ L-1 以下」を引く

最大サイズがちょうど  L であるものを数え上げるためには、「最大サイズが  L 以下である個数」から「最大サイズが  L-1 以下である個数」を引くという定石がある。ここでは「最大サイズが  L 以下」という条件で考えることにする。

dp は、一見次のようにすればよさそうに見える。つまり、dp[ n ][ p ] からの遷移は「残りの  N - n 個から  k ( = 1, 2, \dots, L) 個を選んでパスまたはサイクルを形成する」ことによって実現できそうである。

  • dp[ n + k ][ m + k - 1 ] += dp[ n ][ m ] × C(N - n, k) × (サイズ k のパスの個数)
  • dp[ n + k ][ m + k ] += dp[ n ][ m ] × C(N - n, k) × (サイズ k のサイクルの個数)

しかしこのままではダブルカウントがある。たとえば、 N = 6 のとき

  • 先にサイズ 3 のパス (1-2-3) を作り、その後サイズ 3 のパス (4-5-6) を作る
  • 先にサイズ 3 のパス (4-5-6) を作り、その後サイズ 3 のパス (1-2-3) を作る

という場合をダブルカウントしてしまう。これを防ぐためには


残りの要素から新たにいくつかの要素を選ぶときに、添字が最小のものはかならず選ぶ


というふうにする方法がある (他の方法もある)。よって、DP を次のように修正すれば OK。

  • dp[ n + k ][ m + k - 1 ] += dp[ n ][ m ] × C(N - n - 1, k - 1) × (サイズ k のパスの個数)
  • dp[ n + k ][ m + k ] += dp[ n ][ m ] × C(N - n - 1, k - 1) × (サイズ k のサイクルの個数)

最後に、サイズ k のパスの個数は k! / 2 個 (k = 1 のときは 1 個)、サイズ k のサイクルの個数は (k-1)! / 2 個 (k = 2 のときは 1 個) となる。サイズ k のサイクルの個数を求めるときに、円順列 (k-1)! を 2 で割るのは、ひっくり返して一致するものも同一視するため。

コード

全体として計算量は  O(NML) となる。

#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 = 1000000007;
using mint = Fp<MOD>;

// Binomial coefficient
template<class T> struct BiCoef {
    vector<T> fact_, inv_, finv_;
    constexpr BiCoef() {}
    constexpr BiCoef(int n) noexcept : fact_(n, 1), inv_(n, 1), finv_(n, 1) {
        init(n);
    }
    constexpr void init(int n) noexcept {
        fact_.assign(n, 1), inv_.assign(n, 1), finv_.assign(n, 1);
        int MOD = fact_[0].getmod();
        for(int i = 2; i < n; i++){
            fact_[i] = fact_[i-1] * i;
            inv_[i] = -inv_[MOD%i] * (MOD/i);
            finv_[i] = finv_[i-1] * inv_[i];
        }
    }
    constexpr T com(int n, int k) const noexcept {
        if (n < k || n < 0 || k < 0) return 0;
        return fact_[n] * finv_[k] * finv_[n-k];
    }
    constexpr T fact(int n) const noexcept {
        if (n < 0) return 0;
        return fact_[n];
    }
    constexpr T inv(int n) const noexcept {
        if (n < 0) return 0;
        return inv_[n];
    }
    constexpr T finv(int n) const noexcept {
        if (n < 0) return 0;
        return finv_[n];
    }
};

mint solve(int N, int M, int L, const BiCoef<mint> &bc) { 
    mint ivtwo = modinv(mint(2));
    auto pathnum = [&](int n) { 
        if (n == 1) return mint(1);
        else return bc.fact(n) * ivtwo;
    };
    auto cyclenum = [&](int n) {
        if (n == 2) return mint(1);
        else return bc.fact(n-1) * ivtwo;
    };
    vector<vector<mint>> dp(N+1, vector<mint>(M+1, 0));
    dp[0][0] = 1;
    for (int i = 0; i <= N; ++i) {
        for (int j = 0; j <= M; ++j) {
            if (dp[i][j] == 0) continue;
            // add path
            for (int k = 1; k <= L && i+k <= N && j+k-1 <= M; ++k) {
                dp[i+k][j+k-1] += dp[i][j] * bc.com(N-i-1, k-1) * pathnum(k);
            }
            // add cycle
            for (int k = 2; k <= L && i+k <= N && j+k <= M; ++k) {
                dp[i+k][j+k] += dp[i][j] * bc.com(N-i-1, k-1) * cyclenum(k);
            }
        }
    }
    return dp[N][M];
}

int main() {
    BiCoef<mint> bc(110000);
    int N, M, L;
    cin >> N >> M >> L;
    cout << solve(N, M, L, bc) - solve(N, M, L-1, bc) << endl;
}

 

解法 (2):サイズ L を作ったかどうかフラグを持たせる

連結成分の最大サイズが  L という条件は、

  • どの連結成分のサイズも  L 以下とする
  • すでにサイズ  L の連結成分を作ったかどうかを DP のフラグにもたせる

という方法も考えられる。

 

解法 (3):サイズごとに考えていく

また、「単純に DP をするとダブルカウントが発生する」という部分については、サイズの昇順に DP する手もある。

  • dp[ k ][ n ][ m ] := サイズ k 以下の連結成分のみを用いて、n 個の頂点を採用し、m 辺を形成した状態にする方法の個数

このようにするのは、「 N 人のメンバーをグループ分けするときに、人数が異なるグループ分けは区別できるが、人数が同一になるグループ分けは区別できない」という性質を利用している。一般に  ak 人を  k 人グループ  a 組に分ける方法の個数は  \frac{(ak)!}{(k!)^{a}} \times \frac{1}{a!} 通り (f(k, a) とおく) となる。

  • dp[ k ][ i + ka ][ j + (k-1)a ] += dp[ k - 1 ][ i ][ j ] × C(N - i, ka) × f(k, a) × pow(サイズ k のパスの本数, a)

という感じの遷移をする (サイクルについても)。計算量は一見  O(N^{2}ML) になりそうだけど、調和級数の和のような感じになっていて  O(N \log N ML) となる。

コード

#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 = 1000000007;
//const int MOD = 998244353;
using mint = Fp<MOD>;

// Binomial coefficient
template<class T> struct BiCoef {
    vector<T> fact_, inv_, finv_;
    constexpr BiCoef() {}
    constexpr BiCoef(int n) noexcept : fact_(n, 1), inv_(n, 1), finv_(n, 1) {
        init(n);
    }
    constexpr void init(int n) noexcept {
        fact_.assign(n, 1), inv_.assign(n, 1), finv_.assign(n, 1);
        int MOD = fact_[0].getmod();
        for(int i = 2; i < n; i++){
            fact_[i] = fact_[i-1] * i;
            inv_[i] = -inv_[MOD%i] * (MOD/i);
            finv_[i] = finv_[i-1] * inv_[i];
        }
    }
    constexpr T com(int n, int k) const noexcept {
        if (n < k || n < 0 || k < 0) return 0;
        return fact_[n] * finv_[k] * finv_[n-k];
    }
    constexpr T fact(int n) const noexcept {
        if (n < 0) return 0;
        return fact_[n];
    }
    constexpr T inv(int n) const noexcept {
        if (n < 0) return 0;
        return inv_[n];
    }
    constexpr T finv(int n) const noexcept {
        if (n < 0) return 0;
        return finv_[n];
    }
};

mint solve(int N, int M, int L, const BiCoef<mint> &bc) { 
    mint ivtwo = modinv(mint(2));
    auto pathnum = [&](int n) { 
        if (n == 1) return mint(1);
        else return bc.fact(n) * ivtwo;
    };
    auto cyclenum = [&](int n) {
        if (n == 2) return mint(1);
        else return bc.fact(n-1) * ivtwo;
    };

    mint resL, resL1;
    vector<vector<mint>> dp(N+1, vector<mint>(M+1, 0)), ndp = dp;
    dp[0][0] = 1;
    for (int k = 1; k <= L; ++k) {
        // add path
        ndp.assign(N+1, vector<mint>(M+1, 0));
        for (int i = 0; i <= N; ++i) {
            for (int j = 0; j <= M; ++j) {
                if (dp[i][j] == 0) continue;
                for (int a = 0; i+a*k <= N && j+a*(k-1) <= M; ++a) {
                    ndp[i+a*k][j+a*(k-1)] += dp[i][j] * bc.com(N-i, a*k) 
                    * bc.fact(a*k) * modpow(bc.finv(k) * pathnum(k), a) * bc.finv(a);
                }
            }
        }
        swap(dp, ndp);

        // add cycle
        if (k > 1) {
            ndp.assign(N+1, vector<mint>(M+1, 0));
            for (int i = 0; i <= N; ++i) {
               for (int j = 0; j <= M; ++j) {
                    if (dp[i][j] == 0) continue;
                    for (int a = 0; i+a*k <= N && j+a*k <= M; ++a) {
                        ndp[i+a*k][j+a*k] += dp[i][j] * bc.com(N-i, a*k) 
                        * bc.fact(a*k) * modpow(bc.finv(k) * cyclenum(k), a) * bc.finv(a);
                    }
                }
            }
            swap(dp, ndp);
        }

        if (k == L) resL = dp[N][M];
        else if (k == L-1) resL1 = dp[N][M];
    }
    return resL - resL1;
}

int main() {
    BiCoef<mint> bc(110000);
    int N, M, L;
    cin >> N >> M >> L;
    cout << solve(N, M, L, bc) << endl;
}