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-----
3 1 3 1 2

-----Sample Output-----
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.

Difficulty: competition

Starter code:

Solutions:
def main():
    mod = 10**9+7
    inv_n = [0]*1001
    nCr = [[1]*(i+1) for i in range(1001)]
    for i in range(1001):
        inv_n[i] = pow(i, mod-2, mod)
    for i in range(2, 1001):
        for j in range(1, i):
            nCr[i][j] = (nCr[i-1][j-1]+nCr[i-1][j]) % mod
    n, a, b, c, d = list(map(int, input().split()))
    dp = [0]*(n+1)
    dp[0] = 1

    for A in range(b, a-1, -1):
        dp2 = [i for i in dp]
        for N in range(n-c*A, -1, -1):
            e = dp[N]
            if e:
                temp = 1
                for C in range(1, c):
                    temp = temp*nCr[n-N-(C-1)*A][A]*inv_n[C] % mod
                for C in range(c, min(d, (n-N)//A)+1):
                    temp = temp*nCr[n-N-(C-1)*A][A]*inv_n[C] % mod
                    dp2[N+C*A] = (dp2[N+C*A]+temp*e) % mod
        dp = dp2
    print((dp[-1]))


main()


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(list(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(list(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(list(range(D))):
            invss[i][k] = (invss[i][k + 1] * facts[i]) % MOD
    return invss


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

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

# dp[i][j]: i人以下のグループのみで、j人使っている場合の数
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]))


def main():
    mod = 10**9+7
    inv_n = [0]*1001
    nCr = [[1]*(i+1) for i in range(1001)]
    for i in range(1001):
        inv_n[i] = pow(i, mod-2, mod)
    for i in range(2, 1001):
        for j in range(1, i):
            nCr[i][j] = (nCr[i-1][j-1]+nCr[i-1][j]) % mod
    n, a, b, c, d = map(int, input().split())
    dp = [0]*(n+1)
    dp[0] = 1

    for A in range(b, a-1, -1):
        dp2 = [i for i in dp]
        for N in range(n-c*A, -1, -1):
            e = dp[N]
            if e:
                temp = 1
                for C in range(1, c):
                    temp = temp*nCr[n-N-(C-1)*A][A]*inv_n[C] % mod
                for C in range(c, min(d, (n-N)//A)+1):
                    temp = temp*nCr[n-N-(C-1)*A][A]*inv_n[C] % mod
                    dp2[N+C*A] = (dp2[N+C*A]+temp*e) % mod
        dp = dp2
    print(dp[-1])


main()

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)

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

def prepare(n, MOD):
    f = 1
    factorials = [1]
    for m in range(1, n + 1):
        f = f * m % 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 = inv * m % MOD
        invs[m - 1] = inv

    return factorials, invs


def solve(n, a, b, c, d, MOD):
    facts, invs = prepare(n, MOD)

    pre = {}
    for i in range(a, b + 1):
        iv = invs[i]
        for k in range(c, d + 1):
            if i * k > n:
                break
            pre[i, k] = pow(iv, k, MOD) * invs[k] % MOD

    dp = [0] * (n + 1)
    dp[n] = 1
    for i in range(a, b + 1):
        iv = invs[i]
        for j in range(i * c, n + 1):
            base = dp[j]
            for k in range(c, d + 1):
                ik = i * k
                if ik > j:
                    break
                dp[j - ik] = (dp[j - ik] + base * pre[i, k]) % MOD
    return dp[0] * facts[n] % MOD


n, a, b, c, d = list(map(int, input().split()))
MOD = 10 ** 9 + 7
print((solve(n, a, b, c, d, MOD)))


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] #便宜上0の逆元は0にする
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)

mod=10**9+7
MAX=10**3+100

#g1は階乗、g2はその逆元
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=list(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))



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

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


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

MAX = 10**3+1
inv_t = [0,1]
for i in range(2,MAX+1):
  inv_t.append(inv_t[MOD % i] * (MOD - int(MOD / i)) % MOD)
fact_inv = [1,1]
for i in range(2,MAX+1):
  fact_inv.append(fact_inv[-1]*inv_t[i] % MOD)
  
DP = [[0 for _ in range(N+1)] for _ in range(B+1)]
DP[0][0] = 1

#配るDP
div_fct = fact_inv[C] * pow(fact_inv[A],C,MOD)
DP[A][0] = 1
for c in range(C,D+1):
  if c*A > N: break
  DP[A][c*A] += div_fct
  DP[A][c*A] %= MOD
  div_fct *= fact_inv[A]*inv_t[c+1]
  div_fct %= MOD

for E in range(A,B):
  div_fct_0 = fact_inv[C] * pow(fact_inv[E+1],C,MOD)
  for n in range(N+1):
    now = DP[E][n]
    #print(E,n,now)
    DP[E+1][n] += DP[E][n]
    if now == 0: continue
    div_fct = div_fct_0
    for c in range(C,D+1):#配るDP
      nxt = n + c*(E+1)
      if nxt > N: break
      DP[E+1][nxt] += now * div_fct
      DP[E+1][nxt] %= MOD
      div_fct *= fact_inv[E+1]*inv_t[c+1]
      div_fct %= MOD

fact_N = 1
for i in range(1,N+1):
  fact_N *= i
  fact_N %= MOD
#print(*DP, sep="\n")  
print((DP[B][N]*fact_N % MOD))
  


N,A,B,C,D=list(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))


        
  


    


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):
    # 使うなら、C～D人で使う
    prev = dp
    dp = prev.copy()
    for n in range(C,D+1):
        y = n * x
        if y > N:
            break
        # x人組をn組とる
        # dp[i] += dp[i-y] * comb((N-i+y),y) * (y! / (x!)^n / n!)
        coef = fact[y] * pow(int(fact_inv[x]), n, MOD) % MOD * fact_inv[n] % MOD
        #print(x,n,coef)
        dp[y:] += prev[:-y] * comb[N:y-1:-1,y] % MOD * coef
        dp %= MOD

answer = dp[N]
print(answer)

SIZE=5*10**3; 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
fac[0] = fac[1] = 1
finv[0] = finv[1] = 1
for i in range(2,SIZE):
    fac[i] = fac[i-1]*i%MOD
finv[-1] = pow(fac[-1],MOD-2,MOD)
for i in range(SIZE-1,0,-1):
    finv[i-1] = finv[i]*i%MOD
    inv[i] = finv[i]*fac[i-1]%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

def chofuku(ball,box): # nHk mod MOD の計算
    return choose(box+ball-1,box)

# coding: utf-8
# Your code here!
import sys
read = sys.stdin.read
readline = sys.stdin.readline

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



dp = [0]*(n+1) # i まで見て、x人グループが決まっている
dp[0] = 1


#for i in range(b,a-1,-1):
for i in range(a,b+1):
    invmod = [1]
    for p in range(1001//i+2): invmod.append(invmod[-1]*finv[i]%MOD)
    coeff = [fac[i*p]*invmod[p]%MOD*finv[p]%MOD for p in range(1001//i+2)]    
    
    ndp = dp[:]
    for j in range(n+1):
        for p in range(c,d+1):
            if j+i*p > n: break
            ndp[j+i*p] += coeff[p]*dp[j]%MOD*choose(n-j,i*p)
            ndp[j+i*p] %= MOD

    #print(invmod,i)
    #print(coeff)
    dp = ndp
    #if i < 100:break
    #print(dp,i)    
print((dp[n]%MOD))