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
 
MAXN = N*M + 10
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 nCr(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



ans = 0

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

    tmp *= nCr(N*M-2, K-2)

    tmp %= MOD
    
    ans += tmp

    ans %= MOD


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

    tmp *= nCr(N*M-2, K-2)

    tmp %= MOD
    ans += tmp
    ans %= MOD


print(ans)

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

MOD = 10**9+7
n = N*M-2
fac = [1]*(n+1)
rev = [1]*(n+1)
 
for i in range(1,n+1):
  fac[i] = i*fac[i-1]%MOD
  rev[i] = pow(fac[i], MOD-2, MOD)
 
comb = lambda a,b:(fac[a]*rev[a-b]*rev[b])%MOD

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

for d in range(M):
    ans += d*(M-d)*(N**2)*C
    ans %= MOD

print(ans)

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

MOD = 10**9+7
n = N*M-2
fac = [1]*(n+1)
rev = [1]*(n+1)
 
for i in range(1,n+1):
  fac[i] = i*fac[i-1]%MOD
  rev[i] = pow(fac[i], MOD-2, MOD)
 
comb = lambda a,b:(fac[a]*rev[a-b]*rev[b])%MOD

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

for d in range(M):
    ans += d*(M-d)*(N**2)*C
    ans %= MOD

print(ans)

Solution (Python):
n,m,k=map(int,input().split())
M=10**9+7
c=1
a=0
for i in range(1,k-1):c=(c*(n*m-1-i)*pow(i,M-2,M))%M
for i in range(n):a=(a+i*(n-i)*(m**2)*c)%M
for i in range(m):a=(a+i*(m-i)*(n**2)*c)%M
print(a)

Solution (Python):
import numpy as np
import sys
sys.setrecursionlimit(10 ** 7)



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

MOD = 10 ** 9 + 7


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]


def make_fact(U, 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]
    return fact, fact_inv


U = N * M + 10
fact, fact_inv = make_fact(U, MOD)
fact, fact_inv = fact.tolist(), fact_inv.tolist()


def mod_comb_k(n, k, mod):
    return fact[n] * fact_inv[k] % mod * fact_inv[n - k] % mod



for dx in range(1, N):
    pat = N - dx
    ans += dx * pat * M * M
    ans %= MOD


for dy in range(1, M):
    pat = M - dy
    ans += dy * pat * N * N
    ans %= MOD

ans *= mod_comb_k(N * M - 2, K-2, MOD)
ans %= MOD

print(ans)

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

def comb(n, r, mod=10**9+7):
  if n==0 and r==0:
    return 1
  else:
    a = b = 1
    A, B = n, r
    for i in range(r):
      a = a*A % mod
      b = b*B % mod
      A -= 1
      B -= 1
    return a * pow(b, mod-2, mod) % mod

ans = comb(N*M-2, K-2) * ((M*(N-1)*(N+1) + N*(M-1)*(M+1))*N*M//6)%MOD % MOD
print(ans)

Solution (Python):
class Factorial():
    def __init__(self,n,mod):
        self.mod = mod
        self.factorial = [0 for _ in range(n+1)]
        self.inv = [0 for _ in range(n+1)]
        self.factorial[0] = 1
        self.inv[0] = 1
        for i in range(n):
            self.factorial[i+1] = self.factorial[i]*(i+1)%mod
        self.inv[n] = pow(self.factorial[n],mod-2,mod)
        for i in range(n)[::-1]:
            self.inv[i] = self.inv[i+1]*(i+1)%mod
    def comb(self,m,k): 
        return self.factorial[m]*self.inv[k]*self.inv[m-k]%self.mod

from collections import deque

MOD = 1000000007

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

F = Factorial(N*M,MOD)

print(F.comb(N*M-2,K-2)*(N+M)*N*M*(N*M-1)//6%MOD)

Solution (Python):
N, M, K = 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 
MAX_N = 10**6 
g1 = [1, 1] 
g2 = [1, 1] 
inverse = [0, 1] 

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

k = cmb(M * N - 2,K - 2,mod)



def calc(p, q, d):
    return d * (p - d) * (q ** 2)

ans = 0

for i in range(N):
    ans += calc(N, M, i)
for i in range(M):
    ans += calc(M, N, i)

ans %= mod

ans *= k

ans %= mod

print(ans)

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

def main():
    N, M, K = list(map(int, input().split()))
    MOD = 10**9+7

    class Comb:
        def __init__(self, N):
            self.fac = [1] * (N+1)
            for i in range(2, N+1): self.fac[i] = self.fac[i-1] * i % MOD

        def pow(self, a, b):
            res = 1
            while b:
                if b & 1: res = res * a % MOD
                a = a**2 % MOD
                b >>= 1
            return res

        def comb(self, n, r):
            if r < 0 or r > n: return 0
            return (self.fac[n] * self.pow(self.fac[r], MOD-2)) % MOD * self.pow(self.fac[n-r], MOD-2) % MOD

        def permutation(self, n, r):
            if r == 0: return 1
            return self.fac[n] * self.pow(self.fac[n-r], MOD-2) % MOD

    comb = Comb(N*M-2)

    ans = 0
    for i in range(M):
        ans += i * (M-i) * (N**2)

    for i in range(N):
        ans += i * (N-i) * (M**2)

    print(ans*comb.comb(N*M-2, K-2)%MOD)

if __name__ == '__main__':
    main()

Solution (Python):
N,M,K=map(int,input().split())
def cmb(x,y):
  X=1
  Y=1
  for i in range(1,y+1):
    X= X*(x+1-i)%MOD
    Y= Y*i%MOD
  return (X*pow(Y,MOD-2,MOD))%MOD
MOD=10**9+7
A=cmb(M*N-2, K-2)
ans=0
for d in range(1,N):
  ans+=d*(N-d)*M*M*A%MOD
for d in range(1,M):
  ans+=d*(M-d)*N*N*A%MOD
print(ans%MOD)

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



mod = 10**9 + 7
max = n*m + 10
fac, finv, inv = [0]*max, [0]*max, [0]*max

def comInit(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

comInit(max)


def com(n,k):
    if(n < k):
        return 0
    if( n<0 | k < 0):
        return 0
    return fac[n] * (finv[k] * finv[n-k] % mod) % mod



ans = n*n*(m-1)*m*(m+1)//6 + m*m*(n-1)*n*(n+1)//6
ans = ans % mod
ans = (com(n*m -2, k-2) * ans) % mod

print(ans)

Solution (Python):
def power(a, b):
	if b == 0:
		return 1
	elif b == 1:
		return a % 1000000007
	elif b % 2 == 0:
		return (power(a, b//2) ** 2) % 1000000007
	else:
		return (power(a, b//2) ** 2 * a) % 1000000007
 
def divide(a, b):
	return (a * power(b, 1000000005)) % 1000000007

fac_lim = 200000
fac = [None]*(fac_lim+1)
fac[0] = 1
for i in range(fac_lim):
    fac[i+1] = fac[i] * (i+1)
    fac[i+1] = fac[i + 1] % 1000000007

fac_inv = [None]*(fac_lim+1)
fac_inv[fac_lim] = power(fac[fac_lim], 1000000005)
for i in range(fac_lim, 0, -1):
    fac_inv[i-1] = (fac_inv[i] * i)  % 1000000007

def conv(a, b):
    return (fac[a] * fac_inv[a-b] * fac_inv[b]) % 1000000007

MOD = 1000000007
N, M, K = (int(i) for i in input().split())
midN = 0
midM = 0
if N > 1:
    for i in range(1,N):
        midN += i*(N-i)*2
        midN % MOD
    midN = divide(midN, N**2)
else:
    midN = 0
if M > 1:
    for i in range(1,M):
        midM += i*(M-i)*2
        midM % MOD
    midM = divide(midM, M**2)
else:
    midM = 0
kt = (midN+midM)*(N*M)
kt = divide(kt, N*M-1)
pt = conv(N*M, K)
ans = kt*pt % MOD
ans = ans*K*(K-1) % MOD
ans = divide(ans, 2)
print(ans%MOD)

Solution (Python):
import sys
input = sys.stdin.readline
sys.setrecursionlimit(1000000)
from collections import deque, Counter
def getN():
    return int(input())
def getList():
    return list(map(int, input().split()))
import math
INF = 10 ** 20
MOD = 10**9 + 7
def combmod(n, k, mod=MOD):
    if k == 0:
        return 1
    if k == 1:
        return n
    ret = 1
    for i in range(n - k + 1, n + 1):
        ret *= i
        ret %= mod

    for i in range(1, k + 1):
        ret *= pow(i, mod - 2, mod)
        ret %= mod

    return ret

def main():
    n, m, k = getList()
    ans = 0
    mul = combmod(n*m-2, k-2)
    for i in range(0, n):
        for j in range(0, m):
            tmp = (j+i) * (n-i) * (m-j) * mul
            if i != 0 and j != 0:
                tmp *= 2
            ans += tmp
            ans %= MOD
            

    print(ans)
if __name__ == "__main__":
    main()

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


def prepare(n):
    global MOD
    modFacts = [0] * (n + 1)
    modFacts[0] = 1
    for i in range(n):
        modFacts[i + 1] = (modFacts[i] * (i + 1)) % MOD

    invs = [1] * (n + 1)
    invs[n] = pow(modFacts[n], MOD - 2, MOD)
    for i in range(n, 1, -1):
        invs[i - 1] = (invs[i] * i) % MOD

    return modFacts, invs


def modCombi(n, r):
    global MOD
    return (modFacts[n] * invs[n - r] * invs[r]) % MOD


N, M, K = map(int, input().split())
modFacts, invs = prepare(N * M)
var = modCombi(N * M - 2, K - 2)

ans = 0
if N < M:
    N, M = M, N
L = [0] * (N + 1)
for i in range(1, N + 1):
    L[i] = L[i - 1] + i
for x in range(1, N + 1):
    for y in range(1, M + 1):
        Sx = (M * (L[N - x] + L[x - 1])) % MOD
        Sy = (N * (L[M - y] + L[y - 1])) % MOD
        ans += (Sx + Sy) % MOD

print((ans * invs[2] * var) % MOD)

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

factorial = [1, 1]  
inverse = [1, 1]  
inverse_from = [0, 1]  

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


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


ans = 0

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

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

print(ans % MOD)

Solution (Python):

import sys, math, itertools, collections, bisect
input = lambda: sys.stdin.buffer.readline().rstrip().decode('utf-8')
inf = float('inf') ;mod = 10**9+7
mans = inf ;ans = 0 ;count = 0 ;pro = 1

class InverseTable:
    __slots__ = ['In', 'fac', 'Ifac']

    def __init__(self, mod, maxnum):
        self.In = [0, 1]; self.fac = [1, 1]; self.Ifac = [1, 1]
        for i in range( 2, maxnum + 1 ):
            self.In.append( ( -self.In[mod % i] * (mod//i) ) % mod )
            self.fac.append( ( self.fac[-1] * i ) % mod )
            self.Ifac.append( (self.Ifac[-1] * self.In[-1]) % mod )
        
    def cmb(self, n, r):
        if ( r<0 or r>n ): return 0
        r = min(r, n-r)
        return self.fac[n] * self.Ifac[r] * self.Ifac[n-r] % mod
    
    def per(self, n, r):
        if ( r<0 or r>n ): return 0
        return self.fac[n] * self.Ifac[n-r] % mod

IT = InverseTable(mod, 10**6)



n,m,k=map(int,input().split())
N0=n*m
for i in range(n):
  for j in range(m):
    ans += IT.cmb(N0-2,k-2)*\
      (((n-(i+1))*(n-(i+1)+1)+i*(i+1))//2*m+\
        ((m-(j+1))*(m-(j+1)+1)+j*(j+1))//2*n)
    ans %= mod
print(ans*IT.In[2]%mod)

Solution (Python):

def inv(x):
    global mod
    return pow(x, mod - 2, mod)


def fact(x):
    global mod
    res=1
    for i in range(2 , x + 1):
        res=res * i % mod
    return res


def combi(n,r):
    if r < 0 or r > n:
        return 0
    else:
        return fact(n) * inv(fact(r)) * inv(fact(n - r)) % mod

mod = 10 ** 9 + 7
N, M, K = map(int,input().split())
x = (N - 1) * N * (N + 1) * inv(6) % mod
x = x * M * M % mod
y = (M - 1) * M * (M + 1) * inv(6) % mod
y = y * N * N % mod
z = combi(M * N - 2, K - 2)
ans = (x + y) * z % mod
print(ans)

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


mod=10**9+7
kaidan=[1]*(n*m+1)
gyaku=[1]*(n*m+1)
for i in range(1,n*m+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
    elif n==0 or r==0 or n==r:
        return 1
    else:
        0
ans=0
for d in range(1,n):
    ans+=d*(n-d)*m*m*comb(n*m-2,k-2)
for d in range(1,m):
    ans+=d*(m-d)*n*n*comb(n*m-2,k-2)
print(ans%mod)

Solution (Python):
U = 2*10**5
MOD = 10**9+7
 
fact = [1]*(U+1)
fact_inv = [1]*(U+1)
 
for i in range(1,U+1):
    fact[i] = (fact[i-1]*i)%MOD
fact_inv[U] = pow(fact[U], MOD-2, MOD)
 
for i in range(U,0,-1):
    fact_inv[i-1] = (fact_inv[i]*i)%MOD

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

n, m, k = map(int, input().split())
ans = ((n-1)*n*(n+1)//6*m**2 + (m-1)*m*(m+1)//6*n**2) * comb(n*m-2, k-2)
print(ans%MOD)

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

def main():
    N,M,K = map(int,input().split())
    MOD = 10**9+7
    fac = [0 for _ in range(N*M-1)]
    fac[0],fac[1] = 1,1
    inv = copy.deepcopy(fac)
    invfac = copy.deepcopy(fac)
    
    for i in range(2,N*M-1):
        fac[i] = (fac[i-1]*i)%MOD
        inv[i] = MOD-(MOD//i)*inv[MOD%i]%MOD
        invfac[i] = (invfac[i-1]*inv[i])%MOD
        
    def coef(x,y):
        num = (((fac[x]*invfac[y])%MOD)*invfac[x-y]%MOD)
        return num

    xa = ((M**2)*N*(N-1)*(N+1)//6)%MOD
    ya = ((N**2)*M*(M-1)*(M+1)//6)%MOD
    print((xa+ya)*coef(N*M-2,K-2)%MOD)

if __name__ == "__main__":
    main()

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

g1=[1,1]
g2=[1,1]
for i in range(2,MAX+1):
    num_1=g1[-1]*i%mod
    g1.append(num_1)
    g2.append(pow(num_1,mod-2,mod))
    
def cmb(n,r):
    return g1[n]*g2[r]*g2[n-r]%mod


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

ans=0
for t in range(1,M):
    ans=(ans+t*(N**2)*(M-t))%mod
for t in range(1,N):
    ans=(ans+t*(M**2)*(N-t))%mod

print(ans*cmb(N*M-2,K-2)%mod)

Solution (Python):
mod = 10**9 + 7
n, m, k = map(int, input().split())
nm = n * m
fact = [1] * (nm-1)
inv = [1] * (nm-1)
inv_fact = [1] * (nm-1)
for i in range(2, nm-1):
  fact[i] = fact[i-1] * i % mod
  inv[i] = -(mod // i) * inv[mod % i] % mod
  inv_fact[i] = inv_fact[i-1] * inv[i] % mod
comb = fact[nm-2] * inv_fact[k-2] * inv_fact[nm-k] % mod
cost_x = (n-1) * nm * (n+1) // 6 * m % mod
cost_y = (m-1) * nm * (m+1) // 6 * n % mod
print(comb * (cost_x + cost_y) % mod)

Solution (Python):
import sys

MOD = 10 ** 9 + 7 

def make_table(size=10**6, p=MOD):
    fac = [None] * (size + 1)
    fac[0] = 1
    for i in range(size):
        fac[i+1] = fac[i] * (i + 1) % p
    ifac = [None] * (size + 1)
    ifac[size] = pow(fac[size], p-2, p)
    for i in range(size, 0, -1):
        ifac[i-1] = ifac[i] * i % p
    return fac, ifac

fac, ifac = make_table(2*10**5)

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

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

def main():
    vert = list(range(n+1))
    horiz = list(range(m+1))
    for i in range(1, n):
        vert[i+1] += vert[i]
        vert[i+1] %= MOD
    for i in range(1, m):
        horiz[i+1] += horiz[i]
        horiz[i+1] %= MOD
    
    res = 0
    for i in range(1, n):
        res += (vert[n] - vert[i] - i * (n - i)) * m % MOD * m % MOD
        res %= MOD
    for i in range(1, m):
        res += (horiz[m] - horiz[i] - i * (m - i)) * n % MOD * n % MOD
        res %= MOD
    
    res *= comb(n*m-2, k-2)
    res %= MOD
    return res

if __name__ == '__main__':
    ans = main()
    print(ans)

Solution (Python):


import sys

def input(): return sys.stdin.readline().strip()
def list2d(a, b, c): return [[c] * b for i in range(a)]
def list3d(a, b, c, d): return [[[d] * c for j in range(b)] for i in range(a)]
def list4d(a, b, c, d, e): return [[[[e] * d for j in range(c)] for j in range(b)] for i in range(a)]
def ceil(x, y=1): return int(-(-x // y))
def INT(): return int(input())
def MAP(): return map(int, input().split())
def LIST(N=None): return list(MAP()) if N is None else [INT() for i in range(N)]
def Yes(): print('Yes')
def No(): print('No')
def YES(): print('YES')
def NO(): print('NO')
sys.setrecursionlimit(10 ** 9)
INF = 10 ** 18
MOD = 10 ** 9 + 7

class ModTools:
    """ 階乗たくさん使う時用のテーブル準備 """

    def __init__(self, MAX, MOD):
        """ MAX：階乗に使う数値の最大以上まで作る """
        
        MAX += 1
        self.MAX = MAX
        self.MOD = MOD
        
        factorial = [1] * MAX
        factorial[0] = factorial[1] = 1
        for i in range(2, MAX):
            factorial[i] = factorial[i-1] * i % MOD
        
        inverse = [1] * MAX
        
        inverse[MAX-1] = pow(factorial[MAX-1], MOD-2, MOD)
        for i in range(MAX-2, 0, -1):
            
            inverse[i] = inverse[i+1] * (i+1) % MOD
        self.fact = factorial
        self.inv = inverse
    
    def nCr(self, n, r):
        """ 組み合わせの数 (必要な階乗と逆元のテーブルを事前に作っておく) """

        if n < r: return 0
        
        r = min(r, n-r)
        
        numerator = self.fact[n]
        
        denominator = self.inv[r] * self.inv[n-r] % self.MOD
        return numerator * denominator % self.MOD

H, W, K = MAP()

N = H * W
mt = ModTools(N, MOD)

ansh = 0
for i in range(1, H+1):
    ansh += i * (H-i) * pow(W, 2, MOD)
    ansh %= MOD
ansh *= mt.nCr(N-2, K-2)
ansh %= MOD

answ = 0
for i in range(1, W+1):
    answ += i * (W-i) * pow(H, 2, MOD)
    answ %= MOD
answ *= mt.nCr(N-2, K-2)
answ %= MOD

print((ansh + answ) % MOD)

Solution (Python):
n,m,k=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 = 3*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 )

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

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

MOD = 1000000007
def pow_mod(x, y):
    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
mod_inv = lambda x: pow_mod(x, MOD - 2)

a = 0
for i in range(1, n):
  a += (i * (n - i)) % MOD 
  a %= MOD
a *= (m * m) % MOD
a %= MOD

b = 0
for i in range(1, m):
  b += (i * (m - i)) % MOD 
  b %= MOD
b *= (n * n) % MOD
b %= MOD

w = (m * n - 2) % MOD
u = 1
v = 1
for i in range(1, k-1):
  u *= (w - i + 1)
  u %= MOD
  v *= i
  v %= MOD

ans = ((a + b) % MOD) 
ans *= u
ans %= MOD
ans *= mod_inv(v)
ans %= MOD


print(ans)

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

mod = 1000000007
fac = [1]
for i in range(1, n * m + 1):
    fac.append((fac[-1] * i) % mod)

def pow(x, n):
    ret = 1
    while n > 0:
        if (n & 1) == 1:
            ret = (ret * x) % mod
        x = (x * x) % mod
        n >>= 1
    return ret

def inv(i):
    return pow(fac[i], mod - 2)

def cmb(n, k):
    return (fac[n] * inv(k) * inv(n - k)) % mod

x = sum(d * (m - d) for d in range(m)) % mod
y = sum(d * (n - d) for d in range(n)) % mod

c = cmb(n * m - 2, k - 2)
x = (x * n * n * c) % mod
y = (y * m * m * c) % mod

print((x + y) % mod)

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

mod = 1000000007
def pow(x, n):
    ret = 1
    while n > 0:
        if (n & 1) == 1:
            ret = (ret * x) % mod
        x = (x * x) % mod
        n >>= 1
    return ret

fac = [1]
inv = [1]
for i in range(1, n * m + 1):
    fac.append((fac[-1] * i) % mod)
    inv.append(pow(fac[i], mod - 2))

def cmb(n, k):
    return (fac[n] * inv[k] * inv[n - k]) % mod

def doit(n, m, k):
    ret = 0
    for d in range(m):
        ret = (ret + d * (m - d)) % mod
    return (ret * n * n * cmb(n * m - 2, k - 2)) % mod

print((doit(n, m, k) + doit(m, n, k)) % mod)

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

mod = 1000000007
def pow(x, n):
    ret = 1
    while n > 0:
        if (n & 1) == 1:
            ret = (ret * x) % mod
        x = (x * x) % mod
        n //= 2
    return ret

fac = [1]
inv = [1]
for i in range(1, n * m + 1):
    fac.append((fac[-1] * i) % mod)
    inv.append(pow(fac[i], mod - 2))

def cmb(n, k):
    return ((fac[n] * inv[k]) % mod * inv[n - k]) % mod

def doit(n, m, k):
    ret = 0
    for d in range(m):
        ret = (ret + d * (m - d) * n * n * cmb(n * m - 2, k - 2)) % mod
    return ret

print((doit(n, m, k) + doit(m, n, k)) % mod)

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

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

def f(x):
    return (x * (x + 1) * (x - 1)) // 6 % MOD

tmp = (N * N * f(M)) % MOD + (M * M *f(N)) % MOD

class Factorial:
    def __init__(self, n, mod):
        self.f = [1]
        self.mod = mod
        for j in range(1, n + 1):
            self.f.append(self.f[-1] * j % mod)
        self.i = [pow(self.f[-1], mod - 2, mod)]
        for j in range(n, 0, -1):
            self.i.append(self.i[-1] * j % mod)
        self.i.reverse()
    def factorial(self, j):
        return self.f[j]
    def ifactorial(self, j):
        return self.i[j]
    def comb(self, n, k):
        return self.f[n] * self.i[n - k] % self.mod * self.i[k] % self.mod if n >= k else 0


c = Factorial(N * M - 2, MOD).comb

tmp *= c(N * M - 2, K - 2)
tmp %= MOD

print (tmp)

Solution (Python):
import math
def combinations_count(n, r):
    return math.factorial(n) // (math.factorial(n - r) * math.factorial(r))


n, m, k = map(int, input().split())
cmb = combinations_count(n * m - 2, k - 2)
cmb %= 10 ** 9 + 7
ans = 0
for i in range(1, n):
    ans += i * cmb * m * m * (n - i)
    ans %= 10 ** 9 + 7

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

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

Solution (Python):



def modmal(a,b,mod): 

    return a * b % mod




def moddiv(a,b,mod): 

    return (a * pow(b,mod-2,mod)) % mod




def inverse(a,mod): 
    return pow(a,mod-2,mod)







def modfac(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 modnCr(n,r,mod,fac,inv): 

    return fac[n] * inv[n-r] * inv[r] % mod

"""
とある二つのxのみに着目すると
距離0になるのは N * M * (M-1)通り
距離xになるのは (N-x) * M * M 通り
すべての場合において他の駒のおき方が (N*M-2)! / (N*M-K)! 通りある
駒を区別しないので駒の割り当て方→K!でわる

110 1
101 2
011 1
"""

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

fac,inv = modfac(N*M+1,mod)


ansx = 0

for dis in range(N):

    if dis == 0:
        ansx += N*M*(M-1) * dis
    else:
        ansx += 2 * (N-dis) * M * M * dis

    ansx %= mod

ansx *= fac[N*M-2] * inv[N*M-K]
ansx %= mod
ansx *= modnCr(K,2,mod,fac,inv)
ansx *= inv[K]

ansy = 0
for dis in range(M):

    if dis == 0:
        ansy += M*N*(N-1) * dis
    else:
        ansy += 2 * (M-dis) * N * N * dis

    ansy %= mod

ansy *= fac[N*M-2] * inv[N*M-K]
ansy %= mod
ansy *= modnCr(K,2,mod,fac,inv)
ansy *= inv[K]

print ((ansx+ansy)%mod)

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

mod = 10**9 + 7

A = ((N**2)*(N-1))//2 - (N*(2*N-1)*(N-1))//6
A = A*(M**2)
A %= mod
B = 1
for i in range(1,K-1):
    B *= N*M-1-i
    B %= mod
    B *= pow(i,mod-2,mod)
    B %= mod
ans1 = A*B

N,M = M,N
A = ((N**2)*(N-1))//2 - (N*(2*N-1)*(N-1))//6
A = A*(M**2)
A %= mod
B = 1
for i in range(1,K-1):
    B *= N*M-1-i
    B %= mod
    B *= pow(i,mod-2,mod)
    B %= mod
ans2 = A*B

print((ans1+ans2)%mod)

Solution (Python):
n,m,k = map(int,input().split())
mod = 10**9+7
rng = n*m+10
fctr = [1]
finv = [1]
for i in range(1,rng):
  fctr.append(fctr[-1]*i%mod)
for i in range(1,rng):
  finv.append(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

cost = [[0 for i in range(m+1)] for j in range(n+1)]
cost[1][1] = ((n-1)+(m-1))*n*m//2
for i in range(1,n+1):
  for j in range(1,m+1):
    if j>1:
      cost[i][j] = cost[i][j-1]-(m-2*(j-1))*n
    elif i>1:
      cost[i][j] = cost[i-1][j]-(n-2*(i-1))*m
ans = 0
for i in cost:
  ans = (ans+sum(i))%mod
ans = ans*pow(2,mod-2,mod)%mod
ans = ans*cmb(n*m-2,k-2)%mod
print(ans)

Solution (Python):

import re
import math
import fractions
import random
import time

mod=int(10**9+7)
inf=int(10**20)
fac=[1]*250000
for a in range(2,250000):
    fac[a]=fac[a-1]*a%mod
class union_find():
    def __init__(self,n):
        self.n=n
        self.P=[a for a in range(N)]
        self.rank=[0]*n

    def find(self,x):
        if(x!=self.P[x]):self.P[x]=self.find(self.P[x])
        return self.P[x]

    def same(self,x,y):
        return self.find(x)==self.find(y)

    def link(self,x,y):
        if self.rank[x]<self.rank[y]:
            self.P[x]=y
        elif self.rank[y]<self.rank[x]:
            self.P[y]=x
        else:
            self.P[x]=y
            self.rank[y]+=1

    def unite(self,x,y):
        self.link(self.find(x),self.find(y))

    def size(self):
        S=set()
        for a in range(self.n):
            S.add(self.find(a))
        return len(S)
def bin_(num,size):
    A=[0]*size
    for a in range(size):
        if (num>>(size-a-1))&1==1:
            A[a]=1
        else:
            A[a]=0
    return A
def comb(n,r):return math.factorial(n)//math.factorial(n-r)//math.factorial(r)
def next_comb(num,size):
    x=num&(-num)
    y=num+x
    z=num&(~y)
    z//=x
    z=z>>1
    num=(y|z)
    if(num>=(1<<size)):return False
    else:
        return num



def comb_(n,r):return fac[n]*pow(fac[n-r],mod-2,mod)*pow(fac[r],mod-2,mod)

def solve():
    return


N,M,K=map(int,input().split())
sum_=0
for a in range(1,N):
    sum_+=M*M*abs(N-a)*a*comb_(N*M-2,K-2)
    sum_%=mod
for a in range(1,M):
    sum_+=N*N*abs(M-a)*a*comb_(N*M-2,K-2)
    sum_%=mod

print(sum_)

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


n, m, k = [int(_) for _ in input().split()]

mod = 10**9+7 
N = n * m
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 )

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

Solution (Python):
import math
def com(n, r):
    return math.factorial(n) // (math.factorial(n - r) * math.factorial(r))
N,M,K = map(int,input().split())
mod = 10**9+7

ans=0
for i in range(M):
  l=(0+i)*(i+1)+(M-1-i)*(M-i)
  ans+=l//2*N*N%mod
  
for i in range(N):
  l=(0+i)*(i+1)+(N-1-i)*(N-i)
  ans+=l//2*M*M%mod

         
print(ans*com(N*M-2,K-2)//2%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

mod = 10**9+7 
N = 10**6
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 )


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

    for j in range(1,M):
        ans += j * (M - j) * N * N
        ans %= mod

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



if __name__ == "__main__":
    main()

Solution (Python):
import sys
mod = pow(10, 9) + 7
def C(n, r):
    if ( r<0 or r>n ):
        return 0
    r = min(r, n-r)
    return g1[n] * g2[r] * g2[n-r] % mod
NNN = (10**6)
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 )
n,m,k=map(int,input().split())

print(C(n*m,k)*C(k,2)%mod*(n+m)*inverse[3]%mod)

Solution (Python):
n,m,k=map(int,input().split())
U = 2*10**5
MOD = 10**9+7
 
fact = [1]*(U+1)
fact_inv = [1]*(U+1)
 
for i in range(1,U+1):
  fact[i] = (fact[i-1]*i)%MOD
fact_inv[U] = pow(fact[U],MOD-2,MOD)
 
for i in range(U,0,-1):
	fact_inv[i-1] = (fact_inv[i]*i)%MOD
 
def comb(n,k):
  if k < 0 or k > n:
    return 0
  x = fact[n]
  x *= fact_inv[k]
  x %= MOD
  x *= fact_inv[n-k]
  x %= MOD
  return x

ans=0
for i in range(n):
    for j in range(m):
        if i!=0 and j!=0:
            ans+=2*(n-i)*(m-j)*(i+j)
        else:
            ans+=(n-i)*(m-j)*(i+j)
        ans%=MOD

print(ans*comb(n*m-2,k-2)%MOD)

Solution (Python):
n, m, k = map(int, input().split())
mod = 10**9 + 7
cum_n = 0
t = 0
for i in range(1, n):
  t += i
  cum_n += t

cum_m = 0
t = 0
for i in range(1, m):
  t += i
  cum_m += t

ans = (cum_n * m**2 + cum_m * n**2) % mod

def nCr(n,r,mod = 10**9+7):
    r = min(n-r,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

if n*m-2 == 0:
  p = 1
else:
  p = nCr(n*m-2, k-2)

ans = (ans*p)%mod
print(ans)

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


def nCr(n, r, mod=10**9+7):
    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

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

for x in range(M):
    ans += x * (M-x) * N*N
for y in range(N):
    ans += y * (N-y) * M*M


ans *= nCr(N*M-2, K-2, MOD)
ans %= MOD

print(ans)