Problem:
 Score : 900 points    Problem Statement There is a rectangular sheet of paper with height H+1 and width W+1 . We introduce an xy -coordinate system so that the four corners of the sheet are (0, 0) , (W + 1, 0) , (0, H + 1) and (W + 1, H + 1) .  This sheet can be cut along the lines x = 1,2,...,W and the lines y = 1,2,...,H . Consider a sequence of operations of length K where we choose K of these H + W lines and cut the sheet along those lines one by one in some order.  Let the score of a cut be the number of pieces of paper that exist just after the cut. The score of a sequence of operations is the sum of the scores of all of the K cuts.  Find the sum of the scores of all possible sequences of operations of length K . Since this value can be extremely large, print the number modulo 10^9 + 7 .      Constraints  1 \leq H,W \leq 10^7  1 \leq K \leq H + W  H, W and K are integers.         Input Input is given from Standard Input in the following format:  H  W  K       Output Print the sum of the scores, modulo 10^9 + 7 .        Sample Input 1 2 1 2      Sample Output 1 34  Let x_1 , y_1 and y_2 denote the cuts along the lines x = 1 , y = 1 and y = 2 , respectively. The six possible sequences of operations and the score of each of them are as follows:   y_1, y_2 : 2 + 3 = 5  y_2, y_1 : 2 + 3 = 5  y_1, x_1 : 2 + 4 = 6  y_2, x_1 : 2 + 4 = 6  x_1, y_1 : 2 + 4 = 6  x_1, y_2 : 2 + 4 = 6   The sum of these is 34 .       Sample Input 2 30 40 50      Sample Output 2 616365902  Be sure to print the sum modulo 10^9 + 7 .

Solution (Python):
H,W,K=map(int,input().split());m,c=10**9+7,1;H+=W
for i in range(1,K):c=c*(H-i)%m;
v=pow(H-1,m-2,m);d=K*(K+1)//2;a=c*H*(d+K);a+=c*v*2*(H-W)*W*(K-1)*d//3;print(a%m)

Solution (Python):
H,W,K=map(int,input().split());m=10**9+7;c=1;H+=W
for i in range(H-1,H-K,-1):c=c*i %m;
v=pow(H-1,m-2,m);d=K*(K+1)//2;a=c*H*(d+K);d=(K-1)*d//3;a=(a+d*c*v*2*(H-W)*W)%m;print(a)

Solution (Python):
H,W,K=map(int,input().split());m=10**9+7;a,c=0,1;H+=W;mul=lambda x,y:(x*y)%m
for i in range(H-1,H-K,-1): c=c*i %m;
v=pow(H-1,m-2,m);temp1=mul(H,mul((K*(K+1)//2)%m,c));c=mul(c,v);u=((K*(K**2-1))//6)%m;temp2=mul(2*(H-W)*W,mul(u,c));c=mul(c,H-1);ans=mul(K,mul(H,c));ans=(ans+temp1+temp2)%m;print(ans)

Solution (Python):
mod=10**9+7
def mul(a,b):
    x,y=a,b
    return (x*y)%mod

H,W,K=map(int,input().split())
mod=10**9+7
ans=0
const=1
for i in range(H+W-1,H+W-K,-1):
    const=(const*i)%mod

inv=pow(H+W-1,mod-2,mod)
temp1=mul(H+W,mul((K*(K+1)//2)%mod,const))
const=mul(const,inv)
u=(K*(K**2-1))//6
u%=mod
temp2=mul(2*H*W,mul(u,const))
ans=(temp1+temp2)%mod
const=mul(const,H+W-1)
ans+=mul(K,mul(H+W,const))
ans%=mod
print(ans)

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

f_list = [1] * list_size
f_r_list = [1] * list_size

for i in range(list_size-1):
	f_list[i+1] = (f_list[i] * (i+1)) % MOD

f_r_list[-1] = pow(f_list[-1], MOD - 2, MOD)

for i in range(list_size-2, -1, -1):
	f_r_list[i] = (f_r_list[i+1] * (i+1)) % MOD

h, w, k = map(int, input().split())
ans = k * f_list[h+w] * f_r_list[h+w-k]
ans += (h+w) * k * (k+1) * f_list[h+w-1] * f_r_list[2] * f_r_list[h+w-k]
ans %= MOD
ans += h * w * f_list[h+w-2] * (k-1) * k * (k+1) * pow(3, MOD-2, MOD) * f_r_list[h+w-k]
ans %= MOD
print(ans)

Solution (Python):
k = 72
kk = k // 4
KK = 1<<k
nu = lambda L: int("".join([hex(KK+a)[3:] for a in L[::-1]]), 16)
st = lambda n: hex(n)[2:]
li = lambda s, l, r: [int(a, 16) % P if len(a) else 0 for a in [s[-(i+1)*kk:-i*kk] for i in range(l, r)]]

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)]
    nuf = nu(f)
    a = d * inv[v] % P
    t = [1] * (3*d+3)
    for i in range(1, 3*d+3): t[i] = t[i-1] * (a - d + i - 1) % P
    ti = [1] * (3*d+3)
    ti[-1] = pow(t[-1], P-2, P)
    for i in range(1, 3*d+3)[::-1]: ti[i-1] = ti[i] * (a - d + i - 1) % P
    iv = [1] * (3*d+3)
    for i in range(1, 3*d+3):
        iv[i] = ti[i] * t[i-1] % P

    
    g = [inv[i] for i in range(1, 2*d+2)]
    fg = li(st(nuf * nu(g)), d, d * 2 + 1)
    for i in range(d):
        h[i+d+1] = fg[i] * fa[d+i+1] % P * fainv[i] % P

    
    g = [iv[i] for i in range(1, 2*d+2)]
    fg = li(st(nuf * nu(g)), d, d * 2 + 1)
    for i in range(d+1):
        h[i] = h[i] * (fg[i] * t[d+i+1] % P * ti[i] % P) % P

    
    g = [iv[i] for i in range(d+2, 3*d+3)]
    fg = li(st(nuf * nu(g)), d, d * 2 + 1)
    for i in range(d):
        h[i+d+1] = h[i+d+1] * (fg[i] * t[2*d+i+2] % P * ti[d+i+1] % 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

H, W, K = map(int, input().split())
N = H + W
P = 10**9+7

v = 1 << (N.bit_length() + 1) // 2
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
inv = [0] * (2*v+2)
for i in range(1, 2*v+2):
    inv[i] = fainv[i] * fa[i-1] % P

T = create_table(v)

print((K*(K+3)//2 + K*(K+1)*(K-1)//3*H*W*pow(N * (N-1), P-2, P)) * fact(N, T) * pow(fact(N-K, T), P-2, P) % P)

Solution (Python):
k = 72
kk = k // 4
KK = 1<<k
nu = lambda L: int("".join([hex(KK+a)[3:] for a in L[::-1]]), 16)
st = lambda n: hex(n)[2:]
li = lambda s, l, r: [int(a, 16) % P if len(a) else 0 for a in [s[-(i+1)*kk:-i*kk] for i in range(l, r)]]

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)]
    nuf = nu(f)
    a = d * inv[v] % P
    t = [1] * (3*d+3)
    for i in range(1, 3*d+3): t[i] = t[i-1] * (a - d + i - 1) % P
    ti = [1] * (3*d+3)
    ti[-1] = pow(t[-1], P-2, P)
    for i in range(1, 3*d+3)[::-1]: ti[i-1] = ti[i] * (a - d + i - 1) % P
    iv = [1] * (3*d+3)
    for i in range(1, 3*d+3):
        iv[i] = ti[i] * t[i-1] % P

    
    g = [inv[i] for i in range(1, 2*d+2)]
    fg = li(st(nuf * nu(g)), d, d * 2 + 1)
    for i in range(d):
        h[i+d+1] = fg[i] * fa[d+i+1] % P * fainv[i] % P

    
    g = [iv[i] for i in range(1, 2*d+2)]
    fg = li(st(nuf * nu(g)), d, d * 2 + 1)
    for i in range(d+1):
        h[i] = h[i] * (fg[i] * t[d+i+1] % P * ti[i] % P) % P

    
    g = [iv[i] for i in range(d+2, 3*d+3)]
    fg = li(st(nuf * nu(g)), d, d * 2 + 1)
    for i in range(d):
        h[i+d+1] = h[i+d+1] * (fg[i] * t[2*d+i+2] % P * ti[d+i+1] % 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

H, W, K = map(int, input().split())
N = H + W
P = 10**9+7

v = 1 << (N.bit_length() + 1) // 2
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
inv = [0] * (2*v+2)
for i in range(1, 2*v+2):
    inv[i] = fainv[i] * fa[i-1] % P

T = create_table(v)

print((K*(K+3)//2 + K*(K+1)*(K-1)//3*H*W*pow(N * (N-1), P-2, P)) * fact(N, T) * pow(fact(N-K, T), P-2, P) % P)

Solution (Python):
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines
sys.setrecursionlimit(10 ** 7)

import numpy as np

MOD = 10**9 + 7

H,W,K = map(int,readline().split())

"""
・長方形ピースに対して右上のセルがどれであるかに注目
・各xに対してx回目の操作がそのセルについて加点対象になる確率を求めればよい
"""

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

p1 = K
p2 = (K*(K+1)//2)
p3 = ((K-1)*K*(K+1)//3) * H*W * pow((H+W)*(H+W-1),MOD-2,MOD)
p = (p1 + p2 + p3) % MOD

x = np.arange(H+W,H+W-K,-1,dtype=np.int64)
answer = p * cumprod(x,MOD)[-1] % MOD
print(answer)

Solution (Python):
H,W,K=map(int,input().split())
P,N=10**9+7,H+W
x=K*(K+3)//2+(K**3-K)//3*H*W*pow(N*N-N,P-2,P)
for i in range(K):x=x*(N-i)%P
print(x)

Solution (Python):
def solve(h, w, k):
    n = h + w
    MOD = 10 ** 9 + 7

    if k == 1:
        return 2 * n

    ans = 0
    
    ans += k * n * (n - 1) % MOD
    
    tmp = k * (k + 1) // 2 * (n - 1) % MOD
    ans += tmp * n % MOD
    
    
    tmp = (k - 1) * k * (k + 1) // 3 % MOD
    ans += tmp * h * w % MOD

    for m in range(n - 2, n - k, -1):
        ans = ans * m % MOD
    return ans % MOD


h, w, k = map(int, input().split())
print(solve(h, w, k))

Solution (Python):
h,w,k=map(int,input().split())
p,n=10**9+7,h+w
x=k+k*(k+1)//2+k*(k-1)*(k+1)//3*h*w*pow(n*(n-1),p-2,p)
for i in range(n,n-k,-1):x=x*i%p
print(x)

Solution (Python):
H, W, K = map(int, input().split())
P = 10**9 + 7
x = K + K*(K+1)//2 + K*(K+1)*(K-1)//3*H*W*pow((H+W)*(H+W-1),P-2,P)
for i in range(H+W, H+W-K, -1):
    x *= i
    x %= P
print(x)