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

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

CPSCO2019 Session3 F - Flexible Permutation (600 点設定)

このセットの writer をしていました

問題概要

 1, 2, \dots, N を並び替えてできる数列  P_{1}, P_{2}, \dots, P_{N} N! 通りあるが、そのうち以下の条件を満たすものが何通りあるかを、1000000007 で割ったあまりを求めよ。

  •  P_{i} \gt i を満たすような  i がちょうど  A
  •  P_{i} \lt i を満たすような  i がちょうど  B
  •  P_{i} = i を満たすような  i がちょうど  N - A - B

制約

  •  1 \le N \le 300

 

解法 (1): 挿入 DP (O(N3))

まずは想定解法から。順列の数え上げは

 1, 2, \dots, i の順列に対して、 i+1 を上手に追加することで  1, 2, \dots, i+1 の順列を作る」

という考え方がよくとられる。たとえば「3, 2, 5, 4, 1」という順列に 6 を加えるとき、以下のように 6 通りの方法がある。

  • 3, 2, 5, 4, 1, 6 (末尾に追加)
  • 6, 2, 5, 4, 1, 3 (1 番目の要素と swap)
  • 3, 6, 5, 4, 1, 2 (2 番目の要素と swap)
  • 3, 2, 6, 4, 1, 5 (3 番目の要素と swap)
  • 3, 2, 5, 6, 1, 4 (4 番目の要素と swap)
  • 3, 2, 5, 4, 6, 1 (5 番目の要素と swap)

このような操作を各「 1, 2, \dots, i の順列」に対して実施することで、 1, 2, \dots, i+1 の順列」をすべて生成することができる。

これを踏まえて次のような DP をする。

  •  {\rm dp} \lbrack i \rbrack \lbrack a \rbrack \lbrack b \rbrack =  1, 2, \dots, i の順列のうち、 P_{i} \gt i を満たすものが  a 個で、 P_{i} \lt i を満たすものが  b 個となる場合の数

このとき、遷移は次のように考えられる。

  • 末尾に追加するとき: a, b は変化なし
  •  P_{i} \gt i を満たす  P_{i} ( a 箇所) と swap するとき: b が増加
  •  P_{i} \lt i を満たす  P_{i} ( b 箇所) と swap するとき: a が増加
  •  P_{i} = i を満たす  P_{i} ( i - a - b 箇所) を swap するとき: a, b がともに増加

それぞれ式に表すと次のようになる。

 {\rm dp}\lbrack i+1 \rbrack \lbrack a \rbrack \lbrack b \rbrack +=  {\rm dp} \lbrack i \rbrack \lbrack a \rbrack \lbrack b \rbrack
 {\rm dp} \lbrack i+1 \rbrack \lbrack a \rbrack \lbrack b+1 \rbrack +=  {\rm dp} \lbrack i \rbrack \lbrack a \rbrack \lbrack b \rbrack \times a
 {\rm dp} \lbrack i+1 \rbrack \lbrack a+1 \rbrack \lbrack b \rbrack +=  {\rm dp} \lbrack i \rbrack \lbrack a \rbrack \lbrack b \rbrack \times b
 {\rm dp} \lbrack i+1 \rbrack \lbrack a+1 \rbrack \lbrack b+1 \rbrack +=  {\rm dp} \lbrack i \rbrack \lbrack a \rbrack \lbrack b \rbrack \times (i-a-b)

計算量は  O(N^{3})

#include <iostream>
#include <fstream>
#include <string>
#include <cstring>
using namespace std;

const int MOD = 1000000007;
const int MAX = 310;
long long dp[MAX][MAX][MAX];

void add(long long &a, long long b, long long c) {
    a += b * c % MOD;
    a %= MOD;
}

long long solve(long long N, long long A, long long B) {
    memset(dp, 0, sizeof(dp)); // (1...n) の並び替えのうち、a 個が P_i > i で b 個が P_i < i なもの
    dp[0][0][0] = 1;
    for (int n = 0; n <= N; ++n) {
        for (int a = 0; a <= n; ++a) {
            for (int b = 0; a + b <= n; ++b) {
                // i+1 を i+1 番目に
                add(dp[n+1][a][b], dp[n][a][b], 1);
                
                // i+1 を a 族と swap
                add(dp[n+1][a][b+1], dp[n][a][b], a);
                
                // i+1 を b 族と swap
                add(dp[n+1][a+1][b], dp[n][a][b], b);
                
                // i+1 を c 族と swap
                add(dp[n+1][a+1][b+1], dp[n][a][b], n - a - b);
            }
        }
    }
    return dp[N][A][B];
}

int main() {
    long long N, A, B;
    cin >> N >> A >> B;
    cout << solve(N, A, B) << endl;
}

 

解法 (2):箱根駅伝 DP (O(N3))

順列の数え上げといえば、箱根駅伝のような DP をする方法も代表的。

まず、 P_{i} = i となる箇所は  N - A - B 箇所であることがわかっているので、それらをあらかじめ決め打つことにする。決め打つ方法は  {}_{N}{\rm C}_{N-A-B} 通りある。

これをあとで掛けることにすると、次の問題に帰着される。


 1, 2, \dots, A+B の順列  P のうち、以下の条件を満たすものの個数を求めよ。

  •  P_{i} \gt i を満たす  i A
  •  P_{i} \lt i を満たす  i B

ここから先は、箱根駅伝 DP をする。

drken1215.hatenablog.com

まず順列を「左側が  N 個」「右側が  N 個」という要素の並んだマッチングとみなすことにする。そして、「左側の  n 個」と「右側の  n 個」だけを考えた状態を考えて、次のように DP を定義する。

  •  {\rm dp}\lbrack n \rbrack \lbrack a \rbrack \lbrack b \rbrack := 左側の [tex: n 個の要素と、右側の  n 個の要素の間で、すでにマッチングが  a+b 本成立していて、そのうち左側から右側へと下がっている辺が  a 本あり、左側から右側へと上がっている辺が  b 本あるような場合の数

そうすると、DP 遷移は次のようになる。

  •  {\rm dp}\lbrack n+1 \rbrack \lbrack a+1 \rbrack \lbrack b+1 \rbrack +=  {\rm dp}\lbrack n \rbrack \lbrack a \rbrack \lbrack b \rbrack \times (n-a-b)^{2}
    • 左側の頂点  n+1 も、右側の頂点  n+1 も、自分より上側の頂点とマッチングする場合)
  •  {\rm dp}\lbrack n+1 \rbrack \lbrack a+1 \rbrack \lbrack b \rbrack +=  {\rm dp}\lbrack n \rbrack \lbrack a \rbrack \lbrack b \rbrack \times (n-a-b)
    • 右側の頂点  n+1 のみが、自分より上側の左側頂点とマッチングする場合
  •  {\rm dp}\lbrack n+1 \rbrack \lbrack a \rbrack \lbrack b+1 \rbrack +=  {\rm dp}\lbrack n \rbrack \lbrack a \rbrack \lbrack b \rbrack \times (n-a-b)
    • 左側の頂点  n+1 のみが、自分より上側の右側頂点とマッチングする場合
  •  {\rm dp}\lbrack n+1 \rbrack \lbrack a \rbrack \lbrack b \rbrack +=  {\rm dp}\lbrack n \rbrack \lbrack a \rbrack \lbrack b \rbrack \times (n-a-b)
    • 新たなマッチングを形成しない場合

この場合も計算量は  O(N^{3}) となる。

 

解法 (3):撹乱順列に対する漸化式の応用 (O(N2))

解法 (2) で考えた次の問題を考える。


 1, 2, \dots, A+B の順列  P のうち、以下の条件を満たすものの個数を求めよ。

  •  P_{i} \gt i を満たす  i A
  •  P_{i} \lt i を満たす  i B

この問題は「撹乱順列」の応用と言える。撹乱順列のうち、「= でないところ」の左右へのズレの個数まで指定したものを数え上げよ、というわけだ。でも撹乱順列の個数を求める漸化式解法はよく知られている。今回の問題にそれを応用できそうだ。

mathtrain.jp

  •  {\rm dp}\lbrack n \rbrack \lbrack r \rbrack :=  1, 2, \dots, n の撹乱順列のうち、 P_{i} \gt i を満たす  i r 個であるようなものの個数

とする。集める DP で考えることにする。そして、 P_{a} = N P_{N} = b として、 a, b の様相で場合分けする。

 
 a = b であるとき
このとき、 a, n を除いた  n-2 個の順列の場合に帰着される。 a の選び方が  n-1 通りあることから、次のようになる。

 {\rm dp}\lbrack n \rbrack \lbrack r \rbrack +=  {\rm dp}\lbrack n-2 \rbrack \lbrack r-1 \rbrack \times (n-1)

 
 a \lt b であるとき
この場合は、以下の場合と一対一に対応することが示せる。

  •  1, 2, \dots, n-1 の順列のうち
  •  P_{i} \lt i となる箇所が  r 箇所あるようなもの  P' と、
  • その  r 箇所のうちから 1 つ  i を選ぶ方法のペア  (P', i)

具体的には、 n 個の順列において  n を削除して代わりに  P_{a} = b とすれば、上記の条件を満たす  n-1 個の順列に帰着される。逆も同様。よって、

 {\rm dp}\lbrack n \rbrack \lbrack r \rbrack +=  {\rm dp}\lbrack n-1 \rbrack \lbrack r \rbrack \times r

となる。

 
 a \gt b であるとき
この場合は、以下の場合と一対一に対応することが示せる。

  •  1, 2, \dots, n-1 の順列のうち
  •  P_{i} \lt i となる箇所が  r-1 箇所あるようなもの  P' と、
  • そうなっていない  n-r 箇所のうちから 1 つ  i を選ぶ方法のペア  (P', i)

具体的には、 n 個の順列において  n を削除して代わりに  P_{a} = b とすれば、上記の条件を満たす  n-1 個の順列に帰着される。今回は、 a \gt b であるにもかかわらず  a \lt N であることから、 n-1 個の順列 (上記) と  n 個の順列とで、 P_{i} \gt i を満たす  i の個数が 1 ずれることに注意。逆も同様。よって、

 {\rm dp}\lbrack n \rbrack \lbrack r \rbrack +=  {\rm dp}\lbrack n-1 \rbrack \lbrack r-1 \rbrack \times (n-r)

となる。

コード

以上の解法によって、計算量は  O(N^{2}) に落ちる!!!

#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];
    }
};

int main() {
    int N, A, B;
    cin >> N >> A >> B;
    BiCoef<mint> bc(N+1);
    vector<vector<mint>> dp(A+B+1, vector<mint>(A+B+1, 0));
    dp[0][0] = 1;
    for (int n = 1; n <= A+B; ++n) {
        for (int r = 0; r <= A+B; ++r) {
            if (n >= 2 && r >= 1) dp[n][r] += dp[n-2][r-1] * (n-1);
            dp[n][r] += dp[n-1][r] * r;
            if (r >= 1) dp[n][r] += dp[n-1][r-1] * (n-r);
        }
    }
    cout << dp[A+B][A] * bc.com(N, N-A-B) << endl;
}

 

解法 (4):撹乱順列に対する包除原理の応用 (O(N2))

撹乱順列の数え上げには、漸化式を用いる解法とともに、包除原理を適用する解法もある。とすれば今回の問題にもそのような解法もありそうだ。

アルメリアさんの記事に詳しい解説がある!

betrue12.hateblo.jp

AOJ 2566 Restore Calculation (JAG 模擬地区 2013 B) (350 点)

これが生まれて初めての作問だった!!!
AOJ-ICPC で ☆4 個ついてて嬉しい!

問題概要

虫食算が与えられる。解の個数を 1000000007 で割ったあまりを求めたい。

より正確には長さ  N の 3 個の文字列  A, B, C が与えられる。これらは '?' か '0'〜'9' の文字で構成されている。 A, B, C の '?' に '0'〜'9' の数値を当てはめる方法のうち

  • 完成された  A, B, C を整数値としてみたときに、 A + B = C が成立
  •  A, B, C の先頭は '0' ではない

という条件を満たすものの個数を求めよ。

制約

  •  1 \le N \le 50
  •  A, B, C の先頭の文字は '0' ではない

サンプル

3?4
12?
5?6

これの解は以下の 2 通りある。

  • 384 + 122 = 506
  • 394 + 122 = 516

解法

1 の位から順番に数値を当てはめていく (文字列としては末尾から見ていくことになるので、あらかじめ reverse しておく)。この際に「繰り上がり」の発生についても考慮していく。

さて、d 桁目まで数値を当てはめた状態では、以下のいずれかの状態がありうる。

  • 繰り上がりが 0 (繰り上がりなし)
  • 繰り上がりが 1

これらの情報を持って DP しよう。

  • dp[ d ][ 0 or 1 ] := d 桁目まで数値を当てはめたときに、繰り上がりが (0 or 1) であるようにするときの場合の数

最終的な答えは dp[ N ][ 0 ] となる (最後は繰り上がってはいけない)。

DP

DP の遷移は、一見実装が大変だけど、次のように for 文三重ループを回してしまうと楽できる。下の実装では

  • num0:d 桁目から d+1 桁目にかけて繰り上がりが発生しないように、d 桁目に数値を入れる方法の数
  • num1:d 桁目から d+1 桁目にかけて繰り上がりが発生するように、d 桁目に数値を入れる方法の数

としている。これらの値を用いて、DP 遷移を行う。計算量は定数倍が重めの  O(N) となる。

int num0 = 0, num1 = 0;
for (int a = 0; a <= 9; ++a) {
     for (int b = 0; b <= 9; ++b) {
          for (int c = 0; c <= 9; ++c) {
               if (d 桁目で、A に a, B に b, C に c を当てはめるのが valid) {
                      if (繰り上がりなし) ++num0;
                      else ++num1;
               }
          }
     }
}

コード

#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>;

mint solve(string A, string B, string C) {
    int N = (int)A.size();
    reverse(A.begin(), A.end());
    reverse(B.begin(), B.end());
    reverse(C.begin(), C.end());

    auto isvalid = [&](int a, int b, int c, int d, int kuriagari) {
        // 先頭に 0 はダメ
        if (d == N-1) if (a == 0 || b == 0 || c == 0) return false;

        // すでに入っている数値と矛盾がないか
        if (A[d] != '?' && A[d] != (char)('0'+a)) return false;
        if (B[d] != '?' && B[d] != (char)('0'+b)) return false;
        if (C[d] != '?' && C[d] != (char)('0'+c)) return false;

        // ok
        return true;
    };
    vector<vector<mint>> dp(N+1, vector<mint>(2, 0));
    dp[0][0] = 1;
    for (int d = 0; d < N; ++d) {
        for (int kuriagari = 0; kuriagari <= 1; ++kuriagari) {
            int num0 = 0, num1 = 0;
            for (int a = 0; a <= 9; ++a) {
                for (int b = 0; b <= 9; ++b) {
                    for (int c = 0; c <= 9; ++c) {
                        if (isvalid(a, b, c, d, kuriagari)) {
                            if (a + b + kuriagari == c) ++num0;
                            else if (a + b + kuriagari == c + 10) ++num1;
                        }
                    }
                }
            }
            dp[d+1][0] += dp[d][kuriagari] * num0;
            dp[d+1][1] += dp[d][kuriagari] * num1;
        }
    }
    return dp[N][0];
}

int main() {
    string A, B, C;
    while (cin >> A) {
        if (A == "0") break;
        cin >> B >> C;
        cout << solve(A, B, C) << endl;
    }   
}

AtCoder ARC 106 E - Medals (赤色, 700 点)

高速ゼータ変換を思いつくのに時間かかった

問題概要

あなたは  N 人の従業員を持つ店の店長です。  i 番目の従業員は今日から「 A_{i} 日連続で働いた後  A_{i} 日連続で休む」ことを繰り返します。

あなたは今日から毎日出勤し、その日に出勤している従業員の中から 1 人選び、メダルを 1 枚配ります。(その日に出勤している従業員が 1 人もいなければ何もしません)

全ての従業員に  K 枚以上のメダルを配るためには、最短で何日かかるでしょうか。

制約

  •  1 \le N \le 18
  •  1 \le K, A_{i} \le 10^{5}

考えたこと

最初に、「実は全員が × な日を除けば毎日埋められるんじゃないか」と想像してみた。しかしこれはすぐに反例が見つかった。

  •  A = (1, 1, 1, 1, 1, 6)
  •  K = 1

のとき、上記の仮説が正しければ 6 日間で実現できることになる。しかし実際は  (1, 1, 1, 1, 1) のところで、2 日置きに揃って x になる壁が分厚くて、10 日間を要する。

上記のことから、ボトルネックとなるような「従業員の部分集合」を特定できれば、その従業員全員に行き渡らせるのに必要な人数が答えになるのではないかと考えた!!具体的には、 D 日以内に終えられるかどうかは、次のように判定できそうに思えた。


従業員の部分集合  S に対して、 D 日間のうち「全員が ×」という日が  d(S) 日分であったとする。

このとき、任意の  S に対して  D - d(S) \ge K|S| であるならば、可能である。


そしてこれが正しいことが Hall の定理から示せるのだ。Hall の定理とは、頂点集合  U, V とに分割された二部グラフに対して、以下が同値だと主張している!

  •  U をすべてカバーするマッチングが存在する
  • (Hall の条件)  U の任意の部分集合  A に対して、 \Gamma(A) \ge |A| が成立する ( \Gamma(A) A に含まれるいずれかの頂点に隣接している頂点の個数)

そして従業員の問題に戻り、左側を「従業員 ×  K 日分」として、右側を「日程」とすると、上記の予想は Hall の条件に該当する。

よって結論として、下の条件を満たすような  D の最小値を求めればよい。二分探索がピッタリとなる。なお、 D = 2NK とすると条件を満たすことに注意する。これを二分探索の上限に設定できる。


任意の  S に対して  D - d(S) \ge K|S| が成立する


判定パート

残りは、判定問題を解くこと。

  •  f(S) :=  S に該当する従業員のみ x で、それ以外の従業員は o であるような日数

としたとき、

  •  d(S) :=  \sum_{S \subset T} f(T)

が成立することに注意しよう。よって、 f(S) さえ求めておけば高速ゼータ変換によって  d(S) を求めることができる。 f(S) は、あらかじめ  2NK 日分の「働く従業員の集合」を求めておけば、簡単に求められる。

以上をまとめると、 O(N^{2}K + N2^{N}(\log N + \log K)) で求められることがわかった。

コード

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

const int MAX = 4000000;
int main() {
    int N, K; 
    cin >> N >> K;
    vector<int> A(N);
    for (int i = 0; i < N; ++i) cin >> A[i];
    vector<int> bits(MAX, (1<<N)-1);
    for (int t = 0; t < MAX; ++t) {
        for (int i = 0; i < N; ++i) {
            int r = t % (A[i] * 2);
            if (0 <= r && r < A[i]) bits[t] &= ~(1<<i);
        }
    }
  
    auto biok = [&](int x) {
        vector<int> d(1<<N, 0);
        for (int t = 0; t < x; ++t) d[bits[t]]++;
        for (int i = 0; i < N; ++i) {
            for (int bit = (1<<N)-1; bit >= 0; --bit) {
                if (bit & (1<<i)) continue;
                d[bit] += d[bit | (1<<i)];
            }
        }
        for (int bit = 0; bit < (1<<N); ++bit) {
            int con = 0;
            for (int i = 0; i < N; ++i) if (bit & (1<<i)) ++con;
            int lim = x - con * K;
            if (d[bit] > lim) return false;
        }
        return true;
    };

    int left = 0, right = MAX;
    while (right - left > 1) {
        int mid = (left + right) / 2;
        if (biok(mid)) right = mid;
        else left = mid;
    }
    cout << right << endl;
}

AtCoder ARC 106 F - Figures (橙色, 800 点)

コンテスト本番、こっちをやればよかった...ところで解説が天才すぎる!

問題概要

 N 個の部品と、 N-1 個の接続用部品とがある。これらを用いてフィギュアを作ろうとしている。

 i 番目の部品には  d_{i} 個の穴がついている。接続用部品は、2 個の部品を選んで、それぞれの部品の穴を選択して挿し込むことができるようになっている。

フィギュアは  N 個の部品が連結でなければならない (木になる)。フィギュアの組み立て方が何通りあるか、998244353 で割ったあまりを求めよ。木として同一のものであっても、異なる穴組を接続したものは異なるものとみなす。

制約

  •  2 \le N \le 2 \times 10^{5}
  •  0 \le d_{i} \lt 998244353

考えたこと

まず、次の有名事実がある。


頂点数が  N であるような完全グラフの全域木であって、頂点  v の次数が  d_{v} ( d_{v} の総和が  2N-2) であるものの個数は

 \frac{(N-2)!}{(d_{1}-1)!(d_{2}-1)! \dots (d_{N}-1)!}

で与えられる。


このことの証明は、以下の記事に詳しく書いた。

drken1215.hatenablog.com

これを利用すると、以下を求めればよいことがわかる。

 \sum_{e_{i} \ge 1, e_{1} + \dots + e_{N} = 2N-2} (\frac{(N-2)!}{(e_{1}-1)! \dots (e_{N}-1)!} \times \frac{d_{1}!}{(d_{1} - e_{1})!} \dots \frac{d_{N}!}{(d_{N} - e_{N})!})
 = \sum_{e_{i} \ge 0, e_{1} + \dots + e_{N} = N-2} (\frac{(N-2)!}{e_{1}! \dots e_{N}!}\times \frac{d_{1}!}{(d_{1} - e_{1} - 1)!} \dots \frac{d_{N}!}{(d_{N} - e_{N} - 1)!})
 = (N-2)! d_{1} \dots d_{N} \sum_{e_{i} \ge 0, e_{1} + \dots + e_{N} = N-2} ({}_{d_{1}-1}{\rm C}_{e_{1}} \times \dots \times {}_{d_{N}-1}{\rm C}_{e_{N}})

ここで、 \sum_{e_{i} \ge 0, e_{1} + \dots + e_{N} = N-2} ({}_{d_{1}-1}{\rm C}_{e_{1}} \times \dots \times {}_{d_{N}-1}{\rm C}_{e_{N}}) は、

 f(x) = (1 + x)^{d_{1}-1} \dots (1 + x)^{d_{N}-1} = (1 + x)^{d_{1} + \dots + d_{N} - N}

 e_{1} + \dots + e_{N} ( = N-2) 次の係数になっている。よって、求める値は

 (N-2)! d_{1} \dots d_{N} \lbrack x^{N-2} \rbrack f(x)

と求められる。 f(x) は形式的冪級数の pow を用いることで高速に計算できる。

 

解法 (2):形式的冪級数を使わない

もう少し式変形する。 S = d_{1} + \dots + d_{N} とおく。

 (N-2)! d_{1} \dots d_{N} \lbrack x^{N-2} \rbrack f(x)
 = d_{1} \dots d_{N} \frac{(S-N)!}{(S-2N+2)!}
 = d_{1} \dots d_{N} (S-N)(S-N-1)\dots(S-2N+3)

と求められる。よって  O(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;
    cin >> N;
    mint res = 1;
    long long S = 0;
    vector<long long> d(N);
    for (int i = 0; i < N; ++i) {
        cin >> d[i];
        res *= d[i];
        S += d[i];
    }
    for (int n = N; n <= N*2 - 3; ++n) res *= S - n;
    cout << res << endl;
}

プリューファーコード:全域木の数え上げ (頂点次数制約つき)

ARC 106 F に関連して、頂点次数制約のついた全域木の個数を求める問題がまさにあったので、その解説を。

問題概要 (New Year Contest 2015 E - ひも)

頂点数が  N であるような完全グラフの全域木であって、以下の条件を満たすものが何通りあるか、1000000007 で割ったあまりを求めよ。

  • 各頂点  v における次数が  a_{v} である

制約

  •  2 \le N \le 10^{5}
  •  1 \le a_{i} \le 3

 

解法

次数列の制約のない「全域木の数え上げ」は、Cayley の定理としてよく知られている!!


【Cayley の定理】
頂点数が  N であるような完全グラフの部分グラフであって、全域木をなすものの個数は

 N^{N-2}

で与えられる。


Cayley の定理の証明方法としては、たとえば以下の記事に 6 通りの方法が載っている。

joisino.hatenablog.com

 

次数制約つきの Cayley の定理

次数の指定があっても綺麗になる。


頂点数が  N であるような完全グラフの全域木であって、頂点  v の次数が  d_{v} ( d_{v} の総和が  2N-2) であるものの個数は

 \frac{(N-2)!}{(d_{1}-1)!(d_{2}-1)! \dots (d_{N}-1)!}

で与えられる。


むしろこれを補題として先に示すことで Cayley の定理を証明する方法も考えられる。補題がわかっていれば、Cayley の定理は次のように簡単に導出できる。

 \sum_{d_{i} \ge 1, d_{1} + \dots + d_{N} = 2N-2} \frac{(N-2)!}{(d_{1}-1)! \dots (d_{N}-1)!}
 = \sum_{d_{i} \ge 0, d_{1} + \dots + d_{N} = N-2} \frac{(N-2)!}{d_{1}! \dots d_{N}!}

の値が  N^{N-2} となることを示せばよい。そしてこの各項は、

 (x_{1} + \dots + x_{N})^{N-2}

の次数列が  d_{1}, \dots, d_{N} である項の係数と一致することに注意する。よって、 x_{1} = \dots = x_{N} = 1 とすることで計算できて、

 \sum_{d_{i} \ge 0, d_{1} + \dots + d_{N} = N-2} \frac{(N-2)!}{d_{1}! \dots d_{N}!} = N^{N-2}

となることが示された。

 

補題の証明

補題を示す。数学的帰納法を用いる。また、Cayley の定理の証明方法として有名な「プリューファーコード」による方法を活用する。

ブリューファーコードとは、全域木に対して、「各要素が  1, 2, \dots, N であるような数列  c_{1}, \dots, c_{N-2}」に一対一に対応付けるものである。全域木からプリューファーコードへの対応は次の通り。


  • 全域木の頂点数が 2 になるまで以下の操作を繰り返す (初期状態では数列は空とする)
    • 全域木の「葉」のうち、頂点番号が最小のものに隣接している頂点の番号を数列に末尾に追加する
    • その「葉」を削除する

下図は  N = 7 の木に対して、プリューファーコードを生成する一例を示している。

この「全域木」から「プリューファーコード」への写像は、「逆写像」が自然に作れる (後述) ので、全単射になっていることがわかる。このことから Cayley の定理 (全域木が  N^{N-2} 個) が直接導かれる。

さて、プリューファーコードを真似して、以下の場合に場合分けする

  • プリューファーコードの先頭の要素が頂点 1 のとき
  • プリューファーコードの先頭の要素が頂点 2 のとき
  • ...
  • プリューファーコードの先頭の要素が頂点 N のとき

数学的帰納法の仮定より、プリューファーコードの先頭要素が頂点  1 であるようなものの個数は

 \frac{(N-3)!}{(d_{1}-2)!(d_{2}-1)! \dots (d_{N}-1)!}

となる。一般に、先頭が頂点  v であるものの個数は、上の式の分母の積において、 v 番目のみが  (d_{v}-2)! の形をしていて、それ以外は  (d_{v} - 1)! の形をしているようなものになる。よって、分母を通分して総和をとると、

 \frac{(N-3)!( (d_{1}-1) + \dots + (d_{N}-1) )}{(d_{1}-1)! \dots (d_{N}-1)!}
 = \frac{(N-2)!}{(d_{1}-1)! \dots (d_{N}-1)!}

が導かれる。補題は示された。

 

プリューファーコードから木への対応

最後に、プリューファーコード  c_{1}, \dots, c_{N-2} から、頂点数  N の木を作るのは、次のようにしてできる。ちょうどいい感じの構築問題と言えそう。

  • まず、 c_{1}, c_{2}, \dots, c_{N-1}, c_{N-2} に含まれていない最小の値を  a_{1} として、2 頂点  a_{1}, c_{1} を結ぶ
  • 次に、 a_{1}, c_{2}, \dots, c_{N-1}, c_{N-2} に含まれていない最小の値を  a_{2} として、2 頂点  a_{2}, c_{2} を結ぶ
  • ...
  • 次に、 a_{1}, a_{2}, \dots, a_{N-1}, c_{N-2} に含まれていない最小の値を  a_{N-2} として、2 頂点  a_{N-2}, c_{N-2} を結ぶ
  • 最後に、 a_{1}, a_{2}, \dots, a_{N-2} に含まれていない 2 つの値を  d, e として、2 頂点  d, e を結ぶ

上図のように、プリューファーコードに登場する値が、ちょうど木の「内点」になることがわかる。そして上の手続きは

「指定された内点に対して、まだ使っていない最小番号の頂点を葉にしてくっつける」

という意味を持っていることがわかる。

 

コード

補題が示されれば、それを実装するだけ。ただし

 a_{1} + \dots + a_{N} = 2N-2

とならない場合は 0 通りとすることに注意。

#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];
    }
};

BiCoef<mint> bc;

int main() {
    int N;
    cin >> N;
    bc.init(N);
    long long sum = 0;
    mint res = bc.fact(N-2);
    for (int i = 0; i < N; ++i) {
        int a; cin >> a;
        res *= bc.finv(a-1);
        sum += a-1;
    }
    cout << (sum == N-2 ? res : 0) << endl;
}