AtCoder Regular Contest 096 - Everything on It

In "Takahashi-ya", a ramen restaurant, basically they have one menu: "ramen", but N kinds of toppings are also offered. When a customer orders a bowl of ramen, for each kind of topping, he/she can choose whether to put it on top of his/her ramen or not. There is no limit on the number of toppings, and it is allowed to have all kinds of toppings or no topping at all. That is, considering the combination of the toppings, 2^N types of ramen can be ordered.
Akaki entered Takahashi-ya. She is thinking of ordering some bowls of ramen that satisfy both of the following two conditions:
 - Do not order multiple bowls of ramen with the exactly same set of toppings.
 - Each of the N kinds of toppings is on two or more bowls of ramen ordered.
You are given N and a prime number M. Find the number of the sets of bowls of ramen that satisfy these conditions, disregarding order, modulo M. Since she is in extreme hunger, ordering any number of bowls of ramen is fine.

-----Constraints-----
 - 2 \leq N \leq 3000
 - 10^8 \leq M \leq 10^9 + 9
 - N is an integer.
 - M is a prime number.

-----Subscores-----
 - 600 points will be awarded for passing the test set satisfying N ≤ 50.

-----Input-----
Input is given from Standard Input in the following format:
N M

-----Output-----
Print the number of the sets of bowls of ramen that satisfy the conditions, disregarding order, modulo M.

-----Sample Input-----
2 1000000007

-----Sample Output-----
2

Let the two kinds of toppings be A and B. Four types of ramen can be ordered: "no toppings", "with A", "with B" and "with A, B". There are two sets of ramen that satisfy the conditions:
 - The following three ramen: "with A", "with B", "with A, B".
 - Four ramen, one for each type.

Source: atcoder | Difficulty: VERY_HARD


Solutions:
(N, M) = list(map(int, input().split()))
fact = [1] * (N + 1)
rfact = [1] * (N + 1)
for i in range(1, N + 1):
	fact[i] = r = i * fact[i - 1] % M
	rfact[i] = pow(r, M - 2, M)
memo = {}

def f(N, K):
	if (N, K) in memo:
		return memo[N, K]
	if K == 0:
		return N == 0
	elif K == 1 or N == K:
		return 1
	r = memo[N, K] = (f(N - 1, K - 1) + K * f(N - 1, K)) % M
	return r
S = [1]
rev2 = pow(2, M - 2, M)
base = pow(2, N, M)
ans = 0
S = [1]
for K in range(N + 1):
	r = fact[N] * rfact[K] * rfact[N - K] % M
	r = r * pow(2, pow(2, N - K, M - 1), M) % M
	b = 1
	v = 0
	T = [0] * (K + 2)
	for L in range(K):
		T[L + 1] = s = (S[L] + (L + 1) * S[L + 1]) % M
		v += s * b
		b = b * base % M
	v += b
	T[K + 1] = 1
	S = T
	r = r * v % M
	if K % 2:
		ans -= r
	else:
		ans += r
	ans %= M
	base = base * rev2 % M
print(ans)


import numpy as np

def prepare(n, m):
	f = 1
	facts = [1]
	for i in range(1, n + 1):
		f = f * i % m
		facts.append(f)
	f = pow(f, m - 2, m)
	invs = [1] * (n + 1)
	invs[n] = f
	for i in range(n, 0, -1):
		f = f * i % m
		invs[i - 1] = f
	return (facts, invs)

def nik(n2, k, m):
	ret = 1
	yield ret
	for _ in range(k):
		ret = ret * n2 % m
		yield ret
(n, m) = list(map(int, input().split()))
(facts, invs) = prepare(n, m)
fn = facts[n]
np_range = np.arange(2, n + 4, dtype=np.int64)
dp = np.zeros(n + 2, dtype=np.int64)
dp[0] = 1
nn2 = [2]
for i in range(n):
	nn2.append(nn2[-1] ** 2 % m)
ans = 0
parity = 1
for i in range(n + 1):
	n2 = pow(2, n - i, m)
	coef = np.fromiter(nik(n2, i, m), dtype=np.int64)
	ncr = fn * invs[i] * invs[n - i] % m
	left = (dp[:i + 1] * coef % m).sum() % m
	right = nn2[n - i]
	ans = (ans + ncr * left % m * right % m * parity) % m
	dp[1:i + 2] = (dp[1:i + 2] * np_range[:i + 1] + dp[:i + 1]) % m
	parity *= -1
print(ans)


(n, m) = map(int, input().split())
MOD = m
list_size = n + 2
f_list = [1] * list_size
f_r_list = [1] * list_size
for i in range(list_size - 1):
	f_list[i + 1] = int(f_list[i] * (i + 1) % MOD)
f_r_list = [pow(x, MOD - 2, MOD) for x in f_list]

def comb(N, r):
	if N < r or r < 0:
		return 0
	else:
		return f_list[N] * f_r_list[N - r] % MOD * f_r_list[r] % MOD
st = [[0 for _ in range(n + 1)] for _ in range(n + 1)]
st[0][0] = 1
for i in range(1, n + 1):
	st[i][0] = 1
	for j in range(1, i + 1):
		st[i][j] = (st[i - 1][j - 1] + (j + 1) * st[i - 1][j]) % MOD
ans = 0
for i in range(n + 1):
	res = comb(n, i) * pow(2, pow(2, n - i, MOD - 1), MOD) % MOD
	tmp = 0
	for j in range(i + 1):
		tmp += st[i][j] * pow(2, (n - i) * j, MOD)
		tmp %= MOD
	res *= tmp
	res %= MOD
	if i % 2 == 0:
		ans += res
	else:
		ans -= res
	ans %= MOD
print(ans)


import numpy as np

def prepare(n, m):
	f = 1
	for i in range(1, n + 1):
		f = f * i % m
	fn = f
	inv = [1] * (n + 1)
	f = pow(f, m - 2, m)
	inv[n] = f
	for i in range(n, 0, -1):
		f = f * i % m
		inv[i - 1] = f
	return (fn, inv)

def a_x(a, x, m):
	ret = 1
	yield ret
	for _ in range(x):
		ret = ret * a % m
		yield ret
(N, M) = list(map(int, input().split()))
(fn, inv) = prepare(N, M)
stir2 = np.zeros(N + 2, dtype=np.int64)
stir2[0] = 1
upd = np.arange(2, N + 3, dtype=np.int64)
ex2 = [2]
for i in range(N):
	ex2.append(ex2[-1] ** 2 % M)
ans = 0
si = 1
for i in range(N + 1):
	nCi = fn * inv[i] * inv[N - i] % M
	i_with = np.fromiter(a_x(pow(2, N - i, M), i, M), dtype=np.int64)
	i_on = (stir2[:i + 1] * i_with % M).sum() % M
	ans = (ans + nCi * i_on % M * ex2[N - i] % M * si) % M
	stir2[1:i + 2] = (stir2[1:i + 2] * upd[:i + 1] + stir2[:i + 1]) % M
	si *= -1
print(ans)


from collections import defaultdict

class Combinatorics:

	def __init__(self, N, mod):
		self.mod = mod
		self.fact = {i: None for i in range(N + 1)}
		self.inverse = {i: None for i in range(1, N + 1)}
		self.fact_inverse = {i: None for i in range(N + 1)}
		self.fact[0] = self.fact[1] = 1
		self.fact_inverse[0] = self.fact_inverse[1] = 1
		self.inverse[1] = 1
		for i in range(2, N + 1):
			self.fact[i] = i * self.fact[i - 1] % self.mod
			(q, r) = divmod(self.mod, i)
			self.inverse[i] = -(q % self.mod) * self.inverse[r] % self.mod
			self.fact_inverse[i] = self.inverse[i] * self.fact_inverse[i - 1] % self.mod

	def binom(self, n, r):
		if n < r or n < 0 or r < 0:
			return 0
		else:
			return self.fact[n] * (self.fact_inverse[r] * self.fact_inverse[n - r] % self.mod) % self.mod
(N, M) = map(int, input().split())
com = Combinatorics(N, M)
ans = 0
pow2 = [0] * (N * N // 4 + 1)
pow_pow2 = [0] * (N + 1)
pow2[0] = 1
pow_pow2[0] = 2
for i in range(1, N * N // 4 + 1):
	pow2[i] = pow2[i - 1] * 2 % M
for i in range(1, N + 1):
	pow_pow2[i] = pow(pow_pow2[i - 1], 2, M)
ways2 = {n: defaultdict(int) for n in range(N + 1)}
for n in range(N + 1):
	ways = 0
	temp = 0
	if n >= 1:
		for i in range(n + 1):
			ways2[n][i] = (ways2[n][i] + ways2[n - 1][i]) % M
			if i >= 1:
				ways2[n][i] = (ways2[n][i] + ways2[n - 1][i - 1]) % M
			ways2[n][i] = (ways2[n][i] + i * ways2[n - 1][i] % M) % M
			temp = (temp + ways2[n][i] * pow2[(N - n) * i] % M) % M
	else:
		ways2[n][0] = 1
		temp += ways2[n][0]
	ways = temp * pow_pow2[N - n] % M
	if n % 2 == 0:
		ans = (ans + com.binom(N, n) * ways % M) % M
	else:
		ans = (ans - com.binom(N, n) * ways % M) % M
print(ans)


(N, M) = list(map(int, input().split()))
fact = [1] * (N + 1)
rfact = [1] * (N + 1)
for i in range(1, N + 1):
	fact[i] = r = i * fact[i - 1] % M
	rfact[i] = pow(r, M - 2, M)
S = [1]
rev2 = pow(2, M - 2, M)
base = pow(2, N, M)
ans = 0
S = [1]
for K in range(N + 1):
	res = fact[N] * rfact[K] * rfact[N - K] % M
	res = res * pow(2, pow(2, N - K, M - 1), M) % M
	b = 1
	v = 0
	T = [0] * (K + 2)
	for L in range(K):
		T[L + 1] = s = (S[L] + (L + 1) * S[L + 1]) % M
		v += s * b
		b = b * base % M
	v += b
	T[K + 1] = 1
	S = T
	res = res * v % M
	if K % 2:
		ans -= res
	else:
		ans += res
	ans %= M
	base = base * rev2 % M
print(ans)


import sys
input = sys.stdin.readline
import numpy as np
(N, MOD) = map(int, input().split())
B = np.zeros((N + 1, N + 1), dtype=np.int64)
B[0, 0] = 1
for n in range(1, N + 1):
	B[n, 1:] = B[n - 1, :-1]
	B[n, 1:] += B[n - 1, 1:] * np.arange(1, N + 1) % MOD
	B[n] %= MOD
pow_2 = np.ones((N + 1, N + 1), dtype=np.int64)
for n in range(1, N + 1):
	pow_2[1, n] = 2 * pow_2[1, n - 1] % MOD
for n in range(2, N + 1):
	pow_2[n] = pow_2[n - 1] * pow_2[1] % MOD
pow_pow = np.zeros(N + 1, dtype=np.int64)
pow_pow[0] = 2
for n in range(1, N + 1):
	pow_pow[n] = pow_pow[n - 1] ** 2 % MOD
A = (B[:, 1:] * (pow_2[::-1, 1:] + pow_2[::-1, :-1] * np.arange(1, N + 1) % MOD) % MOD).sum(axis=1)
A %= MOD
A *= pow_pow[::-1]
A %= MOD
fact = [1] * (N + 1)
fact_inv = [1] * (N + 1)
for n in range(1, N + 1):
	fact[n] = fact[n - 1] * n % MOD
fact_inv[N] = pow(fact[N], MOD - 2, MOD)
for n in range(N, 0, -1):
	fact_inv[n - 1] = fact_inv[n] * n % MOD
fact = np.array(fact, dtype=np.int64)
fact_inv = np.array(fact_inv, dtype=np.int64)
comb = fact_inv * fact_inv[::-1] % MOD * fact[N] % MOD
A[::2] *= -1
A *= comb
A %= MOD
answer = pow(2, pow(2, N, MOD - 1), MOD) - A.sum()
answer %= MOD
print(answer)


import sys
input = sys.stdin.readline
import numpy as np
(N, MOD) = list(map(int, input().split()))
B = np.zeros((N + 1, N + 1), dtype=np.int64)
B[0, 0] = 1
for n in range(1, N + 1):
	B[n, 1:] = B[n - 1, :-1]
	B[n, 1:] += B[n - 1, 1:] * np.arange(1, N + 1) % MOD
	B[n] %= MOD
pow_2 = np.ones((N + 1, N + 1), dtype=np.int64)
for n in range(1, N + 1):
	pow_2[1, n] = 2 * pow_2[1, n - 1] % MOD
for n in range(2, N + 1):
	pow_2[n] = pow_2[n - 1] * pow_2[1] % MOD
A = np.zeros(N + 1, dtype=np.int64)
for n in range(N + 1):
	A[n] = (pow(2, pow(2, N - n, MOD - 1), MOD) * B[n, 1:] % MOD * (pow_2[N - n, 1:] + pow_2[N - n, :-1] * np.arange(1, N + 1) % MOD) % MOD).sum() % MOD
comb = np.zeros((N + 1, N + 1), dtype=np.int64)
comb[:, 0] = 1
for n in range(1, N + 1):
	comb[n, 1:] = (comb[n - 1, 1:] + comb[n - 1, :-1]) % MOD
A[::2] *= -1
A *= comb[N]
A %= MOD
answer = pow(2, pow(2, N, MOD - 1), MOD) - A.sum()
answer %= MOD
print(answer)


import sys
sys.setrecursionlimit(10 ** 6)
readline = sys.stdin.readline
(n, MOD) = [int(i) for i in readline().split()]
SIZE = 3001
SIZE += 1
inv = [0] * SIZE
fac = [0] * SIZE
finv = [0] * SIZE
inv[1] = 1
fac[0] = fac[1] = 1
finv[0] = finv[1] = 1
for i in range(2, SIZE):
	inv[i] = MOD - MOD // i * inv[MOD % i] % MOD
	fac[i] = fac[i - 1] * i % MOD
	finv[i] = finv[i - 1] * inv[i] % MOD

def choose(n, r):
	if 0 <= r <= n:
		return fac[n] * finv[r] % MOD * finv[n - r] % MOD
	else:
		return 0
SIZE = n + 1
f = [[0] * SIZE for _ in range(SIZE)]
f[0][0] = 1
for i in range(1, SIZE):
	f[i][0] = 1
	for k in range(1, i + 1):
		f[i][k] = (f[i - 1][k] + f[i - 1][k - 1] + k * f[i - 1][k]) % MOD
pp = [1] * SIZE
for i in range(1, SIZE):
	pp[i] = 2 * pp[i - 1] % (MOD - 1)

def size(s):
	res = 0
	p = 1
	e = pow(2, n - s, MOD)
	for b in range(SIZE):
		res += f[s][b] * p
		p *= e
		p %= MOD
	return res * pow(2, pp[n - s], MOD) % MOD
ans = 0
sgn = 1
for j in range(n + 1):
	ans += sgn * choose(n, j) * size(j)
	ans %= MOD
	sgn *= -1
print(ans)