Problem:
 Score : 600 points    Problem Statement There are N people, conveniently numbered 1 through N .
We want to divide them into some number of groups, under the following two conditions:    Every group contains between A and B people, inclusive.    Let F_i be the number of the groups containing exactly i people. Then, for all i , either F_i=0 or C≤F_i≤D holds.    Find the number of these ways to divide the people into groups.
Here, two ways to divide them into groups is considered different if and only if there exists two people such that they belong to the same group in exactly one of the two ways.
Since the number of these ways can be extremely large, print the count modulo 10^9+7 .      Constraints  1≤N≤10^3  1≤A≤B≤N  1≤C≤D≤N         Input The input is given from Standard Input in the following format:  N  A  B  C  D       Output Print the number of ways to divide the people into groups under the conditions, modulo 10^9+7 .        Sample Input 1 3 1 3 1 2      Sample Output 1 4  There are four ways to divide the people:   (1,2),(3)  (1,3),(2)  (2,3),(1)  (1,2,3)   The following way to divide the people does not count: (1),(2),(3) . This is because it only satisfies the first condition and not the second.       Sample Input 2 7 2 3 1 3      Sample Output 2 105  The only ways to divide the people under the conditions are the ones where there are two groups of two people, and one group of three people.
There are 105 such ways.       Sample Input 3 1000 1 1000 1 1000      Sample Output 3 465231251       Sample Input 4 10 3 4 2 5      Sample Output 4 0  The answer can be 0 .

Solution (Python):
from math import floor
N,A,B,C,D = map(int,input().split())
mod = 10**9+7
factorials = [None for i in range(N+1)]
factorials[0] = 1
for i in range(1, N+1):
    factorials[i] = (factorials[i-1]*i) % mod
def inv(x):
    ret = 1
    k = mod - 2
    y = x
    while k:
        if (k & 1):
            ret = (ret * y) % mod
        y = (y * y) % mod
        k //= 2
    
    return ret

finv = [0]*(N+1)
finv[N] = inv(factorials[N])
for i in range(N, 0, -1):
    finv[i-1] = (finv[i]*i)%mod
def calc(i,j,k):
    tmp = finv[N-j]
    tmp = (tmp*finv[k])%mod
    tmp = (tmp*factorials[N-j+i*k])%mod
    y = finv[i]
    ret = 1
    while k:
        if (k & 1):
            ret = (ret * y) % mod
        y = (y * y) % mod
        k //= 2
    return (ret*tmp) % mod

dp = [[0]*(N+1) for _ in range(N+1)]
for i in range(B+1):
    dp[i][0] = 1
ls = [0]+list(range(C,D+1))
l = len(ls)
for i in range(A,B+1):
    for j in range(1,N+1):
        tmp = 0
        for k in ls:
            if k > j/i:
                break
            tmp = (tmp + dp[i-1][j-i*k]*calc(i,j,k)) % mod
        dp[i][j] = tmp
print(dp[B][N])

Solution (Python):
def main():
    import sys
    input = sys.stdin.readline

    mod = 1000000007
    nmax = 1001  
    fac = [0] * (nmax + 1)
    finv = [0] * (nmax + 1)
    inv = [0] * (nmax + 1)
    fac[0] = 1
    fac[1] = 1
    finv[0] = 1
    finv[1] = 1
    inv[1] = 1
    for i in range(2, nmax):
        fac[i] = fac[i - 1] * i % mod
        inv[i] = mod - inv[mod % i] * (mod // i) % mod
        finv[i] = finv[i - 1] * inv[i] % mod

    N, A, B, C, D = map(int, input().split())

    dp = [[0] * (N+1) for _ in range(B+1)]
    for i in range(B+1):
        dp[i][0] = 1

    for i in range(A, B+1):
        for j in range(N+1):
            dp[i][j] = dp[i-1][j]
        for j in range(N+1):
            if (N-j) // i < C:
                continue
            n_max = min(D, (N-j)//i)
            for n in range(C, n_max+1):
                tmp = (((fac[N-j] * finv[N-j-i*n])%mod * pow(finv[i], n, mod))%mod * finv[n])%mod
                dp[i][j+i*n] = (dp[i][j+i*n] + (dp[i-1][j] * tmp)%mod) % mod
    print(dp[-1][-1])


if __name__ == '__main__':
    main()

Solution (Python):
mod=10**9+7
MAX=10**3+100


g1=[1,1]
g2=[1,1]
for i in range(2,MAX+1):
    num_1=g1[-1]*i%mod
    g1.append(num_1)
    g2.append(pow(num_1,mod-2,mod))
    
def cmb(n,r,MOD):
    return g1[n]*g2[r]*g2[n-r]%MOD



N,A,B,C,D=map(int,input().split())

data=[[] for i in range(B+1)]
for i in range(A,B+1):
    for j in range(C,D+1):
        if i*j>N:
            break
        data[i].append([i*j, pow(g2[i], j, mod)*g2[j]%mod])

dp=[0]*(N+1)
dp[0]=1
for i in range(A,B+1):
    H=dp[:]
    for u,v in data[i]:
        for j in range(u,N+1):
            H[j]=(H[j]+dp[j-u]*v)%mod
    dp=H
print(dp[N]*g1[N]%mod)

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

NN = 10**4 
g1 = [1, 1] 
g2 = [1, 1] 
inverse = [0, 1] 

for i in range( 2, NN + 1 ):
    g1.append( ( g1[-1] * i ) % mod )
    inverse.append( ( -inverse[mod % i] * (mod//i) ) % mod )
    g2.append( (g2[-1] * inverse[-1]) % mod )

N, A, B, C, D = map(int, input().split())

Map = [[] for _ in range(N+1)]
for n in range(1, N+1):
    a = 1
    for _ in range(N+1):
        Map[n].append(a)
        a = a * g2[n] % mod

L = [[] for _ in range(N+1)]
for n in range(A, B+1):
    for c in range(C, D+1):
        if n*c > N:
            continue
        L[n].append(c)

dp = [0]*(N+1)
dp[0] = 1
for n in range(1, N+1):
    for i, c in enumerate(L[n]):
        if i == 0:
            dp2 = dp[:]
        for nn in range(N+1-n*c):
            dp2[nn+n*c] = (dp2[nn+n*c] + dp[nn] * Map[n][c] % mod * Map[c][1] % mod) % mod
        if i == len(L[n])-1:
            dp = dp2

ans = dp[N] * g1[N] % mod
print(ans)

Solution (Python):
import os
import sys



if os.getenv("LOCAL"):
    sys.stdin = open("_in.txt", "r")

sys.setrecursionlimit(10 ** 9)
INF = float("inf")
IINF = 10 ** 18
MOD = 10 ** 9 + 7





def get_factorials(max, mod=None):
    """
    階乗 0!, 1!, 2!, ..., max!
    :param int max:
    :param int mod:
    :return:
    """
    ret = [1]
    n = 1
    if mod:
        for i in range(1, max + 1):
            n *= i
            n %= mod
            ret.append(n)
    else:
        for i in range(1, max + 1):
            n *= i
            ret.append(n)
    return ret


def mod_invs(max, mod):
    """
    逆元のリスト 0 から max まで
    :param int max:
    :param int mod:
    """
    invs = [1] * (max + 1)
    for x in range(2, max + 1):
        invs[x] = (-(mod // x) * invs[mod % x]) % mod
    return invs


def factorial_invs(max, mod):
    """
    階乗 0!, 1!, 2!, ..., max! の逆元
    :param int max:
    :param int mod:
    """
    ret = []
    r = 1
    for inv in mod_invs(max, mod):
        r = r * inv % mod
        ret.append(r)
    return ret


class Combination:
    def __init__(self, max, mod):
        """
        :param int max:
        :param int mod: 3 以上の素数であること
        """
        self._factorials = get_factorials(max, mod)
        self._finvs = factorial_invs(max, mod)
        self._mod = mod

    
    def ncr(self, n, r):
        """
        :param int n:
        :param int r:
        :rtype: int
        """
        if n < r:
            return 0
        return self._factorials[n] * self._finvs[r] % self._mod * self._finvs[n - r] % self._mod


N, A, B, C, D = list(map(int, sys.stdin.buffer.readline().split()))

comb = Combination(max=N + 1, mod=MOD)

dp = [[0] * (N + 1) for _ in range(N + 1)]
for u in range(N + 1):
    dp[u][0] = 1
for unit in range(A, B + 1):
    for j in reversed(range(N + 1)):
        dp[unit][j] = dp[unit - 1][j]
        for cnt in range(C, D + 1):
            if j + unit * cnt > N:
                break
            
            ncr = comb._factorials[N - j] * comb._finvs[N - j - unit * cnt] % MOD * \
                  pow(comb._finvs[unit], cnt, MOD) % MOD * comb._finvs[cnt] % MOD
            dp[unit][j + unit * cnt] += dp[unit - 1][j] * ncr % MOD
            dp[unit][j + unit * cnt] %= MOD

print(dp[B][N])

Solution (Python):

import sys
from collections import deque, defaultdict
import copy
import bisect

import math
import heapq
from itertools import product, permutations,combinations
import fractions

import sys
def input():
	return sys.stdin.readline().strip()


def pow(x, y, mod=1000000007):
	pow_list = [x]
	i = 1
	while 2**i <= y:
		a = pow_list[-1]**2
		if mod != 0:
			a = a % mod
		pow_list.append(a)
		i += 1
	ans = 1
	for bit in range(len(pow_list)):
		if (2**bit) & y != 0:
			ans *= pow_list[bit]
			if mod != 0:
				ans = ans % mod
	return ans

def fact(n, mod=100000007):
	ans = 1
	for i in range(n):
		ans *= i + 1
		if mod != 0 and ans >= mod:
			ans = ans % mod
	return ans

def mod_rev(x, mod):
	"""
	関数powが必要
	"""
	return pow(x, mod - 2, mod)

N, A, B, C, D = list(map(int, input().split()))
mod = 1000000007

mod_sho = [0]
for i in range(N):
	mod_sho.append(mod_rev(i + 1, mod))

comb = [[0]*(N + 1) for _ in range(N + 1)]
comb[0][0] = 1
for n in range(1, N + 1):
	comb[n][0] = 1
	comb[n][n] = 1
	for i in range(1, n):
		comb[n][i] = comb[n - 1][i - 1] + comb[n - 1][i]
		comb[n][i] = comb[n][i] % mod
plus_list = []
for i in range(A, B + 1):
	G = []
	for j in range(C, D + 1):
		if i*j > N:
			break
		else:
			G.append(i*j)
	if len(G) > 0:
		plus_list.append(G)

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

if len(plus_list) == 0:
	print(0)
	exit()
num = A - 1
for number in plus_list:
	
	num += 1
	sum_next = [0] * (N + 1)
	for n in range(0, N + 1):
		if sum[n] == 0:
			continue
		sum_next[n] += sum[n]
		sum_next[n] = sum_next[n] % mod
		times = 1
		rest = N - n
		seki = 1
		for x in range(1, D + 1):

			if n + x*num <= N:
				seki *= comb[rest][num]
				seki = seki*mod_sho[times]
				seki = seki % mod
				if x >= C:
					sum_next[n + x*num] += (sum[n]*seki) % mod
					sum_next[n + x*num] = sum_next[n + x*num] % mod
					
				times += 1
				rest -= num
			else:
				break
	sum = copy.deepcopy(sum_next)


print(sum[N])

Solution (Python):
class Combination:
    """
    O(n)の前計算を1回行うことで，O(1)でnCr mod mを求められる
    n_max = 10**6のとき前処理は約950ms (PyPyなら約340ms, 10**7で約1800ms)
    使用例：
    comb = Combination(1000000)
    print(comb(5, 3))  
    """
    def __init__(self, n_max, mod=10**9+7):
        self.mod = mod
        self.modinv = self.make_modinv_list(n_max)
        self.fac, self.facinv = self.make_factorial_list(n_max)

    def __call__(self, n, r):
        return self.fac[n] * self.facinv[r] % self.mod * self.facinv[n-r] % self.mod

    def make_factorial_list(self, n):
        
        
        fac = [1]
        facinv = [1]
        for i in range(1, n+1):
            fac.append(fac[i-1] * i % self.mod)
            facinv.append(facinv[i-1] * self.modinv[i] % self.mod)
        return fac, facinv

    def make_modinv_list(self, n):
        
        modinv = [0] * (n+1)
        modinv[1] = 1
        for i in range(2, n+1):
            modinv[i] = self.mod - self.mod//i * modinv[self.mod%i] % self.mod
        return modinv

N, A, B, C, D = map(int, input().split())
mod = 10**9+7
comb = Combination(1010)
dp = [0]*(N+1)
dp[0] = 1
for n_person in range(A, B+1):
    dp_new = dp[:]
    for i, v in enumerate(dp):
        c = 1
        rem = N-i
        for n_groups in range(1, C):
            if rem < n_person:
                break
            c = c * comb(rem, n_person) % mod
            rem -= n_person
        else:
            for n_groups in range(C, D+1):
                if rem < n_person:
                    break
                c = c * comb(rem, n_person) % mod
                rem -= n_person
                dp_new[N-rem] = (dp_new[N-rem] + c * comb.facinv[n_groups] * v) % mod
    dp = dp_new
print(dp[N])

Solution (Python):

import sys
import math
from collections import defaultdict


sys.setrecursionlimit(10**7)
def input():
  return sys.stdin.readline()[:-1]


mod = 10**9 + 7


def I(): return int(input())
def II(): return map(int, input().split())
def III(): return list(map(int, input().split()))
def Line(N):
  read_all = [tuple(map(int, input().split())) for _ in range(N)]
  return map(list,zip(*read_all))




def extgcd(a,b):
  r = [1,0,a]
  w = [0,1,b]
  while w[2]!=1:
    q = r[2]//w[2]
    r2 = w
    w2 = [r[0]-q*w[0],r[1]-q*w[1],r[2]-q*w[2]]
    r = r2
    w = w2
  return [w[0],w[1]]


def mod_inv(a,M=mod):
  x = extgcd(a,M)[0]
  return (M+x%M)%M


def mod_fact(n,Mod=mod):
  d = [1]*(n+1)
  for i in range(1,n+1):
    d[i] = d[i-1]*i%Mod
  return d

N,A,B,C,D = II()

value = [0]
for i in range(C,D+1):
  value.append(i)

fact = mod_fact(N)

inv = [0]*(N+1)
inv[N] = mod_inv(fact[N])
for i in range(1,N+1)[::-1]:
  inv[i-1] = inv[i]*i%mod

dp=[[0]*(N+1) for _ in range(B+1)]
for i in range(B):
  dp[i+1][0] = 1
for i in range(A,B+1):
  for j in range(A,N+1):
    if i==A:
      if j%A!=0:
        continue
      else:
        num = j//A
        if C<=num<=D:
          dp[i][j] = (fact[N]*inv[N-j]*inv[num]*pow(inv[A],num,mod))%mod
    else:
      for k in value:
        if k>j//i:
          break
        else:
          dp[i][j] = (dp[i][j] + dp[i-1][j-i*k]*fact[N-j+i*k]*inv[N-j]*pow(inv[i],k,mod)*inv[k])%mod
print(dp[B][N])

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

N,A,B,C,D = map(int,input().split())

MOD = 10 ** 9 + 7

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 = 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

dp = np.zeros(N+1, dtype=np.int64)
dp[0] = 1
for x in range(A,B+1):
    
    prev = dp
    dp = prev.copy()
    for n in range(C,D+1):
        y = n * x
        if y > N:
            break
        
        
        coef = fact[y] * pow(int(fact_inv[x]), n, MOD) % MOD * fact_inv[n] % MOD
        
        dp[y:] += prev[:-y] * comb[N:y-1:-1,y] % MOD * coef
        dp %= MOD

answer = dp[N]
print(answer)

Solution (Python):
N, A, B, C, D = map(int, input().split())
MOD = 10 ** 9 + 7

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

finv = [1] * (N + 1)
finv[-1] = pow(fac[-1], MOD - 2, MOD)
for i in reversed(range(1, N + 1)):
    finv[i - 1] = finv[i] * i % MOD

dp = [[0] * (N + 1) for _ in range(B + 1)]
dp[A - 1][0] = 1

for i in range(A, B + 1):
    for j in range(N + 1):
        dp[i][j] = dp[i - 1][j]
        for k in range(C, min(D, j // i) + 1):
            t = (fac[N - j + i * k] * finv[N - j] \
                      * pow(finv[i], k, MOD) * finv[k]) % MOD
            dp[i][j] += dp[i - 1][j - i * k] * t
            dp[i][j] %= MOD

print(dp[B][N])

Solution (Python):
from itertools import*
import math
from collections import*
from heapq import*
from bisect import bisect_left,bisect_right
from copy import deepcopy
inf = 10**18
mod = 10**9+7
from functools import reduce
import sys
sys.setrecursionlimit(10**7)

n,a,b,c,d = map(int,input().split())
Max = 1005

fac = [0]*(Max)
finv = [0]*(Max)
inv = [0]*(Max)
fac[0]=fac[1]=1
finv[0]=finv[1]=1
inv[1]=1
for i in range(2,Max):
        fac[i] = fac[i-1] * i % mod
        inv[i] = mod - inv[mod%i]*(mod//i)%mod
        finv[i] = finv[i-1]*inv[i]%mod

def C(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

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


dp = [[0]*(n+1) for i in range(n+1)]
for i in range(n+1):
    dp[i][0]=1

for i in range(a,b+1):
    for j in range(1,n+1):
        dp[i][j] = dp[i-1][j]
        k = c
        while k <= j//i and k <= d:
            x = (fac[j] * finv[j - k * i] * finv[k] * pow(finv[i], k, mod)) % mod
            dp[i][j] += (dp[i - 1][j - k * i] * x) % mod
            k += 1
        dp[i][j] %= mod
print(dp[b][n])

Solution (Python):
from itertools import*
import math
from collections import*
from heapq import*
from bisect import bisect_left,bisect_right
from copy import deepcopy
inf = 10**18
mod = 10**9+7
from functools import reduce
import sys
sys.setrecursionlimit(10**7)

n,a,b,c,d = map(int,input().split())
Max = 1005

fac = [0]*(Max)
finv = [0]*(Max)
inv = [0]*(Max)
fac[0]=fac[1]=1
finv[0]=finv[1]=1
inv[1]=1
for i in range(2,Max):
        fac[i] = fac[i-1] * i % mod
        inv[i] = mod - inv[mod%i]*(mod//i)%mod
        finv[i] = finv[i-1]*inv[i]%mod

def C(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

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


dp = [[0]*(n+1) for i in range(n+1)]
for i in range(n+1):
    dp[i][0]=1

for i in range(a,b+1):
    for j in range(1,n+1):
        dpij = dp[i-1][j]
        k = c
        while k <= j//i and k <= d:
            x = (fac[j] * finv[j - k * i] * finv[k] * pow(finv[i], k, mod)) % mod
            dpij += (dp[i - 1][j - k * i] * x) % mod
            k += 1
        dp[i][j] = dpij%mod
print(dp[b][n])

Solution (Python):
import sys
stdin = sys.stdin
 
sys.setrecursionlimit(10**5) 
 
def li(): return map(int, stdin.readline().split())
def li_(): return map(lambda x: int(x)-1, stdin.readline().split())
def lf(): return map(float, stdin.readline().split())
def ls(): return stdin.readline().split()
def ns(): return stdin.readline().rstrip()
def lc(): return list(ns())
def ni(): return int(stdin.readline())
def nf(): return float(stdin.readline())

def fact(n: int, mod: int):
    fac = []
    for i in range(n+1):
        fac.append(1 if i < 2 else i * fac[-1] % mod)
        
    return fac

def fact_inv(fac: list, mod: int):
    return [1 if i < 2 else pow(faci, mod-2, mod) for i, faci in enumerate(fac)]

def nCr(n: int, r: int, mod: int, fac: list, fac_inv: list):
    return fac[n] * fac_inv[r] * fac_inv[n-r] % mod

def nPr(n: int, r: int, mod: int, fac: list, fac_inv: list):
    return fac[n] * fac_inv[n-r] % mod



n,a,b,c,d = li()
MOD = 10**9+7

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


fac = fact(n, MOD)
fac_inv = fact_inv(fac, MOD)

fac_i_k_mod = [[0]*(n+1) for _ in range(n+1)]
fac_inv_i_k_mod = [[0]*(n+1) for _ in range(n+1)]
for i in range(n+1):
    for k in range(n+1):
        fac_i_k_mod[i][k] = pow(fac[i], k, MOD)
        
for i in range(n+1):
    for k in range(n+1):
        fac_inv_i_k_mod[i][k] = pow(fac_i_k_mod[i][k], MOD-2, MOD)        


for i in range(1, n+1):
    for j in range(n+1):
        dp[i][j] += dp[i-1][j]
        if i < a or i > b:
            continue
        
        for k in range(c, min(d, j//i)+1):
            if j-i*k < 0:
                continue
            
            dp[i][j] += dp[i-1][j-i*k] * nPr(j, i*k, MOD, fac, fac_inv) % MOD\
                        * fac_inv_i_k_mod[i][k] % MOD\
                        * fac_inv[k] % MOD
            dp[i][j] %= MOD
            
print(dp[n][n])

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] = int((f_list[i] * (i + 2)) % MOD)

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

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

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

def comb_suc(n, r, k):
	if n == r*k:
		return (f_list[n-1] * pow(f_r_list[r-1], k, MOD) * f_r_list[k-1]) % MOD
	else:
		return (f_list[n-1] * f_r_list[n-(r*k)-1] * pow(f_r_list[r-1], k, MOD) * f_r_list[k-1]) % MOD

import sys
def input():
	return sys.stdin.readline()[:-1]


n, a, b, c, d = map(int, input().split())
dp = [[0 for _ in range(n+1)] for _ in range(b+1)]
dp[a-1][0] = 1
for i in range(a, b+1):
	for j in range(n+1):
		dp[i][j] += dp[i-1][j]
		if j-c*i >= 0:
			for k in range(c, min(d+1, j//i+1)):
				
				dp[i][j] += dp[i-1][j-k*i] * comb_suc(n-(j-k*i), i, k)
		dp[i][j] %= MOD


print(dp[b][n])

Solution (Python):
N, A, B, C, D = map(int, input().split())


MOD = 10 ** 9 + 7
M = 1005
fact = [0] * M
fact[0] = 1
for i in range(1, M):
    fact[i] = fact[i - 1] * i
    fact[i] %= MOD

rfact = [0] * M
rfact[M - 1] = pow(fact[M - 1], MOD - 2, MOD)
for i in range(M - 2, -1, -1):
    rfact[i] = rfact[i + 1] * (i + 1)
    rfact[i] %= MOD

dp = [[0] * (N + 1) for i in range(N + 1)]
dp[A - 1][0] = 1
for i in range(A, B + 1):
    for j in range(N + 1):
        dp[i][j] += dp[i - 1][j]
        k = C
        while k <= D and j - k * i >= 0:
            dp[i][j] += dp[i - 1][j - k * i] * fact[N - (j - k * i)] * pow(rfact[i], k, MOD) * rfact[k] * rfact[N - j]
            dp[i][j] %= MOD
            k += 1

print(dp[B][N])

Solution (Python):
N, A, B, C, D = map(int, input().split())
mod = 7 + 10 ** 9

fact = [1] * (N+1)
frev = [1] * (N+1)
for i in range(1, N+1):
    temp = fact[i] = (fact[i-1] * i) % mod
    frev[i] = pow(temp, mod-2, mod)

def P(n, r):
    return (fact[n] * frev[n-r]) % mod

DP = [[0 for j in range(N+1)] for i in range(N+1)]
for i in range(N+1): DP[i][0] = 1
for i in range(A, B+1):
    revi = frev[i]
    for j in range(1, N+1):
       dpij = 0 + DP[i-1][j]
       for k in range(C, D + 1):
           if j - i * k < 0: break
           temp = (P(N-j+i*k, i*k) * pow(revi, k, mod) * frev[k]) % mod
           dpij += (temp * DP[i-1][j-i*k]) % mod
       DP[i][j] = dpij % mod

print(DP[B][N])

Solution (Python):
import sys
sys.setrecursionlimit(1000000)

MOD = 10**9 + 7

def ex_euclid(a,b,c,d):
    if (b == 1) : return d
    x = a // b
    y = a % b
    return ex_euclid(b,y,d,(c[0]-x*d[0],c[1]-x*d[1]))
    
 
def mod_inv(a,p):
    x = ex_euclid(p,a,(1,0),(0,1))
    return x[1] % p

n,a,b,c,d = map(int,input().split())

gp = b-a+1

dp = [[-1 for j in range(n+1)] for i in range(gp)]

tb1=[0,1]
tb2=[1] 
for i in range(2,n+1):
    tb1.append((tb1[i-1]*i) % MOD)
for i in range(1,n+1):
    tb2.append(mod_inv(tb1[i],MOD))
    

for i in range(gp):dp[i][0] = 1


def dfs(sum,k):
    summ = sum
    if(k < 0):
        if (sum != 0):return 0
        else : return 1
    if dp[k][sum] != -1: return dp[k][sum]
    kk = k+a
    temp = dfs(sum,k-1)
    temp1 = 1
    for i in range(c-1):
        if (sum < 0):break
        temp1 = (temp1 * tb1[sum] * tb2[sum-kk]*tb2[kk]) % MOD 
        sum -= kk
    for i in range(c,d+1):
        if(sum-kk < 0) : break
        temp1 = (temp1 * tb1[sum] * tb2[sum-kk]*tb2[kk]) % MOD
        sum -= kk
        temp = (temp + temp1*dfs(sum,k-1)*tb2[i]) % MOD

    dp[k][summ] = temp
    return temp
    

        
ans = dfs(n,gp-1)



print(ans)

Solution (Python):
from collections import defaultdict
N,A,B,C,D=map(int,input().split())
mod=10**9+7
table=[1]*(N+3)
t=1
for i in range(1,N+3):
    t*=i
    t%=mod
    table[i]=t
rtable=[1]*(N+3)
t=1
for i in range(1,N+3):
    t*=pow(i,mod-2,mod)
    t%=mod
    rtable[i]=t
dd=defaultdict(int)
dd[N]=1
for i in range(A,B+1):
    nd=defaultdict(int)
    for gro,num in dd.items():
        nd[gro]+=num
        t=gro//i
        for k in range(C,min(t,D)+1):
            
            u=(table[gro] * rtable[gro - k * i] * pow(rtable[i],k,mod) * rtable[k] )%mod
            nd[gro - k*i] += num* u
            nd[gro - k * i]%=mod
    dd=nd
    
print(dd[0]%mod)

Solution (Python):
from functools import lru_cache
n,a,b,c,d=map(int,input().split())

mod=10**9+7
fac=[1]*(n+3)
finv=[1]*(n+3)

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

finvp = [finv]
for i in range(n-1):
    p=finvp[-1]
    q=[p[i]*finv[i]%mod for i in range(n+1)]
    finvp.append(q)

dp = [[0 for _ in range(n+1)] for _ in range(b+1)]
dp[a-1][0] = 1

for i in range(a-1,b):
    for j in range(n+1):
        dp[i+1][j]=dp[i][j]
        for k in range(c,1+min(d,j//(i+1))):
            x=j-k*(i+1)
            dp[i+1][j]+=fac[j]*finv[x]*finvp[k-1][i+1]*finv[k]*dp[i][x]%mod
        dp[i+1][j]%=mod
print(dp[b][n])

Solution (Python):
from collections import defaultdict
N,A,B,C,D=map(int,input().split())

mod=10**9+7
table=[1]*(N+3)
t=1
for i in range(1,N+3):
    t*=i
    t%=mod
    table[i]=t
rtable=[1]*(N+3)
t=1
for i in range(1,N+3):
    t*=pow(i,mod-2,mod)
    t%=mod
    rtable[i]=t
dp=defaultdict(int)
dp[N]=1
for l in range(A,B+1):
    ndp=defaultdict(int)
    for mem,num in dp.items():
        ndp[mem]+=num
        for k in range(C,D+1):
            if mem-k*l<0:
                break
            t=(table[mem]*pow(rtable[l],k,mod)*rtable[mem-k*l]*rtable[k])%mod
            ndp[mem-k*l]+=(num*t)%mod
            ndp[mem-k*l]%mod
            

    dp=ndp
    
print(dp[0]%mod)

Solution (Python):
N,A,B,C,D=map(int,input().split())
mod=10**9+7

F=[1]*(N+1)
for i in range(1,N+1):
  F[i]=F[i-1]*(i)%mod
def power(x,y):
  if y==0:
    return 1
  elif y==1:
    return x%mod
  elif y%2==0:
    return power(x,y//2)**2%mod
  else:
    return (power(x,y//2)**2)*x%mod

invF=[1]*(N+1)
invF[N]=power(F[N],mod-2)
for i in range(0,N)[::-1]:
  invF[i]=(invF[i+1]*(i+1))%mod
invp=[[0]*(N+1) for i in range(N+1)]
for i in range(N+1):  
  for k in range(N+1):
    if k==0:
      invp[i][k]=1
    else:
      invp[i][k]=(invp[i][k-1]*invF[i])%mod
dp=[[0]*(N+1) for i in range(B-A+2)]
dp[0][0]=1
for i in range(A,B+1):
  for j in range(N+1):
    dp[i-A+1][j]=dp[i-A][j]
    for k in range(C,min(D,j//i)+1):
      dp[i-A+1][j]+=(dp[i-A][j-k*i]*F[N-j+k*i]*invF[N-j]*invp[i][k]*invF[k])%mod
      dp[i-A+1][j]%=mod
    
print((dp[B-A+1][N])%mod)

Solution (Python):
mod = 10**9+7
def inpl(): return [int(i) for i in input().split()]
N, A, B, C, D = inpl()
fac = [1 for _ in range(N+1)]
for i in range(N):
    fac[i+1] = (i+1)*fac[i] %mod
inv = [1 for _ in range(N+1)]
for i in range(2,N+1):
    inv[i] = (-(mod//i) * inv[mod%i]) %mod    
facinv = [1 for _ in range(N+1)]
for i in range(N):
    facinv[i+1] = inv[i+1]*facinv[i] %mod
facinvp = [facinv]
for i in range(N-1):
    p = facinvp[-1]
    q = [p[i]*facinv[i]%mod for i in range(N+1)]
    facinvp.append(q)
    
dp = [[0 for _ in range(N+1)] for _ in range(B+1)]
dp[A-1][0] = 1

for i in range(A-1,B):
    for j in range(N+1):
        dp[i+1][j] = dp[i][j]
        for k in range(C,1+min(D, j//(i+1))):
            x = j - k*(i+1)
            dp[i+1][j] += fac[j]*facinv[x]*facinvp[k-1][i+1]*facinv[k]*dp[i][x]%mod
        dp[i+1][j] %= mod
print(dp[B][N])

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

def power(x, n, MOD):
    ans = 1
    while n:
        if n % 2 == 1:
            ans = (ans * x) % MOD
        x = (x * x) % MOD
        n //= 2
    return ans

def getFacts(n, MOD):
    facts = [1] + [0] * n
    for x in range(1, n + 1):
        facts[x] = (facts[x - 1] * x) % MOD
    return facts

def getInvFacts(n, factN, MOD):
    invFacts = [0] * (n + 1)
    invFacts[n] = power(factN, MOD - 2, MOD)
    for x in reversed(range(n)):
        invFacts[x] = (invFacts[x + 1] * (x + 1)) % MOD
    return invFacts

def getInvss(B, D, facts, MOD):
    invss = [[0] * (D + 1) for _ in range(B + 1)]
    invss[B][D] = power(power(facts[B], D, MOD), MOD - 2, MOD)
    for i in reversed(range(B)):
        invss[i][D] = (invss[i + 1][D] * power(i + 1, D, MOD)) % MOD
    for i in range(B + 1):
        for k in reversed(range(D)):
            invss[i][k] = (invss[i][k + 1] * facts[i]) % MOD
    return invss


N, A, B, C, D = map(int, input().split())

facts = getFacts(N, MOD)
invFacts = getInvFacts(N, facts[N], MOD)
invss = getInvss(B, D, facts, MOD)


dp = [[0] * (N + 1) for _ in range(B + 1)]
dp[A - 1][0] = 1

for i in range(A, B + 1):
    for j in range(N + 1):
        dp[i][j] = dp[i - 1][j]
        for k in range(C, min(D, j // i) + 1):
            dp[i][j] += facts[N - j + i * k] * invFacts[N - j] * invss[i][k] * invFacts[k] * dp[i - 1][j - i * k]
            dp[i][j] %= MOD

print(dp[B][N])

Solution (Python):
import math
N, A, B, C, D = map(int, input().split())
dp = [[0] * (N + 1) for _ in range(N + 1)]





mod = 10 ** 9 + 7
dp_F = [0] * (2000) 
dp_Finv = [0] * (2000)


def factorial(n):
    if dp_F[n] > 0:
        return dp_F[n]
    elif n == 0:
        dp_F[0] = 1
        return 1
    else:
        dp_F[n] = n * factorial(n - 1) % mod
        return dp_F[n]

def inverse(n):
    
    return pow(n, mod - 2, mod)

def factorial_inv(n):
    if dp_Finv[n] > 0:
        return dp_Finv[n]
    elif n == 0:
        dp_Finv[0] = 1
        return 1
    else:
        dp_Finv[n] = inverse(n) * factorial_inv(n - 1) % mod
        return dp_Finv[n]

def combination(n, k):
    return (factorial(n) * inverse(factorial(n - k)) * inverse(factorial(k))) % mod

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

for i in range(A, B + 1):
    for j in range(1, N + 1):
        dp_ij = 0
        dp_ij += dp[i - 1][j]
        k = C
        while k <= j // i and k <= D:
            x = (factorial(j) * factorial_inv(j - k * i) * factorial_inv(k) * pow(factorial_inv(i), k, mod)) % mod
            dp_ij += (dp[i - 1][j - k * i] * x) % mod
            k += 1
        dp[i][j] = dp_ij % mod
print(dp[B][N])

Solution (Python):
N, A, B, C, D = map(int, input().split())
MOD = 10 ** 9 + 7
fact = [1] * (N + 1)
frev = [1] * (N + 1)

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

for i in range(1, N + 1):
    v = fact[i] = (fact[i - 1] * i) % MOD
    frev[i] = pow(v, MOD - 2, MOD)

Y = [None] * (N + 1)
for g in range(A, min(N // C, B) + 1):
    p = C * g
    q = D * g
    y = frev[g]
    Y[p: min(q, N) + 1: g] = (pow(y, j, MOD) * frev[j]
                              for j in range(C, min(D, N // g) + 1))
    dp[p:] = (dp[i] + frev[N - i] * sum(dp[i - j] * fact[N - i + j] * Y[j]
              for j in range(p, min(q, i) + 1, g)) % MOD
              for i in range(p, N + 1))

print(dp[N] % MOD)

Solution (Python):
N, A, B, C, D = map(int, raw_input().split())
mod = 10 ** 9 + 7


factorial = [1]
inverse = [1]
for i in xrange(1, N + 2):
    factorial.append(factorial[-1] * i % mod)
    inverse.append(pow(factorial[-1], mod - 2, mod))



def nCr(n, r):
    if n < r or r < 0:
        return 0
    elif r == 0:
        return 1
    return factorial[n] * inverse[r] * inverse[n - r] % mod



def f(n, k, x):
    return nCr(n, k * x) * factorial[k * x] * pow(pow(factorial[k], x, mod), mod-2, mod) * inverse[x] % mod


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

for k in xrange(A, B + 1):  
    for n in range(1, N + 1)[::-1]:  
        X = xrange(C, D + 1)
        for x in X:
            if n - k * x < 0:  
                break
            dp[n] += f(n, k, x) * dp[n - k * x] % mod

print(dp[N] % mod)

Solution (Python):
N,A,B,C,D=map(int,raw_input().split());Z=[1]+[0]*N;M=10**9+7;I=[1,1];F=I+[];R=I+[]
for i in range(2,N+1):I+=[(M-M/i)*I[M%i]%M];F+=[i*F[-1]%M];R+=[I[i]*R[-1]%M]
for i in range(A,B+1):Z=[(Z[j]+sum(Z[j-i*k]*F[N-j+i*k]*R[N-j]*pow(R[i],k,M)*R[k]for k in range(C,min(D,j/i)+1)))%M for j in range(N+1)]
print Z[N]

Solution (Python):
N,A,B,C,D=map(int,input().split());Z=[1]+[0]*N;M=10**9+7;I=[1,1];F=I+[];R=I+[]
for i in range(2,N+1):I+=[(M-M//i)*I[M%i]%M];F+=[i*F[-1]%M];R+=[I[i]*R[-1]%M]
for i in range(A,B+1):Z=[(Z[j]+sum(Z[j-i*k]*F[N-j+i*k]*R[N-j]*pow(R[i],k,M)*R[k]for k in range(C,min(D,j//i)+1)))%M for j in range(N+1)]
print(Z[N])

Solution (Python):
n, a, b, c, d = map(int, raw_input().split())
MOD = 10**9+7
fact = [1]*(n+1)
frev = [1]*(n+1)

rcd = range(c, d+1)

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

for i in xrange(1, n+1):
    v = fact[i] = (fact[i-1] * i) % MOD
    frev[i] = pow(v, MOD-2, MOD)
Y = [None]*(n+1)
for g in xrange(a, min(n/c, b)+1):
    p = c*g; q = d*g
    y = frev[g]
    Y[p:min(q,n)+1:g] = (pow(y, j, MOD)*frev[j] for j in xrange(c, min(d, n/g)+1))
    dp[p:] = (dp[i] + frev[n-i]*sum(dp[i-j] * fact[n-i+j] * Y[j] for j in xrange(p, min(q, i)+1, g)) % MOD for i in xrange(p, n+1))
print dp[n]%MOD