これ 81 人も通してるのか...
問題概要
人の子供がいる。
日間にわたって、 人の中から 人を一様ランダムに選んでクッキーを与える。
日分のあらゆるクッキーの配り方を考えたときの「各子供の最終的にもらったクッキーの個数の積」の総和を求めよ。
制約
解法 (1): 公式解法、積和の言い換え
「積」の総和を求める問題って、やれることがとても限られてるイメージはある。手が出せずに困っていた。とりあえず、解説の方針を覚えておこうと思う。積の総和は、各子供へのクッキーの配り方それぞれについて、「その子供のクッキーをもらった日のうちの 1 日」を選ぶ、という方法の総和をとったものととらえることができる。そうすると元の問題は以下のようになる。
- のマス目を「無地」「白」「黒」に塗る
- 行 の「白」と「黒」が合計で 個
- どの列を見ても黒の個数が 1 個
という方法を数え上げよ
- dp[ i ][ j ] := i -1 行目までを決めて、すでに黒く塗られたマスが j 個分あるような場合の数
とすると、i 行目で新たに k 個の黒マスを増やすことにするとき、
- 黒く塗られたことのない列 (N - j 個) から k 個を選ぶ: C(N - j, k) 通り
- 残りの列から a[ i ] - k 個を選ぶ: C(N - k, a[ i ] - k) 通り
ということで、
- dp[ i + 1 ][ j + k ] += dp[ i ][ j ] × C(N - j, k) × C(N - k, a[ i ] - k)
となる。計算量は 。
#include <iostream> #include <vector> using namespace std; // modint: mod 計算を int を扱うように扱える構造体 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() { 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 ostream& operator << (ostream &os, const Fp<MOD>& x) noexcept { return os << x.val; } friend constexpr Fp<MOD> modpow(const Fp<MOD> &a, long long n) noexcept { if (n == 0) return 1; auto t = modpow(a, n / 2); t = t * t; if (n & 1) t = t * a; return t; } }; // 二項係数ライブラリ 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]; } }; const int MOD = 1000000007; using mint = Fp<MOD>; BiCoef<mint> bc; int main() { int N, K; cin >> N >> K; vector<int> a(K); for (int i = 0; i < K; ++i) cin >> a[i]; bc.init(N+1); vector<vector<mint> > dp(K+1, vector<mint>(N+1, 0)); dp[0][0] = 1; for (int i = 0; i < K; ++i) { for (int j = 0; j <= N; ++j) { for (int k = 0; k <= a[i]; ++k) { dp[i+1][j+k] += dp[i][j] * bc.com(N-j, k) * bc.com(N-k, a[i]-k); } } } cout << dp[K][N] << endl; }
解法 (2): 確率変数の積からの発想
maspy さんが多項式で殴ったと聞いて、僕も多項式で殴れるようにしたいと思った。
そもそも問題は の行列を 0 か 1 で塗る方法を考えるような整理の仕方をするとよさそう。「行 には 個の 1 がある」という条件を満たすように塗るわけだ。さて
- マス を 1 にするとき であり、0 にするとき となるような確率変数
を考える。サンプル 1 について考えてみる。, なので、求める期待値は
= + + + + + + +
となる。全部で 個の項にわかれるわけだ。しかしよくよく考えてみると、対称性からまとめられるところがある。これらのうち
- 行目から 3 個、 行目から 0 個選ぶもの: 1 項 (確率は )
- 行目から 2 個、 行目から 1 個選ぶもの: 3 項 (確率は )
- 行目から 1 個、 行目から 3 個選ぶもの: 3 項 (確率は )
- 行目から 0 個、 行目から 3 個選ぶもの: 3 項 (確率は )
という感じになっていて、求める期待値は と求められる。一般には、確率変数として、それぞれの行から 個ずつ選んであげて、その総和が になるようにする。
というわけで dp を書いてみる。
- dp[ ][ ] := 行分までを見て確率変数に選ばれたものが 個である場合の期待値
とする。次の行で新たに 個分を選ぶとき、
- まだ選ばれてない 行から 列分を選ぶ: 通り
- それらがすべて 1 になる確率:
なので、
dp[ ][ ] = dp[ ][ ]
となる。実はさっきと同じ DP になっている。
さらに:FFT
DP の式を式変形する。
- dp2[ ][ ] = dp1[ [ ]
とおくと、
- dp2[ ][ ] += dp2[ ][ ]
となる。 を外だしすると、
の 次の項が答えになると思われる。FFT を用いることで となり、もっというと繰り返し二乗法を用いることで になると思う。