Problem:


Score : 800 points

Problem Statement

There are N children, numbered 1,2,\ldots,N.
In the next K days, we will give them some cookies.
In the i-th day, we will choose a_i children among the N with equal probability, and give one cookie to each child chosen. (We make these K choices independently.)

Let us define the happiness of the children as c_1 \times c_2 \times \ldots \times c_N, where c_i is the number of cookies received by Child i in the K days.
Find the expected happiness multiplied by \binom{N}{a_1} \times \binom{N}{a_2} \times \ldots \times \binom{N}{a_K} (we can show that this value is an integer), modulo (10^{9}+7).

Notes

\binom{n}{k} denotes the number of possible choices of k objects out of given distinct n objects, disregarding order.

Constraints

1 \leq N \leq 1000

1 \leq K \leq 20

1 \leq a_i \leq N

Input

Input is given from Standard Input in the following format:

N K
a_1 a_2 \ldots a_K

Output

Print the answer.

Sample Input 1

3 2
3 2

Sample Output 1

12

On the first day, all the children 1, 2, and 3 receive one cookie.

On the second day, two out of the three children 1, 2, and 3 receive one cookie.

Their happiness is 4 in any case, so the expected happiness is 4. Print this value multiplied by \binom{3}{3} \times \binom{3}{2}, which is 12.

Sample Input 2

856 16
399 263 665 432 206 61 784 548 422 313 848 478 827 26 398 63

Sample Output 2

337587117

Find the expected value multiplied by \binom{N}{a_1} \times \binom{N}{a_2} \times \ldots \times \binom{N}{a_K}, modulo (10^9+7).



Sample input:
3 2
3 2


Sample output:
12

Submission status: Accepted

Solution code:
def inv(a):
    m=10**9+7
    b=m
    (x, lastx) = (0, 1)
    (y, lasty) = (1, 0)
    while b != 0:
        q = a // b
        (a, b) = (b, a % b)
        (x, lastx) = (lastx - q * x, x)
        (y, lasty) = (lasty - q * y, y)
    return lastx % m


def main():
    MOD=10**9+7
    n,k=map(int,input().split())
    a=list(map(int,input().split()))

    frac=[1]
    for i in range(1,n+1):
        frac.append(frac[-1]*i%MOD)
    
    fracinv=[inv(frac[i]) for i in range(n+1)]


    dp=[[0]*(n+1) for _ in range(k)]

    for x in range(a[0]+1):
        dp[0][x]=frac[n]*frac[n-x]*inv(frac[a[0]-x]*frac[n-a[0]]*frac[x])%MOD

    for i in range(k-1):
        ai=a[i+1]
        for j in range(n+1):
            for x in range(ai+1):
                p=j+x
                if p>n: break
                tmp=dp[i][j]
                tmp=tmp*frac[n-x]%MOD
                tmp=tmp*fracinv[ai-x]%MOD
                tmp=tmp*fracinv[n-ai]%MOD
                tmp=tmp*fracinv[x]%MOD
                dp[i+1][p]=(dp[i+1][p]+tmp)%MOD
    print(dp[k-1][n]%MOD)

if __name__ == "__main__":
    main()

Submission status: Accepted

Solution code:
def inv(a):
    m=10**9+7
    b=m
    (x, lastx) = (0, 1)
    (y, lasty) = (1, 0)
    while b != 0:
        q = a // b
        (a, b) = (b, a % b)
        (x, lastx) = (lastx - q * x, x)
        (y, lasty) = (lasty - q * y, y)
    return lastx % m


def main():
    MOD=10**9+7
    n,k=map(int,input().split())
    a=list(map(int,input().split()))

    frac=[1]
    for i in range(1,n+1):
        frac.append(frac[-1]*i%MOD)
    
    fracinv=[inv(frac[i]) for i in range(n+1)]


    dp=[[0]*(n+1) for _ in range(k)]

    for x in range(a[0]+1):
        dp[0][x]=frac[n]*frac[n-x]*inv(frac[a[0]-x]*frac[n-a[0]]*frac[x])%MOD

    for i in range(k-1):
        ai=a[i+1]
        for j in range(n+1):
            for x in range(ai+1):
                p=j+x
                if p>n: break
                tmp=dp[i][j]
                tmp=(tmp*frac[n-x]*fracinv[ai-x])%MOD
                tmp=(tmp*fracinv[n-ai]*fracinv[x])%MOD
                dp[i+1][p]=(dp[i+1][p]+tmp)%MOD
    print(dp[k-1][n]%MOD)

if __name__ == "__main__":
    main()

Submission status: Accepted

Solution code:
import sys
import numpy as np


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


n, k, *aaa = map(int, sys.stdin.buffer.read().split())
MOD = 10 ** 9 + 7
facts, invs = prepare(n, MOD)

dp = np.zeros(n + 1, dtype=np.int64)
dp[0] = facts[n]
for a in aaa:
    pat0 = facts[n] * invs[a] % MOD * invs[n - a] % MOD
    ndp = dp * pat0 % MOD
    for x in range(1, a + 1):
        pat = facts[n - x] * invs[a - x] % MOD * invs[n - a] % MOD * invs[x] % MOD
        ndp[x:] = (ndp[x:] + dp[:-x] * pat) % MOD
    dp = ndp
print(dp[n])

Submission status: Accepted

Solution code:
import sys
readline = sys.stdin.readline

MOD = 10**9+7
def frac(limit):
    frac = [1]*limit
    for i in range(2,limit):
        frac[i] = i * frac[i-1]%MOD
    fraci = [None]*limit
    fraci[-1] = pow(frac[-1], MOD -2, MOD)
    for i in range(-2, -limit-1, -1):
        fraci[i] = fraci[i+1] * (limit + i + 1) % MOD
    return frac, fraci
frac, fraci = frac(1398)
def comb(a, b):
    if not a >= b >= 0:
        return 0
    return frac[a]*fraci[b]*fraci[a-b]%MOD


N, K = map(int, readline().split())
A = list(map(int, readline().split()))

dp = [0]*(N+1)
dp[0] = 1

for i in range(K):
    a = A[i]
    dp2 = [0]*(N+1)
    for j in range(N+1):
        d = dp[j]
        if not d:
            continue
        for k in range(min(a, N-j)+1):
            dp2[j+k] = (dp2[j+k] + d*comb(N-j, k)*comb(N-k, a-k)) %MOD
    dp = dp2[:]
print(dp[N])

Submission status: Accepted

Solution code:
import sys
input = sys.stdin.readline

N,K=map(int,input().split())
A=list(map(int,input().split()))
mod=10**9+7

Combi=[[] for i in range(N+1)]
Combi[0]=[1,0]

for i in range(1,N+1):
    Combi[i].append(1)
    for j in range(i):
        Combi[i].append((Combi[i-1][j]+Combi[i-1][j+1])%mod)
    Combi[i].append(0)

DP=[0]*(N+1)
DP[0]=1

for k in range(K):
    use=A[k]

    NDP=[0]*(N+1)

    for i in range(N+1):
        ANS=0
        for j in range(i+1):
            if use-(i-j)>=0:
                ANS=(ANS+DP[j]*Combi[N-j][i-j]*Combi[N-(i-j)][use-(i-j)])%mod

        NDP[i]=ANS

    #print(DP)

    DP=NDP

print(DP[N])

Submission status: Accepted

Solution code:
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

import numpy as np

N,K,*A = map(int,read().split())

MOD = 10 ** 9 + 7

class ModArray(np.ndarray):
    def __new__(cls, input_array, info=None):
        obj = np.asarray(input_array, dtype=np.int64).view(cls)
        return obj
    
    @classmethod
    def zeros(cls, *args):
        obj = np.zeros(*args, dtype=np.int64)
        return obj.view(cls)
    
    @classmethod
    def ones(cls, *args):
        obj = np.ones(*args, dtype=np.int64)
        return obj.view(cls)
    
    @classmethod
    def arange(cls, *args):
        obj = np.arange(*args, dtype=np.int64) % MOD
        return obj.view(cls)
    
    @classmethod
    def powers(cls, a, N):
        B = N.bit_length()
        x = cls.ones(1<<B)
        for n in range(B):
            x[1<<n:1<<(n+1)] = x[:1<<n] * a
            a *= a; a %= MOD
        return x[:N]
    
    @classmethod
    def _set_constant_array(cls, N):
        x = cls.arange(N); x[0] = 1
        fact = x.cumprod()
        x = cls.arange(N,0,-1); x[0] = pow(int(fact[-1]), MOD-2, MOD)
        fact_inv = x.cumprod()[::-1]
        inv = cls.zeros(N)
        inv[1:N] = fact_inv[1:N] * fact[0:N-1]
        fact.flags.writeable = False
        fact_inv.flags.writeable = False
        inv.flags.writeable = False
        cls._fact = fact; cls._fact_inv = fact_inv; cls._inverse = inv
        
    @classmethod
    def fact(cls, N):
        if (getattr(cls, '_fact', None) is None) or len(cls._fact) < N:
            cls._set_constant_array(max(N, 10 ** 6))
        return cls._fact[:N]

    @classmethod
    def fact_inv(cls, N):
        if (getattr(cls, '_fact', None) is None) or len(cls._fact) < N:
            cls._set_constant_array(max(N, 10 ** 6))
        return cls._fact_inv[:N]

    @classmethod
    def inverse(cls, N):
        if (getattr(cls, '_fact', None) is None) or len(cls._fact) < N:
            cls._set_constant_array(max(N, 10 ** 6))
        return cls._inverse[:N]
    
    @classmethod
    def _convolve_small(cls, f, g):
        Lf = len(f); Lg = len(g); L = Lf + Lg - 1
        if min(Lf, Lg) < 100 or Lf + Lg < 300:
            return (np.convolve(f,g) % MOD).view(cls)
        else:
            fft = np.fft.rfft; ifft = np.fft.irfft; n = 1 << L.bit_length()
            return (np.rint(ifft(fft(f, n) * fft(g, n))[:L]).astype(np.int64) % MOD).view(cls)
    
    @classmethod
    def convolve(cls, f, g):
        if f.max() < (1<<15) and g.max() < (1<<15):
            return cls._convolve_small(f, g)
        Lf = len(f); Lg = len(g); L = Lf + Lg - 1
        fl = f & (1 << 15) - 1; fh = f >> 15
        gl = g & (1 << 15) - 1; gh = g >> 15
        x = cls._convolve_small(fl, gl)
        z = cls._convolve_small(fh, gh)
        y = cls._convolve_small(fl+fh, gl+gh) - x - z
        return x + (y<<15) + (z<<30)
    
    def resize(self, N):
        L = len(self)
        if L >= N:
            return self[:N]
        A = np.resize(self, N).view(self.__class__)
        A[L:] = 0
        return A
    
    def diff(self):
        return np.diff(self) % MOD
               
    def cumprod(self):
        L = len(self); Lsq = int(L**.5+1)
        A = self.resize(Lsq**2).reshape(Lsq, Lsq)
        for n in range(1,Lsq):
            A[:,n] *= A[:,n-1]
        for n in range(1,Lsq):
            A[n] *= A[n-1,-1]
        return A.ravel()[:L]

    def __array_wrap__(self, out_arr, context=None):
        if context and context[0] == np.mod:
             return out_arr
        if out_arr.dtype == np.int64:
            np.mod(out_arr, MOD, out=out_arr)
            return out_arr
        else:
            return np.asarray(out_arr)

U = 10 ** 5
fact = ModArray.fact(U)
fact_inv = ModArray.fact_inv(U)

f = np.zeros(N+1,np.int64)
f[0] = 1
for a in A:
    g = fact[N-a:N+1] * fact_inv[0:a+1] * fact_inv[N-a]
    g = g[::-1]
    g *= fact_inv[:len(g)]
    f = ModArray.convolve(f,g)[:N+1]

f *= fact[:N+1]

answer = f[N]
print(answer)

Submission status: Accepted

Solution code:
N,K = map(int,input().split())
A = list(map(int,input().split()))
s = sum(A)
if s < N:
    print(0)
    exit()

MOD = 10**9+7

MAXN = N+5
fac = [1,1] + [0]*MAXN
finv = [1,1] + [0]*MAXN
inv = [0,1] + [0]*MAXN
for i in range(2,MAXN+2):
    fac[i] = fac[i-1] * i % MOD
    inv[i] = -inv[MOD%i] * (MOD // i) % MOD
    finv[i] = finv[i-1] * inv[i] % MOD

def ncr(n,r):
    if n < r: return 0
    if n < 0 or r < 0: return 0
    return fac[n] * (finv[r] * finv[n-r] % MOD) % MOD

dp = [[0]*(N+1) for i in range(K+1)]
dp[0][0] = 1
for i,a in enumerate(A):
    for j in range(N+1):
        for k in range(j,N+1):
            if k-j > a: break
            dp[i+1][k] += dp[i][j] * ncr(N-j, k-j) * ncr(N-(k-j), a-(k-j))
            dp[i+1][k] %= MOD
print(dp[-1][-1])

Submission status: Accepted

Solution code:
import sys

sys.setrecursionlimit(10 ** 6)
def MI(): return map(int, sys.stdin.readline().split())
def LI(): return list(map(int, sys.stdin.readline().split()))

def main():
    def com(com_n, com_r):
        return fac[com_n] * inv[com_r] * inv[com_n - com_r] % md

    n, k = MI()
    aa = LI()
    # combinationの準備
    md = 10 ** 9 + 7
    n_max = n + 3
    fac = [1] * (n_max + 1)
    inv = [1] * (n_max + 1)
    for i in range(2, n_max + 1): fac[i] = fac[i - 1] * i % md
    inv[n_max] = pow(fac[n_max], md - 2, md)
    for i in range(n_max - 1, 1, -1): inv[i] = inv[i + 1] * (i + 1) % md

    dp = [0] * (n + 1)
    dp[0] = 1
    for a in aa:
        for j in range(n, -1, -1):
            pre = dp[j]
            if pre == 0: continue
            dp[j] = 0
            for dj in range(a + 1):
                nj = j + dj
                if nj > n or a - dj < 0 or n - j < dj: break
                dp[nj] += pre * com(n - j, dj) * com(n - dj, a - dj)
                dp[nj] %= md
        # print(dp)
    print(dp[n])

main()

Submission status: Accepted

Solution code:
# http://drken1215.hatenablog.com/entry/2020/01/12/011700

# 本質的な速度改善の方法は分からない
# dpを二次元配列にするとTLEするので、
# 配列2本で処理した

# 修正点:
# [1]
# enumerate(a) -> a
# 不要なenumerateを消した
# [2]
# counted in range(n + 1) -> enumerate(dp)
# enumerateの方がわずかに速いはず
# [3]
# cmbのif文を消した

def main():
    MOD = 10 ** 9 + 7

    N, K = map(int, input().split())
    *a, = map(int, input().split())

    # -+-+-+-prepare combination-+-+-+-
    # -+-+-+-START-+-+-+-
    fact = [0] * (N + 1)
    fact[0] = 1
    fact[1] = 1
    for x in range(2, N + 1):
        fact[x] = x * fact[x - 1] % MOD

    invs = [1] * (N + 1)
    invs[N] = pow(fact[N], MOD - 2, MOD)
    for x in range(N - 1, 0, -1):
        invs[x] = invs[x + 1] * (x + 1) % MOD

    def cmb(n, r):
        return (((invs[r] * invs[n - r]) % MOD) * fact[n]) % MOD

    # -+-+-+-END-+-+-+-

    dp = [0] * (N + 1)
    dp[0] = 1

    for ai in a:
        ndp = [0] * (N + 1)
        for counted, dp_counted in enumerate(dp):
            for to_count in range(min(N - counted, ai) + 1):
                ndp[counted + to_count] = (
                                                  ndp[counted + to_count]
                                                  + dp_counted
                                                  * cmb(N - counted, to_count)
                                                  * cmb(N - to_count, ai - to_count)
                                          ) % MOD
        dp = ndp

    print(dp[N])


if __name__ == '__main__':
    main()

Submission status: Accepted

Solution code:
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

import numpy as np

N,K,*A = map(int,read().split())

MOD = 10 ** 9 + 7

def cumprod(A, MOD = MOD):
    L = len(A); Lsq = int(L**.5+1)
    A = np.resize(A, Lsq**2).reshape(Lsq,Lsq)
    for n in range(1,Lsq):
        A[:,n] *= A[:,n-1]; A[:,n] %= MOD
    for n in range(1,Lsq):
        A[n] *= A[n-1,-1]; A[n] %= MOD
    return A.ravel()[:L]
def make_fact(U, MOD = MOD):
    x = np.arange(U, dtype = np.int64); x[0] = 1
    fact = cumprod(x, MOD)
    x = np.arange(U, 0, -1, dtype=np.int64); x[0] = pow(int(fact[-1]), MOD-2, MOD)
    fact_inv = cumprod(x, MOD)[::-1]
    return fact,fact_inv

def convolve(f, g, MOD=MOD):
    """
    数列 (多項式) f, g の畳み込みの計算．上下 15 bitずつ分けて計算することで，
    30 bit以下の整数，長さ 250000 程度の数列での計算が正確に行える．
    """
    Lf = len(f); Lg = len(g); L = Lf + Lg - 1
    fl = f & (1 << 15) - 1; fh = f >> 15
    gl = g & (1 << 15) - 1; gh = g >> 15
    if L < (1<<8):
        conv = lambda f,g: \
            np.convolve(f,g) % MOD
    else:
        fft = np.fft.rfft; ifft = np.fft.irfft; fft_len = 1 << L.bit_length()
        conv = lambda f,g: \
            np.rint(ifft(fft(f,fft_len) * fft(g,fft_len))[:L]).astype(np.int64) % MOD
    x = conv(fl, gl)
    z = conv(fh, gh)
    y = conv(fl+fh, gl+gh)
    return (x + ((y - x - z) << 15) + (z << 30)) % MOD

fact, fact_inv = make_fact(10 ** 5)

f = np.zeros(N+1,np.int64)
f[0] = 1
for a in A:
    g = fact[N-a:N+1] * fact_inv[0:a+1] % MOD * fact_inv[N-a] % MOD
    g = g[::-1]
    g *= fact_inv[:len(g)]; g %= MOD
    f = convolve(f,g)[:N+1]

f *= fact[:N+1]; f %= MOD

answer = f[N]
print(answer)

Submission status: Accepted

Solution code:
# http://drken1215.hatenablog.com/entry/2020/01/12/011700

def main():
    MOD = 10 ** 9 + 7

    N, K = map(int, input().split())
    *a, = map(int, input().split())

    # -+-+-+-prepare combination-+-+-+-
    # -+-+-+-START-+-+-+-
    fact = [0] * (N + 1)
    fact[0] = 1
    fact[1] = 1
    for x in range(2, N + 1):
        fact[x] = x * fact[x - 1] % MOD

    invs = [1] * (N + 1)
    invs[N] = pow(fact[N], MOD - 2, MOD)
    for x in range(N - 1, 0, -1):
        invs[x] = invs[x + 1] * (x + 1) % MOD

    def cmb(n, r):
        if n - r < r: return cmb(n, n - r)
        if r < 0: return 0
        if r == 0: return 1
        if r == 1: return n
        return (((fact[n] * invs[r]) % MOD) * invs[n - r]) % MOD

    # -+-+-+-END-+-+-+-

    dp = [0] * (N + 1)
    dp[0] = 1

    for i, ai in enumerate(a):
        ndp = [0] * (N + 1)
        for counted in range(N + 1):
            for to_count in range(min(N - counted, ai) + 1):
                ndp[counted + to_count] = (
                                                  ndp[counted + to_count]
                                                  + dp[counted]
                                                  * cmb(N - counted, to_count)
                                                  * cmb(N - to_count, ai - to_count)
                                          ) % MOD
        dp = ndp

    print(dp[N])


if __name__ == '__main__':
    main()

Submission status: Accepted

Solution code:
P = 10**9+7
fa = [1]
for i in range(1, 1010): fa.append(fa[-1] * i % P)
fainv = [pow(fa[-1], P-2, P)]
for i in range(1, 1010)[::-1]: fainv.append(fainv[-1] * i % P)
fainv = fainv[::-1]

k = 70
kk = 1<<k
nu = lambda L: int("".join([bin(kk+a)[-k:] for a in L[::-1]]), 2)
st = lambda n: bin(n)[2:] + "0"
li = lambda s: [int(a, 2) % P if len(a) else 0 for a in [s[-(i+1)*k-1:-i*k-1] for i in range(N+1)]]

N, K = map(int, input().split())
A = [int(a) for a in input().split()]
X = [fa[N]] + [0] * N
for a in A:
    X = li(st(nu(X) * nu([fainv[i] * fa[N-i] * fainv[a-i] * fainv[N-a] % P if i <= a else 0 for i in range(N+1)])))

print(X[-1])

Submission status: Accepted

Solution code:
SIZE=1001; MOD=10**9+7 #998244353 #ここを変更する

SIZE += 1
inv = [0]*SIZE  # inv[j] = j^{-1} mod MOD
fac = [0]*SIZE  # fac[j] = j! mod MOD
finv = [0]*SIZE # finv[j] = (j!)^{-1} mod MOD
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): # nCk mod MOD の計算
    if 0 <= r <= n:
        return (fac[n]*finv[r]%MOD)*finv[n-r]%MOD
    else:
        return 0

# coding: utf-8
# Your code here!
import sys
sys.setrecursionlimit(10**6)
readline = sys.stdin.readline
read = sys.stdin.read

n,k = [int(i) for i in readline().split()]
a = [int(i) for i in readline().split()]

def product(p,q):
    res = [0]*(n+1)
    for i in range(n+1):
        for j in range(n-i+1):
            res[i+j] += p[i]*q[j]
            res[i+j] %= MOD
    return res
    

ans = [1]+[0]*n

for ai in a:
    q = [0]*(n+1)
    for j in range(n+1):
        q[j] = choose(n-j,n-ai)*finv[j]%MOD
    ans = product(ans,q)

res = ans[n]
for i in range(1,n+1):
    res *= i
    res %= MOD

print(res)

Submission status: Accepted

Solution code:
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

import numpy as np

N,K,*A = map(int,read().split())

MOD = 10 ** 9 + 7

def cumprod(A, MOD = MOD):
    L = len(A); Lsq = int(L**.5+1)
    A = np.resize(A, Lsq**2).reshape(Lsq,Lsq)
    for n in range(1,Lsq):
        A[:,n] *= A[:,n-1]; A[:,n] %= MOD
    for n in range(1,Lsq):
        A[n] *= A[n-1,-1]; A[n] %= MOD
    return A.ravel()[:L]
def make_fact(U, MOD = MOD):
    x = np.arange(U, dtype = np.int64); x[0] = 1
    fact = cumprod(x, MOD)
    x = np.arange(U, 0, -1, dtype=np.int64); x[0] = pow(int(fact[-1]), MOD-2, MOD)
    fact_inv = cumprod(x, MOD)[::-1]
    return fact,fact_inv

def convolve(f, g, MOD=MOD):
    """
    数列 (多項式) f, g の畳み込みの計算．上下 15 bitずつ分けて計算することで，
    30 bit以下の整数，長さ 250000 程度の数列での計算が正確に行える．
    """
    Lf = len(f); Lg = len(g); L = Lf + Lg - 1
    fl = f & (1 << 15) - 1; fh = f >> 15
    gl = g & (1 << 15) - 1; gh = g >> 15
    if L < (1<<8):
        conv = lambda f,g: \
            np.convolve(f,g) % MOD
    else:
        fft = np.fft.rfft; ifft = np.fft.irfft; fft_len = 1 << L.bit_length()
        conv = lambda f,g: \
            np.rint(ifft(fft(f,fft_len) * fft(g,fft_len))[:L]).astype(np.int64) % MOD
    x = conv(fl, gl)
    z = conv(fh, gh)
    y = conv(fl+fh, gl+gh)
    return (x + ((y - x - z) << 15) + (z << 30)) % MOD

fact, fact_inv = make_fact(10 ** 5)

f = np.zeros(N+1,np.int64)
f[0] = 1
for a in A:
    g = fact[N-a:N+1] * fact_inv[0:a+1] % MOD * fact_inv[N-a] % MOD
    g = g[::-1]
    g *= fact_inv[:len(g)]; g %= MOD
    f = convolve(f,g)[:N+1]

f *= fact[:N+1]; f %= MOD

answer = f[N]
print(answer)

Submission status: Runtime Error

Solution code:
M=10**9+7;n,k,*A=map(int,open(0).read().split());C=[[0]*(n+1)for _ in"_"*(n+1)]
for i in range(n+1):
 for j in range(n):C[i][j]=(1,C[i-1][j-1]+C[i-1][j])[j>0]%M
F=[[0]*(n+1)for _ in"_"*n];F[0][0]=1
for i in range(k):
 for j in range(n+1):
  for m in range(j+1):F[i+1][j]+=C[j][m]*F[i][m]*C[n-j+m][A[i]-j+m]
print(F[k][n]%M)

Submission status: Runtime Error

Solution code:
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
//#include <boost/multiprecision/cpp_int.hpp>
//typedef boost::multiprecision::cpp_int ll;
typedef long double dd;

#define i_7 (ll)(1E9+7)
//#define i_7 998244353
#define i_5 i_7-2
ll mod(ll a){
    ll c=a%i_7;
    if(c>=0)return c;
    return c+i_7;
}
 
typedef pair<ll,ll> l_l;
ll inf=(ll)1E16;
#define rep(i,l,r) for(ll i=l;i<=r;i++)
#define pb push_back
ll max(ll a,ll b){if(a<b)return b;else return a;}
ll min(ll a,ll b){if(a>b)return b;else return a;}

void Max(ll &pos,ll val){pos=max(pos,val);}//Max(dp[n],dp[n-1]);
void Min(ll &pos,ll val){pos=min(pos,val);}
void Add(ll &pos,ll val){pos=mod(pos+val);}
dd EPS=1E-9;
#define fastio ios::sync_with_stdio(false); cin.tie(0); cout.tie(0);
#define fi first
#define se second
////////////////////////////

#define N 20004

ll kai[N];
ll kai2[N];
ll po(ll i,ll p){
    if(p==0)return 1;
    else{
        i=mod(i);
        if(p==1)return i;
        if(p%2==0)return po(mod(i*i),p/2);
        return mod(i*po(i,p-1));
    }
}
void calc(){
    kai[0]=1;
    kai2[0]=1;
    rep(i,1,N-1){
        kai[i]=mod(kai[i-1]*i);
    }
    kai2[N-1]=po(kai[N-1],i_5);
    for(ll i=N-2;i>=0;i--){
        kai2[i]=mod(kai2[i+1]*(i+1));
    }
}
ll comb(ll n,ll k){
    if(n<k)return 0;
    if(n==0)return 1;
    return mod(mod(kai[n]*kai2[n-k])*kai2[k]);
}
ll bunbo(ll n){
    return po(n,i_5);
}

int main(){fastio
    calc();
    ll n,k;cin>>n>>k;
    ll a[k];rep(i,0,k-1)cin>>a[i];
    ll dp[k+1][n+1];memset(dp,0,sizeof(dp));
    dp[0][0]=1;
    rep(i,0,k-1){
        rep(j,0,n){
            rep(p,0,a[i]){
                if(j+p<=n){
                    Add(dp[i+1][j+p],mod(dp[i][j]*comb(n-j,p))*comb(n-p,a[i]-p));
                }
            }
        }
    }
    /*
    rep(i,0,k){
        rep(j,0,n){
            cout<<dp[i][j]<<" ";
        }cout<<endl;
    }*/
    cout<<mod(dp[k][n])<<endl;
    
    return 0;
}

Submission status: Runtime Error

Solution code:
from math import factorial
from numpy import prod
def probab(n,a):
    return factorial(n)/(factorial(a)*factorial(n-a))
N,K = map(int,input().split())
ListNum = list(map(int,input().split()))
num = sum(ListNum)
List = [num//N]*N
while(sum(List)<=num):
    index = List.index(min(List))
    List[index]+=1
    num-=1
print(int((prod(List)*prod(list(probab(N,i) for i in ListNum)))%(10**9 +7)))

Submission status: Runtime Error

Solution code:
import operator as op
from functools import reduce

N, K = map(int, input().split())
cookies = list(map(int, input().split()))
happines = (sum(cookies) - N) * K
multiplied = 1

def ncr(n, r):
    r = min(r, n-r)
    numer = reduce(op.mul, range(n, n-r, -1), 1)
    denom = reduce(op.mul, range(1, r+1), 1)
    return numer / denom

for i in cookies:
    print(ncr(N,i))
    multiplied = multiplied * ncr(N, i)
res = int(happines * multiplied)
print(res)

Submission status: Runtime Error

Solution code:
from math import factorial
import random


def get_comb(n, r):
    combination = int(factorial(n) / factorial(r) / factorial(n - r))
    return combination


n, k = map(int, input().split())
a_list = list(map(int, input().split()))
children = [0 for _ in range(n)]
happiness = 1

for a in a_list:
    c = get_comb(n, a)
    random_num = random.sample([i for i in range(n)], a)
    for r in random_num:
        children[r] += 1
    happiness *= c

child_expected = 1
for c in children:
    child_expected *= c
print(happiness*child_expected % (10**9+7))