Problem:
 Score : 800 points    Problem Statement 12:17 (UTC): The sample input 1 and 2 were swapped. The error is now fixed. We are very sorry for your inconvenience.  There are N children in AtCoder Kindergarten, conveniently numbered 1 through N . Mr. Evi will distribute C indistinguishable candies to the children.  If child i is given a candies, the child's happiness will become x_i^a , where x_i is the child's excitement level . The activity level of the kindergarten is the product of the happiness of all the N children.  For each possible way to distribute C candies to the children by giving zero or more candies to each child, calculate the activity level of the kindergarten. Then, calculate the sum over all possible way to distribute C candies. This sum can be seen as a function of the children's excitement levels x_1,..,x_N , thus we call it f(x_1,..,x_N) .  You are given integers A_i,B_i (1≦i≦N) . Find modulo 10^9+7 .      Constraints  1≦N≦400  1≦C≦400  1≦A_i≦B_i≦400 (1≦i≦N)       Partial Score  400 points will be awarded for passing the test set satisfying A_i=B_i (1≦i≦N) .         Input The input is given from Standard Input in the following format:  N  C  A_1  A_2 ... A_N  B_1  B_2 ... B_N       Output Print the value of modulo 10^9+7 .        Sample Input 1 2 3
1 1
1 1      Sample Output 1 4  This case is included in the test set for the partial score, since A_i=B_i .
We only have to consider the sum of the activity level of the kindergarten where the excitement level of both child 1 and child 2 are 1 ( f(1,1) ).   If child 1 is given 0 candy, and child 2 is given 3 candies, the activity level of the kindergarten is 1^0*1^3=1 .  If child 1 is given 1 candy, and child 2 is given 2 candies, the activity level of the kindergarten is 1^1*1^2=1 .  If child 1 is given 2 candies, and child 2 is given 1 candy, the activity level of the kindergarten is 1^2*1^1=1 .  If child 1 is given 3 candies, and child 2 is given 0 candy, the activity level of the kindergarten is 1^3*1^0=1 .   Thus, f(1,1)=1+1+1+1=4 , and the sum over all f is also 4 .       Sample Input 2 1 2
1
3      Sample Output 2 14  Since there is only one child, child 1 's happiness itself will be the activity level of the kindergarten. Since the only possible way to distribute 2 candies is to give both candies to child 1 , the activity level in this case will become the value of f .   When the excitement level of child 1 is 1 , f(1)=1^2=1 .  When the excitement level of child 1 is 2 , f(2)=2^2=4 .  When the excitement level of child 1 is 3 , f(3)=3^2=9 .   Thus, the answer is 1+4+9=14 .       Sample Input 3 2 3
1 1
2 2      Sample Output 3 66  Since it can be seen that f(1,1)=4 , f(1,2)=15 , f(2,1)=15 , f(2,2)=32 , the answer is 4+15+15+32=66 .       Sample Input 4 4 8
3 1 4 1
3 1 4 1      Sample Output 4 421749  This case is included in the test set for the partial score.       Sample Input 5 3 100
7 6 5
9 9 9      Sample Output 5 139123417

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

mod = 10**9 + 7
N, C = map(int, input().split())
A = list(map(int, input().split()))
B = list(map(int, input().split()))
dp = [0] * (C + 1)
dp[0] = 1
for i in range(N):
    X = [0] * (C + 1)
    for j in range(A[i], B[i] + 1):
        x = 1
        for k in range(C+1):
            X[k] = (X[k] + x) % mod
            x = x * j % mod
    ndp = [0] * (C + 1)
    for j in range(C+1):
        for k in range(C+1-j):
            ndp[j + k] = (ndp[j + k] + dp[j] * X[k]) % mod
    dp = ndp
print(dp[C])

Solution (Python):
from itertools import accumulate
mod = 10 ** 9 + 7

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

memo = []
for p in range(C + 1):
    tmp = [pow(i, p, mod) for i in range(401)]
    memo.append(tuple(accumulate(tmp)))


dp = [[0] * (C + 1) for _ in range(N + 1)]
dp[0][0] = 1
for i, (a, b) in enumerate(zip(A, B)):
    for j in range(C + 1):  
        for k in range(j + 1):  
            dp[i + 1][j] += dp[i][k] * (memo[j-k][b] - memo[j-k][a-1])
            dp[i + 1][j] %= mod

print(dp[N][C])

Solution (Python):


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


def rev(X):
    return pow(X, mod-2, mod)


pow_sum = [[B[i]-A[i]+1 for j in range(C+1)] for i in range(N)]
pow_table = [[1 for j in range(401)] for i in range(401)]
for i in range(401):
    for j in range(401):
        pow_table[i][j] = pow(i, j, mod)

for cnt in range(C):
    for i in range(N):
        pow_sum[i][cnt+1] = 0
        for j in range(A[i], B[i]+1):
            pow_sum[i][cnt+1] += pow_table[j][cnt+1]
        pow_sum[i][cnt+1] %= mod

ans_dp = [[0 for i in range(C+1)] for j in range(N+1)]

ans_dp[0][0] = 1
for n in range(1, N+1):
    
    for c in range(C+1):
        
        for k in range(c+1):
            ans_dp[n][c] += ans_dp[n-1][k]*pow_sum[n-1][c-k]
        ans_dp[n][c] %= mod
print(ans_dp[N][C] % mod)

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

ex400 = [[1] * 401 for _ in range(401)]
for i in range(2,401):
    for j in range(1,401):
        ex400[i][j] = (ex400[i][j-1] * i) % mod

cumsum400 = [[1] * 401 for _ in range(401)]
for i in range(1,401):
    for j in range(401):
        cumsum400[i][j] = (cumsum400[i-1][j] + ex400[i][j]) % mod

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


for i in range(1,n+1):
    for j in range(c+1):
        for k in range(j+1):
            dp[i][j] += dp[i-1][j-k] * (cumsum400[b[i-1]][k] - cumsum400[a[i-1]-1][k])
            dp[i][j] %= mod

print(dp[-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
import pprint
sys.setrecursionlimit(10 ** 9)


INF = 10 ** 20
def LI(): return list(map(int, sys.stdin.buffer.readline().split()))
def I(): return int(sys.stdin.buffer.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, c = LI()
A = LI()
B = LI()

pre_compute = [[0] * 401 for _ in range(401)]
for ii in range(401):
    for jj in range(401):
        pre_compute[ii][jj] = pre_compute[ii - 1][jj]
        pre_compute[ii][jj] += pow(ii, jj, mod)
        pre_compute[ii][jj] %= mod


dp = [[0] * (c + 1) for _ in range(n + 1)]
dp[0][0] = 1
for i in range(n):
    for j in range(c + 1):
        for k in range(j + 1):
            dp[i+1][j] += dp[i][k] * (pre_compute[B[i]][j - k] - pre_compute[A[i] - 1][j - k]) % mod
            
            dp[i+1][j] %= mod



print(dp[n][c])

Solution (Python):
'''
研究室PCでの解答
'''
import math

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


ipt = stdin.readline
setrecursionlimit(10**7)
mod = 10**9+7

def main():
    n,c = map(int,ipt().split())
    a = [int(i) for i in ipt().split()]
    b = [int(i) for i in ipt().split()]
    nj = [[] for j in range(c+1)]
    nj[0] = [1]*401
    for j in range(1,c+1):
        for i in range(401):
            nj[j].append(nj[j-1][i]*i%mod)
    snj = [[0] for j in range(c+1)] 
    for j in range(c+1):
        sj = snj[j]
        ssj = nj[j]
        for i in range(1,401):
            sj.append((sj[-1]+ssj[i])%mod)
    dp = [snj[k][b[0]]-snj[k][a[0]-1] for k in range(c+1)]
    for i in range(1,n):

        ai = a[i]
        bi = b[i]
        for j in range(c,-1,-1):
            dp[j] = dp[j]*(bi-ai+1)%mod
            for k in range(1,j+1):
                dp[j] = (dp[j]+dp[j-k]*(snj[k][bi]-snj[k][ai-1]))%mod

    print(dp[c]%mod)


    return None

if __name__ == '__main__':
    main()

Solution (Python):
N,C=map(int,input().split())
A=list(map(int,input().split()))
B=list(map(int,input().split()))
DP=[[0]*(C+1) for i in range(N+1)]
DP[0][0]=1
mod=10**9+7
X=[[0] for i in range(401)]
for i in range(401):
  for j in range(401):
    X[i].append((X[i][j]+pow(j,i,mod))%mod)
for i in range(N):
  for j in range(C+1):
    for k in range(j+1):
      DP[i+1][j]=(DP[i][j-k]*(X[k][B[i]+1]-X[k][A[i]])+DP[i+1][j])%mod
print(DP[N][C])

Solution (Python):


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

n,c = map(int,readline().split())
*a, = map(int,readline().split())
*b, = map(int,readline().split())

M = 401
MOD = 10**9+7

powsum = []
for i in range(M):
    res = [pow(j,i,MOD) for j in range(M)]
    for j in range(1,M):
        res[j] = (res[j]+res[j-1])%MOD
    powsum.append(res)

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

for ai,bi in zip(a,b):
    ndp = [0]*(c+1)
    for k in range(c+1):
        for j in range(k+1):
            ndp[k] += (powsum[j][bi]-powsum[j][ai-1])*dp[k-j]%MOD
        ndp[k] %= MOD

    dp = ndp

print(dp[c]%MOD)

Solution (Python):
import sys
input = sys.stdin.readline
import numpy as np
MOD = 10**9+7

n, c = map(int, input().split())
A = np.array(tuple(map(int, input().split())))
B = np.array(tuple(map(int, input().split())))
E = np.zeros((n, c+1), dtype=np.int64)
for j in range(c+1):
  cum = np.array(tuple(pow(k, j, MOD) for k in range(401))).cumsum()%MOD
  E[:, j] = cum[B] - cum[A-1]

dp = np.zeros((n+1, c+1), dtype=np.int64)
dp[0, 0] = 1
for i, e in enumerate(E):
  for j, f in enumerate(e):
    dp[i+1, j:] += dp[i, :c+1-j]*f
    dp[i+1] %= MOD
ans = dp[n, c]
print(ans)

Solution (Python):
import sys
input = sys.stdin.readline
import numpy as np
MOD = 10**9+7

n, c = map(int, input().split())
A = np.array(tuple(map(int, input().split())))
B = np.array(tuple(map(int, input().split())))
E = np.zeros((n, c+1), dtype=np.int64)
for j in range(c+1):
  cum = np.array(tuple(pow(k, j, MOD) for k in range(401))).cumsum()%MOD
  E[:, j] = cum[B] - cum[A-1]

dp = np.zeros((n+1, c+1), dtype=np.int64)
dp[0, 0] = 1
for i, e in enumerate(E):
  for j, f in enumerate(e):
    dp[i+1, j:] += dp[i, :c+1-j]*f
    dp[i+1, j:] %= MOD
ans = dp[n, c]
print(ans)

Solution (Python):
import sys
input = sys.stdin.readline
import numpy as np
MOD = 10**9+7

n, c = map(int, input().split())
A = np.array(tuple(map(int, input().split())))
B = np.array(tuple(map(int, input().split())))
E = np.zeros((n, c+1), dtype=np.int64)
for j in range(c+1):
  cum = np.array(tuple(pow(k, j, MOD) for k in range(401))).cumsum()%MOD
  E[:, j] = cum[B] - cum[A-1]

dp = np.zeros((n+1, c+1), dtype=np.int64)
dp[0, 0] = 1
for i, e in enumerate(E):
  for j, f in enumerate(e):
    dp[i+1][j:] += dp[i][:c+1-j]*f
    dp[i+1] %= MOD
ans = dp[n, c]
print(ans)

Solution (Python):
import sys
input = sys.stdin.readline
from itertools import accumulate

MOD = 10**9+7

n, c = map(int, input().split())
A = tuple(map(int, input().split()))
B = tuple(map(int, input().split()))
E = [[0]*(c+1) for _ in range(n)]
for j in range(c+1):
  cumsum = tuple(pow(k, j, MOD) for k in range(401))
  cumsum = tuple(i%MOD for i in accumulate(cumsum))
  for i, (a, b) in enumerate(zip(A, B)):
    E[i][j] = cumsum[b] - cumsum[a-1]

dp = [[0]*(c+1) for _ in range(n+1)]
dp[0][0] = 1
for i, e in enumerate(E):
  for j, f in enumerate(e):
    for k in range(c+1):
      if j+k <= c:
        dp[i+1][j+k] += dp[i][k]*f
        dp[i+1][j+k] %= MOD
ans = dp[n][c]
print(ans)

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

table=[[0 for i in range(401)] for j in range(401)]
for i in range(401):
    S=0
    for j in range(1,401):
        S+=pow(j,i,mod)
        S%=mod
        table[i][j]=S

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


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

for i in range(C+1):
    dp[0][i]=table[i][B[0]]-table[i][A[0]-1]
    dp[0][i]%=mod

for i in range(1,N):
    for j in range(C+1):
        dp[i][j]=sum(dp[i-1][j-k]*(table[k][B[i]]-table[k][A[i]-1]) for k in range(j+1))
        dp[i][j]%=mod

print(dp[N-1][C])

Solution (Python):
def main():
    mod = 10**9+7
    n, c = map(int, input().split())
    A = list(map(int, input().split()))
    B = list(map(int, input().split()))
    dp = [0]*(c+1)
    dp[0] = 1
    ans = 0

    p = [[0]*(401) for _ in [0]*401]  
    for i in range(401):
        for j in range(401):
            p[i][j] = pow(j, i, mod)

    for a, b in zip(A, B):
        dp2 = [0]*(c+1)
        q = [0]*(c+1)
        for i in range(c+1):
            q[i] = sum(p[i][a:b+1]) % mod
        for i in range(c+1):
            temp = 0
            for j in range(i+1):
                temp += dp[i-j]*q[j]
            dp2[i] = temp % mod
        dp = dp2
    ans += dp[-1]
    print(ans % mod)


main()

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

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

ru = [[0] * 410 for i in range(410)]
for cnt in range(410):
    for x in range(409):
        ru[cnt][x + 1] += ru[cnt][x] + pow(x, cnt, MOD)
        ru[cnt][x + 1] %= MOD
      
for ni in range(n):
    for cnt in range(c + 1):
        for ci in range(cnt, c + 1):
            dp[ni + 1][ci] += (ru[cnt][b[ni] + 1] - ru[cnt][a[ni]]) * dp[ni][ci - cnt] 
            dp[ni + 1][ci] %= MOD
print(dp[-1][-1])

Solution (Python):
n,c = map(int,input().split())
a = list(map(int,input().split()))
b = list(map(int,input().split()))
ab = list(zip(a,b))
table = [[0 for i in range(c+1)] for j in range(401)]
mod = 10**9+7
for j in range(c+1):
  if j == 0:
    for i in range(401):
      table[i][j] = 1
  else:
    for i in range(401):
      table[i][j] = table[i][j-1]*i%mod
dp = [[0 for i in range(c+1)] for j in range(n+1)]
for i in range(c+1):
  dp[0][i] = 1
for i in range(1,n+1):
  a,b = ab[i-1]
  for j in range(a,b+1):
    for k in range(c+1):
      dp[i][k] = (dp[i][k]+table[j][k])%mod
ans = [[0 for i in range(c+1)] for j in range(n+1)]
ans[0][0] = 1
for i in range(1,n+1):
  for j in range(c+1):
    for k in range(j+1):
      ans[i][j] = (ans[i][j]+dp[i][k]*ans[i-1][j-k])%mod
print(ans[n][c])

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

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

ruiseki = [[0] * (410) for _ in range(c+1)]
for i in range(c + 1):
    for j in range(409):
        ruiseki[i][j + 1] = ruiseki[i][j] + pow(j, i, MOD)
        ruiseki[i][j + 1] %= MOD
  
for i in range(n):
    for all_num in range(c + 1):
        for j in range(all_num + 1):
            
            dp[i+1][all_num] += (ruiseki[j][b[i]+1] - ruiseki[j][a[i]]) * dp[i][all_num - j]
            dp[i+1][all_num] %= MOD
print(dp[-1][-1])

Solution (Python):
mod = 10**9+7
N, C = map(int, input().split())
A = list(map(int, input().split()))
B = list(map(int, input().split()))
dp = [[0]*(C+1) for i in range(N+1)]
dp[-1][0] = 1
S = [[0]*(C+1) for _ in range(N+1)]
for i in range(N):
    power = [1]*(400+1)
    for j in range(C+1):
        for x in range(A[i], B[i]+1):
            S[i][j]+=power[x]
            S[i][j]%=mod
            power[x]*=x
            power[x]%=mod
for i in range(N):
    for j in range(C+1):
        for k in range(j+1):
            dp[i][j] += dp[i-1][j-k]*S[i][k]
            dp[i][j]%=mod
print(dp[-2][-1])

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


def power_func(a,n,mod):
    bi=str(format(n,"b"))
    res=1
    for i in range(len(bi)):
        res=(res*res) %mod
        if bi[i]=="1":
            res=(res*a) %mod
    return res

def main():
    N, C = map(int, input().split())
    A = list(map(int, input().split()))
    B = list(map(int, input().split()))
    mod = 10**9+7
    L = 400

    dp = [[0]*(L+1) for _ in range(C+1)]
    count = B[0] - A[0] + 1
    for c in range(C+1):
        a, b = A[0], B[0]
        dp[c][0] = power_func(a, c, mod)
        for k in range(1, b-a+1):
            dp[c][k] = (dp[c][k-1] + power_func(a+k, c, mod)) % mod
        dp[c][L] = dp[c][b-a]

    for i, (a, b) in enumerate(zip(A, B)):
        if i == 0: continue
        for k in range(b-a+1):
            dp[0][k] = count * (k+1) % mod
        dp[0][L] = dp[0][b-a]
        for c in range(1, C+1):
            dp[c][0] = (dp[c][L] + a*dp[c-1][0]) % mod
            for k in range(1, b-a+1):
                R = dp[c][k-1] + dp[c][L] + (a+k)*(dp[c-1][k]-dp[c-1][k-1])
                if R < 0: R += mod
                dp[c][k] = R%mod
            dp[c][L] = dp[c][b-a]
        count = (count * (b-a+1)) % mod

    print(dp[C][L])


if __name__ == "__main__":
    main()

Solution (Python):
N, C = map(int, input().split())
A = [int(a) for a in input().split()]
B = [int(a) for a in input().split()]
P = 10**9+7
Y = [[pow(i, j, P) for j in range(401)] for i in range(401)]
for i in range(1, 401):
    for j in range(401):
        Y[i][j] = (Y[i][j] + Y[i-1][j]) % P
X = [[0] * (C+1) for _ in range(N+1)]
X[0][0] = 1
for i in range(1, N+1):
    for j in range(C+1):
        X[i][j] = sum([X[i-1][k] * (Y[B[i-1]][j-k] - Y[A[i-1]-1][j-k]) % P for k in range(j+1)]) % P
print(X[N][C])

Solution (Python):
import sys
from itertools import accumulate
readline = sys.stdin.readline
MOD = 10**9+7

N, C = map(int, readline().split())
A = [0] + list(map(int, readline().split()))
B = [0] + list(map(int, readline().split()))

acsq = []
lb = max(B)+1
for j in range(C+1):
    ac = [0]*lb
    ac[0] = 1
    for i in range(1, lb):
        ac[i] = (ac[i-1] + pow(i, j, MOD)) % MOD
    acsq.append(ac)
dp = [[0]*(C+1) for _ in range(N+1)]
dp[0][0] = 1
for i in range(1, N+1):
    a, b = A[i], B[i]
    for c in range(C+1):
        res = 0
        for j in range(c+1):
            res = (res + (acsq[c-j][b] - acsq[c-j][a-1])*dp[i-1][j]) % MOD
        dp[i][c] = res
print(dp[N][C])

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

acc = [[0] for _ in range(401)]
for i in range(401):
    for j in range(1, 401):
        acc[i].append(((acc[i][-1] + pow(j, i, mod)) % mod))

dp = [[0 for _ in range(C+1)] for _ in range(N+1)]
dp[0][0] = 1
for i in range(N):
    for candy_cur in range(C+1):
        for candy_plus in range(C - candy_cur+1):
            dp[i+1][candy_cur + candy_plus] = (dp[i+1][candy_cur + candy_plus] + (dp[i][candy_cur]\
                                               * ((acc[candy_plus][B[i]] - acc[candy_plus][A[i]-1])%mod)%mod)) % mod

print(dp[-1][C])

Solution (Python):
import sys
input = sys.stdin.buffer.readline
N, C = map(int, input().split())
a = list(map(int, input().split()))
b = list(map(int, input().split()))
c = [0] * (N * (C + 1))
mod = 10 ** 9 + 7
for i in range(N):
  c[i * (C + 1)] = b[i] - a[i] + 1
for i in range(N):
  for k in range(a[i], b[i] + 1):
    t = 1
    for j in range(C):
      t *= k
      t %= mod 
      c[i * (C + 1) + j + 1] += t
      c[i * (C + 1) + j + 1] %= mod
dpx = [0] * (C + 1)
dpy = [0] * (C + 1)
dpx[0] = 1
for i in range(N):
  if i % 2:
    dpx = [0] * (C + 1)
  else:
    dpy = [0] * (C + 1)
  for j in range(C + 1):
    for k in range(j, C + 1):
      if i % 2:
        dpx[k] += dpy[j] * c[i * (C + 1) + k - j]
        dpx[k] %= mod
      else:
        dpy[k] += dpx[j] * c[i * (C + 1) + k - j]
        dpy[k] %= mod
if N % 2:
  print(dpy[-1] % mod)
else:
  print(dpx[-1] % mod)

Solution (Python):
import sys
input = sys.stdin.readline
sys.setrecursionlimit(pow(10, 6))

MOD = pow(10, 9) + 7

def main():
    n, c = map(int, input().split())
    a = list(map(int, input().split()))
    b = list(map(int, input().split()))
    
    pow_table = [[0 for _ in range(401)] for _ in range(401)]
    
    for i in range(1, 401):
        pow_table[i][0] = 1
    
    for i in range(1, 401):
        for j in range(1, 401):
            pow_table[i][j] = pow_table[i][j-1] * i % MOD
    for i in range(1, 401):
        for j in range(401):
            pow_table[i][j] += pow_table[i-1][j]
            if pow_table[i][j] >= MOD:
                pow_table[i][j] -= MOD
    
    dp = [[0 for _ in range(c+1)] for _ in range(n+1)]
    dp[0][0] = 1
    for i in range(1, n+1):
        for j in range(c+1):
            for k in range(j+1):
                ab = pow_table[b[i-1]][k] - pow_table[a[i-1]-1][k]
                if ab < 0:
                    ab += MOD
                dp[i][j] += dp[i-1][j-k] * ab % MOD
                if dp[i][j] >= MOD:
                    dp[i][j] -= MOD

    print(dp[n][c])


if __name__ == '__main__':
    main()

Solution (Python):
def main(N,C,A,B):
    import sys
    input = sys.stdin.readline
    
    
    
    mod = 10**9+7
    powA = [[0]*(401) for i in range(401)]
    for i in range(401):
        for j in range(401):
            powA[i][j] = pow(i,j,mod)
    
    S = [[0]*401 for i in range(C+1)]
    for i in range(C+1):
        S[i][0] = 0
    for i in range(C+1):
        for j in range(1,401):
            S[i][j] = (S[i][j-1] + powA[j][i])%mod
    
    dp = [[0]*(C+1) for i in range(N)]
    for i in range(C+1):
        dp[0][i] = S[i][B[0]] - S[i][A[0]-1]
    for i in range(1,N):
        for j in range(C+1):
            tmp = 0
            for k in range(j+1):
                tmp = (tmp + (S[k][B[i]] - S[k][A[i]-1])*dp[i-1][j-k])%mod
            dp[i][j] = tmp
    print(dp[N-1][C])
    
    
def main2(N,C,A,B):
    import sys
    input = sys.stdin.readline
    
    
    
    mod = 10**9+7
    powA = [[0]*(401) for i in range(401)]
    for i in range(401):
        for j in range(401):
            powA[i][j] = pow(i,j,mod)
    
    S = [[0]*400 for i in range(C+1)]
    for i in range(C+1):
        S[i][0] = 1
    for i in range(C+1):
        for j in range(1,400):
            S[i][j] = (S[i][j-1] + powA[j+1][i])%mod
    
    dp = [[0]*(C+1) for i in range(N)]
    for i in range(C+1):
        dp[0][i] = pow(A[0],i,mod)
    for i in range(1,N):
        for j in range(C+1):
            tmp = 0
            for k in range(j+1):
                tmp = (tmp + powA[A[i]][k]*dp[i-1][j-k])%mod
            dp[i][j] = tmp
    print(dp[N-1][C])
    
import sys
input = sys.stdin.readline
N,C = map(int,input().split())
A = list(map(int,input().split()))
B = list(map(int,input().split()))

if A==B:
    main2(N,C,A,B)
else:
    main(N,C,A,B)

Solution (Python):
import sys
input = sys.stdin.readline
sys.setrecursionlimit(pow(10, 6))

MOD = pow(10, 9)+7

def main():
    n, c = map(int, input().split())
    a = list(map(int, input().split()))
    b = list(map(int, input().split()))
    
    
    
    
    
    
    
    
    
    abc = [[0 for _ in range(401)] for _ in range(401)]
    
    for i in range(401):
        abc[0][i] = 1
    abc[0][0] = 0
    for i in range(1, 401):
        for j in range(1, 401):
            abc[j][i] = abc[j-1][i] * i % MOD
    
    for i in range(1, 401):
        for j in range(401):
            abc[j][i] += abc[j][i-1]
            if abc[j][i] >= MOD:
                abc[j][i] -= MOD
    
    dp = [[0 for _ in range(c+1)] for _ in range(n+1)]
    dp[0][0] = 1
    for i in range(1, n+1):
        for j in range(c+1):
            for k in range(j+1):
                accm = abc[k][b[i-1]] - abc[k][a[i-1]-1]
                if accm < 0:
                    accm += MOD
                dp[i][j] += dp[i-1][j-k] * accm % MOD
                if dp[i][j] >= MOD:
                    dp[i][j] -= MOD
    
    print(dp[n][c])


if __name__ == '__main__':
    main()

Solution (Python):
import os
import sys

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

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


def ModInt(mod):
    class _ModInt:
        def __init__(self, value):
            self.value = value % mod

        def __add__(self, other):
            if isinstance(other, _ModInt):
                return _ModInt(self.value + other.value)
            else:
                return _ModInt(self.value + other)

        def __sub__(self, other):
            if isinstance(other, _ModInt):
                return _ModInt(self.value - other.value)
            else:
                return _ModInt(self.value - other)

        def __radd__(self, other):
            return self.__add__(other)

        def __mul__(self, other):
            if isinstance(other, _ModInt):
                return _ModInt(self.value * other.value)
            else:
                return _ModInt(self.value * other)

        def __truediv__(self, other):
            
            raise NotImplementedError()

        def __repr__(self):
            return str(self.value)

    return _ModInt


MI = ModInt(MOD)

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


def solve():
    
    P = [[1] * (C + 1) for _ in range(max(B) + 1)]
    P[0] = [MI(0)] * (C + 1)
    for i in range(1, len(P)):
        for c in range(1, C + 1):
            P[i][c] = P[i][c - 1] * i

    
    cs = [[0] * (C + 1) for _ in range(max(B) + 1)]
    for c in range(C + 1):
        s = 0
        for i in range(len(P)):
            s += P[i][c]
            cs[i][c] = s

    
    S = [[0] * (C + 1) for _ in range(N)]
    for i in range(N):
        for c in range(C + 1):
            S[i][c] = cs[B[i]][c] - cs[A[i] - 1][c]

    
    dp = S[0][:]
    for i in range(1, N):
        for c in reversed(range(C + 1)):
            s = 0
            for j in range(c + 1):
                s += dp[c - j] * S[i][j]
            dp[c] = s
    return dp[-1]






print(solve())

Solution (Python):

from collections import defaultdict,deque
from heapq import heappush, heappop
import sys
import math
import bisect
import random
def LI(): return [int(x) for x in sys.stdin.readline().split()]
def I(): return int(sys.stdin.readline())
def LS():return [list(x) for x in sys.stdin.readline().split()]
def S():
    res = list(sys.stdin.readline())
    if res[-1] == "\n":
        return res[:-1]
    return res
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)]

sys.setrecursionlimit(1000000)
mod = 1000000007


def A():

    return


def B():

    return


def C():
    n,c = LI()
    a = LI()
    b = LI()
    s = [[0]*(c+1) for i in range(n)]
    for i in range(n):
        ai,bi = a[i],b[i]
        for x in range(ai,bi+1):
            e = 1
            for k in range(c+1):
                s[i][k] += e
                if s[i][k] >= mod:
                    s[i][k] %= mod
                e *= x
                if e >= mod:
                    e %= mod
    dp = [[0]*(c+1) for i in range(n)]
    for k in range(c+1):
        dp[0][k] = s[0][k]

    for i in range(1,n):
        for k in range(c+1):
            for j in range(k+1):
                dp[i][k] += dp[i-1][j]*s[i][k-j]
                if dp[i][k] >= mod:
                    dp[i][k] %= mod
    print(dp[-1][c])
    return


def D():

    return


def E():

    return


def F():

    return


if __name__ == "__main__":
    C()

Solution (Python):
MOD = 10**9 + 7
n, m = [int(item) for item in input().split()]
a = [int(item) for item in input().split()]
b = [int(item) for item in input().split()]

cumsum = [[1] * 410 for _ in range(410)]
for order in range(405):
    for i in range(1, 405):
        cumsum[order][i] = pow(i, order, MOD) + cumsum[order][i-1]
        cumsum[order][i] %= MOD

dp = [[0] * (m + 1) for _ in range(n + 1)]
dp[0][0] = 1
for i in range(1, n+1):
    for j in range(0, m+1):
        for k in range(j+1):
            l = a[i - 1]
            r = b[i - 1]
            x = (cumsum[j-k][r] - cumsum[j-k][l-1] + MOD) % MOD
            dp[i][j] += dp[i-1][k] * x 
        dp[i][j] %= MOD

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

Solution (Python):
import numpy as np
N, C = map(int, input().split())
A = list(map(int, input().split()))
B = list(map(int, input().split()))
mod = 10**9+7
P = np.empty((404, 404), dtype=np.int64)
P[0, :] = 1
ar = np.arange(404, dtype=np.int64)
for i in range(1, 404):
    P[i] = P[i-1] * ar % mod


P = P.cumsum(axis=1, dtype=np.int64) % mod  
P = P.T

dp = np.zeros(C+1, dtype=np.int64)
dp[0] = 1
for a, b in zip(A, B):
    dp_new = np.zeros(C+1, dtype=np.int64)
    p = (P[b] - P[a-1]) % mod
    for c in range(C+1):
        dp_new[c] = (dp[:c+1] * p[c::-1] % mod).sum()
    dp = dp_new % mod
print(dp[C])

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

MOD = 10 ** 9 + 7

N,C = map(int,input().split())
A = [int(x) for x in input().split()]
B = [int(x) for x in input().split()]


kth_pow = np.ones((1024, 401), dtype=np.int64)
rng = np.arange(401, dtype=np.int64)
for i in range(1,C+1):
    kth_pow[i] = kth_pow[i-1] * rng % MOD
kth_pow_cum = kth_pow.cumsum(axis = 1) % MOD
kth_pow_cum[C+1:] = 0

def convolve(A,B,n=32):
    if n == 8:
        return np.rint(np.fft.irfft(np.fft.rfft(A) * np.fft.rfft(B))).astype(np.int64)
    n //= 2
    M = 1 << n
    A1 = A // M
    A2 = A - M * A1
    B1 = A // M
    B2 = B - M * B1
    X = convolve(A1,B1,n)
    Y = convolve(A1-A2,B1-B2,n)
    Z = convolve(A2,B2,n)
    return (X * (M * M % MOD) + (X + Z - Y) * M + Z) % MOD

dp = np.zeros(1024, dtype=np.int64) 
dp[0] = 1
for a,b in zip(A,B):
    arr = kth_pow_cum[:,b] - kth_pow_cum[:,a-1]
    dp = convolve(dp,arr)
    dp[C+1:] = 0

answer = dp[C]
print(answer)

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

MOD = 10 ** 9 + 7

N,C = map(int,input().split())
A = [int(x) for x in input().split()]
B = [int(x) for x in input().split()]


kth_pow = np.ones((C+1, 401), dtype=np.int64)
rng = np.arange(401, dtype=np.int64)
for i in range(1,C+1):
    kth_pow[i] = kth_pow[i-1] * rng % MOD
kth_pow_cum = kth_pow.cumsum(axis = 1) % MOD

dp = np.zeros((C+1), dtype=np.int64) 
dp[0] = 1
for a,b in zip(A,B):
    arr = kth_pow_cum[:,b] - kth_pow_cum[:,a-1]
    prev = dp
    dp = np.zeros(C+1, dtype=np.int64)
    for n in range(C+1):
        dp[n:] += arr[n] * prev[:C+1-n] % MOD
    dp %= MOD

answer = dp[C]
print(answer)

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

sig = [[0 for _ in range(max(b)+1)] for _ in range(c+1)]
for i in range(c+1):
	for j in range(1, max(b)+1):
		sig[i][j] = sig[i][j-1] + pow(j, i, MOD)
		sig[i][j] %= MOD

def sigma(C, A, B):
	return (sig[C][B] - sig[C][A-1]) % MOD

dp = [[0 for _ in range(c+1)] for _ in range(n)]
for j in range(c+1):
	dp[0][j] = sigma(j, a[0], b[0])

for i in range(1, n):
	for j in range(c+1):
		for k in range(j+1):
			dp[i][j] += dp[i-1][k] * sigma(j-k, a[i], b[i])
			dp[i][j] %= MOD

print(dp[n-1][c])

Solution (Python):
from itertools import accumulate
N, C = map(int, input().split())
A = list(map(int, input().split()))
B = list(map(int, input().split()))
mod = 10**9 + 7
X = [list(accumulate([pow(x, d, mod) for x in range(401)])) for d in range(C+1)]
dp = [[0]*(C+1) for _ in range(N+1)]
dp[0][0] = 1
for i in range(1, N+1):
    for d in range(C+1):
        dp[i][d] = sum([dp[i-1][d-k]*(X[k][B[i-1]] - X[k][A[i-1]-1]) for k in range(d+1)])%mod
        
print(dp[-1][-1])