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())

ans = 0
mod = 10**9+7

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

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

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

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

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

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

Solution (Python):


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

difx = 0
dify = 0

for i in range(m):
    difx += (i * (m-i) * pow(n, 2, mod))%mod

for j in range(n):
    dify += (j * (n-j) * pow(m, 2, mod))%mod

def comb_mod(n, r, p):
    num = 1
    den = 1
    r = min(r, n-r)
    for i in range(r):
        num = (num * (n-i))%p
        den = (den * (i+1))%p
    return num * pow(den, p-2, p)  

pat = comb_mod(n*m-2, k-2, mod)

ans = ((difx + dify) * pat)%mod

print(ans)

Solution (Python):
N, M, K = map(int, input().split())
mod = int(1e9) + 7
fact = [1 for _ in range(N*M+1)]
invf = [1 for _ in range(N*M+1)]
def inved(a):
  x, y, u, v, k, l = 1, 0, 0, 1, a, mod
  while l != 0:
    x, y, u, v = u, v, x - u * (k // l), y - v * (k // l)
    k, l = l, k % l
  return x % mod
for i in range(N*M):
  fact[i+1] = (fact[i] * (i + 1)) % mod
invf[-1] = inved(fact[-1])
for i in range(N*M, 0, -1):
  invf[i-1] = (invf[i] * i) % mod
Z = (fact[M*N-2] * invf[K-2] * invf[M*N-K]) % mod
X = M * (N + 1) * (N - 1)
X %= mod
Y = N * (M + 1) * (M - 1)
Y %= mod
XY = (inved(6) * M * N * (X + Y)) % mod
print((XY * Z) % mod)

Solution (Python):
n,m,k=map(int,input().split())
mod=10**9+7
f=[1]
for i in range(2*10**5+7):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 i in range(1,n):
  ans+=(i*(n-i)*m**2)*(comb(n*m-2,k-2))
  ans%=mod
for i in range(1,m):
  ans+=(i*(m-i)*n**2)*(comb(n*m-2,k-2))
  ans%=mod
print(ans)

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

def powmod(a,p):
  if p==0:
    return 1
  elif p==1:
    return a
  else:
    pow2=powmod(a,p//2)
    if p%2==0:
      return (pow2**2)%MOD
    else:
      return (a*pow2**2)%MOD

def invmod(a):
  return powmod(a,MOD-2)

def comb_mod(n,r):
  nPr=1
  fact_r=1
  for i in range(r):
    nPr*=n-i
    nPr%=MOD
    fact_r*=r-i
    fact_r%=MOD  
  return (nPr*invmod(fact_r))%MOD

answer_tate=0
for d in range(1,N):
  answer_tate+=d*M**2*(N-d)
  answer_tate%=MOD
  
answer_yoko=0
for d in range(1,M):
  answer_yoko+=d*N**2*(M-d)
  answer_yoko%=MOD

comb=comb_mod(N*M-2,K-2)
print(comb*(answer_tate+answer_yoko)%MOD)

Solution (Python):
import math
N,M,K=map(int,input().split())
mod=10**9+7
ans=0
for i in range(1,N):
  l=(N+1-i)*(N-i)//2%mod 
  ans+=l*M*M
  ans%mod
for i in range(1,M):
  l=(M+1-i)*(M-i)//2%mod 
  ans+=l*N*N
  ans%mod
  
def c_count(n, r):
    return math.factorial(n) // (math.factorial(n - r) * math.factorial(r))
                                 
ans*=c_count(M*N-2,K-2)%mod
print(ans%mod)

Solution (Python):
def comb(n, k, mod=10**9+7):
    num = den = 1
    for i in range(k):
        num = num * (n-i) % mod
        den = den * (i+1) % mod
    return num * pow(den, mod-2, mod) % mod

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

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

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

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

comb = getComb(N*M-2, K-2, MOD)

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

for j in range(M):
    ans += j*(M-j)*N*N
    ans %= MOD

ans *= comb
ans %= MOD

print(ans)

Solution (Python):
import sys

sys.setrecursionlimit(10 ** 7)
rl = sys.stdin.readline


def cmb(n, r, mod=10 ** 9 + 7):
    if r < 0 or n < r:
        return 0
    r = min(r, n - r)
    numerator, denominator = 1, 1
    for i in range(1, r + 1):
        numerator = (numerator * (n + 1 - i)) % mod
        denominator = (denominator * i) % mod
    return numerator * pow(denominator, mod - 2, mod) % mod


def solve():
    N, M, K = map(int, rl().split())
    MOD = 10 ** 9 + 7
    
    ans = 0
    for i in range(N):
        for j in range(M):
            tmp = (N - i) * (M - j) * (i + j) % MOD
            if i != 0 and j != 0:
                tmp *= 2
            ans = (ans + tmp) % MOD
    ans = ans * cmb(N * M - 2, K - 2) % MOD
    print(ans)


if __name__ == '__main__':
    solve()

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

fact = [1]*(N+1)
rfact = [1]*(N+1)
r = 1
for i in range(1, N+1):
  fact[i] = r = r * i % MOD
rfact[N] = r = pow(fact[N], MOD-2, MOD)
for i in range(N, 0, -1):
  rfact[i-1] = r = r * i % MOD


def perm(n, k):
  return fact[n] * rfact[n-k] % MOD


def comb(n, k):
  return fact[n] * rfact[k] * rfact[n-k] % MOD

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

p = comb(n*m-2, k-2)

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

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

print(ans)

Solution (Python):


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


def mod_nCr_init(n, mod):
    fac[0] = fac[1] = 1
    facinv[0] = facinv[1] = 1
    inv[1] = 1
    for i in range(2, n):
        fac[i] = fac[i-1] * i % mod
        inv[i] = -inv[mod % i] * (mod // i) % mod
        facinv[i] = facinv[i-1] * inv[i] % mod


def mod_nCr(n, r, mod):
    if n < r or n < 0 or r < 0:
        return 0
    return fac[n] * (facinv[r] * facinv[n-r] % mod) % mod


N, M, K = map(int, input().split())
fac = [0]*(N*M)
facinv = [0]*(N*M)
inv = [0]*(N*M)

mod = 10**9 + 7
mod_nCr_init(N*M, mod)
C = mod_nCr(N*M-2, K-2, mod)

ans = 0
for _ in [0, 0]:
    for d in range(1, N):
        ans += (N-d)*M*M*d*C

    N, M = M, N

print(ans % mod)

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 combination(self, n, r):
        if n - r < r:
            return self.combination(n, n - 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 main():
    MOD = 10 ** 9 + 7

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

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

    xsum = 0
    for d in range(1, N):
        xsum = (xsum + d * (N - d)) % MOD
    xsum = (xsum * M * M) % MOD

    ysum = 0
    for d in range(1, M):
        ysum = (ysum + d * (M - d)) % MOD
    ysum = (ysum * N * N) % MOD

    ans = (xsum + ysum) * calc.combination(N * M - 2, K - 2) % MOD

    print(ans)


if __name__ == '__main__':
    main()

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

def comb(n, r, fact, revfact):
    return (fact[n] * revfact[n-r] * revfact[r]) % mod
def solve():
    input = sys.stdin.readline
    N, M, K = map(int, input().split())
    fact = [1] * (N * M + 1)
    for i in range(1, N * M + 1): fact[i] = (fact[i-1] * i) % mod
    revfact = [1] * ( N * M + 1)
    revfact[N * M] = pow(fact[N * M], mod-2, mod)
    for i in reversed(range(1, N * M)): revfact[i] = ((i + 1) * revfact[i+1]) % mod

    others = comb(N * M - 2, K - 2, fact, revfact)
    pair_type = 0
    for i in range(N): 
        if i == 0:
            for j in range(1, M): 
                pair_type += (i + j) * (N - i) * (M - j) % mod
                pair_type %= mod
        else:
            for j in range(M):
                if j == 0: pair_type += (i + j) * (N - i) * (M - j) % mod
                else: pair_type += 2 * (i + j) * (N - i) * (M - j) % mod
                pair_type %= mod
    
    print(pair_type * others % mod)
    return 0

if __name__ == "__main__":
    solve()

Solution (Python):
class Combination:
    def __init__(self, size, mod=10**9 + 7):
        self.size = size + 2
        self.mod = mod
        self.fact = [1, 1] + [0] * size
        self.factInv = [1, 1] + [0] * size
        self.inv = [0, 1] + [0] * size

        for i in range(2, self.size):
            self.fact[i] = self.fact[i - 1] * i % self.mod
            self.inv[i] = -self.inv[self.mod % i] * (self.mod // i) % self.mod
            self.factInv[i] = self.factInv[i - 1] * self.inv[i] % self.mod

    def npr(self, n, r):
        if n < r or n < 0 or r < 0:
            return 0
        return self.fact[n] * self.factInv[n - r] % self.mod

    def ncr(self, n, r):
        if n < r or n < 0 or r < 0:
            return 0
        return self.fact[n] * (self.factInv[r] * self.factInv[n - r] % self.mod) % self.mod

    def nhr(self, n, r):  
        return self.ncr(n + r - 1, n - 1)

    def factN(self, n):
        if n < 0:
            return 0
        return self.fact[n]

N, M, K = map(int, input().split())
MOD = 10**9 + 7
comb = Combination(N * M + 100)

S = 0
for x in range(1, N):
    S += (N - x) * (M**2) * x
for y in range(1, M):
    S += (M - y) * (N**2) * y

ans = S * comb.ncr(N * M - 2, K - 2)
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

n,m,k = map(int,input().split())
nm = n*m
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

invt = InverseTable(mod, 10**6)


for i in range(n):
  for j in range(m):
    
    ans += ((i * (i+1) // 2 + (n-1-i) * (n-i) //2) * m + (j * (j+1) // 2 + (m-1-j) * (m-j) //2) * n) * invt.cmb(nm-2,k-2)
    ans %= mod
print(ans*pow(2,mod-2,mod)%mod)

Solution (Python):
import sys
input = sys.stdin.readline
N, M, K = map(int, input().split())
mod = 10 ** 9 + 7
class Factorial:
  def __init__(self, n, mod):
    self.f = [1]
    for i in range(1, n + 1):
      self.f.append(self.f[-1] * i % mod)
    self.i = [pow(self.f[-1], mod - 2, mod)]
    for i in range(1, n + 1)[: : -1]:
      self.i.append(self.i[-1] * i % mod)
    self.i.reverse()
  def factorial(self, i):
    return self.f[i]
  def ifactorial(self, i):
    return self.i[i]
  def combi(self, n, k):
    return self.f[n] * self.i[n - k] % mod * self.i[k] % mod

f = Factorial(N * M, mod)
t = f.combi(N * M - 2, K - 2)

res = 0
for d in range(1, M):
  x = N ** 2 % mod * (M - d) % mod
  res += x * d % mod
  res %= mod
  
for d in range(1, N):
  y = M ** 2 % mod * (N - d) % mod
  res += y * d % mod
  res %= mod

print(res * t % mod)

Solution (Python):


class Cmber:
    def __init__(self,gensize,mod=10**9+7):
        self.mod = mod
        self.g1=[1,1]
        self.g2 = [1, 1] 
        self.inverse = [0, 1] 
        self.getgen(gensize)
    def getgen(self,N):
        for i in range( 2, N + 1 ):
            self.g1.append( ( self.g1[-1] * i ) % self.mod )
            self.inverse.append( ( -self.inverse[self.mod % i] * (self.mod//i) ) % self.mod )
            self.g2.append( (self.g2[-1] * self.inverse[-1]) % self.mod )
 
    def cmb(self,n, r):
        if ( r<0 or r>n ):
            return 0
        r = min(r, n-r)
        return self.g1[n] * self.g2[r] * self.g2[n-r] % self.mod
 
    def cmb_very_high_n(self,n, r):
        if ( r<0 or r>n ):
            return 0
        g1=1
        r = min(r, n-r)
        for i in range(n-r+1,n+1):
            g1*=i
            g1%=self.mod
        return g1 * self.g2[r] % self.mod

    def multiple_cmb(self,vals):
        ans = self.g1[sum(vals)] % self.mod
        for i in vals:
            ans = (ans*self.g2[i])%self.mod
        return ans










def resolve():
    N,M,K = map(int,input().split())
    cmber = Cmber(N*M)
    score = 0
    for dist in range(1,N):
        cmb_of_2 = dist*(((((N-dist)*M)%cmber.mod)*M)% cmber.mod) % cmber.mod

        score = (score+(cmb_of_2*cmber.cmb(N*M-2,K-2))%cmber.mod)%cmber.mod
        
    for dist in range(1,M):
        cmb_of_2 = dist*(((((M-dist)*N)%cmber.mod)*N)% cmber.mod) % cmber.mod

        score = (score+(cmb_of_2*cmber.cmb(N*M-2,K-2))%cmber.mod)%cmber.mod

    print(score)

if __name__ == "__main__":
    resolve()

Solution (Python):
n,m,k=map(int,input().split())
mod=10**9+7
x=1
y=1
for i in range(1,k-1):
    y*=i
    y%=mod
for i in range(k-2):
    x*=n*m-2-i
    x%=mod
x*=pow(y,mod-2,mod)
x%=mod
def f(n,m):
    s=0
    for i in range(1,n):
        s+=i*(n-i)
        s%=mod
    return s*m*m%mod
print((f(n,m)+f(m,n))*x%mod)

Solution (Python):
n,m,k = map(int,input().split())
mod = 10**9+7
def combination(n,r,p):
    if n < r or r < 0:
        return 0
    r = min(n-r,r)
    
    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)
    return fact[n]*factinv[r]*factinv[n-r]%p
temp = combination(n*m-2,k-2,mod)
ans = 0
for d in range(1,n):
    ans += d*((n-d)*(m**2))*temp
for d in range(1,m):
    ans += d*((m-d)*(n**2))*temp
print(ans%mod)

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

mod = 10**9+7 
fact = [1,1]
finv = [1,1]
inv = [0,1]
 
for i in range(2,n*m+1):
    fact.append((fact[-1]*i)%mod)
    inv.append((inv[mod%i]*(mod-mod//i))%mod)
    finv.append((finv[-1]*inv[-1])%mod)
 
def nCr(n,r,mod):
    if r > n:
        return 0
    else: 
        return fact[n]*finv[r]*finv[n-r]%mod    


ans = 0

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

Solution (Python):

mod=10**9+7
def find_power(n,mod):
    
    powlist=[0]*(n+1)
    powlist[0]=1
    powlist[1]=1
    for i in range(2,n+1):
        powlist[i]=powlist[i-1]*i%(mod)
    return powlist
 
def pow_cal(x,n,mod):
    if n==0:
         return 1
    elif n==1:
         return x%mod
    elif n>=2:
        if n%2==0:
             return (pow_cal(x,n//2,mod)**2)%mod
        else:
             return (x*pow_cal(x,n//2,mod)**2)%mod

def find_inv_power(n,mod):
    
    c=1
    uselist=[0 for i in range(n+1)]
    for i in range(1,n+1):
        c*=i
        c%=mod
    first=pow_cal(c,mod-2,mod)
    uselist[n]=first
    for i in range(n,0,-1):
        uselist[i-1]=(uselist[i]*i)%(mod)
    return uselist
    
A=find_power(4*10**5,mod)
B=find_inv_power(4*10**5,mod)

def combi(n,r,mod):
    if n<r:
        return 0
    elif n>=r:
        return (A[n]*B[r]*B[n-r])%(mod)

    
M,N,K=map(int,input().split())
ans=M*N*(M+N)*(M*N-1)//6

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

Solution (Python):

import sys
input = sys.stdin.readline




class Combination:
    def __init__(self, n, mod):
        self.mod = mod
        self.fact = [1]*n
        self.finv = [1]*n
        self.inv  = [1]*n
        for i in range(2,n):
            self.fact[i] = (self.fact[i-1]*i) % self.mod
            self.inv[i]  = self.mod - self.inv[self.mod%i] * (self.mod//i)%self.mod
            self.finv[i] = self.finv[i-1] * self.inv[i] % self.mod
 
    def nCr(self, n, r):
        if n<r:
            return 0
        else:
            return self.fact[n] * (self.finv[r] * self.finv[n-r] % self.mod) % self.mod

mod = 10**9+7
Combination = Combination(200010, mod)





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

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

print(ans)

Solution (Python):
def framod(n, mod, a=1):
    for i in range(1,n+1):
        a=a * i % mod
    return a

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

def comb(n, k, mod):
    a=framod(n, mod)
    b=framod(k, mod)
    c=framod(n-k, mod)
    return (a * power(b, mod-2, mod) * power(c, mod-2, mod)) % mod

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

ans = 0

for i in range(1, M):
    ans += (i * (M-i) * N**2) % (10**9 + 7) 
for i in range(1, N):
    ans += (i * (N-i) * M**2) % (10**9 + 7)

ans *= comb(N*M-2 , K-2, 10**9+7)

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

Solution (Python):

import sys
input = sys.stdin.readline




class Combination:
    def __init__(self, n, mod):
        self.mod = mod
        self.fact = [1]*n
        self.finv = [1]*n
        self.inv  = [1]*n
        for i in range(2,n):
            self.fact[i] = (self.fact[i-1]*i) % self.mod
            self.inv[i]  = self.mod - self.inv[self.mod%i] * (self.mod//i)%self.mod
            self.finv[i] = self.finv[i-1] * self.inv[i] % self.mod
 
    def nCr(self, n, r):
        if n<r:
            return 0
        else:
            return self.fact[n] * (self.finv[r] * self.finv[n-r] % self.mod) % self.mod

mod = 10**9+7
Combination = Combination(200010, mod)





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

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

print(ans)

Solution (Python):

import sys
input = sys.stdin.readline




def init_fact(n,mod):
    fact,finv,inv = [1]*n,[1]*n,[1]*n
    for i in range(2,n):
        fact[i] = (fact[i-1]*i) % mod
        inv[i]  = mod - inv[mod%i] * (mod//i)%mod
        finv[i] = finv[i-1] * inv[i] % mod
    return (fact,finv,inv)
 
def nCr(n,r,mod,fact,finv):
    if n<r:
        return 0
    else:
        return fact[n] * (finv[r] * finv[n-r] % mod) % mod

mod = 10**9+7
fact,finv,inv = init_fact(200010,mod)






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

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

print(ans)

Solution (Python):
import math

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

facts = [1 for i in range(MAX)]
invs = [1 for i in range(MAX)]

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

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

def nCr(n,r):
	return facts[n] * invs[r] * invs[n-r] % mod

def main():
    ans = 0
    complement = nCr((N*M-2)%mod,(K-2)%mod)
    for i in range(1,N):
        ans += i * (N-i) * M**2
    for i in range(1,M):
        ans += i * (M-i) * N**2
    print((ans*complement) % mod)

main()

Solution (Python):
N, M, K = map(int, input().split())
import sys
mod = pow(10, 9) + 7
sys.setrecursionlimit(pow(10, 8))

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

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

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

def cmb(n, r):
    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 )

def f(n, m):
    x = m*m%mod
    y = (n*n*(n-1)//2 - n*(n-1)*(2*n-1)//6) % mod
    return x*y%mod

print(mul(f(N, M) + f(M, N), cmb(N*M-2, K-2)))

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

mod = 10 ** 9 + 7
MAX = n * m

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

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


def comb(n, k):
    return fact[n] * inv[n-k] * inv[k] % mod


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

        ans += add
        ans %= mod

ans *= comb(MAX - 2, k - 2)
ans %= mod
print(ans)

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


class ModInt:
    def __init__(self, num):
        self.num = num

    def __str__(self):
        return str(self.num)

    def __repr__(self):
        return "ModInt(num: {}, mod: {}".format(self.num, mod)

    def __add__(self, other):
        ret = self.num + other.num
        ret %= mod
        return ModInt(ret)

    def __sub__(self, other):
        ret = self.num - other.num
        ret %= mod
        return ModInt(ret)

    def __mul__(self, other):
        ret = self.num * other.num
        ret %= mod
        return ModInt(ret)

    def pow(self, times):
        pw = pow(self.num, times, mod)
        return ModInt(pw)

    def inverse(self):
        return ModInt(inv[self.num])

    def __truediv__(self, other):
        num = self * other.inverse()
        return ModInt(num)


def comb(n, k):
    return fact[n] * inv[n-k] * inv[k]


fact = [None] * (MAX + 1)
fact[0] = ModInt(1)
for i in range(1, MAX + 1):
    fact[i] = fact[i-1] * ModInt(i)

inv = [None] * (MAX + 1)
inv[MAX] = fact[MAX].pow(mod - 2)
for i in range(MAX, 0, -1):
    inv[i-1] = inv[i] * ModInt(i)

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

        ans += add

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

Solution (Python):
class ModInt:
    def __init__(self, num, mod):
        self.num = num
        self.mod = mod

    def __str__(self):
        return str(self.num)

    def __repr__(self):
        return "ModInt(num: {}, mod: {}".format(self.num, self.mod)

    def __add__(self, other):
        ret = self.num + other.num
        ret %= self.mod
        return ModInt(ret, self.mod)

    def __sub__(self, other):
        ret = self.num - other.num
        ret %= self.mod
        return ModInt(ret, self.mod)

    def __mul__(self, other):
        ret = self.num * other.num
        ret %= self.mod
        return ModInt(ret, self.mod)

    def pow(self, times):
        pw = pow(self.num, times, self.mod)
        return ModInt(pw, self.mod)

    def inverse(self):
        return self.pow(self.mod - 2)

    def __truediv__(self, other):
        num = self * other.inverse()
        return ModInt(num, self.mod)


class Combination:
    def __init__(self, n, mod):
        self.mod = mod
        self.fact = [ModInt(1, mod)]
        self.inverse = [ModInt(1, mod)]
        for i in range(1, n + 1):
            self.fact.append(self.fact[-1] * ModInt(i, mod))
            self.inverse.append(self.inverse[-1] * ModInt(i, mod).inverse())

    def comb(self, n, k):
        if k < 0 or n < k:
            return ModInt(0, self.mod)

        return self.fact[n] * self.inverse[k] * self.inverse[n-k]


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

cmb = Combination(2 * 10 ** 5, mod)
ans = ModInt(0, mod)
for i in range(n):
    for j in range(m):
        add = ModInt(n - i, mod) * ModInt(m - j, mod) * ModInt(i + j, mod)
        if i != 0 and j != 0:
            add *= ModInt(2, mod)

        ans += add

ans *= cmb.comb(n * m - 2, k - 2)
print(ans)

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

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

def function(i,k):
    return -((k*(k+1)*(2*k+1))//3)+k*(k+1)*(N+i-M)+2*N*(k+1)*(M-i)

N,M,K=map(int,input().split())
if M>N:
    N,M=M,N
ans=0
for i in range(1,N+M-1):
    if M-1>=i:
        ans+=i*(function(i,i)-N*(M-i)-M*(N-i))
        ans=ans%mod
    if i>=M and N-1>=i:
        ans+=i*(function(i,i)-function(i,i-M)-M*(N-i))
        ans=ans%mod
    if i>=N:
        ans+=i*(function(i,N-1)-function(i,i-M))
        ans=ans%mod

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

Solution (Python):

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

    def cmb(self, n, r):
        if n > 10 ** 6:
            return self.cmbl(n, r)
        return self.cmbr(n, r)

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

    def cmbl(self, n, r):
        t = 1
        r = min(r, n-r)
        for i in range(n - r + 1, n + 1):
            t = t * i % self.mod
        return t * self.g2[r] % self.mod

def main():
    N, M, K = map(int, input().split())
    mod = 10**9 + 7
    c = cmbs(mod)
    r = c.cmb((N*M)-2, K-2)
    r2 = 0
    for i in range(N):
        r2 += ((i + 1) * i + (N - i - 1) * (N-i))
    r2 *= pow(M, 2, mod)
    r2 %= mod
    r3 = 0
    for i in range(M):
        r3 += ((i + 1) * i + (M - i - 1) * (M-i))
    r3 *= pow(N, 2, mod)
    r3 %= mod
    r = r * (r2 + r3) * c.inverse[4]
    return r % mod
print(main())

Solution (Python):
class cmbs(object):
    def __init__(self, mod):
        self.mod = mod
        self.g1 = [1, 1]
        self.g2 = [1, 1]
        inverse = [0, 1]
        for i in range(2, 10 ** 6 + 1):
            self.g1.append((self.g1[-1] * i) % mod)
            inverse.append((-inverse[mod % i] * (mod // i)) % mod)
            self.g2.append((self.g2[-1] * inverse[-1]) % mod)

    def cmb(self, n, r):
        if n > 10 ** 6:
            return self.cmbl(n, r)
        return self.cmbr(n, r)

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

    def cmbl(self, n, r):
        t = 1
        r = min(r, n-r)
        for i in range(n - r + 1, n + 1):
            t = t * i % self.mod
        return t * self.g2[r] % self.mod

def main():
    N, M, K = map(int, input().split())
    mod = 10**9 + 7
    c = cmbs(mod)
    r = c.cmb((N*M)-2, K-2)
    r2 = 0
    for i in range(N):
        r2 += i * (N - i)
        
    r2 *= pow(M, 2, mod)
    r2 %= mod
    r3 = 0
    for i in range(M):
        r3 += i * (M-i)
        
    r3 *= pow(N, 2, mod)
    r3 %= mod
    r = r * (r2 + r3) 
    return r % mod
print(main())

Solution (Python):
from itertools import combinations

def solve(N, M, K):
    
    r = 0
    for x in range(N):
        r += x * (N - x) * M * M
    for y in range(M):
        r += y * (M - y) * N * N
    return r

def main():
    N, M, K = map(int, input().split())
    g1 = [1, 1]
    g2 = [1, 1]
    inverse = [0, 1]
    mod = 10 ** 9 + 7
    for i in range(2, N * M):
        g1.append((g1[-1] * i) % mod)
        inverse.append((-inverse[mod % i] * (mod // i)) % mod)
        g2.append((g2[-1] * inverse[-1]) % mod)
    t = solve(N, M, 2)
    for k in range(2, K):
        t = t * (N * M - k) * inverse[k - 1] % mod
    print(t)

main()

Solution (Python):
"""

AtCoder Beginner Contest 127  E - Cell Distance


x,yは独立して考えられるので分ける。
差分の絶対値の和の問題は
ある項について、＋として働くのが何回か、－として働くのが何回か
を考えてその分を足してあげればいい

例1の 2 x 2で考えると、x=2,y=1のについて、
xが＋に働くのは、x=2より小さい、(x,y)=(1,1),(1,2)
xが±0なのは、xが同じ、(x,y)=(2,2)
xが－に働くのはなし（x=2までなので）
的な。



Kについて考えると、ある駒a,b間のコストが求められる回数は、nCr(N*M-2, K-2)
（全部でK個選ぶときにa,bを先に選んで、残りのN-2からK-2個選ぶから。）

N x Mのマスで、(ax,ay)の物のコストについて、
X：+ax*(ax - 1) * N - ax * (M - ax) * N = ax*N*(ax-1-M+ax) = ax*N*(2ax-1-M)
これがマスの高さの分だけある(ayが1~Nまで)&これがnCr(N*M-2, K-2)回
Yも同様
"""


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)