Problem:


Score: 1000 points

Problem Statement

Find the number of sequences of N non-negative integers A_1, A_2, ..., A_N that satisfy the following conditions:

L \leq A_1 + A_2 + ... + A_N \leq R

When the N elements are sorted in non-increasing order, the M-th and (M+1)-th elements are equal.

Since the answer can be enormous, print it modulo 10^9+7.

Constraints

All values in input are integers.

1 \leq M < N \leq 3 \times 10^5

1 \leq L \leq R \leq 3 \times 10^5

Input

Input is given from Standard Input in the following format:

N M L R

Output

Print the number of sequences of N non-negative integers, modulo 10^9+7.

Sample Input 1

4 2 3 7

Sample Output 1

105

The sequences of non-negative integers that satisfy the conditions are:

\begin{eqnarray} &(1, 1, 1, 0), (1, 1, 1, 1), (2, 1, 1, 0), (2, 1, 1, 1), (2, 2, 2, 0), (2, 2, 2, 1), \\ &(3, 0, 0, 0), (3, 1, 1, 0), (3, 1, 1, 1), (3, 2, 2, 0), (4, 0, 0, 0), (4, 1, 1, 0), \\ &(4, 1, 1, 1), (5, 0, 0, 0), (5, 1, 1, 0), (6, 0, 0, 0), (7, 0, 0, 0)\end{eqnarray}

and their permutations, for a total of 105 sequences.

Sample Input 2

2 1 4 8

Sample Output 2

3

The three sequences that satisfy the conditions are (2, 2), (3, 3), and (4, 4).

Sample Input 3

141592 6535 89793 238462

Sample Output 3

933832916



Sample input:
4 2 3 7


Sample output:
105

Submission status: Wrong Answer

Solution code:
# -*- coding: utf-8 -*-

#############
# Libraries #
#############

import sys
input = sys.stdin.readline

import math
#from math import gcd
import bisect
import heapq
from collections import defaultdict
from collections import deque
from collections import Counter
from functools import lru_cache

#############
# Constants #
#############

MOD = 10**9+7
INF = float('inf')
AZ = "abcdefghijklmnopqrstuvwxyz"

#############
# Functions #
#############

######INPUT######
def I(): return int(input().strip())
def S(): return input().strip()
def IL(): return list(map(int,input().split()))
def SL(): return list(map(str,input().split()))
def ILs(n): return list(int(input()) for _ in range(n))
def SLs(n): return list(input().strip() for _ in range(n))
def ILL(n): return [list(map(int, input().split())) for _ in range(n)]
def SLL(n): return [list(map(str, input().split())) for _ in range(n)]

######OUTPUT######
def P(arg): print(arg); return
def Y(): print("Yes"); return
def N(): print("No"); return
def E(): exit()
def PE(arg): print(arg); exit()
def YE(): print("Yes"); exit()
def NE(): print("No"); exit()

#####Shorten#####
def DD(arg): return defaultdict(arg)

#####Inverse#####
def inv(n): return pow(n, MOD-2, MOD)

######Combination######
kaijo_memo = []
def kaijo(n):
  if(len(kaijo_memo) > n): return kaijo_memo[n]
  if(len(kaijo_memo) == 0): kaijo_memo.append(1)
  while(len(kaijo_memo) <= n): kaijo_memo.append(kaijo_memo[-1] * len(kaijo_memo) % MOD)
  return kaijo_memo[n]

gyaku_kaijo_memo = []
def gyaku_kaijo(n):
  if(len(gyaku_kaijo_memo) > n): return gyaku_kaijo_memo[n]
  if(len(gyaku_kaijo_memo) == 0): gyaku_kaijo_memo.append(1)
  while(len(gyaku_kaijo_memo) <= n): gyaku_kaijo_memo.append(gyaku_kaijo_memo[-1] * pow(len(gyaku_kaijo_memo),MOD-2,MOD) % MOD)
  return gyaku_kaijo_memo[n]

def nCr(n,r):
  if n == r: return 1
  if n < r or r < 0: return 0
  ret = 1
  ret = ret * kaijo(n) % MOD
  ret = ret * gyaku_kaijo(r) % MOD
  ret = ret * gyaku_kaijo(n-r) % MOD
  return ret

######Factorization######
def factorization(n):
  arr = []
  temp = n
  for i in range(2, int(-(-n**0.5//1))+1):
    if temp%i==0:
      cnt=0
      while temp%i==0: 
        cnt+=1 
        temp //= i
      arr.append([i, cnt])
  if temp!=1: arr.append([temp, 1])
  if arr==[]: arr.append([n, 1])
  return arr

#####MakeDivisors######
def make_divisors(n):
  divisors = []
  for i in range(1, int(n**0.5)+1):
    if n % i == 0:
      divisors.append(i)
      if i != n // i: 
        divisors.append(n//i)
  return divisors

#####MakePrimes######
def make_primes(N):
  max = int(math.sqrt(N))
  seachList = [i for i in range(2,N+1)]
  primeNum = []
  while seachList[0] <= max:
    primeNum.append(seachList[0])
    tmp = seachList[0]
    seachList = [i for i in seachList if i % tmp != 0]
  primeNum.extend(seachList)
  return primeNum

#####GCD#####
def gcd(a, b):
    while b: a, b = b, a % b
    return a

#####LCM#####
def lcm(a, b):
    return a * b // gcd (a, b)

#####BitCount#####
def count_bit(n):
  count = 0
  while n:
    n &= n-1
    count += 1
  return count

#####ChangeBase#####
def base_10_to_n(X, n):
  if X//n: return base_10_to_n(X//n, n)+[X%n]
  return [X%n]

def base_n_to_10(X, n):
  return sum(int(str(X)[-i-1])*n**i for i in range(len(str(X))))

def base_10_to_n_without_0(X, n):
  X -= 1
  if X//n: return base_10_to_n_without_0(X//n, n)+[X%n]
  return [X%n]

#####IntLog#####
def int_log(n, a):
  count = 0
  while n>=a:
    n //= a
    count += 1
  return count

#############
# Main Code #
#############
N,M,L,R = IL()
M = max(M,N-M)

@lru_cache(maxsize=None)
def P(t):
  ans = 0
  for k in range(t//M+1):
    u = t-M*(k+1)
    if u>=0:
      if u%(k+1) == 0:
        ans += nCr(N-M,u//(k+1))*(-1)**((u//(k+1))%2)
      if k!=0:
        if u%k == 0:
          ans -= nCr(N-M,u//k)*(-1)**((u//k)%2)
    ans %= MOD
  return ans

@lru_cache(maxsize=None)
def Q(t):
  ans = nCr(t+N,N)
  return ans

def f(MAX):
  ans = 0
  for t in range(MAX+1):
    ans += P(t)*Q(MAX-t)
    ans %= MOD
  return ans

def F(x):
  return (Q(x)-nCr(N,M)*f(x))%MOD
  
print((F(R)-F(L-1))%MOD)

Submission status: Time Limit Exceeded

Solution code:
# -*- coding: utf-8 -*-

#############
# Libraries #
#############

import sys
input = sys.stdin.readline

import math
#from math import gcd
import bisect
import heapq
from collections import defaultdict
from collections import deque
from collections import Counter
from functools import lru_cache

#############
# Constants #
#############

MOD = 10**9+7
INF = float('inf')
AZ = "abcdefghijklmnopqrstuvwxyz"

#############
# Functions #
#############

######INPUT######
def I(): return int(input().strip())
def S(): return input().strip()
def IL(): return list(map(int,input().split()))
def SL(): return list(map(str,input().split()))
def ILs(n): return list(int(input()) for _ in range(n))
def SLs(n): return list(input().strip() for _ in range(n))
def ILL(n): return [list(map(int, input().split())) for _ in range(n)]
def SLL(n): return [list(map(str, input().split())) for _ in range(n)]

######OUTPUT######
def P(arg): print(arg); return
def Y(): print("Yes"); return
def N(): print("No"); return
def E(): exit()
def PE(arg): print(arg); exit()
def YE(): print("Yes"); exit()
def NE(): print("No"); exit()

#####Shorten#####
def DD(arg): return defaultdict(arg)

#####Inverse#####
def inv(n): return pow(n, MOD-2, MOD)

######Combination######
kaijo_memo = []
def kaijo(n):
  if(len(kaijo_memo) > n): return kaijo_memo[n]
  if(len(kaijo_memo) == 0): kaijo_memo.append(1)
  while(len(kaijo_memo) <= n): kaijo_memo.append(kaijo_memo[-1] * len(kaijo_memo) % MOD)
  return kaijo_memo[n]

gyaku_kaijo_memo = []
def gyaku_kaijo(n):
  if(len(gyaku_kaijo_memo) > n): return gyaku_kaijo_memo[n]
  if(len(gyaku_kaijo_memo) == 0): gyaku_kaijo_memo.append(1)
  while(len(gyaku_kaijo_memo) <= n): gyaku_kaijo_memo.append(gyaku_kaijo_memo[-1] * pow(len(gyaku_kaijo_memo),MOD-2,MOD) % MOD)
  return gyaku_kaijo_memo[n]

def nCr(n,r):
  if n == r: return 1
  if n < r or r < 0: return 0
  ret = 1
  ret = ret * kaijo(n) % MOD
  ret = ret * gyaku_kaijo(r) % MOD
  ret = ret * gyaku_kaijo(n-r) % MOD
  return ret

######Factorization######
def factorization(n):
  arr = []
  temp = n
  for i in range(2, int(-(-n**0.5//1))+1):
    if temp%i==0:
      cnt=0
      while temp%i==0: 
        cnt+=1 
        temp //= i
      arr.append([i, cnt])
  if temp!=1: arr.append([temp, 1])
  if arr==[]: arr.append([n, 1])
  return arr

#####MakeDivisors######
def make_divisors(n):
  divisors = []
  for i in range(1, int(n**0.5)+1):
    if n % i == 0:
      divisors.append(i)
      if i != n // i: 
        divisors.append(n//i)
  return divisors

#####MakePrimes######
def make_primes(N):
  max = int(math.sqrt(N))
  seachList = [i for i in range(2,N+1)]
  primeNum = []
  while seachList[0] <= max:
    primeNum.append(seachList[0])
    tmp = seachList[0]
    seachList = [i for i in seachList if i % tmp != 0]
  primeNum.extend(seachList)
  return primeNum

#####GCD#####
def gcd(a, b):
    while b: a, b = b, a % b
    return a

#####LCM#####
def lcm(a, b):
    return a * b // gcd (a, b)

#####BitCount#####
def count_bit(n):
  count = 0
  while n:
    n &= n-1
    count += 1
  return count

#####ChangeBase#####
def base_10_to_n(X, n):
  if X//n: return base_10_to_n(X//n, n)+[X%n]
  return [X%n]

def base_n_to_10(X, n):
  return sum(int(str(X)[-i-1])*n**i for i in range(len(str(X))))

def base_10_to_n_without_0(X, n):
  X -= 1
  if X//n: return base_10_to_n_without_0(X//n, n)+[X%n]
  return [X%n]

#####IntLog#####
def int_log(n, a):
  count = 0
  while n>=a:
    n //= a
    count += 1
  return count

#############
# Main Code #
#############
N,M,L,R = IL()

@lru_cache(maxsize=None)
def P(t):
  ans = 0
  for k in range(t//M+1):
    u = t-M*(k+1)
    if u>=0:
      if u%(k+1) == 0:
        ans += nCr(N-M,u//(k+1))*(-1)**((u//(k+1))%2)
      if k!=0:
        if u%k == 0:
          ans -= nCr(N-M,u//k)*(-1)**((u//k)%2)
    ans %= MOD
  return ans

@lru_cache(maxsize=None)
def Q(t):
  ans = nCr(t+N,N)
  return ans

def f(MAX):
  ans = 0
  for t in range(MAX+1):
    ans += P(t)*Q(MAX-t)
    ans %= MOD
  return ans

def F(x):
  return (Q(x)-nCr(N,M)*f(x))%MOD
  
print((F(R)-F(L-1))%MOD)

Submission status: Time Limit Exceeded

Solution code:
import sys
input = sys.stdin.readline

import numpy as np

N,M,L,R = map(int,input().split())

MOD = 10 ** 9 + 7

X=1000
fact = np.arange(X*X,dtype=np.int64).reshape(X,X)
fact[0,0] = 1
for n in range(1,X):
    fact[:,n] *= fact[:,n-1]
    fact[:,n] %= MOD
for n in range(1,X):
    fact[n] *= fact[n-1,-1]
    fact[n] %= MOD
fact = fact.ravel()

fact_inv = np.arange(1,X*X+1,dtype=np.int64)[::-1].reshape(X,X)
fact_inv[0,0] = pow(int(fact[-1]),MOD-2,MOD)
for n in range(1,X):
    fact_inv[:,n] *= fact_inv[:,n-1]
    fact_inv[:,n] %= MOD
for n in range(1,X):
    fact_inv[n] *= fact_inv[n-1,-1]
    fact_inv[n] %= MOD
fact_inv = fact_inv.ravel()[::-1]

U = N - M
comb = fact[U] * fact_inv[:U+1] % MOD * fact_inv[U::-1] % MOD
comb[1::2] *= (-1)

F = np.zeros(N+M+1,dtype=np.int64)
F[M:N+1] = comb

P = np.zeros(R+1,dtype=np.int64)
for i,x in enumerate(F[1:],1):
    P[::i] += x

P[:M-1:-1] -= P[-M-1::-1]
P %= MOD

# (1-x)^{N+1}で割る
Q = fact[N:N+R+1] * fact_inv[:R+1] % MOD * fact_inv[N] % MOD

def f(R):
    total = fact[R+N] * fact_inv[R] % MOD * fact_inv[N] % MOD
    x = (Q[R::-1] * P[:R+1] % MOD).sum() % MOD
    coef = fact[N] * fact_inv[M] % MOD * fact_inv[N-M] % MOD
    return (total - x*coef)%MOD

answer = (f(R) - f(L-1))%MOD
print(answer)

Submission status: Time Limit Exceeded

Solution code:
import numpy as np

N,M,L,R = map(int,input().split())

MOD = 10 ** 9 + 7

U = 6*10**5 + 100
fact = [1] * (U+1)
for n in range(1,U+1):
    fact[n] = fact[n-1]*n%MOD
fact_inv = [1] * (U+1)
fact_inv[U] = pow(fact[U],MOD-2,MOD)
for n in range(U,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)

U = N - M
comb = fact[U] * fact_inv[:U+1] % MOD * fact_inv[U::-1] % MOD
comb[1::2] *= (-1)

F = np.zeros(N+M+1,dtype=np.int64)
F[M:N+1] = comb

P = np.zeros(R+1,dtype=np.int64)
for i,x in enumerate(F[1:],1):
    P[::i] += x

P[:M-1:-1] -= P[-M-1::-1]
P %= MOD

# (1-x)^{N+1}で割る
Q = fact[N:N+R+1] * fact_inv[:R+1] % MOD * fact_inv[N] % MOD

def f(R):
    total = fact[R+N] * fact_inv[R] % MOD * fact_inv[N] % MOD
    x = (Q[R::-1] * P[:R+1] % MOD).sum() % MOD
    coef = fact[N] * fact_inv[M] % MOD * fact_inv[N-M] % MOD
    return (total - x*coef)%MOD

answer = (f(R) - f(L-1))%MOD
print(answer)

Submission status: Time Limit Exceeded

Solution code:
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10 ** 7)

import numpy as np

N,M,L,R = map(int,input().split())

MOD = 10 ** 9 + 7

U = 6*10**5 + 100
fact = [1] * (U+1)
for n in range(1,U+1):
    fact[n] = fact[n-1]*n%MOD
fact_inv = [1] * (U+1)
fact_inv[U] = pow(fact[U],MOD-2,MOD)
for n in range(U,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)

U = N - M
comb = fact[U] * fact_inv[:U+1] % MOD * fact_inv[U::-1] % MOD
comb[1::2] *= (-1)

F = np.zeros(N+M+1,dtype=np.int64)
F[M:N+1] = comb

P = np.zeros(R+1,dtype=np.int64)
for i,x in enumerate(F[1:],1):
    P[::i] += x

P[M:] -= P.copy()[:-M]
P %= MOD

# (1-x)^{N+1}で割る
Q = fact[N:N+R+1] * fact_inv[:R+1] % MOD * fact_inv[N] % MOD

def f(R):
    total = fact[R+N] * fact_inv[R] % MOD * fact_inv[N] % MOD
    x = (Q[R::-1] * P[:R+1] % MOD).sum() % MOD
    coef = fact[N] * fact_inv[M] % MOD * fact_inv[N-M] % MOD
    return (total - x*coef)%MOD

answer = (f(R) - f(L-1))%MOD
print(answer)

Submission status: Time Limit Exceeded

Solution code:
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10 ** 7)

import numpy as np

N,M,L,R = map(int,input().split())

MOD = 10 ** 9 + 7

U = 6*10**5 + 100
fact = [1] * (U+1)
for n in range(1,U+1):
    fact[n] = fact[n-1]*n%MOD
fact_inv = [1] * (U+1)
fact_inv[U] = pow(fact[U],MOD-2,MOD)
for n in range(U,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)

U = N - M
comb = fact[U] * fact_inv[:U+1] % MOD * fact_inv[U::-1] % MOD
comb[1::2] *= (-1)

F = np.zeros(N+M+1,dtype=np.int64)
F[M:N+1] = comb

P = np.zeros(R+1,dtype=np.int64)
for a in range(1,R+1):
    n = min(len(F),(R+a)//a)
    P[:a*n:a] += F[:n]

P[M:] -= P.copy()[:-M]
P %= MOD

# (1-x)^{N+1}で割る
Q = fact[N:N+R+1] * fact_inv[:R+1] % MOD * fact_inv[N] % MOD

def f(R):
    total = fact[R+N] * fact_inv[R] % MOD * fact_inv[N] % MOD
    x = (Q[R::-1] * P[:R+1] % MOD).sum() % MOD
    coef = fact[N] * fact_inv[M] % MOD * fact_inv[N-M] % MOD
    return (total - x*coef)%MOD

answer = (f(R) - f(L-1))%MOD
print(answer)

Submission status: Accepted

Solution code:
MOD = 10**9+7

kaijo_memo = []
def kaijo(n):
  if(len(kaijo_memo) > n): return kaijo_memo[n]
  if(len(kaijo_memo) == 0): kaijo_memo.append(1)
  while(len(kaijo_memo) <= n): kaijo_memo.append(kaijo_memo[-1] * len(kaijo_memo) % MOD)
  return kaijo_memo[n]

gyaku_kaijo_memo = []
def gyaku_kaijo(n):
  if(len(gyaku_kaijo_memo) > n): return gyaku_kaijo_memo[n]
  if(len(gyaku_kaijo_memo) == 0): gyaku_kaijo_memo.append(1)
  while(len(gyaku_kaijo_memo) <= n): gyaku_kaijo_memo.append(gyaku_kaijo_memo[-1] * pow(len(gyaku_kaijo_memo),MOD-2,MOD) % MOD)
  return gyaku_kaijo_memo[n]

def nCr(n,r):
  if n == r: return 1
  if n < r or r < 0: return 0
  ret = 1
  ret = ret * kaijo(n) % MOD
  ret = ret * gyaku_kaijo(r) % MOD
  ret = ret * gyaku_kaijo(n-r) % MOD
  return ret

N,M,L,R = map(int,input().split())

P = [0 for i in range(R+1)]
Q = [0 for i in range(R+1)]
Z = [0 for i in range(R+1)]
for k in range(1,R+1):
  u = (R-M*k)//k
  for j in range(u+1):
    P[k*j+M*k] += nCr(N-M,j)*(-1)**(j%2)
    P[k*j+M*k] %= MOD
for k in range(R+1):
  Q[k] += nCr(k+N,N)
  Q[k] %= MOD
  Z[k] += nCr(k+N,N)
  Z[k] %= MOD
for k in range(R-M+1):
  Q[k+M] -= nCr(k+N,N)
  Q[k+M] %= MOD

def f(MAX):
  ans = 0
  for t in range(MAX+1):
    ans += P[t]*Q[MAX-t]
    ans %= MOD
  return ans

def F(x):
  return (Z[x]-nCr(N,M)*f(x))%MOD

print((F(R)-F(L-1))%MOD)

Submission status: Accepted

Solution code:
MOD = 10**9+7

kaijo_memo = []
def kaijo(n):
  if(len(kaijo_memo) > n): return kaijo_memo[n]
  if(len(kaijo_memo) == 0): kaijo_memo.append(1)
  while(len(kaijo_memo) <= n): kaijo_memo.append(kaijo_memo[-1] * len(kaijo_memo) % MOD)
  return kaijo_memo[n]

gyaku_kaijo_memo = []
def gyaku_kaijo(n):
  if(len(gyaku_kaijo_memo) > n): return gyaku_kaijo_memo[n]
  if(len(gyaku_kaijo_memo) == 0): gyaku_kaijo_memo.append(1)
  while(len(gyaku_kaijo_memo) <= n): gyaku_kaijo_memo.append(gyaku_kaijo_memo[-1] * pow(len(gyaku_kaijo_memo),MOD-2,MOD) % MOD)
  return gyaku_kaijo_memo[n]

def nCr(n,r):
  if n == r: return 1
  if n < r or r < 0: return 0
  ret = 1
  ret = ret * kaijo(n) % MOD
  ret = ret * gyaku_kaijo(r) % MOD
  ret = ret * gyaku_kaijo(n-r) % MOD
  return ret

N,M,L,R = map(int,input().split())

P = [0 for i in range(R+1)]
for k in range(R+1):
  u = (R-M*(k+1))//(k+1)
  for j in range(u+1):
    P[(k+1)*j+M*(k+1)] += nCr(N-M,j)*(-1)**(j%2)
    P[(k+1)*j+M*(k+1)] %= MOD
for k in range(1,R+1):
  u = (R-M*(k+1))//k
  for j in range(u+1):
    P[k*j+M*(k+1)] -= nCr(N-M,j)*(-1)**(j%2)
    P[k*j+M*(k+1)] %= MOD

Q = [nCr(i+N,N) for i in range(R+1)]

def f(MAX):
  ans = 0
  for t in range(MAX+1):
    ans += P[t]*Q[MAX-t]
    ans %= MOD
  return ans

def F(x):
  return (Q[x]-nCr(N,M)*f(x))%MOD

print((F(R)-F(L-1))%MOD)

Submission status: Accepted

Solution code:
# -*- coding: utf-8 -*-

#############
# Libraries #
#############

import sys
input = sys.stdin.readline

import math
#from math import gcd
import bisect
import heapq
from collections import defaultdict
from collections import deque
from collections import Counter
from functools import lru_cache

#############
# Constants #
#############

MOD = 10**9+7
INF = float('inf')
AZ = "abcdefghijklmnopqrstuvwxyz"

#############
# Functions #
#############

######INPUT######
def I(): return int(input().strip())
def S(): return input().strip()
def IL(): return list(map(int,input().split()))
def SL(): return list(map(str,input().split()))
def ILs(n): return list(int(input()) for _ in range(n))
def SLs(n): return list(input().strip() for _ in range(n))
def ILL(n): return [list(map(int, input().split())) for _ in range(n)]
def SLL(n): return [list(map(str, input().split())) for _ in range(n)]

######OUTPUT######
def P(arg): print(arg); return
def Y(): print("Yes"); return
def N(): print("No"); return
def E(): exit()
def PE(arg): print(arg); exit()
def YE(): print("Yes"); exit()
def NE(): print("No"); exit()

#####Shorten#####
def DD(arg): return defaultdict(arg)

#####Inverse#####
def inv(n): return pow(n, MOD-2, MOD)

######Combination######
kaijo_memo = []
def kaijo(n):
  if(len(kaijo_memo) > n): return kaijo_memo[n]
  if(len(kaijo_memo) == 0): kaijo_memo.append(1)
  while(len(kaijo_memo) <= n): kaijo_memo.append(kaijo_memo[-1] * len(kaijo_memo) % MOD)
  return kaijo_memo[n]

gyaku_kaijo_memo = []
def gyaku_kaijo(n):
  if(len(gyaku_kaijo_memo) > n): return gyaku_kaijo_memo[n]
  if(len(gyaku_kaijo_memo) == 0): gyaku_kaijo_memo.append(1)
  while(len(gyaku_kaijo_memo) <= n): gyaku_kaijo_memo.append(gyaku_kaijo_memo[-1] * pow(len(gyaku_kaijo_memo),MOD-2,MOD) % MOD)
  return gyaku_kaijo_memo[n]

def nCr(n,r):
  if n == r: return 1
  if n < r or r < 0: return 0
  ret = 1
  ret = ret * kaijo(n) % MOD
  ret = ret * gyaku_kaijo(r) % MOD
  ret = ret * gyaku_kaijo(n-r) % MOD
  return ret

######Factorization######
def factorization(n):
  arr = []
  temp = n
  for i in range(2, int(-(-n**0.5//1))+1):
    if temp%i==0:
      cnt=0
      while temp%i==0: 
        cnt+=1 
        temp //= i
      arr.append([i, cnt])
  if temp!=1: arr.append([temp, 1])
  if arr==[]: arr.append([n, 1])
  return arr

#####MakeDivisors######
def make_divisors(n):
  divisors = []
  for i in range(1, int(n**0.5)+1):
    if n % i == 0:
      divisors.append(i)
      if i != n // i: 
        divisors.append(n//i)
  return divisors

#####MakePrimes######
def make_primes(N):
  max = int(math.sqrt(N))
  seachList = [i for i in range(2,N+1)]
  primeNum = []
  while seachList[0] <= max:
    primeNum.append(seachList[0])
    tmp = seachList[0]
    seachList = [i for i in seachList if i % tmp != 0]
  primeNum.extend(seachList)
  return primeNum

#####GCD#####
def gcd(a, b):
    while b: a, b = b, a % b
    return a

#####LCM#####
def lcm(a, b):
    return a * b // gcd (a, b)

#####BitCount#####
def count_bit(n):
  count = 0
  while n:
    n &= n-1
    count += 1
  return count

#####ChangeBase#####
def base_10_to_n(X, n):
  if X//n: return base_10_to_n(X//n, n)+[X%n]
  return [X%n]

def base_n_to_10(X, n):
  return sum(int(str(X)[-i-1])*n**i for i in range(len(str(X))))

def base_10_to_n_without_0(X, n):
  X -= 1
  if X//n: return base_10_to_n_without_0(X//n, n)+[X%n]
  return [X%n]

#####IntLog#####
def int_log(n, a):
  count = 0
  while n>=a:
    n //= a
    count += 1
  return count

#############
# Main Code #
#############
N,M,L,R = IL()

P = [0 for i in range(R+1)]

for k in range(R+1):
  u = (R-M*(k+1))//(k+1)
  for j in range(u+1):
    P[(k+1)*j+M*(k+1)] += nCr(N-M,j)*(-1)**(j%2)
    P[(k+1)*j+M*(k+1)] %= MOD
for k in range(1,R+1):
  u = (R-M*(k+1))//k
  for j in range(u+1):
    P[k*j+M*(k+1)] -= nCr(N-M,j)*(-1)**(j%2)
    P[k*j+M*(k+1)] %= MOD
    
"""
@lru_cache(maxsize=None)
def P(t):
  ans = 0
  for k in range(t+1):
    u = t-M*(k+1)
    if u>=0:
      if u%(k+1) == 0:
        ans += nCr(N-M,u//(k+1))*(-1)**((u//(k+1))%2)
    ans %= MOD
  for k in range(t+1):
    u = t-M*(k+1)
    if u>=0:
      if k!=0:
        if u%k == 0:
          ans -= nCr(N-M,u//k)*(-1)**((u//k)%2)
    ans %= MOD
  return ans
"""

@lru_cache(maxsize=None)
def Q(t):
  ans = nCr(t+N,N)
  return ans

def f(MAX):
  ans = 0
  for t in range(MAX+1):
    ans += P[t]*Q(MAX-t)
    ans %= MOD
  return ans

def F(x):
  return (Q(x)-nCr(N,M)*f(x))%MOD

print((F(R)-F(L-1))%MOD)

Submission status: Accepted

Solution code:
import sys
input = sys.stdin.readline

import numpy as np

N,M,L,R = map(int,input().split())

MOD = 10 ** 9 + 7

X=1000
fact = np.arange(X*X,dtype=np.int64).reshape(X,X)
fact[0,0] = 1
for n in range(1,X):
    fact[:,n] *= fact[:,n-1]
    fact[:,n] %= MOD
for n in range(1,X):
    fact[n] *= fact[n-1,-1]
    fact[n] %= MOD
fact = fact.ravel()

fact_inv = np.arange(1,X*X+1,dtype=np.int64)[::-1].reshape(X,X)
fact_inv[0,0] = pow(int(fact[-1]),MOD-2,MOD)
for n in range(1,X):
    fact_inv[:,n] *= fact_inv[:,n-1]
    fact_inv[:,n] %= MOD
for n in range(1,X):
    fact_inv[n] *= fact_inv[n-1,-1]
    fact_inv[n] %= MOD
fact_inv = fact_inv.ravel()[::-1]

U = N - M
comb = fact[U] * fact_inv[:U+1] % MOD * fact_inv[U::-1] % MOD
comb[1::2] *= (-1)

P = np.zeros(R+1,dtype=np.int64)
for i,x in enumerate(comb,M):
    P[::i] += x

P[M:] -= P.copy()[:-M]
P %= MOD

# (1-x)^{N+1}で割る
Q = fact[N:N+R+1] * fact_inv[:R+1] % MOD * fact_inv[N] % MOD

def f(R):
    total = fact[R+N] * fact_inv[R] % MOD * fact_inv[N] % MOD
    x = (Q[R::-1] * P[:R+1] % MOD).sum() % MOD
    coef = fact[N] * fact_inv[M] % MOD * fact_inv[N-M] % MOD
    return (total - x*coef)%MOD

answer = (f(R) - f(L-1))%MOD
print(answer)

Submission status: Accepted

Solution code:
import sys
input = sys.stdin.readline

import numpy as np

N,M,L,R = map(int,input().split())

MOD = 10 ** 9 + 7

X=1000
fact = np.arange(X*X,dtype=np.int64).reshape(X,X)
fact[0,0] = 1
for n in range(1,X):
    fact[:,n] *= fact[:,n-1]
    fact[:,n] %= MOD
for n in range(1,X):
    fact[n] *= fact[n-1,-1]
    fact[n] %= MOD
fact = fact.ravel()

fact_inv = np.arange(1,X*X+1,dtype=np.int64)[::-1].reshape(X,X)
fact_inv[0,0] = pow(int(fact[-1]),MOD-2,MOD)
for n in range(1,X):
    fact_inv[:,n] *= fact_inv[:,n-1]
    fact_inv[:,n] %= MOD
for n in range(1,X):
    fact_inv[n] *= fact_inv[n-1,-1]
    fact_inv[n] %= MOD
fact_inv = fact_inv.ravel()[::-1]

U = N - M
comb = fact[U] * fact_inv[:U+1] % MOD * fact_inv[U::-1] % MOD
comb[1::2] *= (-1)

P = np.zeros(R+1,dtype=np.int64)
for i,x in enumerate(comb,M):
    P[::i] += x

P[:M-1:-1] -= P[-M-1::-1]
P %= MOD

# (1-x)^{N+1}で割る
Q = fact[N:N+R+1] * fact_inv[:R+1] % MOD * fact_inv[N] % MOD

def f(R):
    total = fact[R+N] * fact_inv[R] % MOD * fact_inv[N] % MOD
    x = (Q[R::-1] * P[:R+1] % MOD).sum() % MOD
    coef = fact[N] * fact_inv[M] % MOD * fact_inv[N-M] % MOD
    return (total - x*coef)%MOD

answer = (f(R) - f(L-1))%MOD
print(answer)

Submission status: Accepted

Solution code:
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10 ** 7)

import numpy as np

N,M,L,R = map(int,input().split())

MOD = 10 ** 9 + 7

U = 6*10**5 + 100
fact = [1] * (U+1)
for n in range(1,U+1):
    fact[n] = fact[n-1]*n%MOD
fact_inv = [1] * (U+1)
fact_inv[U] = pow(fact[U],MOD-2,MOD)
for n in range(U,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)


U = N - M
comb = fact[U] * fact_inv[:U+1] % MOD * fact_inv[U::-1] % MOD
comb[1::2] *= (-1)

P = np.zeros(R+1,dtype=np.int64)
for i,x in enumerate(comb,M):
    P[::i] += x


P[M:] -= P.copy()[:-M]
P %= MOD

# (1-x)^{N+1}で割る
Q = fact[N:N+R+1] * fact_inv[:R+1] % MOD * fact_inv[N] % MOD

def f(R):
    total = fact[R+N] * fact_inv[R] % MOD * fact_inv[N] % MOD
    x = (Q[R::-1] * P[:R+1] % MOD).sum() % MOD
    coef = fact[N] * fact_inv[M] % MOD * fact_inv[N-M] % MOD
    return (total - x*coef)%MOD

answer = (f(R) - f(L-1))%MOD
print(answer)

Submission status: Runtime Error

Solution code:
#include <bits/stdc++.h>
using namespace std;

typedef pair<int, int> pii;
typedef long long ll;
template<class T> using V = vector<T>;
template<class T> using VV = V<V<T>>;

#define pb push_back
#define eb emplace_back
#define mp make_pair
#define fi first
#define se second
#define rep(i,n) rep2(i,0,n)
#define rep2(i,m,n) for(int i=m;i<(n);i++)
#define ALL(c) (c).begin(),(c).end()
#define dump(x) cerr << #x << " = " << (x) << endl
constexpr ll TEN(int n) { return (n == 0) ? 1 : 10 * TEN(n-1); }

template<class T, class U> void chmin(T& t, const U& u) { if (t > u) t = u; }
template<class T, class U> void chmax(T& t, const U& u) { if (t < u) t = u; }

template<class T, class U>
ostream& operator<<(ostream& os, const pair<T, U>& p) {
	os<<"("<<p.first<<","<<p.second<<")";
	return os;
}

template<class T>
ostream& operator<<(ostream& os, const vector<T>& v) {
	os<<"{";
	rep(i, v.size()) {
		if (i) os<<",";
		os<<v[i];
	}
	os<<"}";
	return os;
}

const ll MOD = TEN(9) + 7;
const int MX = TEN(6);

ll inv[MX], fact[MX], ifact[MX];

void init() {
    inv[1] = 1;
    for (int i = 2; i < MX; ++i) {
        inv[i] = inv[MOD % i] * (MOD - MOD / i) % MOD;
    }
    fact[0] = ifact[0] = 1;
    for (int i = 1; i < MX; ++i) {
        fact[i] = fact[i-1] * i % MOD;
        ifact[i] = ifact[i-1] * inv[i] % MOD;
    }
}

ll comb(int n, int r) {
    if (n < 0 || r < 0 || r > n) return 0;
    return fact[n] * ifact[r] % MOD * ifact[n - r] % MOD;
}

ll solve_sub(int N, int M, int S, int r, int l) {
	ll v = 0;
	for (int i = 0; i <= N - M; ++i) {
		ll rem = S - (ll)r * M - (ll)(l + 1) * i;
		if (rem < 0) break;
		ll t = comb(rem + N, N);
		if (i % 2 == 1) {
			t *= -1;
		}

		t = t * comb(N - M, i) % MOD;
		if (t < 0) t += MOD;
		v = (v + t) % MOD;
	}
	v = v * comb(N, M) % MOD;
	return v;
}

ll solve(int N, int M, int S) {
	ll all = comb(S + N, N);
	ll ng = 0;

	for (int x = 1; (ll)x * M <= S; ++x) { //minimum of bigger part
		ll a = solve_sub(N, M, S, x, x-1) - solve_sub(N, M, S, x+1, x-1);
		ng = (ng + a) % MOD;
		if (ng < 0) ng += MOD;
	}

	return (all - ng + MOD) % MOD;
}

int main() {
	int N, M, L, R;
	cin >> N >> M >> L >> R;

	init();

	ll ans = (solve(N, M, R) - solve(N, M, L-1) + MOD) % MOD;
	cout << ans << endl;
	return 0;
}