Problem:
 Score : 500 points    Problem Statement We have a grid of squares with N rows and M columns. Let (i, j) denote the square at the i -th row from the top and j -th column from the left. We will choose K of the squares and put a piece on each of them.  If we place the K pieces on squares (x_1, y_1) , (x_2, y_2) , ..., and (x_K, y_K) , the cost of this arrangement is computed as:  \sum_{i=1}^{K-1} \sum_{j=i+1}^K (|x_i - x_j| + |y_i - y_j|)  Find the sum of the costs of all possible arrangements of the pieces. Since this value can be tremendous, print it modulo 10^9+7 .  We consider two arrangements of the pieces different if and only if there is a square that contains a piece in one of the arrangements but not in the other.      Constraints  2 \leq N \times M \leq 2 \times 10^5  2 \leq K \leq N \times M  All values in input are integers.         Input Input is given from Standard Input in the following format:  N  M  K       Output Print the sum of the costs of all possible arrangements of the pieces, modulo 10^9+7 .        Sample Input 1 2 2 2      Sample Output 1 8  There are six possible arrangements of the pieces, as follows:   ((1,1),(1,2)) , with the cost |1-1|+|1-2| = 1  ((1,1),(2,1)) , with the cost |1-2|+|1-1| = 1  ((1,1),(2,2)) , with the cost |1-2|+|1-2| = 2  ((1,2),(2,1)) , with the cost |1-2|+|2-1| = 2  ((1,2),(2,2)) , with the cost |1-2|+|2-2| = 1  ((2,1),(2,2)) , with the cost |2-2|+|1-2| = 1   The sum of these costs is 8 .       Sample Input 2 4 5 4      Sample Output 2 87210       Sample Input 3 100 100 5000      Sample Output 3 817260251  Be sure to print the sum modulo 10^9+7 .

Solution (Python):
n, m, k = map(int, input().split())

mod = 10**9+7
N = 2*10**5+50
fac = [1]*(N+1)
finv = [1]*(N+1)
for i in range(N):
    fac[i+1] = fac[i] * (i+1) % mod
finv[-1] = pow(fac[-1], mod-2, mod)
for i in reversed(range(N)):
    finv[i] = finv[i+1] * (i+1) % mod

def cmb1(n, r, mod):
    if r <0 or r > n:
        return 0
    r = min(r, n-r)
    return fac[n] * finv[r] * finv[n-r] % mod

ans = 0
for i in range(n*m):
    x, y = divmod(i, m)
    
    ans += x*m*cmb1(n*m-2, k-2, mod)*x
    ans -= (n-1-x)*m*cmb1(n*m-2, k-2, mod)*x
    
    ans += y*n*cmb1(n*m-2, k-2, mod)*y
    ans -= (m-1-y)*n*cmb1(n*m-2, k-2, mod)*y

ans %= mod
print(ans)

Solution (Python):
n, m, k = map(int, input().split())

mod = 10**9+7
N = 2*10**5+50
fac = [1]*(N+1)
finv = [1]*(N+1)
for i in range(N):
    fac[i+1] = fac[i] * (i+1) % mod
finv[-1] = pow(fac[-1], mod-2, mod)
for i in reversed(range(N)):
    finv[i] = finv[i+1] * (i+1) % mod

def cmb1(n, r, mod):
    if r <0 or r > n:
        return 0
    r = min(r, n-r)
    return fac[n] * finv[r] * finv[n-r] % mod

ans = 0
for i in range(n*m):
    x, y = divmod(i, m)
    
    ans += x*m*cmb1(n*m-2, k-2, mod)*x
    ans -= (n-1-x)*m*cmb1(n*m-2, k-2, mod)*x
    
    ans += y*n*cmb1(n*m-2, k-2, mod)*y
    ans -= (m-1-y)*n*cmb1(n*m-2, k-2, mod)*y

ans %= mod
print(ans)

Solution (Python):
N,M,K = map(int,input().split())
MOD = 10**9+7

MAXN = N*M+5
fac = [1,1] + [0]*MAXN
finv = [1,1] + [0]*MAXN
inv = [0,1] + [0]*MAXN
for i in range(2,MAXN+2):
    fac[i] = fac[i-1] * i % MOD
    inv[i] = -inv[MOD%i] * (MOD // i) % MOD
    finv[i] = finv[i-1] * inv[i] % MOD

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

a = 0
for d in range(1,N):
    a += d*(N-d)*M*M
b = 0
for d in range(1,M):
    b += d*(M-d)*N*N

ans = (a+b) * comb(N*M-2, K-2)
ans %= MOD
print(ans)

Solution (Python):
class Combination:
    def __init__(self, n_max=10**6, mod=10**9+7):
        
        self._fac, self._finv, self._inv = [0]*n_max, [0]*n_max, [0]*n_max
        self._fac[0], self._fac[1] = 1, 1
        self._finv[0], self._finv[1] = 1, 1
        self._inv[1] = 1
        for i in range(2, n_max):
            self._fac[i] = self._fac[i - 1] * i % MOD
            self._inv[i] = MOD - self._inv[MOD%i] * (MOD // i) % MOD
            self._finv[i] = self._finv[i - 1] * self._inv[i] % MOD
    def com(self, n, r):
        if n < r: return 0
        if n < 0 or r < 0: return 0
        return self._fac[n] * (self._finv[r] * self._finv[n - r] % MOD) % MOD


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


MOD = 10**9+7
comb = Combination(10**6, MOD)
n,m,k = map(int, input().split())

ans = 0
for i in range(m):
    other_sum = i*(i+1)//2 + (m-1-i)*(m-i)//2
    
    pattern = n*n*comb.com(n*m-2,k-2)
    ans += pattern*other_sum
    ans%=MOD

for i in range(n):
    other_sum = i*(i+1)//2 + (n-1-i)*(n-i)//2
    
    pattern = m*m*comb.com(n*m-2,k-2)
    ans += pattern*other_sum
    ans%=MOD

a = modinv(2,MOD)
ans = (ans*a)%MOD
print(ans)

Solution (Python):
n,m,k=map(int,input().split())

mod = 10**9+7
f = [1]
for i in range(n*m):
  f += [f[-1]*(i+1) % mod]
def comb(a, b):
  return f[a]*pow(f[b], mod-2, mod)*pow(f[a-b], mod-2, mod) % mod

ans=0
for dist in range(1,n):
  c=n-dist
  p=(c*m*m)%mod
  ans+=p*dist
for dist in range(1,m):
  c=m-dist
  p=(c*n*n)%mod
  ans+=p*dist
ans%=mod
print(ans*comb(n*m-2,k-2)%mod)

Solution (Python):
import sys
input = sys.stdin.readline

def I(): return int(input())
def MI(): return map(int, input().split())
def LI(): return list(map(int, input().split()))

def main():
    mod=10**9+7
    """
    まず，xとyを分離する=>xだけを考えるので1行とみなす．
    いっぺんにいくつも考えられないので，コマを2つだけ置いた時を考える．
    2この組み合わせで距離が1=>M-1通り，2=>M-2通り,...,M-1=>1通り
    (以下のように考えてΣしても良いがやや面倒
        (j列目にコマを固定すると，相手との距離はj,j-1,...,1,0,1,2,...,M-j))
    実際にはN列あるので，選び方がN^2
        同じ行・列を選んだ場合は不可能だけど，どうせ距離も0なので関係ない
        
    ある2つの組みに注目した時，それを何回数えるかは，外野の並べ方に依存
    (N*M-2)C(K-2)通り
    """
    N,M,K=MI()
    
    
    def cmb(n, r, mod):
        if (r < 0) or (n < r):
            return 0
        r = min(r, n - r)
        return (fact[n] * factinv[r] * factinv[n-r])%mod

    fact=[1,1]
    factinv=[1,1]
    inv=[0,1]
    
    for i in range(2, N*M+5+ 1):
        fact.append((fact[-1] * i) % mod)
        inv.append((-inv[mod % i] * (mod // i)) % mod)
        factinv.append((factinv[-1] * inv[-1]) % mod)
    
    tempj=0
    for j in range(M):
        tempj+=(M-j)*j
        
    tempi=0
    for i in range(N):
        tempi+=(N-i)*i

        
    ans=(tempj*N*N)%mod
    ans=(ans+tempi*M*M)%mod
    ans=(ans*cmb(N*M-2,K-2,mod))%mod
    
    print(ans)
    

main()

Solution (Python):
import math
def choose(n,k):
    return math.factorial(n)//(math.factorial(n-k)*math.factorial(k))
mod = 10**9 +7
N,M,K = map(int,input().split())
ans = 0
memo = choose(N*M-2,K-2)
memo %= mod

d = 1
while d < M:
    ans += d*memo*(M-d)*N*N
    d += 1
    ans %= mod


d = 1
while d < N:
    ans += d*memo*(N-d)*M*M
    d += 1
    ans %= mod

print(ans)

Solution (Python):
import sys
readline = sys.stdin.readline
sys.setrecursionlimit(10**8)
mod = 10**9+7

INF = 10**18
eps = 10**-7

m,n,k = map(int,readline().split())

def comb(n, r, mod):
    r = min(r, n-r)
    mol = 1
    deno = 1
    for i in range(1, r+1):
        mol = mol * (n-r+i) % mod
        deno = deno * i % mod
    ret = mol * pow(deno, mod-2, mod) % mod
    return ret

ans = m*n*(m+n)*(m*n-1)//3

ans = (ans * comb(m*n-2,k-2,mod))%mod
ans = ans*pow(2,mod-2,mod)%mod

print(ans)

Solution (Python):
import sys
readline = sys.stdin.readline
sys.setrecursionlimit(10**8)
mod = 10**9+7

INF = 10**18
eps = 10**-7

m,n,k = map(int,readline().split())

def comb(n, r, mod):
    r = min(r, n-r)
    mol = 1
    deno = 1
    for i in range(1, r+1):
        mol = mol * (n-r+i) % mod
        deno = deno * i % mod
    ret = mol * pow(deno, mod-2, mod) % mod
    return ret

def f1(x,y,a,b):
    return (a+b-x-y)*(a-x+1)*(b-y+1)//2

def f2(x,a):
    return (a-x)*(a-x+1)//2

ans = 0

for x in range(1,n+1):
    for y in range(1,m+1):
        a = n-x+1
        b = m-y+1
        c = m*(a*(a-1)+x*(x-1))
        d = n*(b*(b-1)+y*(y-1))
        ans += (c+d)//2
        ans %= mod

ans = (ans * comb(m*n-2,k-2,mod))%mod
ans = ans*pow(2,mod-2,mod)%mod

print(ans)

Solution (Python):
import sys
import itertools

import time
import math
import heapq
from collections import defaultdict
from collections import Counter
sys.setrecursionlimit(10 ** 7)
 
INF = 10 ** 18
MOD = 10 ** 9 + 7
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines



N, M, K = map(int, input().split())
MAX = 2 * 10 ** 5 + 5
MOD = 10 ** 9 + 7

fac = [0 for i in range(MAX)]
finv = [0 for i in range(MAX)]
inv = [0 for i in range(MAX)]

def comInit(mod):
    fac[0], fac[1] = 1, 1
    finv[0], finv[1] = 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 com(n, r, mod):
    if n < r:
        return 0
    if n < 0 or r < 0:
        return 0
    return fac[n] * (finv[r] * finv[n - r] % mod) % mod

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

xs = 0
for i in range(M):
    xs += i * (M - i)
    xs %= MOD
xs *= N ** 2
ys = 0
for i in range(N):
    ys += i * (N - i)
    ys %= MOD
ys *= M ** 2
ys %= MOD
print(((xs + ys) % MOD) * com(N * M - 2, K - 2, MOD) % MOD)

Solution (Python):







class ModInt:
    def __init__(self, x, MOD=10 ** 9 + 7):
        '''
        pypyじゃないと地獄みたいな遅さなので注意
        '''
        self.mod = MOD
        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, self.mod)
        else:
            return ModInt(self.x + other, self.mod)

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

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

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

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

    __radd__ = __add__

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

    __rmul__ = __mul__

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

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


def combination(n, r):
    
    if r > n:
        return 0  
    r = min(r, n - r)
    nf = rf = ModInt(1)
    for i in range(r):
        nf = nf * (n - i)
        rf = rf * (i + 1)
    return nf / rf


N, M, K = map(int, input().split())
C = combination(N * M - 2, K - 2)
x_term = ModInt(0)
for dx in range(1, N):
    x_term += (N - dx) * dx
x_term *= M**2

y_term = ModInt(0)
for dy in range(1, M):
    y_term += (M - dy) * dy
y_term *= N**2

print(C * (x_term + y_term))

Solution (Python):
N, M, K = map(int, input().split())

def cmb(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
mod = 10 ** 9 + 7
n = 3 * 10 ** 5  
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)










ans1 = 0
for i in range(1, M):
  ans1 += i * (M - i) % mod
  
ans1 = (ans1 * N * N) % mod

ans2 = 0
for i in range(1, N):
  ans2 += i * (N - i) % mod
ans2 = (ans2 * M * M) % mod



ans = ans1 + ans2
print(ans * cmb(M * N - 2, K - 2, p) % mod)

Solution (Python):

def main():
    N, M, K = map(int, input().split())
    m = 10**9+7
    
    
    ans = 0
    for d in range(1,M):
        ans += (d*(M-d)%m) * pow(N,2,m) % m
        ans %= m

    for d in range(1,N):
        ans += (d*(N-d)%m) * pow(M,2,m) % m
        ans %= m

    
    
    for r in range(K-2):
        ans *= ((N*M%m-2)%m-r)%m
        ans %= m
        ans *= pow(r+1,m-2,m)
        ans %= m
    print(ans)

if __name__ == "__main__":
    main()

Solution (Python):
def cmb(n, r, p):
  if (r < 0) or (n < r):
    return 0
  r = min(r, n - r)
  return fac[n]*finv[r]*finv[n-r]%p

n,m,K = map(int,input().split())
MOD = pow(10,9)+7
N = pow(10,6)
fac = [-1]*(N+1); fac[0] = 1; fac[1] = 1 
finv = [-1]*(N+1); finv[0] = 1; finv[1] = 1 
inv = [-1]*(N+1); inv[0] = 0; inv[1] = 1 
for i in range(2,N+1):
  fac[i] = fac[i-1]*i%MOD
  inv[i] = MOD - inv[MOD%i]*(MOD//i)%MOD
  finv[i] = finv[i-1]*inv[i]%MOD

factor = cmb(n*m-2,K-2,MOD)
X = factor*(m-1)*m*(m+1)//6*n*n%MOD
Y = factor*(n-1)*n*(n+1)//6*m*m%MOD

ans = (X+Y)%MOD
print(ans)

Solution (Python):
N,M,K = map(int,input().split())
MOD = 10**9+7

MAXN = N*M+5
fac = [1,1] + [0]*MAXN
finv = [1,1] + [0]*MAXN
inv = [0,1] + [0]*MAXN
for i in range(2,MAXN+2):
    fac[i] = fac[i-1] * i % MOD
    inv[i] = -inv[MOD%i] * (MOD // i) % MOD
    finv[i] = finv[i-1] * inv[i] % MOD

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

a = 0
for d in range(1,N):
    a += d*(N-d) * M*M
b = 0
for d in range(1,M):
    b += d*(M-d) * N*N
ans = (a+b) * comb(N*M-2, K-2)
ans %= MOD
print(ans)

Solution (Python):
H, W, K = map(int, input().split())
MOD = 10 ** 9 + 7

factorial = [1, 1]
inverse = [1, 1]
invere_base = [0, 1]
for i in range(2, H * W + 2):
    factorial.append((factorial[-1] * i) % MOD)
    invere_base.append((-invere_base[MOD % i] * (MOD // i)) % MOD)
    inverse.append((inverse[-1] * invere_base[-1]) % MOD)


def nCr(n, r):
    if not 0 <= r <= n:
        return 0
    return factorial[n] * inverse[r] * inverse[n - r] % MOD


ans = 0
for dh in range(1, H):
    ans += dh * (H - dh) * pow(W, 2, MOD) * nCr(H * W - 2, K - 2) % MOD
    ans %= MOD
for dw in range(1, W):
    ans += dw * (W - dw) * pow(H, 2, MOD) * nCr(H * W - 2, K - 2) % MOD
    ans %= MOD

print(ans)

Solution (Python):
MAX = 2*10**5+100
MOD = 10**9+7
fact = [0]*MAX 
inv = [0]*MAX 
finv = [0]*MAX 
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

N, M, K = map(int, input().split())
ans = 0

for i in range(1, N):
    ans += i*(N-i)*M*M*C(N*M-2, K-2)
    ans %= MOD

for i in range(1, M):
    ans += i*(M-i)*N*N*C(N*M-2, K-2)
    ans %= MOD

print(ans)

Solution (Python):
import sys
readline = sys.stdin.readline
sys.setrecursionlimit(10**8)
mod = 10**9+7

INF = 10**18
eps = 10**-7

m,n,k = map(int,readline().split())

def comb(n, r, mod):
    r = min(r, n-r)
    mol = 1
    deno = 1
    for i in range(1, r+1):
        mol = mol * (n-r+i) % mod
        deno = deno * i % mod
    ret = mol * pow(deno, mod-2, mod) % mod
    return ret

def f1(x,y,a,b):
    return (a+b-x-y)*(a-x+1)*(b-y+1)//2

def f2(x,a):
    return (a-x)*(a-x+1)//2

ans = 0

for x in range(1,n+1):
    for y in range(1,m+1):
        ans += f1(x,y,n,m)+f1(1,1,x,y)+f1(1,y,x,m)+f1(x,1,n,y)
        ans -= f2(x,n) + f2(1,x) + f2(1,y) + f2(y,m)
        ans %= mod

ans = (ans * comb(m*n-2,k-2,mod))%mod
ans = ans*pow(2,mod-2,mod)%mod

print(ans)

Solution (Python):
N,M,K =map(int, input().split())
mod = 10**9 + 7






import numpy as np
 
def prepare(n, MOD):
    nrt = int(n ** 0.5) + 1
    nsq = nrt * nrt
    facts = np.arange(nsq, dtype=np.int64).reshape(nrt, nrt)
    facts[0, 0] = 1
    for i in range(1, nrt):
        facts[:, i] = facts[:, i] * facts[:, i - 1] % MOD
    for i in range(1, nrt):
        facts[i] = facts[i] * facts[i - 1, -1] % MOD
    facts = facts.ravel().tolist()
 
    invs = np.arange(1, nsq + 1, dtype=np.int64).reshape(nrt, nrt)
    invs[-1, -1] = pow(facts[-1], MOD - 2, MOD)
    for i in range(nrt - 2, -1, -1):
        invs[:, i] = invs[:, i] * invs[:, i + 1] % MOD
    for i in range(nrt - 2, -1, -1):
        invs[i] = invs[i] * invs[i + 1, 0] % MOD
    invs = invs.ravel().tolist()
 
    return facts, invs

facts, invs = prepare(N*M+10,mod)
def cmb(n,r,MOD):
  return (((facts[n] * invs[n-r]) % mod) * invs[r]) % mod







ans = 0

for dis in range(M):
  ans += N**2 * (M-dis) * dis * cmb(N*M-2, K-2, mod)
  ans %= mod
for dis in range(N):
  ans += M**2 * (N-dis) * dis * cmb(N*M-2, K-2, mod)
  ans %= mod

print(ans)

Solution (Python):
def sum_distance(a, b):
    x = 0
    for i in range(a):
        x += ((i + 1) * (a - i))
        x %= mod
    x *= pow(b, 2)
    return x % mod

def comb(n, r):
    x, y = 1, 1
    for i in range(n, n - r, -1):
        x *= i
        y *= i + r - n
        x %= mod
        y %= mod
    return pow(y, mod - 2, mod) * x % mod

n, m, k = map(int, input().split())
mod = pow(10, 9) + 7
ans = 0
ans += sum_distance(n - 1, m)
ans += sum_distance(m - 1, n)
ans %= mod
ans *= comb(n * m - 2, k - 2)
ans %= mod
print(ans)

Solution (Python):
class Combination:
    def __init__(self, mod, max_n):
        self.MOD = mod
        self.MAX_N = max_n

        self.f = self.factorial(self.MAX_N)
        self.f_inv = [self.inv(x) for x in self.f]

    def inv(self,x):
        return pow(x, self.MOD-2, self.MOD)
    
    def factorial(self, n):
        res = [1]
        for i in range(1,n+1):
            res.append(res[-1] * i % self.MOD)
        return res
    
    def comb(self, n, r):
        return (self.f[n] * self.f_inv[r] % self.MOD) * self.f_inv[n-r] % self.MOD

N, M, K = map(int,input().split())
CB = Combination(10**9+7, 200010)

ans = 0

for d in range(N):
    ans += (((CB.comb(N*M-2, K-2) * (N-d) % CB.MOD ) * M % CB.MOD ) * M % CB.MOD) * d % CB.MOD
    ans %= CB.MOD

for d in range(M):
    ans += (((CB.comb(N*M-2, K-2) * (M-d) % CB.MOD ) * N % CB.MOD ) * N % CB.MOD) * d % CB.MOD
    ans %= CB.MOD

print(ans)

Solution (Python):
from collections import defaultdict,deque
import math
def gcd(a,b):
    return b if not a%b else gcd(b,a%b)
def lcm(a,b):
    return (a*b)//gcd(a,b)
def prepare(n, MOD):
    
    f = 1
    factorials = [1]  
    for m in range(1, n + 1):
        f *= m
        f %= MOD
        factorials.append(f)
    
    inv = pow(f, MOD - 2, MOD)
    
    invs = [1] * (n + 1)
    invs[n] = inv
    for m in range(n, 1, -1):
        inv *= m
        inv %= MOD
        invs[m - 1] = inv

    return factorials, invs
    
def main():
    n,m,k=map(int, input().split())
    MOD = 10**9 + 7
    
    facs,invs = prepare(max(k,n*m),MOD)
    res = 0
    tmp = min(k,n)
    for i in range(1,n):
        res = res + ((n-i)*i*m*m)*(facs[n*m-2]%MOD * invs[k-2]%MOD * invs[n*m-2-(k-2)]%MOD)%MOD
        res%=MOD
    tmp = min(k,m)
    for i in range(1,m):
        res = res + ((m-i)*i*n*n)*(facs[n*m-2]%MOD * invs[k-2]%MOD * invs[n*m-2-(k-2)]%MOD)%MOD
        res%=MOD
    print(res)


if __name__ == '__main__':
    main()

Solution (Python):
import sys
def MI(): return map(int,sys.stdin.readline().rstrip().split())

N,M,K = MI()
mod = 10**9+7


def kaijou(n):
    res = 1
    for i in range(1,n+1):
        res *= i
        res %= mod
    return res

def nCr(n,r):
    return (kaijou(n)*pow(kaijou(r),mod-2,mod)*pow(kaijou(n-r),mod-2,mod)) % mod


a = N*(N+1)*(N-1)//6 * M**2 + M*(M+1)*(M-1)//6 * N**2
a %= mod
print((nCr(N*M-2,K-2) * a) % mod)

Solution (Python):
def inverse(a, mod):
    while(a < 0):
        a += mod
    a %= mod
    p = mod
    ya = 0
    yp = 1
    while (a > 0):
        n = p // a
        r = p % a
        yr = ya - n * yp
        ya = yp
        yp = yr
        p = a
        a = r
    return ya % mod if ya >= 0 else ya + mod

mod = 1000000007

n, m, k = map(int, input().split())

inv = 1
for i in range(1, k-1):
    inv = inv * i % mod
inv = inverse(inv, mod)

combination = inv
for i in range(n*m-k+1, n*m-1):
    combination = combination * i % mod


xn = n * (n + 1) * (n - 1) // 6 % mod * m % mod * m % mod
yn = m * (m + 1) * (m - 1) // 6 % mod * n % mod * n % mod
ans = (xn + yn) * combination % mod
print(ans)

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

    
def comb_mod(N,K,modp):

    K_fact=1
    NK_fact=1
    N_fact=1

    for i in range(1,K+1):
        K_fact*=i
        K_fact=K_fact%modp
    for i in range(1,N-K+1):
        NK_fact*=i
        NK_fact=NK_fact%modp
    for i in range(1,N+1):
        N_fact*=i
        N_fact=N_fact%modp

    K_inv=pow(K_fact, -1, modp)
    NK_inv=pow(NK_fact, -1, modp)
    pat=(N_fact*K_inv)%modp*NK_inv%modp
    return pat

pat=comb_mod(p*q-2,r-2,10**9+7)

ans=0
for i in range(0,q+1):
    ans+=i*(q-i)*p*p
    ans=ans%(10**9+7)
    
for k in range(0,p+1):
    ans+=k*(p-k)*q*q%(10**9+7)
    ans=ans%(10**9+7)

print(ans*pat%(10**9+7))

Solution (Python):

N_MAX = 10**6
MOD = 10**9 + 7
fac, finv, inv = [0]*N_MAX ,[0]*N_MAX, [0]*N_MAX
def com_init():
    fac[0], fac[1] = 1, 1
    finv[0], finv[1] = 1, 1
    inv[1] = 1
    for i in range(2, N_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 com(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


n,m,k = map(int, input().split())
com_init()

ans = 0
for dist in range(1,m):
    comb = (m-dist)*n*n
    other_comb = com(n*m-2,k-2)
    ans += comb*other_comb*dist
    ans %= MOD

for dist in range(1,n):
    comb = (n-dist)*m*m
    other_comb = com(n*m-2,k-2)
    ans += comb*other_comb*dist
    ans %= MOD

print(ans)

Solution (Python):
MOD = 10**9 + 7
N,M,K= map(int,input().split())

MOD = 10**9 + 7
num = 2000020
fact = [1] * (num+1)
ifact = [1] * (num+1)

for i in range(1,num+1):
    fact[i] = (fact[i-1] * i) % MOD
ifact[-1] = pow(fact[-1],MOD-2,MOD)
for i in range(1,num+1)[::-1]:
    ifact[i-1] = (ifact[i] * i)% MOD

memo = {}
def nCr2(n,r):
    if (n,r) in memo:
        return memo[(n,r)]
    else:
        memo[(n,r)] = (fact[n] * ifact[r] * ifact[n-r]) % MOD
        return memo[(n,r)]

def nCr(n,r):
    if r > n:
        return 0
    return (fact[n] * ifact[r] * ifact[n-r]) % MOD

ans = 0
ans1 = 0
ans2 = 0

for j in range(2,K+1):
    ans1 += M**2 * nCr(2*(M-1),j-2)* nCr((N-2)*M, K-j) 
    ans1 %= MOD
ans1 *= sum((i*(N-i) for i in range(1,N)))

for j in range(2,K+1):
    ans2 += N**2 * nCr(2*(N-1),j-2) * nCr((M-2)*N, K-j)
    ans2 %= MOD
ans2 *= sum((i*(M-i) for i in range(1,M)))

ans = (ans1 + ans2) % MOD
'''

for i in range(1,N):
    print(i)
    for j in range(2,K+1):
        for k in range(1,j):
            ans += i * (N-i) *k*(j-k)*nCr(M,k) * nCr(M,j-k)* nCr((N-2)*M, K-j) 
            ans %= MOD

for i in range(1,M):
    print(i)
    for j in range(2,min(K,2*N)+1):
        for k in range(1,j):
            ans += i * (M-i) * k* (j-k) * nCr(N,k) * nCr(N,j-k) * nCr((M-2)*N, K-j)
            ans %= MOD
'''

print(ans)

Solution (Python):
import sys
import heapq
import math
import fractions
import bisect
import itertools
from collections import Counter
from collections import deque
from operator import itemgetter
def input(): return sys.stdin.readline().strip()
def mp(): return map(int,input().split())
def lmp(): return list(map(int,input().split()))

mod = 10**9 + 7


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



def modinv(x):
    return power(x, mod-2)

def binomial_coefficients(n, k):
    numera = 1  
    denomi = 1  

    for i in range(k):
        numera *= n-i
        numera %= mod
        denomi *= i+1
        denomi %= mod
    return numera * modinv(denomi) % mod

ans=0
n,m,k=mp()
for i in range(n):
    for j in range(m):
        if i!=0 or j!=0:
            if i==0 or j==0:
                ans+=(n-i)*(m-j)*(i+j)%mod
            else:
                ans+=2*(n-i)*(m-j)*(i+j)%mod
            ans%=mod
print(ans*binomial_coefficients(n*m-2, k-2)%mod)

Solution (Python):
def ext_euc(a, b):
    
    if b == 0:
        return 1, 0, a
    y, x, v = ext_euc(b, a % b)
    y -= (a // b) * x
    return x, y, v


def mod_inv(a, mod):
    x, _, _ = ext_euc(a, mod)
    return x % mod


def comb(n, k, mod):
    if k >= n // 2:
        return comb(n, n - k, mod)
    ret = 1
    for i in range(k):
        ret *= n - i
        ret *= mod_inv(k - i, mod)
        ret %= mod
    return ret


def main():
    MOD = 10 ** 9 + 7
    N, M, K = list(map(int, input().split(' ')))
    ans = comb(N * M - 2, K - 2, MOD)
    ans *= ((N ** 2) * (M - 1) * M * (M + 1) + (M ** 2) * (N - 1) * N * (N + 1)) // 6
    ans %= MOD
    print(ans)


if __name__ == '__main__':
    main()

Solution (Python):
import sys
input = sys.stdin.readline

def I(): return int(input())
def MI(): return map(int, input().split())
def LI(): return list(map(int, input().split()))

def main():
    mod=10**9+7
    N,M,K=MI()
    
    
    def sum2(x):
        return ((x*(x+1))//2)%mod
    
    
    def calc(i,j):
        temp=0
        temp=sum2(M-j-1)
        temp=(temp+((sum2(M-j-1)+sum2(j))*(N-i-1)))%mod
        temp=(temp+sum2(N-i-1)*M)%mod
        return temp
        
    S=0
    for i in range(N):
        for j in range(M):
            S=(S+calc(i,j))%mod
            
    rem=(N*M)-2
    
    def cmb(n, r, mod):
        if (r < 0) or (n < r):
            return 0
        r = min(r, n - r)
        return (fact[n] * factinv[r] * factinv[n-r])%mod

    fact=[1,1]
    factinv=[1,1]
    inv=[0,1]
    
    for i in range(2, rem + 1):
        fact.append((fact[-1] * i) % mod)
        inv.append((-inv[mod % i] * (mod // i)) % mod)
        factinv.append((factinv[-1] * inv[-1]) % mod)  
        
    
    S=(S*cmb(rem,K-2,mod))%mod
            
    print(S)

main()

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**5)
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 )


N,M,K = map(int,input().split())
a=cmb(N*M-2, K-2, mod)
b=((N**2)*M*(M+1)*(M-1)//6)%(10**9+7)
c=((M**2)*N*(N+1)*(N-1)//6)%(10**9+7)
print((a*(b+c))%(10**9+7))

Solution (Python):
def cmb(a,b,c):
    b = min(b,a-b)
    num = 1
    for i in range(b):
        num = num*(a-i) % c
    den = 1
    for i in range(b):
        den = den*(i+1) % c
    return num * pow(den,c-2,c) % c
n,m,k = map(int,input().split())
mod = 10**9+7
ans = 0
C = cmb(n*m-2,k-2,mod)
for i in range(1,n):
    ans += i*(n-i)*m**2*C
    ans %= mod
for i in range(1,m):
    ans += i*(m-i)*n**2*C
    ans %= mod
print(ans)

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

fact = [1] * (MAX + 1)
for i in range(1, MAX + 1):
    fact[i] = (fact[i-1] * i) % mod

inv = [1] * (MAX + 1)
for i in range(2, MAX + 1):
    inv[i] = inv[mod % i] * (mod - mod // i) % mod

fact_inv = [1] * (MAX + 1)
for i in range(1, MAX + 1):
    fact_inv[i] = fact_inv[i-1] * inv[i] % mod


def comb(n, k):
    if n < k:
        return 0
    return fact[n] * fact_inv[n-k] * fact_inv[k] % mod


n, m, k = map(int, input().split())

ans = 0
cmb = comb(m * n - 2, k - 2)
for i in range(1, m):
    ans += i * (m - i) * n * n * cmb
    ans %= mod

for i in range(1, n):
    ans += i * m * m * (n - i) * cmb
    ans %= mod

print(ans)

Solution (Python):
n, m, k = map(int, input().split())
mod = 10 ** 9 + 7

res_n = 0
for i in range(1, n+1):
    res_n += (n + 1 - i) * (n - i) // 2
    res_n %= mod
res_n = (res_n * (m ** 2)) % mod

res_m = 0
for i in range(1, m+1):
    res_m += (m + 1 - i) * (m - i) // 2
    res_m %= mod
res_m = (res_m * (n ** 2)) % mod

res = res_n + res_m


f = [1 for _ in range(n*m)]
for i in range(1, n*m):
    f[i] = f[i-1] * i % mod

res = res * f[n*m-2] * pow(f[k-2], mod-2, mod) * pow(f[n*m-k], mod-2, mod) % mod

print(res)

Solution (Python):
import sys
sys.setrecursionlimit(10000000)
MOD = 10 ** 9 + 7
INF = 10 ** 15

def main():
    N,M,K = map(int,input().split())

    MAXN = N*M
    fact = [1]
    for i in range(1,MAXN + 1):
        fact.append(fact[-1]*i%MOD)
    inv_fact = [-1] * (MAXN + 1)
    inv_fact[-1] = pow(fact[-1],MOD - 2,MOD)
    for i in range(MAXN - 1,-1,-1):
        inv_fact[i] = inv_fact[i + 1]*(i + 1)%MOD
    
    nck = lambda N,K: 0 if K > N or K < 0 else fact[N]*inv_fact[N - K]*inv_fact[K]%MOD

    ans = 0
    const = nck(N*M - 2,K - 2)
    for i in range(1,N):
        ans += (N - i) * M*M * const * i % MOD
        ans %= MOD
    
    for i in range(1,M):
        ans += (M - i) * N*N * const * i % MOD
        ans %= MOD
    print(ans)
if __name__ == '__main__':
    main()

Solution (Python):
class Calc:
    def __init__(self, max_value, mod):
        """combination(max_value, all)"""
        fact = [-1] * (max_value + 1)
        fact[0] = 1
        fact[1] = 1
        for x in range(2, max_value + 1):
            fact[x] = x * fact[x - 1] % mod

        invs = [1] * (max_value + 1)
        invs[max_value] = pow(fact[max_value], mod - 2, mod)
        for x in range(max_value - 1, 0, -1):
            invs[x] = invs[x + 1] * (x + 1) % mod

        self.fact = fact
        self.invs = invs
        self.mod = mod

    def nCr(self, n, r):
        r = min(n - r, r)
        if r < 0: return 0
        if r == 0: return 1
        if r == 1: return n
        return self.fact[n] * self.invs[r] * self.invs[n - r] % self.mod

    def nHr(self, n, r):
        return self.nCr(n - 1 + r, r)


def main():
    MOD = 10 ** 9 + 7

    N, M, K = map(int, input().split())

    calc = Calc(max_value=N * M - 2, mod=MOD)

    t = 0
    for d in range(1, N):
        t = (t + (N - d) * d) % MOD
    t = t * M * M * calc.nCr(N * M - 2, K - 2) % MOD

    u = 0
    for d in range(1, M):
        u = (u + (M - d) * d) % MOD
    u = u * N * N * calc.nCr(N * M - 2, K - 2) % MOD

    ans = (t + u) % MOD

    print(ans)


if __name__ == '__main__':
    main()