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
readline = sys.stdin.buffer.readline
r1,c1,r2,c2 = map(int,readline().split())
mod = 10**9+7
def pow(n,p,mod=10**9+7): 
    res = 1
    while p > 0:
        if p % 2 == 0:
            n = n ** 2 % mod
            p //= 2
        else:
            res = res * n % mod
            p -= 1
    return res % mod
def factrial_memo(n=3*10**6+1,mod=10**9+7):
    fact = [1, 1]
    for i in range(2, n + 1):
        fact.append((fact[-1] * i) % mod)
    return fact
def fermat_cmb(n, r, mod=10**9+7): 
    return fact[n] * pow(fact[r],mod-2) * pow(fact[n-r],mod-2) %mod
fact = factrial_memo()
r1,c1,r2,c2 = r1+1,c1+1,r2+1,c2+1
al = fermat_cmb(r2+c2,r2)
sub_r = fermat_cmb(r1-1+c2,r1-1)
sub_c = fermat_cmb(r2+c1-1,c1-1)
add_rc = fermat_cmb(r1-1+c1-1,r1-1)

al -= sub_r+sub_c
al += add_rc

print(al%mod)

Solution (Python):
r,c,R,C=map(int,input().split())
m=10**9+7
max_n=2*10**6+1
fac=[1]*(max_n+1)
inv=[1]*(max_n+1)
ifac=[1]*(max_n+1)
for n in range(2,max_n+1):
    fac[n]=(fac[n-1]*n)%m
    inv[n]=m-inv[m%n]*(m//n)%m
    ifac[n]=(ifac[n-1]*inv[n])%m
def comb(n,k):
    if n<k or n<0 or k<0:
        return 0
    return (fac[n]*ifac[k]*ifac[n-k])%m
def f(r,c):
    return sum([comb(i+1+c,c) for i in range(r+1)])%m
print((f(R,C)-f(R,c-1)-f(r-1,C)+f(r-1,c-1))%m)

Solution (Python):
import sys
import time

def getN():
    return int(input())
def getNM():
    return map(int, input().split())
def getlist():
    return list(map(int, input().split()))
import math
import heapq
import bisect
from collections import defaultdict, Counter, deque
MOD = 10**9 + 7
INF = 10**21

def routes_from_zero(r,c, factorial, gyaku):
    ret = 0
    for i in range(r + 1):
        tmp = 1
        tmp *= factorial[i+c+1]
        
        tmp *= gyaku[c]
        
        tmp *= gyaku[i+1]
        

        ret += tmp
        ret %= MOD

    return ret

def main():
    r1, c1, r2, c2 = getlist()
    
    
    
    
    
    
    
    
    
    
    SIZE = 2 * 10 ** 6 + 10 ** 5
    SIZE += 1
    inv = [0] * SIZE  
    fac = [0] * SIZE  
    finv = [0] * SIZE  
    inv[1] = 1
    fac[0] = fac[1] = 1
    finv[0] = finv[1] = 1
    for i in range(2, SIZE):
        inv[i] = MOD - (MOD // i) * inv[MOD % i] % MOD
        fac[i] = fac[i - 1] * i % MOD
        finv[i] = finv[i - 1] * inv[i] % MOD
    factorial = fac
    gyaku = finv
    
    
    
    
    
    
    
    ans = routes_from_zero(r2, c2, factorial, gyaku) - routes_from_zero(r1 - 1, c2, factorial, gyaku)\
          - routes_from_zero(r2, c1 - 1, factorial, gyaku) + routes_from_zero(r1 - 1, c1 - 1, factorial, gyaku)
    print(ans % MOD)

if __name__ == '__main__':
    main()

Solution (Python):
class nCrMod():
    def __init__(self, mod):
        self.mod = mod
        self.fac = [1, 1]
        self.finv = [1, 1]
        self.inv = [0, 1]

    def prep(self, n):
        mod = self.mod
        f, fi = self.fac[-1], self.finv[-1]
        for i in range(len(self.fac), n + 1):
            fn = f * i % mod
            v = -self.inv[mod % i] * (mod // i) % mod
            fin = fi * v % mod
            f, fi = fn, fin
            self.fac.append(f)
            self.finv.append(fi)
            self.inv.append(v)

    def __call__(self, n, r):
        if len(self.fac) <= n:
            self.prep(n)
        return self.fac[n] * self.finv[r] * self.finv[n - r] % self.mod


def main():
    mod = 10 ** 9 + 7
    ncr = nCrMod(mod)

    r1, c1, r2, c2 = map(int, input().split())
    r1, c1 = r1 - 1, c1 - 1
    t = ncr(r2 + c2 + 2, c2 + 1) + ncr(r1 + c1 + 2, c1 + 1) - ncr(r1 + c2 + 2, c2 + 1) - ncr(r2 + c1 + 2, c1 + 1)
    print(t % mod)

main()

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


def get_euclidian (A, B):
  if B == 1:
    return (1)
  else:
    return (int((1 - A * get_euclidian (B, A % B)) / B))


F = [1]
for i in range(1, r2 + c2 + 3):
  F.append((F[i - 1] * i) % L)


def get_combi (n, r):

  Euc = get_euclidian(L, (F[r] * F[n - r]) % L)
  return ((F[n] * Euc) % L)

combi1 = get_combi(r2 + c2 + 2, r2 + 1)
combi2 = get_combi(r2 + c1 + 1, c1)
combi3 = get_combi(r1 + c2 + 1, r1)
combi4 = get_combi(r1 + c1, r1)

print ((int(combi1 - combi2 - combi3 + combi4) + L) % L)

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

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

def inv(x):
  m = mod
  u = 1
  v = 0
  while(m):
    t = x // m
    x -= t * m
    x, m = m, x
    u -= t * v
    u, v = v, u
  u %= mod
  if u < 0:
    u += mod
  return u

fac = [1]
for i in range(1, r2 + c2 + 3):
  fac.append(fac[i - 1] * i % mod)

def paths(r, c):
  return fac[r + c + 2] * inv(fac[r + 1] * fac[c + 1]) % mod

f = paths(r2, c2) - paths(r2, c1 - 1) - paths(r1 - 1, c2) + paths(r1 - 1, c1 - 1)
f %= mod
if f < 0:
   f += mod
print(f)

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

MOD = 10**9 + 7
MAX = r2 + c2 + 10
fact = [1] * (MAX + 1)  
finv = [1] * (MAX + 1)  
iinv = [1] * (MAX + 1)  
for i in range(2, MAX + 1):
    fact[i] = fact[i - 1] * i % MOD
    iinv[i] = MOD - iinv[MOD % i] * (MOD // i) % MOD
    finv[i] = finv[i - 1] * iinv[i] % MOD


def comb(n: int, k: int) -> int:
    if n < k or n < 0 or k < 0:
        return 0
    return (fact[n] * finv[k] % MOD) * finv[n - k] % MOD


def func(a, b):
    res = 0
    for i in range(a + 1):
        res += comb(i + b + 1, i + 1)
        res %= MOD
    return res


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

Solution (Python):
class Combination:
    """
    SIZEが10^6程度以下の二項係数を何回も呼び出したいときに使う
    使い方:
    comb = Combination(SIZE, MOD)
    comb(10, 3) => 120
    """

    def __init__(self, N, MOD=10 ** 9 + 7):
        self.MOD = MOD
        self.fact, self.inv = self._make_factorial_list(N)

    def __call__(self, n, k):
        if k < 0 or k > n:
            return 0
        res = self.fact[n] * self.inv[k] % self.MOD
        res = res * self.inv[n - k] % self.MOD
        return res

    def _make_factorial_list(self, N):
        fact = [1] * (N + 1)
        inv = [1] * (N + 1)
        MOD = self.MOD
        for i in range(1, N + 1):
            fact[i] = (fact[i - 1] * i) % MOD
        inv[N] = pow(fact[N], MOD - 2, MOD)
        for i in range(N, 0, -1):
            inv[i - 1] = (inv[i] * i) % MOD
        return fact, inv


if __name__ == "__main__":
    r1, c1, r2, c2 = map(int, input().split())
    mod = 10 ** 9 + 7
    comb = Combination(2 * 10 ** 6 + 10, mod)

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

    CD = (f(r2, c2) - f(r2, c1-1) + mod) % mod
    C = (f(r1-1, c2) - f(r1-1, c1-1) + mod) % mod
    ans = (CD - C + mod) % mod
    print(ans)

Solution (Python):
def getComb(n, k, MOD):
    if n < k:
        return 0
    if n-k < k:
        k = n-k
    
    comb = 1
    for x in range(n-k+1, n+1):
        comb = (comb * x) % MOD
    
    d = 1
    for x in range(1, k+1):
        d = (d * x) % MOD
    
    comb *= pow(d, MOD-2, MOD)
    return comb % MOD

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

ans += getComb(r2+c2+2, r2+1, MOD)
ans -= getComb(r2+c1+1, r2+1, MOD)
ans -= getComb(r1+c2+1, r1, MOD)
ans += getComb(r1+c1, r1, MOD)

print(ans % MOD)

Solution (Python):
import sys

def f(r, c, fact, revfact, mod):
    return (fact[r + c] * revfact[r] * revfact[c]) % mod

def g(r, c, fact, revfact, mod):
    if r == c == 0: return 1
    elif r == 0: return c + 1
    elif c == 0: return r + 1
    ans = 0
    for i in range(1, c + 2):
        ans += f(r, i, fact, revfact, mod)
        ans %= mod
    return ans

def solve():
    input = sys.stdin.readline
    r1, c1, r2, c2 = map(int, input().split())
    mod = 7 + 10 ** 9
    fact = [1] * (5 + r2 + c2)
    revfact = [1] * (5 + r2 + c2)
    for i in range(1, 5 + r2 + c2): fact[i] = (i * fact[i-1]) % mod
    revfact[4 + r2 + c2] = pow(fact[4 + r2 + c2], mod - 2, mod)
    for i in reversed(range(1, 4 + r2 + c2)): revfact[i] = ((i + 1) * revfact[i + 1]) % mod

    ans = g(r2, c2, fact, revfact, mod) - g(r2, c1 - 1, fact, revfact, mod) - g(r1 - 1, c2, fact, revfact, mod) + g(r1 - 1, c1 - 1, fact, revfact, mod)
    print((ans + mod) % mod)

    return 0

if __name__ == "__main__":
    solve()

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):
    return ((r+2)*C(c+r+2, r+2)-c-1)*inv[c+1]

r1, c1, r2, c2 = map(int, input().split())
MOD = 10**9+7
C_init()
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():
    r1, c1, r2, c2 = (int(i) for i in input().split())
    m = r2 + c2 + 3
    fac = [0] * m
    finv = [0] * m
    inv = [0] * m
    MOD = 10**9 + 7

    def COMBinitialize(m):
        fac[0] = 1
        finv[0] = 1
        if m > 1:
            fac[1] = 1
            finv[1] = 1
            inv[1] = 1
            for i in range(2, m):
                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 COMB(n, k):
        if n < k:
            return 0
        if n < 0 or k < 0:
            return 0
        return fac[n] * (finv[k] * finv[n - k] % MOD) % MOD

    COMBinitialize(m)
    ans = COMB(r2+c2+2, r2+1)
    ans -= COMB(r1+c2+1, r1) + COMB(r2+c1+1, c1)
    ans %= MOD
    ans += COMB(r1+c1, r1)
    print(ans % MOD)


if __name__ == '__main__':
    main()

Solution (Python):











import sys
sys.setrecursionlimit(10**6)
input=sys.stdin.readline
from math import floor,sqrt,factorial,hypot,log 
from heapq import heappop, heappush, heappushpop
from collections import Counter,defaultdict,deque
from itertools import accumulate,permutations,combinations,product,combinations_with_replacement
from bisect import bisect_left,bisect_right
from copy import deepcopy
from fractions import gcd
from random import randint
def ceil(a,b): return (a+b-1)//b
inf=float('inf')
mod = 10**9+7
def pprint(*A): 
    for a in A:     print(*a,sep='\n')
def INT_(n): return int(n)-1
def MI(): return map(int,input().split())
def MF(): return map(float, input().split())
def MI_(): return map(INT_,input().split())
def LI(): return list(MI())
def LI_(): return [int(x) - 1 for x in input().split()]
def LF(): return list(MF())
def LIN(n:int): return [I() for _ in range(n)]
def LLIN(n: int): return [LI() for _ in range(n)]
def LLIN_(n: int): return [LI_() for _ in range(n)]
def LLI(): return [list(map(int, l.split() )) for l in input()]
def I(): return int(input())
def F(): return float(input())
def ST(): return input().replace('\n', '')

class ModInt:
    def __init__(self, x):
        self.x = x % mod
    
    def __str__(self):
        return str(self.x)
    
    __repr__ = __str__

    def __add__(self, other):
        if isinstance(other, ModInt):
            return ModInt(self.x + other.x)
        else:
            return ModInt(self.x + other)

    __radd__ = __add__
        
    def __sub__(self, other):
        if isinstance(other, ModInt):
            return ModInt(self.x - other.x)
        else:
            return ModInt(self.x - other)

    def __rsub__(self, other):
        if isinstance(other, ModInt):
            return ModInt(other.x - self.x)
        else:
            return ModInt(other - self.x)

    def __mul__(self, other):
        if isinstance(other, ModInt):
            return ModInt(self.x * other.x)
        else:
            return ModInt(self.x * other)

    __rmul__ = __mul__

    def __truediv__(self, other):
        if isinstance(other, ModInt):
            return ModInt(self.x * pow(other.x, mod-2,mod))
        else:
            return ModInt(self.x * pow(other, mod - 2, mod))
            
    def __rtruediv(self, other):
        if isinstance(other, self):
            return ModInt(other * pow(self.x, mod - 2, mod))
        else:
            return ModInt(other.x * pow(self.x, mod - 2, mod))
            

    def __pow__(self, other):
        if isinstance(other, ModInt):
            return ModInt(pow(self.x, other.x, mod))
        else:
            return ModInt(pow(self.x, other, mod))
            

    def __rpow__(self, other):
        if isinstance(other, ModInt):
            return ModInt(pow(other.x, self.x, mod))
        else:
            return ModInt(pow(other, self.x, mod))

class BigCombination(object):

    def __init__(self, mod: int=10**9+7, max_n: int=3*(10**6)):
        fac, inv=[1], []
        fac_append, inv_append=fac.append, inv.append

        for i in range(1, max_n+1):
            fac_append(fac[-1] * i % mod)

        inv_append(pow(fac[-1], mod-2, mod))

        for i in range(max_n, 0, -1):
            inv_append(inv[-1] * i % mod)

        self.mod, self.factorial, self.inverse=mod, fac, inv[::-1]

    def get_combination(self, n, r):
        if n < r:
            return 0
        return self.factorial[n] * self.inverse[r] * self.inverse[n-r] % self.mod

    def get_permutation(self, n, r):
        if n < r:
            return 0
        return self.factorial[n] * self.inverse[n-r] % self.mod

def main():
    R1,C1,R2,C2 = MI()
    ans = ModInt(0)
    B = BigCombination()
    for r in range(R1,R2+1):
        ans += B.get_combination(C2+r+1,C2) - B.get_combination(C1-1+r+1,C1-1)
    print(ans)

if __name__ == '__main__':
    main()

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

MOD = 10**9+7
MAXN = r2+c2+5
fac = [1,1] + [0]*MAXN
for i in range(2,MAXN+2):
    fac[i] = fac[i-1] * i % MOD

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

def comb(n,r):
    if n < r: return 0
    if n < 0 or r < 0: return 0
    return fac[n] * inv(fac[r]) * inv(fac[n-r]) % MOD

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

Solution (Python):
r1,c1,r2,c2=map(int,input().split())
mod=10**9+7
n=[1]*(r2+c2+3)
for i in range(r2+c2+2):
  n[i+1]=(n[i]*(i+1))%mod
def g(a,b):
  x,y,z=n[a+b],n[a],n[b]
  return (((x*pow(y,mod-2,mod))%mod)*pow(z,mod-2,mod))%mod
print((g(r2+1,c2+1)-g(r2+1,c1)-g(r1,c2+1)+g(r1,c1))%mod)

Solution (Python):
def paths(r,c):
  return facts[r+c]*ifacts[r]*ifacts[c]

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

Solution (Python):
p = 10 ** 9 + 7 
N = (10 ** 6) * 2 + 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)

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

a = list(map(int,input().split()))
ans = cmb(a[2]+1,a[2]+a[3]+2,p) - cmb(a[1],a[2]+a[1]+1,p) - cmb(a[0],a[3]+a[0]+1,p) + cmb(a[0],a[0]+a[1],p)
print(ans%p)

Solution (Python):
k = 72
K = 1<<k
nu = lambda L: int("".join([bin(K+a)[-k:] for a in L[::-1]]), 2)
st = lambda n: bin(n)[2:] + "0"
li = lambda s, l: [int(a, 2) if len(a) else 0 for a in [s[-(i+1)*k-1:-i*k-1] for i in range(l)]]

def grow(d, v, h):
    h += [0] * d
    f = [(-1 if (i+d) % 2 else 1) * fainv[i] * fainv[d-i] % P * h[i] % P for i in range(d+1)]
    
    for idx, a in enumerate([d+1, d * fa[v-1] * fainv[v] % P, (d * fa[v-1] * fainv[v] + d + 1) % P]):
        t = [1] * (2*d+2)
        for i in range(1, 2*d+2): t[i] = t[i-1] * (a - d + i - 1) % P
        ti = [1] * (2*d+2)
        ti[-1] = pow(t[-1], P-2, P)
        for i in range(1, 2*d+2)[::-1]: ti[i-1] = ti[i] * (a - d + i - 1) % P
        
        g = [ti[i+1] * t[i] % P for i in range(2*d+1)]
        fg = li(st(nu(f) * nu(g)), d * 2 + 2)
            
        if idx == 1:
            for i in range(d+1):
                h[i] = h[i] * (fg[d+i] * t[d+i+1] % P * ti[i] % P) % P
        elif idx == 0:
            for i in range(d):
                h[i+d+1] = fg[d+i] * t[d+i+1] % P * ti[i] % P
        elif idx == 2:
            for i in range(d):
                h[i+d+1] = h[i+d+1] * (fg[d+i] * t[d+i+1] % P * ti[i] % P) % P
    return h


def create_table(v):
    s = 1
    X = [1, v+1]
    while s < v:
        X = grow(s, v, X)
        s *= 2
    table = [1]
    for x in X:
        table.append(table[-1] * x % P)
    return table

def fact(i, table):
    a = table[i//v]
    for j in range(i//v*v+1, i+1):
        a = a * j % P
    return a

P = 10**9+7
v = 1 << 11
fa = [1] * (2*v+2)
fainv = [1] * (2*v+2)
for i in range(2*v+1):
    fa[i+1] = fa[i] * (i+1) % P
fainv[-1] = pow(fa[-1], P-2, P)
for i in range(2*v+1)[::-1]:
    fainv[i] = fainv[i+1] * (i+1) % P

T = create_table(v)

f = lambda a, b: fact(a+b+2, T) * pow(fact(a+1, T) * fact(b+1, T), P-2, P) % P
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] + [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):
k = 72
K = 1<<k
nu = lambda L: int("".join([bin(K+a)[-k:] for a in L[::-1]]), 2)
st = lambda n: bin(n)[2:] + "0"
li = lambda s, l: [int(a, 2) if len(a) else 0 for a in [s[-(i+1)*k-1:-i*k-1] for i in range(l)]]

def grow(d, v, h):
    h += [0] * d
    f = [(-1 if (i+d) % 2 else 1) * fainv[i] * fainv[d-i] % P * h[i] % P for i in range(d+1)]
    
    for idx, a in enumerate([d+1, d * fa[v-1] * fainv[v] % P, (d * fa[v-1] * fainv[v] + d + 1) % P]):
        g = [pow(a - d + i - 1, P-2, P) if i else 0 for i in range(2*d+2)]
        fg = li(st(nu(f) * nu(g)), d * 8 - 1)
        p = 1
        for i in range(d+1):
            p = p * (a-i) % P
        for i in range(d+1):
            fg[d+i+1] = fg[d+i+1] * p % P
            p = p * (a+i+1) % P * pow(a-d+i, P-2, P) % P
        if idx == 1:
            for i in range(d+1):
                h[i] = h[i] * fg[d+i+1] % P
        elif idx == 0:
            for i in range(d):
                h[i+d+1] = fg[d+i+1]
        elif idx == 2:
            for i in range(d):
                h[i+d+1] = h[i+d+1] * fg[d+i+1] % P
    return h


def create_table(v):
    s = 1
    X = [1, v+1]
    while s < v:
        X = grow(s, v, X)
        s *= 2
    table = [1]
    for x in X:
        table.append(table[-1] * x % P)
    return table

def fact(i, table):
    a = table[i//v]
    for j in range(i//v*v+1, i+1):
        a = a * j % P
    return a

P = 10**9+7
v = 1 << 11
fa = [1] * (2*v+2)
fainv = [1] * (2*v+2)
for i in range(2*v+1):
    fa[i+1] = fa[i] * (i+1) % P
fainv[-1] = pow(fa[-1], P-2, P)
for i in range(2*v+1)[::-1]:
    fainv[i] = fainv[i+1] * (i+1) % P

T = create_table(v)

f = lambda a, b: fact(a+b+2, T) * pow(fact(a+1, T) * fact(b+1, T), P-2, P) % P
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 convolve(a, b):
    def fft():
        d = n // 2
        v1 = w1
        v2 = w2
        v3 = w3
        while d >= 1:
            u1 = 1
            u2 = 1
            u3 = 1
            for i in range(d):
                for j in range(i, n, 2*d):
                    f1[j], f1[j+d] = (f1[j] + f1[j+d]) % p1, u1 * (f1[j] - f1[j+d]) % p1
                    f2[j], f2[j+d] = (f2[j] + f2[j+d]) % p2, u2 * (f2[j] - f2[j+d]) % p2
                    f3[j], f3[j+d] = (f3[j] + f3[j+d]) % p3, u3 * (f3[j] - f3[j+d]) % p3
                    g1[j], g1[j+d] = (g1[j] + g1[j+d]) % p1, u1 * (g1[j] - g1[j+d]) % p1
                    g2[j], g2[j+d] = (g2[j] + g2[j+d]) % p2, u2 * (g2[j] - g2[j+d]) % p2
                    g3[j], g3[j+d] = (g3[j] + g3[j+d]) % p3, u3 * (g3[j] - g3[j+d]) % p3
                u1 = u1 * v1 % p1
                u2 = u2 * v2 % p2
                u3 = u3 * v3 % p3
            v1 = v1 * v1 % p1
            v2 = v2 * v2 % p2
            v3 = v3 * v3 % p3
            d //= 2
 
    def ifft():
        d = 1
        while d < n:
            v1 = pow(invw1, n // (2 * d), p1)
            v2 = pow(invw2, n // (2 * d), p2)
            v3 = pow(invw3, n // (2 * d), p3)
            u1 = 1
            u2 = 1
            u3 = 1
            for i in range(d):
                for j in range(i, n, 2*d):
                    f1[j+d] *= u1
                    f2[j+d] *= u2
                    f3[j+d] *= u3
                    f1[j], f1[j+d] = (f1[j] + f1[j+d]) % p1, (f1[j] - f1[j+d]) % p1
                    f2[j], f2[j+d] = (f2[j] + f2[j+d]) % p2, (f2[j] - f2[j+d]) % p2
                    f3[j], f3[j+d] = (f3[j] + f3[j+d]) % p3, (f3[j] - f3[j+d]) % p3
                u1 = u1 * v1 % p1
                u2 = u2 * v2 % p2
                u3 = u3 * v3 % p3
            d *= 2
 
    p1, g1 = (11 << 21) + 1, 3
    p2, g2 = (25 << 20) + 1, 3
    p3, g3 = (27 << 20) + 1, 5
    
    n0, n1 = len(a), len(b)
    n = 1 << (max(n0, n1) - 1).bit_length() + 1
    a = a + [0] * (n-n0)
    b = b + [0] * (n-n1)
    w1 = pow(g1, (p1 - 1) // n, p1)
    w2 = pow(g2, (p2 - 1) // n, p2)
    w3 = pow(g3, (p3 - 1) // n, p3)
    invw1 = pow(w1, p1-2, p1)
    invw2 = pow(w2, p2-2, p2)
    invw3 = pow(w3, p3-2, p3)
    
    f1 = [aa % p1 for aa in a]
    f2 = [aa % p2 for aa in a]
    f3 = [aa % p3 for aa in a]
    
    g1 = [aa % p1 for aa in b]
    g2 = [aa % p2 for aa in b]
    g3 = [aa % p3 for aa in b]
    
    fft()
    for i in range(n):
        f1[i] = f1[i] * g1[i] % p1
        f2[i] = f2[i] * g2[i] % p2
        f3[i] = f3[i] * g3[i] % p3
    ifft()
    invn1 = pow(n, p1 - 2, p1)
    invn2 = pow(n, p2 - 2, p2)
    invn3 = pow(n, p3 - 2, p3)
    k1 = pow(n * p2 * p3 % p1, p1-2, p1) * p2 * p3
    k2 = pow(n * p1 * p3 % p2, p2-2, p2) * p1 * p3
    k3 = pow(n * p1 * p2 % p3, p3-2, p3) * p1 * p2
    q = p1 * p2 * p3
    return [(f1[i] * k1 + f2[i] * k2 + f3[i] * k3) % q % P for i in range(n0 + n1 - 1)]

def grow(d, v, h):
    h += [0] * d
    f = [(-1 if (i+d) % 2 else 1) * fainv[i] * fainv[d-i] % P * h[i] % P for i in range(d+1)]
    
    for idx, a in enumerate([d+1, d * fa[v-1] * fainv[v] % P, (d * fa[v-1] * fainv[v] + d + 1) % P]):
        g = [pow(a - d + i - 1, P-2, P) if i else 0 for i in range(2*d+2)]
        fg = convolve(f, g)
        p = 1
        for i in range(d+1):
            p = p * (a-i) % P
        for i in range(d+1):
            fg[d+i+1] = fg[d+i+1] * p % P
            p = p * (a+i+1) % P * pow(a-d+i, P-2, P) % P
        if idx == 1:
            for i in range(d+1):
                h[i] = h[i] * fg[d+i+1] % P
        elif idx == 0:
            for i in range(d):
                h[i+d+1] = fg[d+i+1]
        elif idx == 2:
            for i in range(d):
                h[i+d+1] = h[i+d+1] * fg[d+i+1] % P
    return h


def create_table(v):
    s = 1
    X = [1, v+1]
    while s < v:
        X = grow(s, v, X)
        s *= 2
    table = [1]
    for x in X:
        table.append(table[-1] * x % P)
    return table

def fact(i, table):
    a = table[i//v]
    for j in range(i//v*v+1, i+1):
        a = a * j % P
    return a

P = 10**9+7
v = 1 << 11
fa = [1] * (2*v+2)
fainv = [1] * (2*v+2)
for i in range(2*v+1):
    fa[i+1] = fa[i] * (i+1) % P
fainv[-1] = pow(fa[-1], P-2, P)
for i in range(2*v+1)[::-1]:
    fainv[i] = fainv[i+1] * (i+1) % P

T = create_table(v)

f = lambda a, b: fact(a+b+2, T) * pow(fact(a+1, T) * fact(b+1, T), P-2, P) % P
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 convolve(a, b):
    def fft(f):
        d = n // 2
        v = w
        while d >= 1:
            u = 1
            for i in range(d):
                for j in range(i, n, 2*d):
                    f[j], f[j+d] = (f[j] + f[j+d]) % p, u * (f[j] - f[j+d]) % p
                u = u * v % p
            v = v * v % p
            d //= 2
 
    def ifft(f):
        d = 1
        while d < n:
            v = pow(invw, n // (2 * d), p)
            u = 1
            for i in range(d):
                for j in range(i, n, 2*d):
                    f[j+d] *= u
                    f[j], f[j+d] = (f[j] + f[j+d]) % p, (f[j] - f[j+d]) % p
                u = u * v % p
            d *= 2
 
    p, g = (33 << 66) + 1, 5
    n0, n1 = len(a), len(b)
    n = 1 << (max(n0, n1) - 1).bit_length() + 1
    a = a + [0] * (n-n0)
    b = b + [0] * (n-n1)
    w = pow(g, (p - 1) // n, p)
    invw = pow(w, p-2, p)
    fft(a), fft(b)
    for i in range(n):
        a[i] = a[i] * b[i] % p
    ifft(a)
    invn = pow(n, p - 2, p)
    return [a[i] * invn % p % P for i in range(n0 + n1 - 1)]

def grow(d, v, h):
    h += [0] * d
    f = [(-1 if (i+d) % 2 else 1) * fainv[i] * fainv[d-i] % P * h[i] % P for i in range(d+1)]
    
    for idx, a in enumerate([d+1, d * fa[v-1] * fainv[v] % P, (d * fa[v-1] * fainv[v] + d + 1) % P]):
        g = [pow(a - d + i - 1, P-2, P) if i else 0 for i in range(2*d+2)]
        fg = convolve(f, g)
        p = 1
        for i in range(d+1):
            p = p * (a-i) % P
        for i in range(d+1):
            fg[d+i+1] = fg[d+i+1] * p % P
            p = p * (a+i+1) % P * pow(a-d+i, P-2, P) % P
        if idx == 1:
            for i in range(d+1):
                h[i] = h[i] * fg[d+i+1] % P
        elif idx == 0:
            for i in range(d):
                h[i+d+1] = fg[d+i+1]
        elif idx == 2:
            for i in range(d):
                h[i+d+1] = h[i+d+1] * fg[d+i+1] % P
    return h


def create_table(v):
    s = 1
    X = [1, v+1]
    while s < v:
        X = grow(s, v, X)
        s *= 2
    table = [1]
    for x in X:
        table.append(table[-1] * x % P)
    return table

def fact(i, table):
    a = table[i//v]
    for j in range(i//v*v+1, i+1):
        a = a * j % P
    return a

P = 10**9+7
v = 1 << 11
fa = [1] * (2*v+2)
fainv = [1] * (2*v+2)
for i in range(2*v+1):
    fa[i+1] = fa[i] * (i+1) % P
fainv[-1] = pow(fa[-1], P-2, P)
for i in range(2*v+1)[::-1]:
    fainv[i] = fainv[i+1] * (i+1) % P

T = create_table(v)

f = lambda a, b: fact(a+b+2, T) * pow(fact(a+1, T) * fact(b+1, T), P-2, P) % P
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):
class Combination():
    def __init__(self, n, mod):
        self.n = n
        self.mod = mod
        self.fact = self.make_fact(n)
        self.fact_inv = self.make_fact_inv(n)
    def make_fact(self, n):
        res = [1]*(n+1)
        for i in range(1, n+1):
            res[i] = res[i-1]*i%self.mod
        return res
    def make_fact_inv(self, n):
        fact_inv = [1]*(n+1)
        fact_inv[n] = pow(self.fact[n], self.mod-2, self.mod)
        for i in range(n, 0, -1):
            fact_inv[i-1] = fact_inv[i]*i%self.mod
        return fact_inv
    def comb(self, m, k):
        return self.fact[m]*self.fact_inv[k]*self.fact_inv[m-k]%self.mod
mod = 10**9+7
r1, c1, r2, c2 = map(int, input().split())
comb = Combination(r2+c2+2, mod)
print((comb.comb(r2+c2+2, r2+1)-comb.comb(r2+c1+1, r2+1)-comb.comb(r1+c2+1, r1)+comb.comb(r1+c1, r1))%mod)

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

n = [1]*(r2+c2+3)
mod = 10**9 + 7
for i in range(r2+c2+2):
    n[i+1] = (n[i]*((i+1)%mod))%mod

def calc(r,c):
    x = n[r+c]%mod
    y = (n[r]*n[c])%mod
    return (x*(pow(y,mod-2,mod)%mod))%mod

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

print(ans)

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


def nCr(n, r, mod=MOD):
    r = min(r, n-r)
    numer = denom = 1
    for i in range(1, r+1):
        numer = numer * (n+1-i) % mod
        denom = denom * i % mod
    return numer * pow(denom, mod-2, mod) % mod

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

r1, c1, r2, c2 = map(int, input().split())
r1 -= 1; c1 -= 1
ans = g(r2, c2)
ans -= g(r1, c2)
ans -= g(r2, c1)
ans += g(r1, c1)
ans += MOD * 2
ans %= MOD
print(ans)

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



















MAX = 10 ** 6 * 2 + 5
def com_init(MAX, mod):
    fac = [0] * MAX
    finv = [0] * MAX
    inv = [0] * MAX
    fac[0] = fac[1] = 1
    finv[0] = finv[1] = 1
    inv[1] = 1
    for i in range(2, MAX):
        fac[i] = fac[i - 1] * i % mod
        inv[i] = mod - inv[mod % i] * (mod // i) % mod
        finv[i] = finv[i - 1] * inv[i] % mod
    return fac, finv

fac, finv = com_init(MAX, mod)

def nCr(n,r):
    return fac[n] * (finv[r] * finv[n - r] % mod) % mod 

def myfunc(r, c):
    return nCr(r + c + 2, r + 1) - 1


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

Solution (Python):
import numpy as np
import sys
read = sys.stdin.read
r1, c1, r2, c2 = map(int, read().split())
MOD = 10 ** 9 + 7


def cumprod(A, MOD=MOD):
    L = len(A)
    Lsq = int(L**.5 + 1)
    A = np.resize(A, Lsq**2).reshape(Lsq, Lsq)
    for n in range(1, Lsq):
        A[:, n] *= A[:, n - 1]
        A[:, n] %= MOD
    for n in range(1, Lsq):
        A[n] *= A[n - 1, -1]
        A[n] %= MOD
    return A.ravel()[:L]


def make_fact(U, MOD=MOD):
    x = np.arange(U, dtype=np.int64)
    x[0] = 1
    fact = cumprod(x, MOD)
    x = np.arange(U, 0, -1, dtype=np.int64)
    x[0] = pow(int(fact[-1]), MOD - 2, MOD)
    fact_inv = cumprod(x, MOD)[::-1]
    fact.flags.writeable = False
    fact_inv.flags.writeable = False
    return fact, fact_inv


fac, finv = make_fact(2 * 10**6 + 10, MOD)


def comb(n, k, mod=MOD, fac=fac, finv=finv):
    if n < k:
        return 0
    if n < 0 or k < 0:
        return 0
    return fac[n] * (finv[k] * finv[n - k] % mod) % mod


answer = comb(r2 + c2 + 2, r2 + 1) - comb(r1 + c2 + 1, r1) - \
    comb(c1 + r2 + 1, c1) + comb(r1 + c1, r1)
answer %= MOD
print(answer)