ARC096 E - Everything on It

AtCoder Regular Contest

春に解説pdfを見ても全く意味がわからなくてずっと放置してて、さっき解説放送見たらやっと理解できた
700点とかの問題と比べると恐ろしいほどむずいけど勉強になる問題だった

問題

あるラーメン屋はn種類のトッピングを提供している
各トッピングは乗せるか乗せないかを選べて、2^n通りのラーメンが注文できる
ここで、
1) トッピングの組み合わせが全く同じラーメンを複数注文しない
2) n種類のトッピングそれぞれが注文したラーメンのうち2杯以上に乗っている
という2つの条件を満たすようにラーメンを何杯か注文したい
そのような注文の仕方の総数を求めよ

2 <= n <= 3000

解法

求めたいものは
「"n個の要素からなる集合S" の部分集合2^n個からなる集合Tの部分集合Uのうち、Sの全ての要素がUの2個以上の互いに異なる要素に属するようなものの個数」
と言い換えられる
(Uは 2^(2^n) 個存在する)
2個以上のグループに属するというのはあまりにも求めづらいので、
f(k): Sのうち1個以下のグループに属するような要素がk個以上存在するようなグループ分けの仕方の個数 (残りn-k個は何個のグループに属してもいい)
とし、包除原理を用いると、答えは
Σ nCk * f(k) * (-1)^k (0 <= k <= n)
となる

次にf(k)を求めたい
n-k個は何個のグループに属しててもいいから、グループへの入れ方は 2^(2^(n-k)) 通りある
グループがx個あるとして、
dp[k][x]: k個の要素を(各要素が1個以下のグループに属するように)x個のグループに分ける場合の数
とすると、x個のグループそれぞれがn-k個の要素のうちどれを含むかというのが 2^(n-k)x 通りあるから
f(k) = Σ dp[k][x] * 2^(n-k)x * 2^(2^(n-k)) (0 <= x <= k)
となる

今度はdp[k][x]を求めたい
これは第二種スターリング数ってのとめっちゃ似てて、↓の遷移で求められる
dp[k][x] = dp[k-1][x-1] + dp[k-1][x] + dp[k-1][x]*x
第1項) 新しいグループを作ってそこにk個目の要素を入れる場合
第2項) k個目の要素をどのグループにも入れない場合
第3項) 既にあるx個のグループのどれか1つにk個目の要素を入れる場合

全部トータルでO(n^2)

実装

フェルマーの小定理 2^(mod-1) % mod = 1 より、
2 ^ (2^(n-k))
= 2 ^ (mod-1) * 2 ^ (mod-1) * ... * 2 ^ (2^(n-k) % (mod-1))
= 2 ^ (2^(n-k) % (mod-1)) (mod mod)
と変形できる

#define int long long

int n, mod, dp[3030][3030], p[10010010], res;

ll mop(ll a,ll b,ll m=mod) {ll r=1;a%=m;for(;b;b>>=1){if(b&1)r=r*a%m;a=a*a%m;}return r;}

const int MAXS = 10010;
ll fact[MAXS+1], factr[MAXS+1], inv[MAXS+1];
ll comb(ll n, ll r) {
	if (fact[0]==0) {
		fact[0] = factr[0] = inv[1] = 1;
		for (int i=2; i<=MAXS; i++) inv[i] = inv[mod%i] * (mod-mod/i) %mod;
		for (int i=1; i<=MAXS; i++) fact[i] = fact[i-1]*i %mod, factr[i] = factr[i-1]*inv[i] %mod;
	}
	if (r<0 || n<r) return 0;
	return fact[n]*factr[r] %mod *factr[n-r] %mod;
}

int f(int k) {
	int ret = 0;
	rep(x,k+1) (ret += dp[k][x] * p[(n-k)*x] %mod) %= mod;
	return ret * mop(2,mop(2,n-k,mod-1)) %mod;
}

signed main() {
	cin >> n >> mod;
	dp[0][0] = 1;
	reps(i,1,n+1) rep(j,n+1) {
			if (j) (dp[i][j] += dp[i-1][j-1]) %= mod;
			(dp[i][j] += dp[i-1][j]*(j+1)) %= mod;
	}
	p[0] = 1;
	reps(i,1,10010010) p[i] = p[i-1]*2 %mod;
	rep(k,n+1) (res += comb(n,k) * f(k) * (k%2 ? -1 : 1) %mod + mod) %= mod;
	cout << res << endl;
}