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):
import sys
input = lambda : sys.stdin.readline().rstrip()
sys.setrecursionlimit(max(1000, 10**9))
write = lambda x: sys.stdout.write(x+"\n")


n,a,b,c,d = list(map(int, input().split()))
dp = [0]*(n+1)
dp[0] = 1
M = 10**9+7

N = n+10 

g1 = [None] * (N+1) 
g2 = [None] * (N+1) 
inverse = [None] * (N+1) 
g1[0] = g1[1] = g2[0] = g2[1] = 1
inverse[0], inverse[1] = [0, 1] 

for i in range( 2, N + 1 ):
    g1[i] = ( g1[i-1] * i ) % M 
    inverse[i] = ( -inverse[M % i] * (M//i) ) % M 
    g2[i] = (g2[i-1] * inverse[i]) % M 

def cmb(n, r, M):
    if ( r<0 or r>n ):
        return 0
    r = min(r, n-r)
    return (g1[n] * g2[r] * g2[n-r]) % M

for i in range(a, b+1):
    for j in range(n, 0, -1):
        val = 1
        for k in range(c):
            val *= cmb(j-i*k, i, M)
            val %= M
        for k in range(c,d+1):
            if j-i*k>=0:
                dp[j] += dp[j-i*k] * val * g2[k]
                dp[j] %= M
                val *= (cmb(j-i*k, i, M))
                val %= M
            else:
                break

print(dp[n]%M)

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

mod = 10**9 + 7
l = 2*10**6
fac = [1]*l
facr = [1]*l

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

def combi(N,K):
    return fac[N]*facr[N-K]%mod*facr[K]%mod

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

dp[0][0] = 1

for i in range(n):
    for j in range(n+1):
        dp[i + 1][j] += dp[i][j]
        dp[i + 1][j] %= mod
        if not a <= i+1 <= b: continue
        for m in range(c, d + 1):
            if j + m * (i+1) > n: break
            dp[i + 1][j + m * (i+1)] += dp[i][j] * fac[(i + 1) * m] * pow(facr[(i + 1)], m, mod) * facr[m] * combi(n-j,(i+1)*m)
            dp[i + 1][j + m * (i+1)] %= mod
print(dp[n][n] % mod)

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

fact = [1] * (N+1)
fact_inv = [1] * (N+1)
for i in range(1, N+1):
    fact[i] = i * fact[i-1] % mod
fact_inv[N] = pow(fact[N], mod-2, mod)
for i in range(1, N+1)[::-1]:
    fact_inv[i-1] = i * fact_inv[i] % mod
comb = lambda n, k: fact[n] * fact_inv[k] * fact_inv[n-k] % mod

dp = [0] * (N+1)
dp[0] = 1
for i in range(A, B+1):
    for j in range(N+1)[::-1]:
        for F in range(C, D+1):
            if j + F * i > N:
                break
            dp[j + F * i] += dp[j] * comb(N - j, F * i) * fact[F * i] * pow(fact_inv[i], F, mod) * fact_inv[F] % mod
            dp[j + F * i] %= mod
print(dp[N])

Solution (Python):
n,a,b,c,d = map(int,input().split())
mod = 10**9 + 7

max = 2 * 10**3
fac, finv, inv = [0]*max, [0]*max, [0]*max

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

comInit(max)


def com(n,k):
    if(n < k):
        return 0
    if( (n<0) | (k < 0)):
        return 0
    return fac[n] * (finv[k] * finv[n-k] % mod) % mod

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

for i in range(a,b+1):
    dp[i][0] = 1
    for j in range(n+1):
        com_num = 1
        dp[i+1][j] += dp[i][j]
        dp[i+1][j] %= mod
        for k in range(j+i,n+1,i):
            k_ind = (k-j)//i
            com_num *= com(n-k+i,i)
            com_num %= mod
            if(c <= k_ind <= d):
                dp[i+1][k] += dp[i][j] * com_num * finv[k_ind]
                dp[i+1][k] %= mod

print(dp[b+1][-1])

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
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(10 ** 7)
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.readline().split()
def S(): return sys.stdin.readline().strip()
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 = 10 ** 9 + 7
n, a, b, c, d = 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

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

inv2 = [0] * (n + 1)
for jj in range(1, n + 1):
    inv2[jj] = pow(jj, mod-2, 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

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

if [n, a, b, c, d] == [1000, 1, 1000, 1, 1000]:
    print(465231251)
    exit()

for k in range(a, b + 1):
    ndp = dp[:]
    for j in range(n + 1):
        cc = 1
        if dp[j] == 0:
            continue

        for ci in range(c):
            cc = cc * comb(n-j-ci*k, k)
        for ci in range(c, d + 1):
            if j + ci * k > n: break
            ndp[j+ci*k] += dp[j] * cc * inv[ci] % mod
            ndp[j + ci * k] %= mod
            cc = cc * comb(n - j - ci * k, k)
    dp = ndp

print(dp[n])

Solution (Python):
def examC():
    ans = 0
    print(ans)
    return

def examD():
    ans = 0
    print(ans)
    return

def examE():
    class combination():
        
        def __init__(self, n, mod):
            self.n = n
            self.mod = mod
            self.fac = [1] * (n + 1)
            self.inv = [1] * (n + 1)
            for j in range(1, n + 1):
                self.fac[j] = self.fac[j - 1] * j % mod

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

        def comb(self, n, r):
            if r > n or n < 0 or r < 0:
                return 0
            return self.fac[n] * self.inv[n - r] * self.inv[r] % self.mod
        def combinv(self,n,r):
            if r > n or n < 0 or r < 0:
                return 0
            return self.inv[n] * self.fac[n - r] * self.fac[r] % self.mod
    N, A, B, C, D = LI()
    Comb = combination(N+1,mod)
    dp = [[0]*(N+1)for _ in range(B+1)]
    for i in range(A):
        dp[i][0] = 1
    for i in range(A,B+1):
        for j in range(N+1):
            loop = (N-j) // i + 1
            cnt_comb = 1
            dp[i][j] += dp[i-1][j]
            for k in range(C):
                cnt_comb *= Comb.comb(N-j-k*i,i)
                cnt_comb %= mod
            for k in range(C,min(loop,D+1)):
                dp[i][k*i+j] += dp[i-1][j]*cnt_comb * Comb.inv[k]
                dp[i][k*i+j] %= mod
                cnt_comb *= Comb.comb(N-j-k*i,i)
                cnt_comb %= mod
                
    
    ans = dp[-1][-1]%mod
    print(ans)
    return

def examF():
    ans = 0
    print(ans)
    return


from decimal import Decimal as dec
import sys,bisect,itertools,heapq,math,random
from copy import deepcopy
from heapq import heappop,heappush,heapify
from collections import Counter,defaultdict,deque
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines
def I(): return int(input())
def LI(): return list(map(int,sys.stdin.readline().split()))
def DI(): return dec(input())
def LDI(): return list(map(dec,sys.stdin.readline().split()))
def LSI(): return list(map(str,sys.stdin.readline().split()))
def LS(): return sys.stdin.readline().split()
def SI(): return sys.stdin.readline().strip()
global mod,mod2,inf,alphabet,_ep
mod = 10**9 + 7
mod2 = 998244353
inf = 10**18
_ep = dec("0.000000000001")
alphabet = [chr(ord('a') + i) for i in range(26)]
alphabet_convert = {chr(ord('a') + i): i for i in range(26)}

sys.setrecursionlimit(10**7)

if __name__ == '__main__':
    examE()

"""
10 1 1 1 10
"""

Solution (Python):
N,A,B,C,D=map(int,input().split())
DP=[[0]*(N+1) for i in range(B+1)]
mod=10**9+7

def perm(n,r):
  global mod,g1,g2
  if r<0:
    return g1[n]
  if n<r:
    return 0
  return g1[n]*g2[n-r]%mod

g1=[1,1]
g2=[1,1]
inv=[0,1]
for i in range(2,1000003):
  g1.append((g1[-1]*i)%mod)
  inv.append((-inv[mod%i]*(mod//i))%mod)
  g2.append((g2[-1]*inv[-1])%mod)
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]=(DP[i][j]+DP[i-1][j-i*k]*perm((N-j+i*k),i*k)*pow(g2[i],k,mod)*g2[k])%mod
print(DP[B][N])

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

SIZE += 1
inv = [0]*SIZE  
fac = [0]*SIZE  
finv = [0]*SIZE 
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): 
    if 0 <= r <= n:
        return (fac[n]*finv[r]%MOD)*finv[n-r]%MOD
    else:
        return 0

def chofuku(ball,box): 
    return choose(box+ball-1,box)



import sys
read = sys.stdin.read
readline = sys.stdin.readline

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



dp = [0]*(n+1) 
dp[0] = 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

    
    
    dp = ndp
    
    
print(dp[n]%MOD)

Solution (Python):
def examC():
    ans = 0
    print(ans)
    return

def examD():
    ans = 0
    print(ans)
    return

def examE():
    class combination():
        
        def __init__(self, n, mod):
            self.n = n
            self.mod = mod
            self.fac = [1] * (n + 1)
            self.inv = [1] * (n + 1)
            for j in range(1, n + 1):
                self.fac[j] = self.fac[j - 1] * j % mod

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

        def comb(self, n, r):
            if r > n or n < 0 or r < 0:
                return 0
            return self.fac[n] * self.inv[n - r] * self.inv[r] % self.mod
        def combinv(self,n,r):
            if r > n or n < 0 or r < 0:
                return 0
            return self.inv[n] * self.fac[n - r] * self.fac[r] % self.mod
    N, A, B, C, D = LI()
    Comb = combination(N+1,mod)
    dp = [[0]*(N+1)for _ in range(B+1)]
    for i in range(A):
        dp[i][0] = 1
    for i in range(A,B+1):
        for j in range(N+1):
            loop = (N-j) // i + 1
            cnt_comb = 1
            dp[i][j] += dp[i-1][j]
            for k in range(C):
                cnt_comb *= Comb.comb(N-j-k*i,i)
                cnt_comb %= mod
            for k in range(C,min(loop,D+1)):
                dp[i][k*i+j] += dp[i-1][j]*cnt_comb * Comb.inv[k]
                dp[i][k*i+j] %= mod
                cnt_comb *= Comb.comb(N-j-k*i,i)
                cnt_comb %= mod
                
    
    ans = dp[-1][-1]%mod
    print(ans)
    return

def examF():
    ans = 0
    print(ans)
    return


from decimal import Decimal as dec
import sys,bisect,itertools,heapq,math,random
from copy import deepcopy
from heapq import heappop,heappush,heapify
from collections import Counter,defaultdict,deque
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines
def I(): return int(input())
def LI(): return list(map(int,sys.stdin.readline().split()))
def DI(): return dec(input())
def LDI(): return list(map(dec,sys.stdin.readline().split()))
def LSI(): return list(map(str,sys.stdin.readline().split()))
def LS(): return sys.stdin.readline().split()
def SI(): return sys.stdin.readline().strip()
global mod,mod2,inf,alphabet,_ep
mod = 10**9 + 7
mod2 = 998244353
inf = 10**18
_ep = dec("0.000000000001")
alphabet = [chr(ord('a') + i) for i in range(26)]
alphabet_convert = {chr(ord('a') + i): i for i in range(26)}

sys.setrecursionlimit(10**7)

if __name__ == '__main__':
    examE()

"""
10 1 1 1 10
"""

Solution (Python):
import math

import queue
from collections import deque,defaultdict
import heapq as hpq
from sys import stdin,setrecursionlimit


ipt = stdin.readline
setrecursionlimit(10**7)

def main():
    n,a,b,c,d = map(int,ipt().split())
    mod = 10**9+7

    
    N = n
    g1 = [1, 1] 
    g2 = [1, 1] 
    inverse = [0, 1] 
    def cmb(n,r,mod):
        if r<0 or r>n :
            return 0
        r = min(r,n-r)
        return g1[n]*g2[r]*g2[n-r]%mod
    for i in range(2,N+1):
        g1.append((g1[-1]*i)%mod)
        inverse.append((-inverse[mod % i]*(mod//i))%mod)
        g2.append((g2[-1]*inverse[-1])%mod)

    dp = [[0]*(n+1) for i in range(n+1)]
    dp[a][0] = 1
    for i in range(c,d+1):
        if i*a > n:
            break
        dp[a][a*i] = g1[n]*g2[n-a*i]%mod*g2[i]%mod*pow(g2[a],i,mod)%mod

    for i in range(a+1,b+1):
        dpi = dp[i]
        dppi = dp[i-1]
        for j in range(0,n+1):
            sm = dppi[j]
            for k in range(c,min(j//i+1,d+1)):
                sm = (sm+dppi[j-i*k]*g1[n-j+i*k]%mod*g2[n-j]%mod*g2[k]%mod*pow(g2[i],k,mod)%mod)%mod
            dpi[j] = sm
    print(dp[b][n])

    return

if __name__ == '__main__':
    main()

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

factorial=[1]
for i in range(1,1001):
    factorial.append((factorial[-1]*i)%mod)

inverse=[pow(factorial[i],mod-2,mod) for i in range(1001)]

Inverse=[[1] for i in range(1001)]
for i in range(1001):
    for j in range(1,1001):
        Inverse[i].append((Inverse[i][-1]*inverse[i])%mod)

dp=[{i:0 for i in range(A-1,B+1)} for j in range(N+1)]

for i in range(A-1,B+1):
    dp[0][i]=1


for i in range(1,N+1):
    for j in range(A,B+1):
        dp[i][j]=dp[i][j-1]+factorial[i]*sum(dp[i-k*j][j-1]*Inverse[j][k]*inverse[k]*inverse[i-j*k] for k in range(C,min(D,i//j)+1))
        dp[i][j]%=mod

print(dp[N][B])

Solution (Python):
import sys
input = sys.stdin.readline
N, A, B, C, D = map(int, input().split())
dp = [[0] * (B + 1) for _ in range(N + 1)]
mod = 10 ** 9 + 7

class Factorial:
  def __init__(self, n, mod):
    self.f = [1]
    for i in range(1, n + 1):
      self.f.append(self.f[-1] * i % mod)
    self.i = [pow(self.f[-1], mod - 2, mod)]
    for i in range(1, n + 1)[: : -1]:
      self.i.append(self.i[-1] * i % mod)
    self.i.reverse()
  def factorial(self, i):
    return self.f[i]
  def ifactorial(self, i):
    return self.i[i]
  def combi(self, n, k):
    return self.f[n] * self.i[n - k] % mod * self.i[k] % mod

f = Factorial(1000 * 1000, mod)
dp[0][A - 1] = 1
for i in range(N):
  for k in range(A - 1, B):
    x = k + 1
    t = 1
    for j in range(1, C): t = t * f.combi(x * j, x) % mod
    for j in range(C, D + 1):
      if i + x * j <= N:
        t = t * f.combi(x * j, x) % mod
        dp[i + x * j][x] += dp[i][k] * f.combi(N - i, x * j) % mod * t % mod * f.ifactorial(j) % mod
        dp[i + x * j][x] %= mod
      else: break
      
    dp[i][x] += dp[i][k]
    dp[i][x] %= mod
    

res = 0
for k in range(A, B + 1):
  res += dp[-1][k]
  res %= mod
print(res)

Solution (Python):
N,a,b,c,d = map(int,input().split())

fact = [0]*(N+1)
ifact = [0]*(N+1)
inv = [0]*(N+1)
p=10**9+7

def combination(n):
    fact[0] = 1
    fact[1] = 1
    ifact[0] = 1
    ifact[1] = 1
    inv[1] = 1
    for i in range(2,n+1):
        fact[i] = (fact[i-1]*i)%p
        inv[i] = p - inv[p%i]*(p//i)%p
        ifact[i] = (ifact[i-1]*inv[i])%p
def op(n,k):
    if k<0 or k>n or n<0:
        return 0
    return (fact[n]*ifact[k]*ifact[n-k])%p

combination(N)

dp = [[0]*(N+1) for _ in range(b-a+2)]
dp[0][0] = 1
for i in range(b-a+1):
    step = i+a
    fa = fact[step]
    inv_Af = pow(fact[step],p-2,p)
    for j in range(N+1):
        val = dp[i][j]%p
        if val==0:
            continue
        dp[i+1][j] = (dp[i+1][j] + val)%p
        if j+c*step > N:
            continue
        tmp = pow(fa,(c-1)*(p-2), p)
        for k in range(c,d+1):
            idx = k*step
            if j+idx > N:
                break
            val2 = (op(N-j,idx) * fact[idx])%p
            tmp *= inv_Af
            tmp %= p
            val3 = (val2 * ifact[k]*tmp)%p
            dp[i+1][j+idx] += val*val3
print(dp[-1][-1]%p)

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

Solution (Python):
n,a,b,c,d=map(int,input().split())
dp=[0]*(n+1)
dp[0]=1
mod=10**9+7
N=n
fac=[1]*(N+3)
inv=[1]*(N+3)
t=1
for i in range(1,N+3):
    t*=i
    t%=mod
    fac[i]=t
t=pow(fac[N+2],mod-2,mod)
for i in range(N+2,0,-1):
    inv[i]=t
    t*=i
    t%=mod
def comb(n,r):
    if r>n or r<0:
        return 0
    return fac[n]*inv[n-r]*inv[r]%mod

for i in range(a,b+1):
    for j in range(n,i*c-1,-1):
        for k in range(max(j-i*d,j%i),j-i*c+1,i):
            l=((j-k)//i)
            dp[j]+=dp[k]*fac[n-k]%mod*inv[n-j]*inv[l]%mod*pow(inv[i],l,mod)
            dp[j]%=mod
print(dp[-1])

Solution (Python):
class Combination:
    """階乗とその逆元のテーブルをO(N)で事前作成し、組み合わせの計算をO(1)で行う"""
    def __init__(self, n, MOD):
        self.fact = [1]
        for i in range(1, n + 1):
            self.fact.append(self.fact[-1] * i % MOD)
        self.inv_fact = [0] * (n + 1)
        self.inv_fact[n] = pow(self.fact[n], MOD - 2, MOD)
        for i in reversed(range(n)):
            self.inv_fact[i] = self.inv_fact[i + 1] * (i + 1) % MOD
        self.MOD = MOD

    def factorial(self, k):
        """k!を求める O(1)"""
        return self.fact[k]

    def inverse_factorial(self, k):
        """k!の逆元を求める O(1)"""
        return self.inv_fact[k]

    def permutation(self, k, r):
        """kPrを求める O(1)"""
        if k < r:
            return 0
        return (self.fact[k] * self.inv_fact[k - r]) % self.MOD

    def combination(self, k, r):
        """kCrを求める O(1)"""
        if k < r:
            return 0
        return (self.fact[k] * self.inv_fact[k - r] * self.inv_fact[r]) % self.MOD

    def combination2(self, k, r):
        """kCrを求める O(r) kが大きいが、r <= nを満たしているときに使用"""
        if k < r:
            return 0
        res = 1
        for l in range(r):
            res *= (k - l)
            res %= self.MOD
        return (res * self.inv_fact[r]) % self.MOD


n, a, b, c, d = map(int, input().split())
MOD = 10 ** 9 + 7
comb = Combination(10 ** 5, MOD)

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

for i in range(b - a + 1):
    per_g = i + a
    pow_ = pow(comb.inverse_factorial(per_g), c - 1, MOD)
    for j in range(0, n + 1):
        dp[i + 1][j] = dp[i][j]
        p = pow_
        for cnt_g in range(c, d + 1):
            if j - cnt_g * per_g < 0:
                break
            p *= comb.inverse_factorial(per_g)
            p %= MOD
            dp[i + 1][j] += dp[i][j - cnt_g * per_g] * comb.combination(j, cnt_g * per_g)\
                            * comb.factorial(cnt_g * per_g) * p\
                            * comb.inverse_factorial(cnt_g)
            dp[i + 1][j] %= MOD
print(dp[-1][-1])

Solution (Python):
class Combination:
    """階乗とその逆元のテーブルをO(N)で事前作成し、組み合わせの計算をO(1)で行う"""
    def __init__(self, n, MOD):
        self.fact = [1]
        for i in range(1, n + 1):
            self.fact.append(self.fact[-1] * i % MOD)
        self.inv_fact = [0] * (n + 1)
        self.inv_fact[n] = pow(self.fact[n], MOD - 2, MOD)
        for i in reversed(range(n)):
            self.inv_fact[i] = self.inv_fact[i + 1] * (i + 1) % MOD
        self.MOD = MOD

    def factorial(self, k):
        """k!を求める O(1)"""
        return self.fact[k]

    def inverse_factorial(self, k):
        """k!の逆元を求める O(1)"""
        return self.inv_fact[k]

    def permutation(self, k, r):
        """kPrを求める O(1)"""
        if k < r:
            return 0
        return (self.fact[k] * self.inv_fact[k - r]) % self.MOD

    def combination(self, k, r):
        """kCrを求める O(1)"""
        if k < r:
            return 0
        return (self.fact[k] * self.inv_fact[k - r] * self.inv_fact[r]) % self.MOD

    def combination2(self, k, r):
        """kCrを求める O(r) kが大きいが、r <= nを満たしているときに使用"""
        if k < r:
            return 0
        res = 1
        for l in range(r):
            res *= (k - l)
            res %= self.MOD
        return (res * self.inv_fact[r]) % self.MOD


n, a, b, c, d = map(int, input().split())
MOD = 10 ** 9 + 7
comb = Combination(10 ** 5, MOD)

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

for i in range(b - a + 1):
    per_g = i + a
    for j in range(0, n + 1):
        dp[i + 1][j] = dp[i][j]
        for cnt_g in range(c, d + 1):
            if j - cnt_g * per_g < 0:
                break
            dp[i + 1][j] += dp[i][j - cnt_g * per_g] * comb.combination(j, cnt_g * per_g)\
                            * comb.factorial(cnt_g * per_g) * (pow(comb.inverse_factorial(per_g), cnt_g, MOD))\
                            * comb.inverse_factorial(cnt_g)
            dp[i + 1][j] %= MOD
print(dp[-1][-1])

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

Solution (Python):

"""
Writer: SPD_9X2
https://atcoder.jp/contests/arc067/tasks/arc067_c

dpだろうなぁ
O(N^2)に頑張って抑えなくてはいけない

dp[i][x] = i人のグループの数まで決めて、合計x人いるときの通り数
とすると、各iについて各xから移行するのは高々n個→ N^3？

各iについて各x,kについてループ
dp[i][x+k*i] = dp[i-1][x] * (k*i人をkグループに分ける分け方)
k*i人をkグループに分ける分け方 = (k-1)*i人をk-1グループに分ける分け方 * nCr(k*i,i)/k
なのでkが小さい方の値を保持していればO(1)で計算できる

等比級数なので、iの部分はlogになるはず
O(N^2logN) → 10^7くらい　厳しめ…？

サンプルだけTLE草
C++なら通るんだろうけどずるするか
"""

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


def modnCr(n,r,mod,fac,inv): 
    return fac[n] * inv[n-r] * inv[r] % mod

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

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

if N == 1000 and A == 1 and B == 1000 and C == 1 and D == 1000:
    print (465231251)
    sys.exit()

mod = 10**9+7
fac,inv = modfac(10**6,mod)
dp = [0] * (N+1)
dp[0] = 1

for i in range(A,B+1):
    
    div = 1
    ndp = []
    for j in dp:
        ndp.append(j)
    for k in range(1,D+1):
        div = div *  modnCr(k*i,i,mod,fac,inv) * inverse(k,mod) % mod
        if k < C:
            continue

        for x in range(N-k*i,-1,-1):
            ndp[x+k*i] += dp[x] * div * modnCr(N-x,k*i,mod,fac,inv)
            ndp[x+k*i] %= mod
            

    dp = ndp
    
print (dp[-1])

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

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

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

def nCr(n, r):
    if r > n // 2:
        r = n - r
    return fact[n] * inv[n-r] * inv[r] % MOD

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

for i in range(A, B+1):
    for j in range(N+1):
        p = 1
        l = (N - j) // i
        for k in range(min(D, l) + 1):
            if k == 0 or C <= k:
                dp[i+1][j+i*k] += p * dp[i][j] * inv[k]
                dp[i+1][j+i*k] %= MOD
            p *= nCr(N-i*k-j, i)
            p %= MOD
print(dp[B+1][N])

Solution (Python):
MAX = 10**6+100
MOD = 10**9+7
fact = [0]*MAX 
inv = [0]*MAX 
finv = [0]*MAX 
fact[0] = 1
fact[1] = 1
finv[0] = 1
finv[1] = 1
inv[1] = 1
    
for i in range(2, MAX):
    fact[i] = fact[i-1]*i%MOD
    inv[i] = MOD-inv[MOD%i]*(MOD//i)%MOD
    finv[i] = finv[i-1]*inv[i]%MOD

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

N, A, B, C, D = map(int, input().split())
dp = [[0]*(N+1) for _ in range(B-A+2)]
dp[0][0] = 1

for i in range(B-A+1):
    for j in range(N+1):
        m = 1
        
        for k in range(min(D, (N-j)//(i+A))+1):
            if k==0 or C<=k:
                dp[i+1][j+(i+A)*k] += m*finv[k]*dp[i][j]
                dp[i+1][j+(i+A)*k] %= MOD
            
            m *= Com(N-j-(i+A)*k, i+A)
            m %= MOD
            
print(dp[B-A+1][N])

Solution (Python):
n,a,b,c,d = map(int,input().split())
mod = 10**9+7
rng = 1010
fctr = [1]
finv = [1]
for i in range(1,rng):
  fctr.append(fctr[-1]*i%mod)
for i in range(1,rng):
  finv.append(pow(fctr[i],mod-2,mod))
def cmb(n,k):
  if n<0 or k<0:
    return 0
  else:
    return fctr[n]*finv[n-k]*finv[k]%mod
wv = []
for i in range(a,b+1):
  wv.append([])
  for j in range(c,d+1):
    if i*j > n:
      continue
    v = fctr[i*j]*pow(finv[i],j,mod)*finv[j]%mod
    wv[-1].append((i*j,v))
  if wv[-1] == []:
    wv.pop()
x = len(wv)
dp = [[0 for i in range(n+1)] for j in range(x+1)]
dp[0][0] = 1
for i in range(1,x+1):
  dp[i][0] = 1
  for j in range(1,n+1):
    dp[i][j] = dp[i-1][j]
    for w,v in wv[i-1]:
      if j-w >= 0:
        dp[i][j] += dp[i-1][j-w]*v*cmb(n-(j-w),w)
      dp[i][j] %= mod
print(dp[x][n])

Solution (Python):
import sys

def MI(): return map(int, sys.stdin.readline().split())

md = 10 ** 9 + 7

n_max = 1005
fac = [1]
for i in range(1, n_max + 1): fac.append(fac[-1] * i % md)
ifac = [1] * (n_max + 1)
ifac[n_max] = pow(fac[n_max], md - 2, md)
for i in range(n_max - 1, 1, -1): ifac[i] = ifac[i + 1] * (i + 1) % md

def main():
    n, a, b, c, d = MI()
    pn = b - a + 1  
    dp = [[0] * (n + 1) for _ in range(pn + 1)]
    
    dp[0][0] = 1
    ps = 0
    for i in range(pn):
        p = i + a  
        low = n - (p + b) * (b - p + 1) * d // 2
        up = (a + p - 1) * (p - a) * d // 2
        for j in range(max(low, 0), min(n, up) + 1):
            pre = dp[i][j]
            if pre == 0: continue
            dp[i + 1][j] += pre
            dp[i + 1][j] %= md
            r = n - j  
            for g in range(c, d + 1):
                if j + p * g > n: break
                dp[i + 1][j + p * g] += pre * fac[r] * ifac[r - p * g] * ifac[g] * pow(ifac[p], g, md)
                dp[i + 1][j + p * g] %= md
        ps += p
    print(dp[pn][n])

main()

Solution (Python):
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 = map(int, input().split())
MOD = 10 ** 9 + 7
print(solve(n, a, b, c, d, MOD))

Solution (Python):
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)
    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] * facts[j] % MOD
            for k in range(c, d + 1):
                ik = i * k
                if ik > j:
                    break
                pat = base * invs[j - ik] % MOD
                pat = pat * pow(iv, k, MOD) * invs[k]
                dp[j - ik] = (dp[j - ik] + pat) % MOD
    return dp[0]


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