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-----
2 2 2

-----Sample Output-----
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.

Difficulty: interview

Starter code:

Solutions:
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 = list(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)


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

def __starting_point():
    ans = main()
    print(ans)
__starting_point()

# python 3.4.3
import sys
input = sys.stdin.readline

# -------------------------------------------------------------
# library
# -------------------------------------------------------------
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)


# -------------------------------------------------------------
# main
# -------------------------------------------------------------

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)

n,m,k=list(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))



mod = 10 ** 9 + 7
class nCr():
    def __init__(self, n):
        self.n = n
        self.fa = [1] * (self.n + 1)
        self.fi = [1] * (self.n + 1)
        for i in range(1, self.n + 1):
            self.fa[i] = self.fa[i - 1] * i % mod
            self.fi[i] = pow(self.fa[i], mod - 2, mod)
    def comb(self, n, r):
        if n < r:return 0
        if n < 0 or r < 0:return 0
        return self.fa[n] * self.fi[r] % mod * self.fi[n - r] % mod
n, m, k = map(int, input().split())
ans = 0
c = nCr(n * m)
for i in range(n):
    ans += i * (n - i) * m * m
    ans %= mod
for i in range(m):
    ans += i * (m - i) * n * n
    ans %= mod
print(ans * c.comb(n * m - 2, k - 2) % mod)

MOD = 10**9+7

def comb_mod(n,r,mod):
    a,b=1,1
    r = min(r,n-r)
    for i in range(r):
        a = a*(n-i) % mod
        b = b*(r-i) % mod
    return a * pow(b, mod-2, mod) % mod
    

N,M,K = map(int,input().split())
    
ans = 0
for dx in range(1,M):
    ans += N*N % MOD * (M-dx) % MOD * dx % MOD
    ans %= MOD
for dy in range(1,N):
    ans += M*M % MOD * (N-dy) % MOD * dy % MOD
    ans %= MOD
    
ans *= comb_mod(N*M-2, K-2, MOD)
ans %= MOD
print(ans)

n, m, k = list(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(list(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)
    # x
    ans += x*m*cmb1(n*m-2, k-2, mod)*x
    ans -= (n-1-x)*m*cmb1(n*m-2, k-2, mod)*x
    # y
    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)


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)

import sys
input = lambda : sys.stdin.readline().rstrip()
mod = 10**9 + 7

def make_tables(n, mod = 10 ** 9 + 7):
    fac = [1, 1] # 階乗テーブル
    finv = [1, 1] #逆元の階乗テーブル
    inv = [0, 1] #逆元テーブル

    for i in range(2, n + 1):
        fac.append((fac[-1] * i) % mod)
        inv.append((mod -inv[mod % i] * (mod // i)) % mod)
        finv.append((finv[-1] * inv[-1]) % mod)
    return fac, finv

def nCk(n, k, mod = 10 ** 9 + 7):
    k = min(k, n-k)
    return fac[n] * finv[k] * finv[n-k] % mod

##########################

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

fac, finv = make_tables(n*m - 2)

ans = 0

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

print(ans)







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)

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


def __starting_point():
    main()

__starting_point()

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


mod = 10 ** 9 + 7
N = 10 ** 6

#逆元テーブル
inv_t = [0]+[1]
for i in range(2, N):
  inv_t += [inv_t[mod % i] * (mod - int(mod / i)) % mod]

#階乗計算
kai = [1, 1]
rev_kai = [1, inv_t[1]]
for i in range(2, N):
	kai.append(kai[-1] * i % mod)
	rev_kai.append(rev_kai[-1] * inv_t[i] % mod)

# コンビネーション計算
def cmb(n,r):
    res = 1
    for i in range(r):
        res *= n * inv_t[r]
        res %= mod
        n -= 1
        r -= 1
    return res

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

ans = 0
ans += m * m * cmb(n*m -2, k-2) * calc(n)
ans %= mod

ans += n * n * cmb(n*m -2, k-2) * calc(m)
ans %= mod

print(ans)


n,m,k = map(int, input().split())
n_ = n*m
mod = 10**9+7
fun = [1]*(n_+1)
for i in range(1,n_+1):
    fun[i] = fun[i-1]*i%mod
rev = [1]*(n_+1)
rev[n_] = pow(fun[n_],mod-2,mod)
for i in range(n_-1,0,-1):
    rev[i] = rev[i+1]*(i+1)%mod

def nCr(n,r):
    if r > n:
        return 0
    return fun[n]*rev[r]%mod*rev[n-r]%mod
  
x = sum(i*(n-i) for i in range(n))%mod
y = sum(i*(m-i) for i in range(m))%mod
ans = (x*m**2 + y*n**2)%mod*nCr(n*m-2,k-2)%mod
print(ans)

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)

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)

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

fac = [1] * n*m
finv = [1] * n*m
inv = [1] * n*m

ans = 0
def COMinit():
    for i in range(2, n*m):
        fac[i] = fac[i - 1] * i % mod
        inv[i] = mod - inv[mod % i] * (mod // i) % mod
        finv[i] = finv[i - 1] * inv[i] % mod


def COM(n, k):
    if n < 0 or k < 0 or n < k:
        return 0
    return fac[n] * (finv[k] * finv[n - k] % mod) % mod
COMinit()
for i in range(1, n):
    ans += (COM(n*m - 2, k-2)) * ((n-i)*i * m**2)%mod
for j in range(1, m):
    ans += (COM(n*m - 2, k-2)) * ((m-j)*j * n**2)%mod
print(ans%mod)

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

# cmbが10**10くらいだけど求められるか？って感じ
# 問題読み違えていた。。。N*M <= 2*10**5だ。。。まあ普通だ。
# もうライブラリ使おう。昔nCrのrが小さい時の工夫とかあったけど今回は大丈夫だ。

# https://ikatakos.com/pot/programming_algorithm/number_theory/mod_combination
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

## シグマを分離すること。２点i,jの選び方数、かけることのi,jの距離とする。あとはこれに２点以外の選び方のパターン数をかければ良い
# X,Yを独立して考える
# 例えばXについて考えると、使う列の選び方がN*2, 使う行の選び方はM-dis(1ならM-1だし、M-1離れてるのは１個しか選べない）
# disは距離が答えに寄与するから
# cmb()はその２点以外の選び方の個数の話。
#  もちろん他のi,jについて見たら同じ盤面何度も見てるけど、今のi,jに関する値しか計算してないからOK
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)

# python 3.4.3
import sys
input = sys.stdin.readline

# -------------------------------------------------------------
# library
# -------------------------------------------------------------
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)

# -------------------------------------------------------------
# main
# -------------------------------------------------------------

N,M,K = list(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)


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

def make_array_for_comb(N, mod=10**9+7):
    fact = [1,1]
    fact_inv = [1,1]
    inv = [0,1]
    for i in range(2, N+1):
        fact.append((fact[-1]*i) % mod)
        inv.append((-inv[mod%i] * (mod//i)) % mod)
        fact_inv.append((fact_inv[-1]*inv[i]) % mod)
    return fact, fact_inv

def comb(n, r, mod=10**9+7):
  if (r < 0) or (n < r):
      return 0
  r = min(r, n - r)
  return fact[n] * fact_inv[r] * fact_inv[n-r] % mod

fact, fact_inv = make_array_for_comb(n*m,mod)

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

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

print(ans*l%mod)

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)

N,M,K=map(int,input().split())
MOD=10**9+7
ans=0
x,y=1,1
k=min(K-2,N*M-K)
for i in range(k):
    x*=N*M-2-i
    y*=k-i
    x%=MOD
    y%=MOD
comb=x*pow(y,MOD-2,MOD)
comb%=MOD
for i in range(1,N):
    a=i*comb*(N-i)*M**2
    ans+=a
    ans%=MOD
for i in range(1,M):
    a=i*comb*(M-i)*N**2
    ans+=a
    ans%=MOD
print(ans%MOD)

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)

n, m, k = list(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)


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)

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
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)
ans=0
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*=cmb(N*M-2,K-2,mod)
print(ans%mod)

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


# E - Cell Distance

n, m, k = list(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)


# 1 つの組を固定したとき、その他の置き方は (N*M-2)C(K-2) 通り
# x 座標が h 離れる組み合わせは、(N-h)*M*M 通り (h!=0)

def combination(n, r, mod=10**9+7):
    # nCr mod m
    # rがn/2に近いと非常に重くなる
    n1, r = n+1, min(r, n-r)
    numer = denom = 1
    for i in range(1, r+1):
        numer = numer * (n1-i) % mod
        denom = denom * i % mod
    return numer * pow(denom, mod-2, mod) % mod

def main():
    N, M, K = list(map(int, input().split()))
    mod = 10**9 + 7
    ans1 = 0
    for h in range(1, N):
        ans1 += (N - h) * h
    ans1 = ans1 * M * M % mod
    ans2 = 0
    for w in range(1, M):
        ans2 += (M - w) * w
    ans2 = ans2 * N * N % mod
    ans = (ans1 + ans2) * combination(N*M-2, K-2) % mod
    print(ans)

main()


MOD = 10 ** 9 + 7

N, M, K = list(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)


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


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

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


def __starting_point():
    main()

__starting_point()

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

# 数え上げ

N, M, K = list(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


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

# Y
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)


import sys
input = sys.stdin.readline

class Comb:
    def __init__(self, n, mod=pow(10, 9)+7, f=True, i=False):
        self.n = n
        self.mod = mod
        self.factorial_table = None
        self.inverse_table = None
        if f: self.get_factorial()
        if i: self.get_inverse()
        
    def get_factorial(self):
        self.factorial_table = [1]*(self.n+1)
        for i in range(1, self.n+1):
            self.factorial_table[i] = self.factorial_table[i-1]*i%self.mod
        self.inverse_table = None
    
    def get_inverse(self):
        mod = self.mod
        self.inverse_table = [1]*(self.n+1)
        for i in range(self.n+1):
            self.inverse_table[i] = pow(self.factorial_table[i], -1, mod)
    
    def factorial(self, n):
        if self.factorial_table == None:
            return None
        return self.factorial_table[n]
    
    def inverse(self, n):
        if self.inverse_table == None:
            return None
        return self.inverse_table[n]   
    
    def permutation(self, n, k):
        if self.factorial_table == None:
            return None
        mod = self.mod
        if self.inverse_table == None:
            return self.factorial_table[n]*pow(self.factorial_table[n-k], -1, mod)%mod
        return self.factorial_table[n]*self.inverse_table[n-k]%mod
            
    def comb(self, n, k):
        mod = self.mod
        if self.factorial_table == None:
            res = 1
            for i in range(k):
                res *= n-i
                res *= i+1
                res %= mod
            return res
        if self.inverse_table == None:
            res = self.factorial_table[n]*pow(self.factorial_table[n-k], -1, mod)%mod\
                                    *pow(self.factorial_table[k], -1, mod)%mod
        else:
            res = self.factorial_table[n]*self.inverse_table[n-k]%mod\
                                    *self.inverse_table[k]%mod
        return res
    
    def recomb(self, n, k):
        if n+k-1 > self.n:
            return None
        return self.comb(n+k-1, k)
        


def main():
    n, m, k = list(map(int, input().split()))
    mod = pow(10, 9) + 7
    
    comb = Comb(n*m, i=True)
    
    ans = 0
    
    for i in range(n):
        ans += (i) * (n-i) * m**2
        ans %= mod
        
    
    for i in range(m):
        ans += (i) * (m-i) * n**2
        ans %= mod
    
    
    
    ans *= comb.comb(n*m-2, k-2)
    ans %= mod
    
    print(ans)
    
    
def __starting_point():
    main()


__starting_point()

def ext_euc(a, b):
    # return x, y, gcd(a, b) s.t. ax + by = gcd(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)


def __starting_point():
    main()
__starting_point()

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)

n,m,k=map(int,input().split())
#n*mからk個選ぶ
#x方向
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)

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)

MOD = 10 ** 9 + 7
n, m, k = list(map(int, input().split()))
def f(i, j):
    return i*(i-1)*(i+1)//6 * j ** 2
def c(a, b):
    r = 1
    for i in range(a, a - b, -1):
        r = r * i % MOD
    t = 1
    for i in range(1, b + 1):
        t = t * i % MOD
    return r * pow(t, MOD - 2, MOD)
print(((f(n, m) + f(m, n)) * c(n*m - 2, k - 2) % MOD))


N, M, K = list(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  # N は必要分だけ用意する
fact = [1, 1]  # fact[n] = (n! mod p)
factinv = [1, 1]  # factinv[n] = ((n!)^(-1) mod p)
inv = [0, 1]  # factinv 計算用
 
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)
#一つの配置でK C 2本の線
#配置はMN C K 通り
#よって、MNCK * KC2 / MNC2 = (MN - 2)C(K - 2)回
#各2頂点間のコストは足される。
#答えは全ての2頂点間のコストの和 * (MN - 2)C(K - 2)

#全ての2頂点間のコストの和を求める
#距離 i 離れた二点の個数は (1 + 2 + ... N - i) * N
#つまり、N * (N - i + 1) * (N - i) / 2
#縦横あるので、 それぞれ足す
ans1 = 0
for i in range(1, M):
  ans1 += i * (M - i) % mod
  #print(ans1)
ans1 = (ans1 * N * N) % mod

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

#print(ans1, ans2)

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



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)

"""

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 = list(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


#X：+ax*(ax - 1) * N - ax * (M - ax) * N = ax*N*(ax-1-M+ax) = ax*N*(2ax-1-M)がN個＆nCr(N*M-2, K-2))回
ans = 0

for i in range(1,M+1):
    tmp = (((((i * N) % MOD) * (2*i - 1 - M)) % MOD) * N) % MOD
    
    #print("~~~", i*N, 2*i - 1 - M,tmp, nCr(N*M-2, K-2))

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

    tmp %= MOD
    #print("pppp", tmp)
    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)


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)

ma = lambda :map(int,input().split())
lma = lambda :list(map(int,input().split()))
tma = lambda :tuple(map(int,input().split()))
ni = lambda:int(input())
yn = lambda fl:print("Yes") if fl else print("No")
import collections
import math
import itertools
import heapq as hq
import sys
import numpy as np
mod=10**9+7
def ncr_1(n,r,p):##　n>>r の場合のncr　O(r)
    x = 1
    y = 1
    for i in range(r):
        x*=(n-i)
        y*=(i+1)
        x%=p
        y%=p
    return x*modinv(y,p) %p
def modinv(a,p):
    b,u,v = p,1,0
    while b:
        t = a//b
        a -= t*b
        a,b = b,a
        u -=t*v
        u,v = v,u
    u %=p
    return u
n,m,k = ma()
ans=0
for i in range(1,n):
    ans+=i*(n-i)*m**2
    ans%=mod
for i in range(1,m):
    ans+=i*(m-i)*n**2
    ans%=mod
print((ans*ncr_1(n*m-2,k-2,mod))%mod)


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

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

# N行, M列 i-j
# 1-1, 1-2, 1-3, 1-4, ... , 1-M
# 2-1, 2-2, 2-3, 2-4, ... , 2-M
# 3-1. 3-2. 3-3, 3-4, ... , 3-M
# ...
# N-1, N-2, N-3, N-4, ... , N-M

# K個の駒が上記の座表上に置かれた時の配置全種類のコストの総和を調べたい

factorial = [1 for i in range(N*M + 1)]
for i in range(1, N * M + 1):
	if i==1:
		factorial[i] = 1
	else:
		factorial[i] = factorial[i-1]*i % mod

def comb(n, k):
	return (factorial[n]*pow(factorial[n-k]*factorial[k], -1, mod)) % mod

xscore = 0
for i in range(1, N):
	xscore += i * (N-i)
xscore %= mod
xscore *= M**2 * comb(N*M-2, K-2)
xscore %= mod

yscore = 0
for i in range(1, M):
	yscore += i * (M-i)
yscore %= mod
yscore *= N**2 * comb(N*M-2, K-2)
yscore %= mod

print(((xscore + yscore) % mod))


# python 3.4.3
import sys
input = sys.stdin.readline

# -------------------------------------------------------------
# library
# -------------------------------------------------------------
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)

# -------------------------------------------------------------
# main
# -------------------------------------------------------------

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)

import sys
readline = sys.stdin.readline
sys.setrecursionlimit(10**8)
mod = 10**9+7
#mod = 998244353
INF = 10**18
eps = 10**-7

m,n,k = list(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)



H, W, K = list(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)


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)
def __starting_point():
    main()
__starting_point()

import sys
readline = sys.stdin.readline
sys.setrecursionlimit(10**8)
mod = 10**9+7
#mod = 998244353
INF = 10**18
eps = 10**-7

m,n,k = list(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)



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)

n, m, k = [int(i) for i in input().split()]
mod = 10**9+7


def kai(x, p=mod):
    a = 1
    for i in range(1, x + 1):
        a *= i
        a %= p
    return a


def comb(a, b, p=mod):
    if a < 0 or b < 0:
        return 0
    elif a < b:
        return 0
    c = 1
    c *= kai(a, p)
    c *= pow(kai(b, p), p - 2, p)
    c *= pow(kai(a - b, p), p - 2, p)
    return c % p


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

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


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



def main():
    import sys
    input = sys.stdin.readline
    sys.setrecursionlimit(10**7)
    from collections import Counter, deque
    from itertools import combinations, permutations, accumulate, groupby, product
    from bisect import bisect_left,bisect_right
    from heapq import heapify, heappop, heappush
    import math
    #from math import gcd

    #inf = 10**17
    mod = 10**9 + 7

    # max_nCrまで調べられる
    max_n = 4*10**5
    fac, inv = [1]*(max_n+1), [0]*(max_n+1)
    for i in range(2, max_n+1):
        fac[i] = fac[i-1] * i % mod
    inv[-1] = pow(fac[-1], mod-2, mod)
    for i in range(max_n, 0, -1):
        inv[i-1] = inv[i] * i % mod
    
    # nCrを求める
    def ncr(n, r):
        return fac[n]*inv[r]*inv[n-r]%mod

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

    base = ncr(n*m-2, k-2)

    res = 0
    # x方向
    for i in range(1, m):
        res += n**2 * i * (m-i)
        res %= mod

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

    print(res*base%mod)

def __starting_point():
    main()
__starting_point()

n, m, k = list(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)


import math

mod = 10**9+7
N,M,K = list(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()


MOD = 10 ** 9 + 7

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

MOD = 10 ** 9 + 7
N = n * m + 10
fact = [0 for _ in range(N)]
invfact = [0 for _ in range(N)]
fact[0] = 1
for i in range(1, N):
    fact[i] = fact[i - 1] * i % MOD

invfact[N - 1] = pow(fact[N - 1], MOD - 2, MOD)

for i in range(N - 2, -1, -1):
    invfact[i] = invfact[i + 1] * (i + 1) % MOD
def nCk(n, k):
    if k < 0 or n < k:
        return 0
    else:
        return fact[n] * invfact[k] * invfact[n - k] % MOD
        
ans = 0
for i in range(n-1, 0, -1):
    ans += m * m * (i + 1) * i // 2 
    ans %= MOD
for i in range(m - 1, 0, -1):
    ans += n * n * (i + 1) * i // 2
    ans %= MOD

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

import sys
readline = sys.stdin.readline
sys.setrecursionlimit(10**8)
mod = 10**9+7
#mod = 998244353
INF = 10**18
eps = 10**-7

m,n,k = list(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)


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)

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

mod = 10**9 + 7

# (n * m - 2)C(k - 2)
# 距離xの組 -  n - x + 1
# 横 n - x + 1
# 縦 m * m
ans = 0
for y in range(1, n):
    ans += y * (n - y) * m * m % mod
for x in range(1, m):
    ans += x * (m - x) * n * n % mod


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


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

print(ans)


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)

# E - Cell Distance
# https://atcoder.jp/contests/abc127/tasks/abc127_e

def comb(n,k,mod):
    if n<k:
        return 0
    if n<0 or k<0:
        return 0
    k=min(n-k,k)
    ans=1
    inv=[1]*(k+1)
    if k>=1:
        ans*=(n-k+1)%mod
    for i in range(2,k+1):
        inv[i]=mod-inv[mod%i]*(mod//i)%mod
        ans=ans*(n-k+i)*inv[i]%mod
    return ans

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

MOD = 1000000007

distX = 0
for d in range(n):
    distX += d * (n - d) * m ** 2

distY = 0
for d in range(m):
    distY += d * (m - d) * n ** 2

ans = ((distX + distY) * comb(n * m - 2, k - 2, MOD)) % MOD

print(ans)

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)

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)

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

factorial=[1 for i in range(N*M+1)]
for i in range(1,N*M+1):
    if i==1:factorial[i]=1
    else:factorial[i] = factorial[i-1]*i % mod

def comb(n,k):
    return (factorial[n]*pow(factorial[n-k]*factorial[k], -1, mod))%mod
  
xscore = 0
for i in range(1,N):
  xscore += i * (N-i)
xscore %= mod
xscore *= M**2 * comb(N*M-2,K-2)
xscore %= mod

yscore = 0
for i in range(1,M):
  yscore += i * (M-i)
yscore %= mod
yscore *= N**2 * comb(N*M-2,K-2)
yscore %= mod

print((xscore + yscore) %mod)

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)

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)

from functools import reduce
n, m, k = list(map(int, input().split()))
mod = 10**9 + 7


def nCr(n, r, MOD=10**9 + 7):
    from functools import reduce
    if r == 0:
        return 1
    num = reduce(lambda x, y: x * y % MOD, list(range(n, n - r, -1)))
    den = reduce(lambda x, y: x * y % MOD, list(range(1, r + 1)))
    return num * pow(den, MOD - 2, MOD) % MOD


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


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


MOD = 10**9 + 7

N, M, K = list(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)


# inverse x^(-1)
def inv(x):
    nonlocal mod
    return pow(x, mod - 2, mod)

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

# combination nCr
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)

N,M,K=map(int,input().split())
c=[0 for i in range(N*M)]
c[1]=1
c[0]=1
mod=10**9+7
for i in range(1,N*M-1):
  c[i+1]=(i+1)*c[i]%mod
def comb(a,b):
  return c[a]*pow(c[b],mod-2,mod)*pow(c[a-b],mod-2,mod)%mod
ans=0
for i in range(N):
  ans+=i*(N-i)*pow(M,2)
  ans%=mod
for i in range(M):
  ans+=i*(M-i)*pow(N,2)%mod
  ans%=mod
ans*=comb(N*M-2,K-2)

ans%=mod
print(ans)

# -*- coding: utf-8 -*-

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 list(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
        # powに第三引数入れると冪乗のmod付計算を高速にやってくれる
        inverse[MAX-1] = pow(factorial[MAX-1], MOD-2, MOD)
        for i in range(MAX-2, 0, -1):
            # 最後から戻っていくこのループならMAX回powするより処理が速い
            inverse[i] = inverse[i+1] * (i+1) % MOD
        self.fact = factorial
        self.inv = inverse
    
    def nCr(self, n, r):
        """ 組み合わせの数 (必要な階乗と逆元のテーブルを事前に作っておく) """

        if n < r: return 0
        # 10C7 = 10C3
        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))


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)


def xgcd(a, b):
    x0, y0, x1, y1 = 1, 0, 0, 1
    while b != 0:
        q, a, b = a // b, b, a % b
        x0, x1 = x1, x0 - q * x1
        y0, y1 = y1, y0 - q * y1
    return x0

def modinv(a, m):
    x = xgcd(a, m)
    return x % m


max_num = 10 ** 9 + 7
res = 0
n, m, k = list(map(int, input().strip().split()))

for d in range(1, n):
    res += (n - d) * d * m * m
    res %= max_num
for d in range(1, m):
    res += (m - d) * d * n * n
    res %= max_num

r = min(n*m-2, k-2)
for i in range(1, r+1):
    res *= n * m - 1 - i
    res %= max_num
for i in range(1, r+1):
    res *= modinv(i, max_num)
    res %= max_num

print(res)


#!/usr/bin/env python3
def main():
    N, M, K = list(map(int, input().split()))
    m = 10**9+7
    # 方針: ２つの駒の距離だけ考えれば良い、それ以外K-2個の駒の位置が(N*M-2)C(K-2)通りある
    # マンハッタン距離はx方向とy方向を独立に扱える
    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

    # 答えに(NM-2)C(K-2)かける
    # K回しか処理しないから間に合う
    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)

def __starting_point():
    main()

__starting_point()

mod = 10 ** 9 + 7


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

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

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

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


def __starting_point():
    N, M, K = map(int, input().split())
    comb = Combination(N * M + 10, mod)

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

    for d in range(1, M):
        tmp = (M - d) * N * N % mod
        ans += d * tmp * comb(N * M - 2, K - 2)
        ans %= mod

    print(ans)
__starting_point()

def combination_mod(n,r,MOD=10**9+7):
    n_=1
    for i in range(1,n+1):
        n_=(n_*i)%MOD
    r_=1
    for i in range(1,r+1):
        r_=(r_*i)%MOD
    nr_=1
    for i in range(1,n-r+1):
        nr_=(nr_*i)%MOD

    power_r=pow(r_,MOD-2,MOD)
    power_nr=pow(nr_,MOD-2,MOD)

    return (n_*power_r*power_nr)%MOD


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

    MOD=10**9+7
    num=combination_mod(N*M-2,K-2)

    col_ans=0
    row_ans=0

    for i in range(1,M):
        col_ans+=pow(N,2)*(M-i)*num*i

    for i in range(1,N):
        row_ans+=pow(M,2)*(N-i)*num*i

    ans=(col_ans+row_ans)%MOD

    print(ans)

def __starting_point():
    main()


__starting_point()

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)

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 = list(map(int, input().split()))
ans = 0
# 2点の差がxになる組み合わせは同じ行でM-x通り、2点の行の組み合わせはN*N通り
for x in range(M):
    ans += x * (M-x) * N*N
for y in range(N):
    ans += y * (N-y) * M*M

# 2点以外(K-2)をN*M-2のマスに配置する組み合わせはnCr(N*M-2, K-2)通り
ans *= nCr(N*M-2, K-2, MOD)
ans %= MOD

print(ans)


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)

def __starting_point():
    main()
__starting_point()

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)

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 # N >= n
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)

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
#print(res)

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)


N,M,K = list(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


#X：+ax*(ax - 1) * N - ax * (M - ax) * N = ax*N*(ax-1-M+ax) = ax*N*(2ax-1-M)
ans = 0

for i in range(1,M+1):
    tmp = (((((i * N) % MOD) * (2*i - 1 - M)) % MOD) * N) % MOD
    
    #print("~~~", i*N, 2*i - 1 - M,tmp, nCr(N*M-2, K-2))

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

    tmp %= MOD
    #print("pppp", tmp)
    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)


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


# コンビネーション、さらに高速。あらかじめO(N)の計算をすることでのちの計算が早くなる
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
g1=[1,1] # g1[i]=i! % mod　:階乗
g2=[1,1] # g2[i]=(i!)^(-1) % mod　:階乗の逆元
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)

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

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)

def comb(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
N = 10 ** 6  # N は必要分だけ用意する
fact = [1, 1]  # fact[n] = (n! mod p)
factinv = [1, 1]  # factinv[n] = ((n!)^(-1) mod p)
inv = [0, 1]  # factinv 計算用
 
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)
############################################################

N, M, K = map(int, input().split())
MOD = 10**9+7
sumx = 0
sumy = 0

for d in range(1, N):
    cnt = (N-d)*M**2
    cnt %= MOD
    cnt *= comb(N*M-2, K-2, MOD)
    sumx += d*cnt
    sumx %= MOD

for d in range(1, M):
    cnt = (M-d)*N**2
    cnt %= MOD
    cnt *= comb(N*M-2, K-2, MOD)
    sumy += d*cnt
    sumy %= MOD
print((sumx + sumy)%MOD)



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



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

def nCk(n, k, mod):
    nu = 1
    de = 1 
    for i in range(1, k+1):
        nu *= (n-i+1)
        de *= i
        nu = nu % mod
        de = de % mod
    ans = nu * pow(de, mod-2, mod) % mod
    return ans

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

##mod,nについての制約は自分でいかに記入する
mod=10**9+7
def find_power(n,mod):
    # 0!から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
 #あるかずxのn乗をmodで割ったあまりを返す関数pow_cal(x,n,mod)
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):
    #0!から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=list(map(int,input().split()))
ans=M*N*(M+N)*(M*N-1)//6

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



MOD = 10**9 + 7

def MOD_perm(n, r):
    rtn = 1
    for _ in range(r):
        rtn *= n
        rtn %= MOD
        n -= 1
    return rtn

def MOD_inv(a):
    b = MOD
    u = 1
    v = 0
    while b > 0:
        t = a // b
        a -= t * b
        a, b = b, a
        u -= t * v
        u, v = v, u
    u = u % MOD
    if u < 0:
        u += MOD
    return u

def MOD_COM2(n, r):
    return MOD_perm(n, r) * MOD_inv(MOD_perm(r, r)) % MOD

def main():
    n, m, k = map(int, input().split())
    comb = MOD_COM2(n*m-2, k-2)
    ans = 0
    for i in range(1, n+1):
        U = (i-1)*i//2
        U *= m*m
        U *= comb
        U %= MOD
        ans += U
        ans %= MOD
    for j in range(1, m+1):
        L = (j-1)*j//2
        L *= n*n
        L *= comb
        L %= MOD
        ans += L
        ans %= MOD
    print(ans%MOD)

def __starting_point():
    main()
__starting_point()

H, W, K = list(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
# Hについて
for d in range(1, H):
    ans += d * (H - d) * pow(W, 2, MOD) * nCr(H * W - 2, K - 2) % MOD
    ans %= MOD
# Wについて
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))


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): #縦の距離がi
        if i == 0:
            for j in range(1, M): #横の距離がj
                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

def __starting_point():
    solve()
__starting_point()

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)

import sys
input = sys.stdin.readline
import copy

def main():
    N,M,K = list(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))

def __starting_point():
    main()

__starting_point()

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


def __starting_point():
    main()

# import sys
#
# sys.setrecursionlimit(10 ** 7)
#
# input = sys.stdin.readline
# rstrip()
# int(input())
# map(int, input().split())

__starting_point()

#!python3

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

# input
N, M, K = LI()

MOD = 10 ** 9 + 7
MAX = 10 ** 5 * 2
fac, finv, inv = [None] * MAX, [None] * MAX, [None] * MAX


def comb_init():
    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] = -inv[MOD%i] * int(MOD / i) % MOD
        finv[i] = finv[i - 1] * inv[i] % MOD


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


def s(x):
    return x * (x + 1) // 2


def cal(a, b):
    ret = a * s(b - 1) + b * s(a - 1)
    return ret % MOD


def main():
    comb_init()
    ans = 0
    for i in range(N):
        for j in range(M):
            v = cal(N - i, M - j) + cal(N - i, j + 1) - (s(j) + s(N - i - 1))
            ans = (ans + v) % MOD

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


def __starting_point():
    main()

__starting_point()