Problem:
 Score : 900 points    Problem Statement 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 1 2 1000000007      Sample Output 1 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.        Sample Input 2 3 1000000009      Sample Output 2 118  Let the three kinds of toppings be A, B and C. In addition to the four types of ramen above, four more types of ramen can be ordered, where C is added to the above four. There are 118 sets of ramen that satisfy the conditions, and here are some of them:   The following three ramen: "with A, B", "with A, C", "with B, C".  The following five ramen: "no toppings", "with A", "with A, B", "with B, C", "with A, B, C".  Eight ramen, one for each type.   Note that the set of the following three does not satisfy the condition: "'with A', 'with B', 'with A, B'", because C is not on any of them.       Sample Input 3 50 111111113      Sample Output 3 1456748  Remember to print the number of the sets modulo M . Note that these three sample inputs above are included in the test set for the partial score.       Sample Input 4 3000 123456791      Sample Output 4 16369789  This sample input is not included in the test set for the partial score.

Solution (Python):
import numpy as np


def solve(n, m):
    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

    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

    return ans


N, M = map(int, input().split())
print(solve(N, M))

Solution (Python):
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 = 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)

Solution (Python):
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)

Solution (Python):
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)

Solution (Python):



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


"""
make the table of Sterling numbers of the second kind
Sterling[ball][box]

SIZE = n
Sterling2 = [[0]*(SIZE+1) for _ in range(SIZE+1)]
Sterling2[0][0] = Sterling2[1][1] = 1
for i in range(2,SIZE+1):
    for k in range(1,i+1):
        Sterling2[i][k] = (Sterling2[i-1][k-1] + k*Sterling2[i-1][k])%MOD
"""





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):
    """
    for b in range(n+1):
        res += f[s][b]*pow(2,b*(n-s),MOD)
    """
    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)

Solution (Python):
import sys
input = sys.stdin.readline
import numpy as np

N,MOD = map(int,input().split())

"""
余事象を調べる。包除の原理を使う。
A[n] = （1,2,...,n）が1杯以下、他は何でも良い
B[n,l] : (1,2,...,n) をlグループに分ける方法の個数
A[n]
・0杯のグループあり
・なし
"""

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)

Solution (Python):
import sys
input = sys.stdin.readline
import numpy as np

N,MOD = map(int,input().split())

"""
余事象を調べる。包除の原理を使う。
A[n] = （1,2,...,n）が1杯以下、他は何でも良い
B[n,l] : (1,2,...,n) をlグループに分ける方法の個数
A[n]
・0杯のグループあり
・なし
"""

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)

Solution (Python):
import numpy as np


def solve(n, m):
    def prepare(n, m):
        f = 1
        for i in range(1, n + 1):
            f = f * i % m
        fn = 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 fn, invs

    def nik(n2, k, m):
        ret = 1
        yield ret
        for _ in range(k):
            ret = ret * n2 % m
            yield ret

    fn, invs = prepare(n, m)
    np_range = np.arange(2, n + 3, 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):
        coef = np.fromiter(nik(pow(2, n - i, m), i, m), dtype=np.int64)
        ncr = fn * invs[i] * invs[n - i] % m
        left = (dp[:i + 1] * coef % m).sum() % m
        ans = (ans + ncr * left % m * nn2[n - i] % m * parity) % m
        dp[1:i + 2] = (dp[1:i + 2] * np_range[:i + 1] + dp[:i + 1]) % m
        parity *= -1

    return ans


n, m = map(int, input().split())
print(solve(n, m))

Solution (Python):
import numpy as np


def prepare(n, m):
    f = 1
    for i in range(1, n + 1):
        f = f * i % m
    fn = 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 fn, invs


def nik(n2, k, m):
    ret = 1
    yield ret
    for _ in range(k):
        ret = ret * n2 % m
        yield ret


n, m = map(int, input().split())

fn, invs = prepare(n, m)
np_range = np.arange(2, n + 3, 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
    ans = (ans + ncr * left % m * nn2[n - i] % m * parity) % m
    dp[1:i + 2] = (dp[1:i + 2] * np_range[:i + 1] + dp[:i + 1]) % m
    parity *= -1

print(ans)

Solution (Python):
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 = 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)

Solution (Python):
from collections import defaultdict

class Combinatorics:
    def __init__(self, N, mod):
        '''
        Preprocess for calculating binomial coefficients nCr (0 <= r <= n, 0 <= n <= N)
        over the finite field Z/(mod)Z.
        Input:
            N (int): maximum n
            mod (int): a prime number. The order of the field Z/(mod)Z over which nCr is calculated.
        '''
        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):
        '''
        Calculate nCr = n! /(r! (n-r)!) % mod
        '''
        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)

Solution (Python):
N, M = 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)

Solution (Python):
N, M = 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)

Solution (Python):
N, M = 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)