Problem:
 Score : 600 points    Problem Statement Snuke is standing on a two-dimensional plane. In one operation, he can move by 1 in the positive x -direction, or move by 1 in the positive y -direction.  Let us define a function f(r, c) as follows:   f(r,c) := (The number of paths from the point (0, 0) to the point (r, c) that Snuke can trace by repeating the operation above)   Given are integers r_1 , r_2 , c_1 , and c_2 .
Find the sum of f(i, j) over all pair of integers (i, j) such that r_1 ≤ i ≤ r_2 and c_1 ≤ j ≤ c_2 , and compute this value modulo (10^9+7) .      Constraints  1 ≤ r_1 ≤ r_2 ≤ 10^6  1 ≤ c_1 ≤ c_2 ≤ 10^6  All values in input are integers.         Input Input is given from Standard Input in the following format:  r_1  c_1  r_2  c_2       Output Print the sum of f(i, j) modulo (10^9+7) .        Sample Input 1 1 1 2 2      Sample Output 1 14  For example, there are two paths from the point (0, 0) to the point (1, 1) : (0,0) → (0,1) → (1,1) and (0,0) → (1,0) → (1,1) , so f(1,1)=2 .  Similarly, f(1,2)=3 , f(2,1)=3 , and f(2,2)=6 . Thus, the sum is 14 .       Sample Input 2 314 159 2653 589      Sample Output 2 602215194

Solution (Python):
import sys
mod = pow(10, 9) + 7
sys.setrecursionlimit(pow(10, 8))

def power(x, y):
    if   y == 0: return 1
    elif y == 1     : return x % mod
    elif y % 2 == 0 : return power(x, y//2)**2 % mod
    else: return power(x, (y-1)//2)**2 * x % mod
    
def mul(a, b):
    return ((a % mod) * (b % mod)) % mod

def div(a, b):
    return mul(a, power(b, mod-2))
def div2(a, b):
    return mul(a, modinv(b))

def modinv(a):
    b, u, v = mod, 1, 0
    while b:
        t = a//b
        a, u = a-t*b, u-t*v
        a, b, u, v = b, a, v, u
    u %= mod
    return u

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

NNN = (2*10**6+1)
g1 = [1, 1]
g2 = [1, 1]
inverse = [0, 1]

for i in range( 2, NNN + 1 ):
    g1.append( ( g1[-1] * i ) % mod )
    inverse.append( ( -inverse[mod % i] * (mod//i) ) % mod )
    g2.append( (g2[-1] * inverse[-1]) % mod )

r1, c1, r2, c2 = map(int, input().split())
r = 0
for j in range(c1, c2+1):
    r = (r + cmb(j+r2+1, r2, mod) - cmb(j+r1, r1-1, mod) ) %mod
print(r)

Solution (Python):
a,p,b,q=map(int,input().split())

n=2*(10**6)+3
mod=10**9+7
kaidan=[1]*(n+1)
gyaku=[1]*(n+1)
for i in range(1,n+1):
    kaidan[i]=i*kaidan[i-1]%mod
    gyaku[i]=pow(kaidan[i],mod-2,mod)
def comb(n,r):
    if n-r>=0 and n>=0 and r>=0:
        return (kaidan[n]*gyaku[n-r]*gyaku[r])%mod
    else:
        return 0

def ans(a,b):
    return comb(a+b+2,a+1) -(a+b+2) 

print((ans(b,q)-ans(b,p-1)-ans(a-1,q)+ans(a-1,p-1))%mod)

Solution (Python):
class comb():
  F = [1, 1]
  Fi = [1, 1]
  I = [0, 1]
  def __init__(self, num, mod):
    self.MOD = mod
    for i in range(2, num + 1):
      self.F.append((self.F[-1] * i) % mod)
      self.I.append(mod - self.I[mod % i] * (mod // i) % mod)
      self.Fi.append(self.Fi[-1] * self.I[i] % mod)
  def com(self, n, k):
    if n < k: return 0
    if n < 0 or k < 0: return 0
    return self.F[n] * (self.Fi[k] * self.Fi[n - k] % self.MOD) % self.MOD


r1, c1, r2, c2 = list(map(int, input().split()))
MOD = 10 ** 9 + 7

def nasu(n, m):
  if n < 0 or m < 0: return 0
  ans = (com.com(n + 2, m + 1) - 1)
  return ans
N = r2 + c2
com = comb(N + 3, MOD)

ans = nasu(N, r2)
ans = (ans + nasu(c1 + r1 - 2, min(c1, r1) - 1)) % MOD
ans = (ans - nasu(r2 + c1 - 1, c1 - 1))
ans = (ans - nasu(c2 + r1 - 1, r1 - 1))
while ans < 0: ans += MOD

print(ans)

Solution (Python):
r1, c1, r2, c2 = map(int, input().split())
MOD = 10**9 + 7

def factorial(n, mod=10**9+7):
    a = 1
    for i in range(1,n+1):
        a = a * i % mod
    return a


def power(n, r, mod=10**9+7):
    if r == 0: return 1
    if r%2 == 0:
        return power(n*n % mod, r//2, mod) % mod
    if r%2 == 1:
        return n * power(n, r-1, mod) % mod


def comb(n, k, mod=10**9+7):
    if n < k or k < 0:
        result = 0
    else:
        a = factorial(n, mod=mod)
        b = factorial(k, mod=mod)
        c = factorial(n-k, mod=mod)
        result = (a * power(b, mod-2, mod=mod) * power(c, mod-2, mod=mod)) % mod
    return result


def g(r, c):
  return comb(c+r+2, r+1) - 1

ans = g(r2, c2) - g(r2, c1-1) - g(r1-1, c2) + g(r1-1, c1-1)
ans %= MOD
print(ans)

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

mod = 10**9+7 
N = 2*10**6+2
g1 = [1, 1] 
g2 = [1, 1] 
inverse = [0, 1] 

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 )

r1,c1,r2,c2=map(int,input().split())

print((cmb(r2+c2+2,c2+1,mod)+cmb(r1+c1,r1,mod)-cmb(r1+c2+1,c2+1,mod)-cmb(r2+c1+1,c1,mod))%mod)

Solution (Python):
r1,c1,r2,c2=map(int,input().split())
mod=10**9+7
A=2*(10**6)+3
P=[0 for i in range(A+1)]
P[1]=1
for i in range(2,A):
    P[i]=P[i-1]*i%mod

dr2=pow(P[r2],mod-2,mod)
dr1=pow(P[r1-1],mod-2,mod)
def combr2(x):
    y=r2
    if y>x or x<=0:
        return(0)
    elif y==0:
        return(1)
    else:
        return(P[x]*pow(P[x-y],mod-2,mod)%mod*dr2%mod)

def combr1(x):
    y=r1-1
    if y>x or x<=0:
        return(0)
    elif y==0:
        return(1)
    else:
        return(P[x]*pow(P[x-y],mod-2,mod)%mod*dr1%mod)
count=0
for i in range(c1+1,c2+2):
    count=(count+combr2(i+r2)-combr1(i+r1-1))%mod
print(count%mod)

Solution (Python):
r1,c1,r2,c2=map(int,input().split())
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

mod = 10**9+7 
N = 2*10**6+10
g1 = [1, 1] 
g2 = [1, 1] 
inverse = [0, 1] 

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 )
print((cmb(c2+r2+2,c2+1,mod)-cmb(c2+r1+1,c2+1,mod)-cmb(c1+r2+1,c1,mod)+cmb(c1+r1,c1,mod)+mod)%mod)

Solution (Python):
nn = 2002002
P = 10**9+7
fa = [1] + [0] * nn
for i in range(nn): fa[i+1] = fa[i] * (i+1) % P
f = lambda a, b: fa[a+b+2] * pow(fa[a+1] * fa[b+1], P-2, P) % P - 1
r1, c1, r2, c2 = map(int, input().split())
print((f(r2, c2) - f(r2, c1-1) - f(r1-1, c2) + f(r1-1, c1-1)) % P)

Solution (Python):
r1, c1, r2, c2 = map(int, input().split())

g1 = [1, 1]
g2 = [1, 1]
inverse = [0, 1]
MOD = 10**9+7
for i in range(2, 2*10**6+3):
    g1.append((g1[-1] * i) % MOD)
    inverse.append((-inverse[MOD % i] * (MOD // i)) % MOD)
    g2.append((g2[-1] * inverse[-1]) % MOD)

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

def g(r, c):
  return memo_comb(r+c+2, r+1) - 1

ans = g(r2, c2) - g(r2, c1-1) - g(r1-1, c2) + g(r1-1, c1-1)
ans %= MOD
print(ans)

Solution (Python):
def main():
  mod = 10**9+7
  rng = 2000100
  fctr = [1]
  finv = [1]
  fcappend = fctr.append
  fiappend = finv.append
  for i in range(1,rng):
    fcappend(fctr[-1]*i%mod)
  for i in range(1,rng):
    fiappend(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

  x1,y1,x2,y2 = map(int,input().split())
  print((cmb(x2+y2+2,x2+1)-cmb(x2+y1+1,y1)-cmb(x1+y2+1,x1)+cmb(x1+y1,x1))%mod)
main()

Solution (Python):
MAX = 2*10**6+100
fact = [0]*MAX
inv = [0]*MAX
finv = [0]*MAX

def C_init():
    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 C(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
    
def g(r, c):
    res = 0
    
    for i in range(1, r+2):
        res += C(i+c, c)
        res %= MOD
    
    return res
 
r1, c1, r2, c2 = map(int, input().split())
MOD = 10**9+7
C_init()

print((g(r2, c2)-g(r2, c1-1)-g(r1-1, c2)+g(r1-1, c1-1))%MOD)

Solution (Python):
MOD = 10**9+7
r1, c1, r2, c2 = map(int, input().split())
n = (2 * 10 ** 6) + 2
from math import factorial

class Facts():
    
    
    def __init__(self, max_num=10**5, p=10**9 + 7):
        self.p = p
        self.max_num = max_num
        self.fact = [1] * (self.max_num + 1)
        for i in range(1, self.max_num + 1):
            self.fact[i] = self.fact[i-1] * i
            self.fact[i] %= self.p

    def comb(self, n, k):
        
        if n < 0 or k < 0 or n < k:
            return 0
        if n == 0 or k == 0:
            return 1
        a = self.fact[n]
        b = self.fact[k]
        c = self.fact[n-k]
        return (a*self.power_func(b, self.p-2) *
                self.power_func(c, self.p-2)) % self.p

    def perm(self, n, k):
        
        if n < 0 or k < 0 or n < k:
            return 0
        if n == 0 or k == 0:
            return 1
        a = self.fact[n]
        b = self.fact[n-k]
        return (a * self.power_func(b, self.p-2)) % self.p


    def power_func(self, a, b):
        
        if b == 0:
            return 1
        if b % 2 == 0:
            d = self.power_func(a, b//2)
            return d*d % self.p
        if b % 2 == 1:
            return (a*self.power_func(a, b-1)) % self.p

class FLT:
    
    
    def __init__(self, mod=10**9+7):
        self.mod = mod

    def rep_sqr(self, base, k):
        if k == 0:
            return 1
        elif k % 2 == 0:
            return (self.rep_sqr(base, k // 2) ** 2) % self.mod
        else:
            return (self.rep_sqr(base, k - 1) * base) % self.mod

    def inv(self, a):
        
        return self.rep_sqr(a, self.mod-2)

fact = Facts(n)

def f(x, y):
    return fact.comb(x+y, y)

def g(x, y):
    return f(x+1, y+1)

print((g(r2, c2) - g(r2, c1-1) - g(r1 - 1, c2) + g(r1 - 1, c1 - 1)) % MOD)

Solution (Python):
r1, c1, r2, c2 = map(int, input().split())

g1 = [1, 1]
g2 = [1, 1]
inverse = [0, 1]
MOD = 10**9+7
for i in range(2, 2*10**6+2):
    g1.append((g1[-1] * i) % MOD)
    inverse.append((-inverse[MOD % i] * (MOD // i)) % MOD)
    g2.append((g2[-1] * inverse[-1]) % MOD)

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

def g(r, c):
  res = 0
  for i in range(r+1):
    res += memo_comb(i+c+1, i+1)
    res %= MOD
  
  return res

ans = g(r2, c2) - g(r2, c1-1) - g(r1-1, c2) + g(r1-1, c1-1)
ans %= MOD
print(ans)

Solution (Python):
def C(n, r):
    return fact[n]*inv_table[r]*inv_table[n-r]
    
def g(r, c):
    res = 0
    
    for i in range(1, r+2):
        res += C(i+c, c)
        res %= MOD
    
    return res

r1, c1, r2, c2 = map(int, input().split())
MOD = 10**9+7
fact = [1]

for i in range(1, 2*10**6+100):
    fact.append(fact[-1]*i%MOD)
    
inv_table = [pow(fact[-1], MOD-2, MOD)]

for i in range(1, 2*10**6+100)[::-1]:
    inv_table.append(inv_table[-1]*i%MOD)

inv_table.reverse()

print((g(r2, c2)-g(r2, c1-1)-g(r1-1, c2)+g(r1-1, c1-1))%MOD)

Solution (Python):
class ModComb:

    def __init__(self, bound, mod):
        self.bound, self.mod = bound, mod
        self.g1, self.g2, self.inv = [1, 1], [1, 1], [0, 1]
        for i in range(2, self.bound + 1):
            self.g1.append(self.g1[-1] * i % self.mod)
            self.inv.append(-self.inv[self.mod % i]
                            * (self.mod // i) % self.mod)
            self.g2.append(self.g2[-1] * self.inv[-1] % self.mod)

    def f(self, r, c):
        if r == 0 or c == 0:
            return 0
        tmp = self.g1[r + c + 2] * self.g2[r + 1] * self.g2[c + 1] % self.mod
        return (tmp - (r + c + 2)) % self.mod


r1, c1, r2, c2 = map(int, input().split())
modcomb = ModComb(2 * 10 ** 6 + 2, 10 ** 9 + 7)
print((modcomb.f(r2, c2) - modcomb.f(r1 - 1, c2)
       - modcomb.f(r2, c1 - 1) + modcomb.f(r1 - 1, c1 - 1)) % (10 ** 9 + 7))

Solution (Python):
r1,c1,r2,c2=map(int,input().split())
mod=10**9+7
n_max=2*(10**6+1)
F,FI=[0]*(n_max+1),[0]*(n_max+1)
F[0],FI[0]=1,1
for i in range(n_max):
  F[i+1]=(F[i]*(i+1))%mod
FI[n_max-1]=pow(F[n_max-1],mod-2,mod)
for i in reversed(range(n_max-1)):
  FI[i]=(FI[i+1]*(i+1))%mod
def comb(x,y):
  return (F[x]*FI[x-y]*FI[y])%mod

def f(x,y):
  return comb(x+y+2,y+1)-1

print((f(r2,c2)-f(r2,c1-1)-f(r1-1,c2)+f(r1-1,c1-1))%mod)

Solution (Python):
MOD = 10**9+7
class mint:
    def __init__(self, x):
        self.x = x % MOD if isinstance(x, int) else int(x) % MOD
    def __str__(self):
        return str(self.x)
    __repr__ = __str__
    def __iadd__(self, other):
        self.x += other.x if isinstance(other, mint) else other
        self.x -= MOD if self.x >= MOD else 0
        return self
    def __isub__(self, other):
        self.x += MOD-other.x if isinstance(other, mint) else MOD-other.x
        self.x -= MOD if self.x >= MOD else 0
        return self
    def __imul__(self, other):
        self.x *= other.x if isinstance(other, mint) else other.x
        self.x %= MOD
        return self
    def __imod__(self, other):
        self.x %= other.x if isinstance(other, mint) else other.x
        self.x %= MOD
        return self
    def __add__(self, other):
        return (
            mint(self.x + other.x) if isinstance(other, mint) else
            mint(self.x + other)
        )
    def __sub__(self, other):
        return (
            mint(self.x - other.x) if isinstance(other, mint) else
            mint(self.x - other)
        )
    def __mul__(self, other):
        return (
            mint(self.x * other.x) if isinstance(other, mint) else
            mint(self.x * other)
        )
    def __mod__(self, other):
        return (
            mint(self.x % other.x) if isinstance(other, mint) else
            mint(self.x % other)
        )
    def __truediv__(self, other):
        return (
            mint(
                self.x * pow(other.x, MOD - 2, MOD)
            ) if isinstance(other, mint) else
            mint(self.x * pow(other, MOD - 2, MOD))
        )
    def __pow__(self, other):
        return (
            mint(pow(self.x, other.x, MOD)) if isinstance(other, mint) else
            mint(pow(self.x, other, MOD))
        )
    __radd__ = __add__
    def __rsub__(self, other):
        return (
            mint(other.x - self.x) if isinstance(other, mint) else
            mint(other - self.x)
        )
    __rmul__ = __mul__
    def __rmod__(self, other):
        return (
            mint(other.x % self.x) if isinstance(other, mint) else
            mint(other % self.x)
        )
    def __rtruediv__(self, other):
        return (
            mint(
                other.x * pow(self.x, MOD - 2, MOD)
            ) if isinstance(other, mint) else
            mint(other * pow(self.x, MOD - 2, MOD))
        )
    def __rpow__(self, other):
        return (
            mint(pow(other.x, self.x, MOD)) if isinstance(other, mint) else
            mint(pow(other, self.x, MOD))
        )

class comb():
    
    def __init__(self, n):
        self.n = n
        fact = [0]*(n+1) 
        ifact = [0]*(n+1) 
        fact[0] = mint(1)
        for i in range(1, n+1):
            fact[i] = fact[i-1]*i
        ifact[-1] = pow(fact[-1], MOD-2)
        for i in range(n, 0, -1):
            ifact[i-1] = ifact[i]*i
        self.fact = fact
        self.ifact = ifact
    def calc(self, n, k):
        if k<0 or k>n:
            return 0
        return self.fact[n]*self.ifact[k]*self.ifact[n-k]

def g(r,c,cb):
    res=0
    for i in range(r+1):
        res+=cb.calc(c+i+1,i+1)
    return res

r1,c1,r2,c2 = map(int, input().split())
ans=0
cb=comb(r2+c2+1)

ans+=g(r2,c2,cb)
ans-=g(r1-1,c2,cb)
ans-=g(r2,c1-1,cb)
ans+=g(r1-1,c1-1,cb)
print(ans)

Solution (Python):
r1, c1, r2, c2 = map(int, input().split())

MAX = 2 * 10 ** 6 + 2
MOD = 10 ** 9 + 7


fac = [0] * (MAX + 1)
fac[0] = 1
fac[1] = 1
for i in range(2, MAX + 1):
    fac[i] = fac[i - 1] * i % MOD
    

finv = [0] * (MAX + 1)
finv[MAX] = pow(fac[MAX], MOD - 2, MOD)
for i in reversed(range(1, MAX + 1)):
    finv[i - 1] = finv[i] * i % MOD

def f(r, c):
    return fac[r + c] * finv[r] * finv[c]

ans = 0
for i in range(c1, c2 + 1):
    ans += f(r2, i + 1) - f(r1 - 1, i + 1)
    ans %= MOD
print(ans)

Solution (Python):
MOD = 10**9+7
class mint:
    def __init__(self, x):
        self.x = x % MOD if isinstance(x, int) else int(x) % MOD
    def __str__(self):
        return str(self.x)
    __repr__ = __str__
    def __iadd__(self, other):
        self.x += other.x if isinstance(other, mint) else other
        self.x -= MOD if self.x >= MOD else 0
        return self
    def __isub__(self, other):
        self.x -= MOD-other.x if isinstance(other, mint) else MOD-other.x
        self.x -= MOD if self.x >= MOD else 0
        return self
    def __imul__(self, other):
        self.x *= other.x if isinstance(other, mint) else other.x
        self.x %= MOD
        return self
    def __imod__(self, other):
        self.x %= other.x if isinstance(other, mint) else other.x
        self.x %= MOD
        return self
    def __add__(self, other):
        return (
            mint(self.x + other.x) if isinstance(other, mint) else
            mint(self.x + other)
        )
    def __sub__(self, other):
        return (
            mint(self.x - other.x) if isinstance(other, mint) else
            mint(self.x - other)
        )
    def __mul__(self, other):
        return (
            mint(self.x * other.x) if isinstance(other, mint) else
            mint(self.x * other)
        )
    def __mod__(self, other):
        return (
            mint(self.x % other.x) if isinstance(other, mint) else
            mint(self.x % other)
        )
    def __truediv__(self, other):
        return (
            mint(
                self.x * pow(other.x, MOD - 2, MOD)
            ) if isinstance(other, mint) else
            mint(self.x * pow(other, MOD - 2, MOD))
        )
    def __pow__(self, other):
        return (
            mint(pow(self.x, other.x, MOD)) if isinstance(other, mint) else
            mint(pow(self.x, other, MOD))
        )
    __radd__ = __add__
    def __rsub__(self, other):
        return (
            mint(other.x - self.x) if isinstance(other, mint) else
            mint(other - self.x)
        )
    __rmul__ = __mul__
    def __rmod__(self, other):
        return (
            mint(other.x % self.x) if isinstance(other, mint) else
            mint(other % self.x)
        )
    def __rtruediv__(self, other):
        return (
            mint(
                other.x * pow(self.x, MOD - 2, MOD)
            ) if isinstance(other, mint) else
            mint(other * pow(self.x, MOD - 2, MOD))
        )
    def __rpow__(self, other):
        return (
            mint(pow(other.x, self.x, MOD)) if isinstance(other, mint) else
            mint(pow(other, self.x, MOD))
        )

class comb():
    
    def __init__(self, n):
        self.n = n
        fact = [0]*(n+1) 
        ifact = [0]*(n+1) 
        fact[0] = mint(1)
        for i in range(1, n+1):
            fact[i] = fact[i-1]*i
        ifact[-1] = pow(fact[-1], MOD-2)
        for i in range(n, 0, -1):
            ifact[i-1] = ifact[i]*i
        self.fact = fact
        self.ifact = ifact
    def calc(self, n, k):
        if k<0 or k>n:
            return 0
        return self.fact[n]*self.ifact[k]*self.ifact[n-k]

r1,c1,r2,c2 = map(int, input().split())
ans=0
cb=comb(r2+c2)
acr=[1]
acc=[1]
for i in range(1,r2-r1+1):
    acr.append(acr[-1]+cb.calc(i+c2-c1-1,i))

for i in range(1,c2-c1+1):
    acc.append(acc[-1]+cb.calc(i+r2-r1-1,i))

for r in range(r1+1,r2+1):
    ans+=acr[r2-r1-(r-r1-1)]*cb.calc(c1+r, r)

for c in range(c1+1,c2+1):
    ans+=acc[c2-c1-(c-c1-1)]*cb.calc(c+r1, c)

ans+=cb.calc(c1+r1, c1)
print(ans)

Solution (Python):
r1,c1,r2,c2 = map(int, input().split())

def nCr_mod_p(n, r, p):
    if (r < 0) or (n < r):
        return 0
    r = min(r, n - r)
    return fact[n] * factinv[r] * factinv[n-r] % p

p = 10 ** 9 + 7
N = 2*10 ** 6 + 10  
fact = [1, 1]  
factinv = [1, 1]  
inv = [0, 1]  

for i in range(2, N + 1):
    fact.append((fact[-1] * i) % p)
    inv.append((-inv[p % i] * (p // i)) % p)
    factinv.append((factinv[-1] * inv[-1]) % p)

sum22 = nCr_mod_p(r2+1+c2+1,r2+1,p) - 1
sum21 = nCr_mod_p(r2+1+c1,r2+1,p) - 1
sum12 = nCr_mod_p(r1+c2+1,c2+1,p) - 1
sum11 = nCr_mod_p(r1+c1,c1,p) - 1

ans = sum22 - sum21 - sum12 + sum11
ans %= p

print(ans)

Solution (Python):
r1,c1,r2,c2=map(int,input().split())
mod=10**9+7
n_max=2*(10**6+1)
F,FI=[0]*(n_max+1),[0]*(n_max+1)
F[0],FI[0]=1,1
Inv=[0]*(n_max+1)
Inv[1]=1
for i in range(2,n_max):
  Inv[i]=pow(i,mod-2,mod)


for i in range(n_max):
  F[i+1]=(F[i]*(i+1))%mod
  FI[i+1]=(FI[i]*Inv[i+1])%mod
def comb(x,y):
  return (F[x]*FI[x-y]*FI[y])%mod

def f(x,y):
  return comb(x+y+2,y+1)-1

print((f(r2,c2)-f(r2,c1-1)-f(r1-1,c2)+f(r1-1,c1-1))%mod)

Solution (Python):
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
 
r1,c1,r2,c2=map(int,input().split())
mod=10**9+7

g1=[1,1]
g2=[1,1]
inverse=[0,1]
 
for i in range(2,2*10**6+2):
  g1.append((g1[-1]*i)%mod)
  inverse.append((-inverse[mod%i]*(mod//i))%mod)
  g2.append((g2[-1]*inverse[-1])%mod)

def g(r,c):
  ans=0
  for i in range(1,c+1):
    ans+=cmb(r+i-1,i,mod)
  return ans

print((g(r2+1,c2+1)-g(r2+1,c1)-g(r1,c2+1)+g(r1,c1)+mod)%mod)

Solution (Python):
r1,c1,r2,c2=map(int,input().split())
mod=10**9+7
f=[1]
for i in range(2*10**6+7):f+=[f[-1]*(i+1)%mod]
ans=0
for i in range(c1,c2+1):
  d=((r2+1)*f[i+r2+1]%mod)*f[r1]
  d%=mod
  d=(d+(mod-((r1*f[r2+1]%mod)*f[i+r1])%mod))%mod
  d*=pow((((i+1)*f[i]%mod)*f[r1]%mod)*f[r2+1]%mod,mod-2,mod)
  d%=mod
  ans+=d
  ans%=mod
print(ans)

Solution (Python):
m=10**9+7
def f(r,c):
 s=t=1
 for i in range(c):s=s*(r+c-i)%m;t=t*-~i%m
 return s*pow(t,m-2,m)%m
a,b,c,d=map(int,input().split())
print((f(a,b)-f(c+1,b)-f(a,d+1)+f(c+1,d+1))%m)

Solution (Python):
import sys


TEST_INPUT = [
    """
    1 1 2 2
    """,
    """
    314 159 2653 589
    """
]
ANSWER = [
    "14",
    "602215194"
]


class InputHandler:

    def __init__(self, text_lines="", is_test=False):
        self.data = list(text_lines.split("\n"))
        self.index = 0
        self.is_test = is_test

    def input(self):
        if self.is_test:
            self.index += 1
            return self.data[self.index].strip()
        else:
            return sys.stdin.readline().rstrip()


LARGE = 10 ** 9 + 7


def solve(r1, c1, r2, c2):

    fact = [0] * (r2 + c2 + 3)
    fact[0] = 1
    for i in range(1, r2 + c2 + 3):
        fact[i] = fact[i - 1] * i % LARGE

    fact_inv = [0] * (r2 + c2 + 3)
    fact_inv[-1] = pow(fact[-1], LARGE - 2, LARGE)
    for i in range(r2 + c2 + 1, -1, -1):
        fact_inv[i] = fact_inv[i + 1] * (i + 1) % LARGE

    def solve_special(r, c):
        return (fact[r + c + 2] * fact_inv[r + 1] * fact_inv[c + 1] - 1) % LARGE

    res = solve_special(r2, c2) - solve_special(r1 - 1, c2) - solve_special(r2, c1 - 1) + solve_special(r1 - 1, c1 - 1)

    return res % LARGE


def input_and_solve(ih):
    r1, c1, r2, c2 = map(int, ih.input().split())
    res = solve(r1, c1, r2, c2)
    return res


def main():
    ih = InputHandler()
    res = input_and_solve(ih)
    print(res)


def test():
    for test_input, ans in zip(TEST_INPUT, ANSWER):
        ih = InputHandler(test_input, True)
        res = input_and_solve(ih)
        
        assert str(res) == str(ans)


if __name__ == "__main__":
    
    main()

Solution (Python):
nn = 2002002
P = 10**9+7
fa = [1]
for i in range(1, nn): fa.append(fa[-1] * i % P)
f = lambda a, b: fa[a+b+2] * pow(fa[a+1] * fa[b+1], P-2, P) % P - 1
r1, c1, r2, c2 = map(int, input().split())
print((f(r2, c2) - f(r2, c1-1) - f(r1-1, c2) + f(r1-1, c1-1)) % P)

Solution (Python):
nn = 2002002
P = 10**9+7
fa = [1]
for i in range(1, nn): fa.append(fa[-1] * i % P)
f = lambda a, b: fa[a+b+2] * pow(fa[a+1] * fa[b+1], P-2, P) % P - 1
r1, c1, r2, c2 = map(int, input().split())
print((f(r2, c2) - f(r2, c1-1) - f(r1-1, c2) + f(r1-1, c1-1)) % P)

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


def inv(a, mod):
    r = [1, 0, a]
    w = [0, 1, mod]
    while w[2] != 1:
        q = r[2]//w[2]
        r_new = [r[0]-q*w[0], r[1]-q*w[1], r[2]-q*w[2]]
        r = w
        w = r_new
    x, y = w[0], w[1]    
    return (mod+x % mod) % mod


max_num = 2*10**6+100
fact = [0 for _ in range(max_num)]
ifact = [0 for _ in range(max_num)]

fact[0] = fact[1] = 1
ifact[0] = ifact[1] = 1

for i in range(2, max_num):
    fact[i] = fact[i-1] * i % mod
    ifact[i] = ifact[i-1] * inv(i, mod) % mod


def cmb(x, y):
    return fact[x] * ifact[y] * ifact[x-y] % mod


def solve(r, c):
    ans = 0
    for i in range(r+1):
        ans += cmb(c+i+1, c)
        ans %= mod
    return ans


r1, c1, r2, c2 = map(int, input().split())

ans = solve(r2, c2) - solve(r2, c1-1) - solve(r1-1, c2) + solve(r1-1, c1-1)
ans %= mod

print(ans)

Solution (Python):
r1,c1,r2,c2=map(int,input().split())
a=r2-r1+1
b=c2-c1+1
m=max(a,b)
n=min(a,b)
mod=10**9+7
def framod(n, mod, a=1):
    for i in range(1,n+1):
        a=a * i % mod
    return a
def comb(n, k, mod):
    a=framod(n, mod)
    b=framod(k, mod)
    c=framod(n-k, mod)
    return (a * pow(b, mod-2, mod) * pow(c, mod-2, mod)) % mod
def f(x,y):
  return comb(x+y+2,y+1,mod)-1
 
print((f(r2,c2)-f(r2,c1-1)-f(r1-1,c2)+f(r1-1,c1-1))%mod)

Solution (Python):
MOD = 10**9+7
r1,c1,r2,c2 = map(int,input().split())

size = 1010000*2
f_list = [1]*(size)
for i in range(1,size):
    f_list[i] = (f_list[i-1]*i)%MOD

def func(r,c):
    return (f_list[r+c]*pow(f_list[r],MOD-2,MOD)*pow(f_list[c],MOD-2,MOD))

ans = func(r2+1,c2+1)-1
ans += func(r1,c1)-1
ans -= func(r2+1,c1)-1
ans -= func(r1,c2+1)-1
print(ans%MOD)

Solution (Python):
nn = 2002002
P = 10**9+7
fa = [1]
for i in range(1, nn): fa.append(fa[-1] * i % P)
fainv = [pow(fa[-1], P-2, P)]
for i in range(1, nn)[::-1]: fainv.append(fainv[-1] * i % P)
fainv = fainv[::-1]
f = lambda a, b: fa[a+b+2] * fainv[a+1] * fainv[b+1] % P - 1
r1, c1, r2, c2 = map(int, input().split())
print((f(r2, c2) - f(r2, c1-1) - f(r1-1, c2) + f(r1-1, c1-1)) % P)

Solution (Python):
nn = 2002002
P = 10**9+7
fa = [1]
for i in range(1, nn): fa.append(fa[-1] * i % P)
fainv = [pow(fa[-1], P-2, P)]
for i in range(1, nn)[::-1]: fainv.append(fainv[-1] * i % P)
fainv = fainv[::-1]
f = lambda a, b: fa[a+b+2] * fainv[a+1] * fainv[b+1] % P - 1
r1, c1, r2, c2 = map(int, input().split())
print((f(r2, c2) - f(r2, c1-1) - f(r1-1, c2) + f(r1-1, c1-1)) % P)

Solution (Python):


def pow_mod(x, y):
    MOD = 1000000007
    if y == 0: return 1
    ans = 1
    while y > 1:
        if y % 2 != 0: 
            ans *= x
            ans %= MOD
        x *= x
        x %= MOD
        y //= 2
    return ans * x % MOD

def main():
  MOD = 1000000007
  r1, c1, r2, c2 = map(int, input().split()) 

  points = {r2 + c2 + 2, r2 + 1, c2 + 1, r2 + c1 + 1, c1, r1 + c2 + 1, r1, r1 + c1}

  factorial_map = {0: 1}
  for i in range(1, r2 + c2 + 3):
    factorial_map[i] = (factorial_map[i - 1] * i) % MOD
    if i > 1 and (i - 2) not in points:
      del factorial_map[i-2]

  a = 1
  a *= factorial_map[r2 + c2 + 2]
  a %= MOD
  a *= pow_mod(factorial_map[r2 + 1], MOD - 2)
  a %= MOD
  a *= pow_mod(factorial_map[c2 + 1], MOD - 2)
  a %= MOD

  b = 1
  b *= factorial_map[r2 + c1 + 1]
  b %= MOD
  b *= pow_mod(factorial_map[r2 + 1], MOD - 2)
  b %= MOD
  b *= pow_mod(factorial_map[c1], MOD - 2)
  b %= MOD

  c = 1
  c *= factorial_map[r1 + c2 + 1]
  c %= MOD
  c *= pow_mod(factorial_map[r1], MOD - 2)
  c %= MOD
  c *= pow_mod(factorial_map[c2 + 1], MOD - 2)
  c %= MOD

  d = 1
  d *= factorial_map[r1 + c1]
  d %= MOD
  d *= pow_mod(factorial_map[r1], MOD - 2)
  d %= MOD
  d *= pow_mod(factorial_map[c1], MOD - 2)
  d %= MOD

  ans = a
  ans -= b
  ans %= MOD
  ans -= c 
  ans %= MOD
  ans += d
  ans %= MOD

  print(ans)


if __name__ == "__main__":
  main()

Solution (Python):
import fractions
from functools import reduce

MOD = 10**9 + 7


def lcm_base(x, y):
    return (x * y) // fractions.gcd(x, y)


def lcm(*numbers):
    return reduce(lcm_base, numbers, 1)


def lcm_list(numbers):
    return reduce(lcm_base, numbers, 1)


def xgcd(a, b):
    x0, y0, x1, y1 = 1, 0, 0, 1
    while b != 0:
        q, a, b = a // b, b, a % b
        x0, x1 = x1, x0 - q * x1
        y0, y1 = y1, y0 - q * y1
    return a, x0, y0


def modinv(a, m):
    g, x, y = xgcd(a, m)
    if g != 1:
        raise Exception('modular inverse does not exist')
    else:
        return x % m


def pathes(x, y):
    xf = None
    yf = None
    xyf = None
    temp = 1
    for i in range(1, x + y + 1):
        temp = (temp*i) % MOD
        if i == x:
            xf = temp
        if i == y:
            yf = temp
        if i == x + y:
            xyf = temp

    return xyf * modinv(xf * yf % MOD, MOD)


def main():
    r1, c1, r2, c2 = map(int, input().split())

    a = (pathes(r1, c1) - 1) % MOD
    b = (pathes(r1, c2 + 1) - 1) % MOD
    c = (pathes(r2 + 1, c1) - 1) % MOD
    d = (pathes(r2 + 1, c2 + 1) - 1) % MOD

    print((d - c - b + a) % MOD)


main()

Solution (Python):
a,x,b,y=map(int,input().split());m=10**9+7;N=2000010;F=[1 for i in range(N+1)]
def v(x):return pow(x,m-2,m)
for i in range(N):F[i+1]=-~i*F[i]%m
def f(r,c):return F[r+c]*v(F[r])*v(F[c])%m
print((f(b+1,y+1)-f(b+1,x)-f(a,y+1)+f(a,x))%m)

Solution (Python):
import sys
input = sys.stdin.readline
lr, lc, rr, rc = map(int, input().split())
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(rr + rc + 1, mod)
res = 0
for i in range(lr, rr + 1):
  res += f.combi(i + rc + 1, i + 1)
  res %= mod
for i in range(lr, rr + 1):
  res -= f.combi(i + lc, i + 1)
  res %= mod
print(res)

Solution (Python):

r1,c1,r2,c2 = map(int,input().split())

mod = 10 ** 9 + 7
MAX = (10**6)*2 + 100
fac = [0] * MAX
finv = [0] * MAX
inv = [0] * MAX

def cominit():
    fac[0] = 1
    fac[1] = 1
    finv[0] = 1
    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



def cmb(n,r):
    if n < 0 or r < 0 or r > n:return 0
    if r > n/2: r = n-r        
    return fac[n] * (finv[r] * finv[n-r] % mod) % mod
    
cominit()

ans = 0
for c in range(c1, c2+1):
    ans += ((r2+1)*fac[r1]*fac[c+r2+1]-r1*fac[r2+1]*fac[c+r1])*finv[c]*finv[r1]*finv[r2+1]*inv[c+1]
    ans %= mod
    
print(ans)

Solution (Python):
r1,c1,r2,c2=map(int,input().split());mod=10**9+7;N=2000010
def inv(x):return pow(x,mod-2,mod)
Fact=[0 for i in range(N+1)];Fact[0]=1
for i in range(N):Fact[i+1]=((i+1)*Fact[i])%mod
def f(r,c):return (Fact[r+c]*inv(Fact[r])*inv(Fact[c]))%mod
print((f(r2+1,c2+1)-f(r2+1,c1)-f(r1,c2+1)+f(r1,c1))%mod)