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

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

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

Solution (Python):






import sys
input = sys.stdin.readline

import math

import bisect
import heapq
from collections import defaultdict
from collections import deque
from collections import Counter
from functools import lru_cache





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






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


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


def DD(arg): return defaultdict(arg)


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


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


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


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


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


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


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


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


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]


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




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)

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


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)

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


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)

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


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)