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) .

Solution (Python):
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import math
import bisect
import random
from itertools import permutations, accumulate, combinations, product
import sys
from pprint import pprint
from copy import deepcopy
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor
from operator import mul
from functools import reduce
from pprint import pprint


sys.setrecursionlimit(2147483647)
INF = 10 ** 20
def LI(): return list(map(int, sys.stdin.readline().split()))
def I(): return int(sys.stdin.readline())
def LS(): return sys.stdin.buffer.readline().rstrip().decode('utf-8').split()
def S(): return sys.stdin.buffer.readline().rstrip().decode('utf-8')
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]
mod = 1000000007


n, k = LI()

fac = [1] * (n + 1)
inv = [1] * (n + 1)
for j in range(1, n + 1):
    fac[j] = fac[j-1] * j % mod


inv[n] = pow(fac[n], mod-2, mod)
for j in range(n-1, -1, -1):
    inv[j] = inv[j+1] * (j+1) % mod

def comb(n, r):
    if r > n or n < 0 or r < 0:
        return 0
    return fac[n] * inv[n - r] * inv[r] % mod




A = LI()
dp = [[0] * (n + 1) for _ in range(k + 1)]
dp[0][0] = 1
for i in range(k):
    for j in range(n + 1):
        for l in range(min(n - j, A[i]) + 1):
            dp[i + 1][j + l] = (dp[i + 1][j + l] + dp[i][j] * comb(n - j, l) % mod
                                * comb(n - l, A[i] - l)) % mod

print(dp[k][n])

Solution (Python):
mod = 10**9+7
M = 1000

fact = [1]*(M+1)
factinv = [1]*(M+1)
inv = [1]*(M+1)
for n in range(2,M+1):
    fact[n] = n*fact[n-1]%mod
    inv[n] = -inv[mod%n]*(mod//n)%mod
    factinv[n] = inv[n]*factinv[n-1]%mod

def comb(n,k):
    return fact[n]*factinv[k]*factinv[n-k]%mod

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

dp = [[0]*(N+1) for i in range(K+1)]

for k in range(K):
    dp[k][0] = 1
    for n in range(1,N+1):
        for i in range(min(a[k]+1,n+1)):
            dp[k+1][n] += fact[a[k]]*factinv[a[k]-i]*fact[N-i]*factinv[N]*comb(n,i)*dp[k][n-i]
            dp[k+1][n] %= mod
ans = dp[-1][-1]
for i in range(K):
    ans *= comb(N,a[i])
    ans %= mod

print(ans)

Solution (Python):
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(i+1):C[i][j]=(1,C[i-1][j-1]+C[i-1][j])[j>0]%M
F=[[0]*(n+1)for _ in"_"*(k+1)];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)

Solution (Python):
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(i+1):C[i][j]=(1,C[i-1][j-1]+C[i-1][j])[j!=0]%M
F=[[0]*(n+1)for _ in"_"*(k+1)];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)

Solution (Python):
MOD = 10**9 + 7
list_size = 1001

f_list = [1] * list_size
f_r_list = [1] * list_size

for i in range(list_size-1):
	f_list[i+1] = (f_list[i] * (i+1)) % MOD

f_r_list[-1] = pow(f_list[-1], MOD - 2, MOD)

for i in range(list_size-2, -1, -1):
	f_r_list[i] = (f_r_list[i+1] * (i+1)) % MOD

def comb(n, r):
	if n < r or r < 0:
		return 0
	elif n == 0 or r == 0 or n == r:
		return 1
	else:
		return (f_list[n] * f_r_list[n-r] * f_r_list[r]) % MOD 

n, k = map(int, input().split())
a = list(map(int, input().split()))

dp = [[0 for _ in range(n+1)] for _ in range(k+1)]
dp[0][0] = f_list[n]
for i in range(1, k+1):
	for j in range(n+1):
		for x in range(j+1):
			dp[i][j] += dp[i-1][j-x] * comb(n-x, a[i-1]-x) * f_r_list[x]
			dp[i][j] %= MOD

print(dp[k][n])

Solution (Python):
import sys
input = sys.stdin.readline
U = 1000
MOD = 10**9+7
 
fact = [1]*(U+1)
fact_inv = [1]*(U+1)
 
for i in range(1,U+1):
    fact[i] = (fact[i-1]*i)%MOD
fact_inv[U] = pow(fact[U], MOD-2, MOD)
 
for i in range(U,0,-1):
    fact_inv[i-1] = (fact_inv[i]*i)%MOD

def comb(n, k):
    if k < 0 or k > n:
        return 0
    z = fact[n]
    z *= fact_inv[k]
    z %= MOD
    z *= fact_inv[n-k]
    z %= MOD
    return z

n, k = map(int, input().split())
A = list(map(int, input().split()))
F = [[0]*(n+1) for _ in range(k+1)]
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] += comb(j, m) * F[i][m] * comb(n-(j-m), A[i]-(j-m))
    F[i+1][j] %= MOD
ans = F[k][n]
print(ans)

Solution (Python):
n,k = map(int,input().split())
a = tuple(map(int,input().split()))
mod = 10**9+7
rng = 1001
fctr = [1]
finv = [1]
for i in range(1,rng):
  fctr.append(fctr[-1]*i%mod)
  finv.append(pow(fctr[i],mod-2,mod))
dp = [[0 for i in range(n+1)] for j in range(k+1)]
dp[0][0] = 1
for i in range(1,k+1):
  x = a[i-1]
  for j in range(n+1):
    if dp[i-1][j] == 0:
      continue
    y = dp[i-1][j]*fctr[n-j]*finv[n-x]%mod
    for t in range(min(n-j,x)+1):
      dp[i][j+t]=(dp[i][j+t]+y*fctr[n-t]*finv[n-j-t]*finv[t]*finv[x-t])%mod
print(dp[k][n])

Solution (Python):
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]*fracinv[a[0]-x]%MOD*fracinv[n-a[0]]*fracinv[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]*frac[n-x]%MOD*fracinv[ai-x]*fracinv[n-ai]%MOD*fracinv[x]
                dp[i+1][p]=(dp[i+1][p]+tmp)%MOD
    print(dp[k-1][n])

if __name__ == "__main__":
    main()

Solution (Python):
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]*fracinv[a[0]-x]%MOD*fracinv[n-a[0]]*fracinv[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]*frac[n-x]%MOD*fracinv[ai-x]%MOD*fracinv[n-ai]%MOD*fracinv[x]%MOD
                dp[i+1][p]=(dp[i+1][p]+tmp)%MOD
    print(dp[k-1][n]%MOD)

if __name__ == "__main__":
    main()

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

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

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

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

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

    

    DP=NDP

print(DP[N])

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

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

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

main()

Solution (Python):
















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

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

    
    
    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

    

    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()

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

Solution (Python):


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

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

    
    
    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

    

    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()

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

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

Solution (Python):
nnnnn = 1010
P = 10 ** 9 + 7

fa = [1]
for i in range(1, nnnnn):
    fa.append(fa[-1] * i % P)
fainv = [pow(fa[-1], P-2, P)]
for i in range(1, nnnnn)[::-1]:
    fainv.append(fainv[-1] * i % P)
fainv = fainv[::-1]

N, K = map(int, input().split())
A = [int(a) for a in input().split()]

def calc(k):
    return [fainv[i] * fa[N-i] * fainv[k-i] * fainv[N-k] % P if i <= k else 0 for i in range(N+1)]

def conv(a, b):
    L = [0] * (N+1)
    for i in range(N+1):
        for j in range(N+1):
            if i+j <= N: L[i+j] += a[i] * b[j]
    return [l % P for l in L]

X = [1] + [0] * N
for a in A:
    X = conv(X, calc(a))

print(X[-1] * fa[N] % P)

Solution (Python):
SIZE=1001; MOD=10**9+7 

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



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)

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