AtCoder Beginner Contest 127 - Cell Distance

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: MEDIUM

Starter code:

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


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)


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)


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


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


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


(N, M, K) = map(int, input().split())
cost = 0
mod = 10 ** 9 + 7
NM = N * M
kaijo = [0] * (NM - 1)
kaijo[0] = 1
for i in range(1, NM - 1):
	kaijo[i] = kaijo[i - 1] * i % mod
gyaku = [0] * (NM - 1)
gyaku[NM - 2] = pow(kaijo[NM - 2], mod - 2, mod)
for i in range(NM - 2, 0, -1):
	gyaku[i - 1] = gyaku[i] * i % mod

def conb(a, b):
	return kaijo[a] * gyaku[b] * gyaku[a - b] % mod
x = conb(NM - 2, K - 2)
for d in range(1, N):
	cost = (cost + x * d * (N - d) * M ** 2) % mod
for d in range(1, M):
	cost = (cost + x * d * (M - d) * N ** 2) % mod
print(cost)


(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
NM = N * M
kaijo = [0] * (NM + 1)
gyaku = [0] * (NM + 1)
kaijo[0] = 1
for i in range(1, NM + 1):
	kaijo[i] = kaijo[i - 1] * i % mod
gyaku[NM] = pow(kaijo[NM], mod - 2, mod)
for i in range(NM, 0, -1):
	gyaku[i - 1] = gyaku[i] * i % mod
C = kaijo[NM - 2] * gyaku[K - 2] * gyaku[NM - K] % mod
ans = 0
for k in range(1, M):
	x = k * (M - k) * C * N * N % mod
	ans = (ans + x) % mod
for k in range(1, N):
	x = k * (N - k) * C * M * M % mod
	ans = (ans + x) % mod
print(ans)


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


p = 10 ** 9 + 7

def comb(n, k, p):
	from math import factorial
	if n < 0 or k < 0 or n < k:
		return 0
	if n == 0 or k == 0:
		return 1
	a = factorial(n) % p
	b = factorial(k) % p
	c = factorial(n - k) % p
	return a * power_func(b, p - 2, p) * power_func(c, p - 2, p) % p

def power_func(a, b, p):
	if b == 0:
		return 1
	if b % 2 == 0:
		d = power_func(a, b // 2, p)
		return d * d % p
	if b % 2 == 1:
		return a * power_func(a, b - 1, p) % p
(N, M, K) = map(int, input().split())

def keisan(a, b):
	foo = a * (a ** 2 + 3 * a + 2) // 6
	foo *= b ** 2
	foo = foo % p
	return foo
ans = (keisan(N - 1, M) + keisan(M - 1, N)) % p
foo = comb(N * M - 2, K - 2, p)
bar = 9 * 17
bar = bar % p
ans = ans * foo % p
print(ans)


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)


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


(N, M, K) = map(int, input().split())
L = N * M + 10
mod = 10 ** 9 + 7
table = [1] * (L + 3)
t = 1
for i in range(1, L + 3):
	t *= i
	t %= mod
	table[i] = t
rtable = [1] * (L + 3)
t = 1
for i in range(1, L + 3):
	t *= pow(i, mod - 2, mod)
	t %= mod
	rtable[i] = t
ans = 0
t = table[N * M - 2] * rtable[K - 2] * rtable[N * M - K]
for i in range(1, M):
	ans += N * N * (M - i) * i
	ans %= mod
for i in range(1, N):
	ans += M * M * (N - i) * i
	ans %= mod
print(ans * t % mod)


(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())
P = 1000000007
import math
print(M * N * (M + N) * (M * N - 1) // 6 * math.factorial(N * M - 2) * pow(math.factorial(K - 2), P - 2, P) * pow(math.factorial(N * M - K), P - 2, P) % P)


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


(N, M, K) = map(int, input().split())
P = 1000000007
MPL = 1
f = 1
for i in range(1, N * M - 1):
	MPL = MPL * i % P
	MPL %= P
	if i == K - 2:
		f *= pow(MPL, P - 2, P)
	if i == N * M - K:
		f *= pow(MPL, P - 2, P)
MPL = MPL * f % P
ans = 0
for A in [N, M]:
	width = N ^ M ^ A
	for i in range(1, A):
		ans = (ans + width * width * i * (i + 1) // 2) % P
print(ans * MPL % P)


(H, W, K) = map(int, input().split())
MOD = 10 ** 9 + 7
Hd = sum([d * (H - d) * W ** 2 for d in range(H)])
Wd = sum([d * (W - d) * H ** 2 for d in range(W)])
factorial = [1]
inverse = [1]
for i in range(1, H * W + 2):
	factorial.append(factorial[-1] * i % MOD)
	inverse.append(pow(factorial[-1], MOD - 2, MOD))
ans = (Hd + Wd) * factorial[H * W - 2] * inverse[K - 2] * inverse[H * W - K] % MOD
print(ans)


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

def cmb(n, r):
	r = min(r, n - r)
	if r == 0:
		return 1
	if r == 1:
		return n
	f = 1
	for k in range(n, n - r, -1):
		f = f * k % MOD
	d = 1
	for k in range(2, r + 1):
		d = d * k % MOD
	return f * inv(d) % MOD
MOD = 10 ** 9 + 7
(N, M, K) = map(int, input().split())
ans = 0
for i in range(M):
	ans += i * (M - i) * N ** 2
for i in range(N):
	ans += i * (N - i) * M ** 2
pat = cmb(N * M - 2, K - 2)
ans = ans * pat % MOD
print(ans)


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

def cmb(n, r, mod):
	if r < 0 or r > n:
		return 0
	r = min(r, n - r)
	return g1[n] * g2[r] * g2[n - r] % mod
mod = 10 ** 9 + 7
MAX_N = 10 ** 6
g1 = [1, 1]
g2 = [1, 1]
inverse = [0, 1]
for i in range(2, MAX_N + 1):
	g1.append(g1[-1] * i % mod)
	inverse.append(-inverse[mod % i] * (mod // i) % mod)
	g2.append(g2[-1] * inverse[-1] % mod)
k = cmb(M * N - 2, K - 2, mod)

def calc(p, q, d):
	return d * (p - d) * q ** 2
ans = 0
for i in range(N):
	ans += calc(N, M, i)
for i in range(M):
	ans += calc(M, N, i)
ans %= mod
ans *= k
ans %= mod
print(ans)


mod = 10 ** 9 + 7

def mod_combination(n, k, mod):

	def ext_gcd(a, b):
		if b == 0:
			return (a, 1, 0)
		else:
			(d, x, y) = ext_gcd(b, a % b)
			x -= a // b * y
		return (d, y, x)
	(p, q) = (1, 1)
	for i in range(n - k + 1, n + 1):
		p = p * i % mod
	for i in range(2, k + 1):
		q = q * i % mod
	return int(p * (ext_gcd(q, mod)[1] % mod) % mod)
(n, m, k) = map(int, input().split())
print((pow(m, 2, mod) * mod_combination(n * m - 2, k - 2, mod) * sum(((n - i) * i % mod for i in range(1, n))) % mod + pow(n, 2, mod) * mod_combination(n * m - 2, k - 2, mod) * sum(((m - i) * i % mod for i in range(1, m))) % mod) % mod)


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


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


def combnk_mod(n, k, p):
	if 2 * k > n:
		k = n - k
	x = 1
	y = 1
	for i in range(1, k + 1):
		x *= pow(i, p - 2, p)
		y *= n - i + 1
		x %= p
		y %= p
	return x * y % p
(n, m, k) = map(int, input().split())
MOD = 10 ** 9 + 7
p = combnk_mod(n * m - 2, k - 2, MOD)
(x, y) = (0, 0)
for i in range(2):
	for d in range(1, n):
		x += d * (n - d) * m ** 2
		x %= MOD
	(x, y) = (y, x)
	(n, m) = (m, n)
print(p * (x + y) % MOD)


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 += 1
		self.MAX = MAX
		self.MOD = MOD
		factorial = [1] * MAX
		factorial[0] = factorial[1] = 1
		for i in range(2, MAX):
			factorial[i] = factorial[i - 1] * i % MOD
		inverse = [1] * MAX
		inverse[MAX - 1] = pow(factorial[MAX - 1], MOD - 2, MOD)
		for i in range(MAX - 2, 0, -1):
			inverse[i] = inverse[i + 1] * (i + 1) % MOD
		self.fact = factorial
		self.inv = inverse

	def nCr(self, n, r):
		if n < r:
			return 0
		r = min(r, n - r)
		numerator = self.fact[n]
		denominator = self.inv[r] * self.inv[n - r] % self.MOD
		return numerator * denominator % self.MOD
(H, W, K) = MAP()
N = H * W
mt = ModTools(N, MOD)
ansh = 0
for i in range(1, H + 1):
	ansh += i * (H - i) * pow(W, 2, MOD)
	ansh %= MOD
ansh *= mt.nCr(N - 2, K - 2)
ansh %= MOD
answ = 0
for i in range(1, W + 1):
	answ += i * (W - i) * pow(H, 2, MOD)
	answ %= MOD
answ *= mt.nCr(N - 2, K - 2)
answ %= MOD
print((ansh + answ) % MOD)


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


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

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

def combi(n, r):
	if r < 0 or r > n:
		return 0
	else:
		return fact(n) * inv(fact(r)) * inv(fact(n - r)) % mod
mod = 10 ** 9 + 7
(N, M, K) = map(int, input().split())
x = (N - 1) * N * (N + 1) * inv(6) % mod
x = x * M * M % mod
y = (M - 1) * M * (M + 1) * inv(6) % mod
y = y * N * N % mod
z = combi(M * N - 2, K - 2)
ans = (x + y) * z % mod
print(ans)


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)


MOD = 10 ** 9 + 7
(N, M, K) = map(int, input().split())
NM = N * M
U = 2 * 10 ** 5
fact = [1] * (U + 1)
fact_inv = [1] * (U + 1)
for i in range(1, U + 1):
	fact[i] = fact[i - 1] * i % MOD

def comb(n, k):
	x = fact[n]
	x *= pow(fact[n - k], MOD - 2, MOD)
	x %= MOD
	x *= pow(fact[k], MOD - 2, MOD)
	x %= MOD
	return x

def F(N):
	x = N * (N + 1) // 2 * N
	y = N * (N + 1) * (2 * N + 1) // 6
	return (x - y) % MOD
MM = M * M % MOD
NN = N * N % MOD
answer = MM * F(N) % MOD + NN * F(M) % MOD
answer %= MOD
answer *= comb(N * M - 2, K - 2)
answer %= MOD
print(answer)


import sys
readline = sys.stdin.readline
sys.setrecursionlimit(10 ** 8)
mod = 10 ** 9 + 7
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)


MOD = 10 ** 9 + 7
(N, M, K) = map(int, input().split())
U = N * M
fact = [1] * (U + 1)
for i in range(1, U + 1):
	fact[i] = fact[i - 1] * i % MOD

def comb(n, k):
	x = fact[n]
	x *= pow(fact[n - k], MOD - 2, MOD)
	x %= MOD
	x *= pow(fact[k], MOD - 2, MOD)
	return x
res = 0
for i in range(1, N):
	res += i * M * M * (N - i) % MOD
	res %= MOD
for i in range(1, M):
	res += i * N * N * (M - i) % MOD
	res %= MOD
res *= comb(N * M - 2, K - 2)
res %= MOD
print(res)


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

def f(x):
	return math.factorial(x)

def c(x, y):
	return f(x) // (f(x - y) * f(y))
print((n ** 2 * (m - 1) * m * (m + 1) // 6 + m ** 2 * (n - 1) * n * (n + 1) // 6) * c(m * n - 2, k - 2) % (10 ** 9 + 7))


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


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)


def init():
	global factorial
	MOD = 10 ** 9 + 7
	factorial = [0] * 200001
	factorial[0] = 1
	for i in range(1, 200001):
		factorial[i] = factorial[i - 1] * i % MOD

def main():
	init()
	(n, m, k) = map(int, input().split())
	s = 0
	MOD = 10 ** 9 + 7
	for i in range(n):
		s += (n - 1 - i) * (n - i) // 2 % MOD * (m * m % MOD) % MOD
	for i in range(m):
		s += (m - 1 - i) * (m - i) // 2 % MOD * (n * n % MOD) % MOD
	print(s * (factorial[n * m - 2] * pow(factorial[k - 2], MOD - 2, MOD) % MOD) % MOD * pow(factorial[n * m - k], MOD - 2, MOD) % MOD)
main()


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)


(n, m, k) = map(int, input().split())
ans = 0
mod = 10 ** 9 + 7
max = 2 * 10 ** 5
g1 = [1, 1]
g2 = [1, 1]
inverse = [0, 1]
for i in range(2, max + 1):
	g1.append(g1[-1] * i % mod)
	inverse.append(-inverse[mod % i] * (mod // i) % mod)
	g2.append(g2[-1] * inverse[-1] % mod)

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
C = cmb(m * n - 2, k - 2, mod)
for i in range(1, n):
	ans += i * (n - i) * m ** 2 * C
for i in range(1, m):
	ans += i * (m - i) * n ** 2 * C
print(ans % mod)


mod = 10 ** 9 + 7

def comb(n, k):
	(x, y) = (1, 1)
	k = min(k, n - k)
	for i in range(k):
		x *= n - i
		x %= mod
		y *= i + 1
		y %= mod
	return x * pow(y, mod - 2, mod) % mod
(N, M, K) = map(int, input().split())
c = comb(N * M - 2, K - 2)
ans = c * (M ** 2 * (N ** 2 * (N + 1) // 2 - N * (N + 1) * (2 * N + 1) // 6) + N ** 2 * (M ** 2 * (M + 1) // 2 - M * (M + 1) * (2 * M + 1) // 6)) % mod
print(ans)


(n, m, k) = [int(i) for i in input().split()]
p = 10 ** 9 + 7
S = 0
for d in range(n):
	S += m ** 2 * d * (n - d)
	S %= p
for d in range(m):
	S += n ** 2 * d * (m - d)
	S %= p

def fact(n, p=10 ** 9 + 7):
	f = [1]
	for i in range(1, n + 1):
		f.append(f[-1] * i % p)
	return f

def invfact(n, f, p=10 ** 9 + 7):
	inv = [pow(f[n], p - 2, p)]
	for i in range(n, 0, -1):
		inv.append(inv[-1] * i % p)
	return inv[::-1]
f = fact(n * m - 2)
invf = invfact(n * m - 2, f)

def comb(a, b):
	if a < b:
		return 0
	if a < 0 or b < 0:
		return 0
	return f[a] * invf[b] * invf[a - b] % p
S *= comb(n * m - 2, k - 2)
print(S % p)


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


import math

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


MOD = 10 ** 9 + 7

def mod_combination(n, k):

	def extended_gcd(a, b):
		if b == 0:
			return (a, 1, 0)
		(d, x, y) = extended_gcd(b, a % b)
		return (d, y, x - a // b * y)
	(p, q) = (1, 1)
	for i in range(n - k + 1, n + 1):
		p = p * i % MOD
	for i in range(2, k + 1):
		q = q * i % MOD
	return p * (extended_gcd(q, MOD)[1] % MOD) % MOD

def mod_pow(a, b):
	return pow(a, b, MOD)
(N, M, K) = map(int, input().split())
X = mod_pow(M, 2) * (mod_pow(N, 3) - N) % MOD
Y = mod_pow(N, 2) * (mod_pow(M, 3) - M) % MOD
comb = mod_combination(N * M - 2, K - 2)
ans = (X + Y) * mod_pow(6, MOD - 2) * comb % MOD
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()


(N, M, K) = map(int, input().split())
MOD = 10 ** 9 + 7
R = 1
for i in range(K - 2):
	R = R * (N * M - 2 - i) * pow(i + 1, MOD - 2, MOD) % MOD
ans = 0
for i in range(1, N + 1):
	for j in range(1, M + 1):
		u = N - i + 1
		d = i
		l = j
		r = M - j + 1
		D = 0
		D += u * l * (u + l - 2) // 2
		D += u * r * (u + r - 2) // 2
		D += d * l * (d + l - 2) // 2
		D += d * r * (d + r - 2) // 2
		D -= u * (u - 1) // 2
		D -= d * (d - 1) // 2
		D -= l * (l - 1) // 2
		D -= r * (r - 1) // 2
		ans += D
ans = ans * R * pow(2, MOD - 2, MOD) % MOD
print(ans)


MOD = 10 ** 9 + 7
(n, m, k) = 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)


def getInv(N):
	inv = [0] * (N + 1)
	inv[0] = 1
	inv[1] = 1
	for i in range(2, N + 1):
		inv[i] = -(Q // i) * inv[Q % i] % Q
	return inv
Q = 10 ** 9 + 7
modinv = getInv(4 * 10 ** 5)
(N, M, K) = map(int, input().split())
T = N * M
comb = 1
for i in range(K - 2):
	comb *= T - 2 - i
	comb %= Q
	comb *= modinv[K - 2 - i]
	comb %= Q
ans = 0
for i in range(1, M):
	ans += N ** 2 * (M - i) * i
	ans %= Q
for i in range(1, N):
	ans += M ** 2 * (N - i) * i
	ans %= Q
ans %= Q
print(ans * comb % Q)


(H, W, K) = map(int, input().split())
MOD = 10 ** 9 + 7
factorial = [1, 1]
inverse = [1, 1]
invere_base = [0, 1]
for i in range(2, H * W + 2):
	factorial.append(factorial[-1] * i % MOD)
	invere_base.append(-invere_base[MOD % i] * (MOD // i) % MOD)
	inverse.append(inverse[-1] * invere_base[-1] % MOD)

def nCr(n, r):
	if not 0 <= r <= n:
		return 0
	return factorial[n] * inverse[r] * inverse[n - r] % MOD
ans = 0
for dh in range(1, H):
	ans += dh * (H - dh) * pow(W, 2, MOD) * nCr(H * W - 2, K - 2) % MOD
	ans %= MOD
for dw in range(1, W):
	ans += dw * (W - dw) * pow(H, 2, MOD) * nCr(H * W - 2, K - 2) % MOD
	ans %= MOD
print(ans)


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


def dist(arr, n):
	res = 0
	su = 0
	for i in range(n):
		res = (res + arr[i] * i - su) % mod
		su = (su + arr[i]) % mod
	return res % mod
(n, m, k) = map(int, input().split())
mod = int(1000000000.0 + 7)
fact = [1]
for i in range(1, n * m + 1):
	fact.append(fact[-1] * i % mod)
x = []
y = []
for i in range(1, n + 1):
	for j in range(1, m + 1):
		x.append(i)
		y.append(j)
y.sort()
tot = (dist(x, n * m) + dist(y, n * m)) % mod
N = n * m - 2
Need = k - 2
num = fact[N]
den = fact[Need] * fact[N - Need] % mod
ncr = num * pow(den, mod - 2, mod) % mod
ans = tot * ncr % mod
print(ans)


mod = 10 ** 9 + 7

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

def solve(string):
	(n, m, k) = map(int, string.split())
	base = n * m * (n + m) * (n * m - 1) // 6 % mod
	return str(base * comb(n * m - 2, k - 2) % mod)
print(solve(input()))


(N, M, K) = list(map(int, input().split()))
mod = 10 ** 9 + 7
A = N ** 2 * (N - 1) // 2 - N * (2 * N - 1) * (N - 1) // 6
A = A * M ** 2
A %= mod
B = 1
for i in range(1, K - 1):
	B *= N * M - 1 - i
	B %= mod
	B *= pow(i, mod - 2, mod)
	B %= mod
ans1 = A * B
(N, M) = (M, N)
A = N ** 2 * (N - 1) // 2 - N * (2 * N - 1) * (N - 1) // 6
A = A * M ** 2
A %= mod
B = 1
for i in range(1, K - 1):
	B *= N * M - 1 - i
	B %= mod
	B *= pow(i, mod - 2, mod)
	B %= mod
ans2 = A * B
print((ans1 + ans2) % mod)


from functools import reduce
from operator import add

def mod_mul(a, b):
	return a * b % MOD

def cmb(n, r):
	r = min(r, n - r)
	if r == 0:
		return 1
	if r == 1:
		return n
	f = reduce(mod_mul, range(n - r + 1, n + 1))
	d = reduce(mod_mul, range(2, r + 1))
	return f * pow(d, MOD - 2, MOD) % MOD
MOD = 10 ** 9 + 7
(N, M, K) = map(int, input().split())
ans_x = reduce(add, (i * (M - i) for i in range(M))) * N ** 2
ans_y = reduce(add, (i * (N - i) for i in range(N))) * M ** 2
ans = (ans_x + ans_y) % MOD
pat = cmb(N * M - 2, K - 2)
ans = ans * pat % MOD
print(ans)


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

def chs(a, b, md):
	res = 1
	for i in range(a - b + 1, a + 1):
		res *= i
		res %= md
	rev = 1
	for i in range(1, b + 1):
		rev *= i
		rev %= md
	rev2 = pow(rev, md - 2, md)
	res *= rev2
	res %= md
	return res
ans = 0
c = chs(n * m - 2, k - 2, mod)
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 *= c
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)


import math

def combinations_count(n, r):
	return math.factorial(n) // (math.factorial(n - r) * math.factorial(r))
(N, M, K) = list(map(int, input().split(' ')))
ans = 0
for i in range(0, N):
	for j in range(0, M):
		mul = 2
		if i == 0 or j == 0:
			mul = 1
		ans += mul * (N - i) * (M - j) * (i + j)
ans *= combinations_count(N * M - 2, K - 2)
ans %= 10 ** 9 + 7
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 Calc:

	def __init__(self, max_value, mod):
		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
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)


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


from math import factorial
(N, M, K) = map(int, input().split())
MOD = 10 ** 9 + 7
fact = [1] * (N * M + 1)
for i in range(2, N * M + 1):
	fact[i] = fact[i - 1] * i % MOD

def comb(n, r):
	r = n - r if n - r < r else r
	return fact[n] * pow(fact[r], MOD - 2, MOD) * pow(fact[n - r], MOD - 2, MOD) % MOD
p = comb(N * M - 2, K - 2)
ans = 0
for i in range(N):
	for j in range(M):
		if i == 0 and j == 0:
			continue
		d = i + j
		cnt = (N - i) * (M - j)
		if i != 0 and j != 0:
			cnt *= 2
		ans += d * cnt % MOD
ans = ans * p % MOD
print(ans)


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


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

def comb(n, r):
	if n < r:
		return 0
	return factorial[n] * inverse[r] * inverse[n - r] % mod
c = comb(n * m - 2, k - 2)
ans = 0
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)


from math import factorial

def comb(n, k, mod):
	if n < 0 or k < 0 or n < k:
		return 0
	if n == 0 or k == 0:
		return 1
	a = factorial(n) % mod
	b = factorial(k) % mod
	c = factorial(n - k) % mod
	return a * pow(b, mod - 2, mod) * pow(c, mod - 2, mod) % mod
(N, M, K) = map(int, input().split())
mod = 10 ** 9 + 7
x = N * N * (M * M * (M + 1) // 2 - M * (M + 1) * (2 * M + 1) // 6)
y = M * M * (N * N * (N + 1) // 2 - N * (N + 1) * (2 * N + 1) // 6)
ans = comb(N * M - 2, K - 2, mod) * (x + y) % mod
print(ans)


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

def comb_mod(n, k, mod):
	fact_n = factorial_mod(n, mod)
	fact_k = factorial_mod(k, mod)
	fact_n_k = factorial_mod(n - k, mod)
	return fact_n * pow(fact_k, mod - 2, mod) * pow(fact_n_k, mod - 2, mod) % mod
(n, m, k) = map(int, input().split())
s = 0
mod = 10 ** 9 + 7
for x in range(m):
	s += (m - x - 1) * (m - x) // 2 * n ** 2
for y in range(n):
	s += (n - y - 1) * (n - y) // 2 * m ** 2
print(s * comb_mod(n * m - 2, k - 2, mod) % mod)


import sys
readline = sys.stdin.readline
sys.setrecursionlimit(10 ** 8)
mod = 10 ** 9 + 7
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)


Mod = 10 ** 9 + 7
(n, m, k) = map(int, input().split())
inv = [0] * (k + 3)
inv[1] = 1
for i in range(2, k + 2):
	inv[i] = (Mod - Mod // i) * inv[Mod % i] % Mod
s = (n * n * (m + 1) * m * (m - 1) + m * m * (n + 1) * n * (n - 1)) // 6
for i in range(k - 2):
	s = s * (n * m - 2 - i) * inv[i + 1] % Mod
print(s)


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

def combi(x, y, p):
	c = 1
	for i in range(1, y + 1):
		c *= (x - i + 1) * pow(i, p - 2, p)
		c %= p
	return c
ans = 0
C = combi(N * M - 2, K - 2, ppp)
for i in range(1, N):
	ans += C * i * (N - i) * M ** 2
for j in range(1, M):
	ans += C * j * (M - j) * N ** 2
print(ans % ppp)


mod = 10 ** 9 + 7
(n, m, k) = map(int, input().split())
x = k
fact = [1] * (x + 1)
for i in range(1, x + 1):
	fact[i] = fact[i - 1] * i % mod

def mod_comb_k(n, k, mod):
	ret = 1
	if n - k < k:
		k = n - k
	for i in range(n, n - k, -1):
		ret = ret * i % mod
	return ret % mod * pow(fact[k], mod - 2, mod) % mod
res = 0
for d in range(1, n):
	res += m * m * (n - d) * d
	res = res % mod
(n, m) = (m, n)
for d in range(1, n):
	res += m * m * (n - d) * d
	res = res % mod
res *= mod_comb_k(n * m - 2, k - 2, mod)
print(res % mod)


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) = 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
MAX = int(200000.0)

def div(a, b):
	return a * pow(b, MOD - 2, MOD) % MOD
FACT = [1] * (MAX + 1)
for i in range(1, MAX + 1):
	FACT[i] = i * FACT[i - 1] % MOD
INV = [1] * (MAX + 1)
INV[MAX] = div(1, FACT[MAX])
for i in range(MAX, 0, -1):
	INV[i - 1] = INV[i] * i % MOD

def combination(r, c):
	return FACT[r] * INV[r - c] * INV[c]

def main():
	(N, M, K) = map(int, input().split())
	ans1 = 0
	for d in range(1, N):
		ans1 = (ans1 + d * (N - d) * M ** 2) % MOD
	ans2 = 0
	for d in range(1, M):
		ans2 = (ans2 + d * (M - d) * N ** 2) % MOD
	ans = (ans1 + ans2) * combination(M * N - 2, K - 2) % MOD
	print(ans)
main()


(n, m, k) = map(int, input().split())
mod = 10 ** 9 + 7
res_n = 0
for i in range(1, n + 1):
	res_n += (n + 1 - i) * (n - i) // 2
	res_n %= mod
res_n = res_n * m ** 2 % mod
res_m = 0
for i in range(1, m + 1):
	res_m += (m + 1 - i) * (m - i) // 2
	res_m %= mod
res_m = res_m * n ** 2 % mod
res = res_n + res_m
f = [1 for _ in range(n * m)]
for i in range(1, n * m):
	f[i] = f[i - 1] * i % mod
res = res * f[n * m - 2] * pow(f[k - 2], mod - 2, mod) * pow(f[n * m - k], mod - 2, mod) % mod
print(res)


def pmod(a, b, m):
	ret = 1
	while b > 0:
		if b & 1:
			ret = ret * a % m
		a = a * a % m
		b = b >> 1
	return ret

def fmod(n, m):
	ret = 1
	for i in range(n, 0, -1):
		ret = ret * i % m
	return ret

def cmod(n, r, m):
	ret = 1
	for i in range(r):
		ret = ret * (n - i) % m
	return ret * pmod(fmod(r, m), m - 2, m) % m
mod = 10 ** 9 + 7
(n, m, k) = map(int, input().split())
ans = 0
c = cmod(n * m - 2, k - 2, mod)
for i in range(1, n):
	ans += i * c * m ** 2 * (n - i)
	ans = ans % mod
for i in range(1, m):
	ans += i * c * n ** 2 * (m - i)
	ans = ans % mod
print(ans % mod)


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


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

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

def cmb(n, r, MOD):
	return facts[n] * invs[n - r] % mod * invs[r] % mod
ans = 0
for dis in range(M):
	ans += N ** 2 * (M - dis) * dis * cmb(N * M - 2, K - 2, mod)
	ans %= mod
for dis in range(N):
	ans += M ** 2 * (N - dis) * dis * cmb(N * M - 2, K - 2, mod)
	ans %= mod
print(ans)


mod = 10 ** 9 + 7
(n, m, k) = map(int, input().split())
nm = n * m
fact = [1] * (nm + 1)
for i in range(1, nm + 1):
	fact[i] = fact[i - 1] * i % mod

def C(n, k):
	x = fact[n]
	x *= pow(fact[n - k], mod - 2, mod)
	x %= mod
	x *= pow(fact[k], mod - 2, mod)
	x %= mod
	return x

def F(x):
	s1 = x * x * (x + 1) // 2
	s2 = x * (x + 1) * (2 * x + 1) // 6
	return (s1 - s2) % mod
res = (m ** 2 % mod * F(n) % mod + n ** 2 * F(m) % mod) % mod
res *= C(n * m - 2, k - 2)
res %= mod
print(res)


from itertools import accumulate

def inpl():
	return list(map(int, input().split()))
MOD = 10 ** 9 + 7
(N, M, K) = inpl()
X = sum([n * (n + 1) // 2 for n in range(1, N)]) * M * M % MOD
Y = sum([m * (m + 1) // 2 for m in range(1, M)]) * N * N % MOD

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]
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)
print((X + Y) * cmb(N * M - 2, K - 2, MOD) % MOD)


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

def pow(x, n):
	ret = 1
	while n > 0:
		if n & 1 == 1:
			ret = ret * x % mod
		x = x * x % mod
		n //= 2
	return ret
fac = [1]
inv = [1]
for i in range(1, n * m + 1):
	fac.append(fac[-1] * i % mod)
	inv.append(pow(fac[i], mod - 2))

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

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


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

def inv(x):
	y = 1
	while x != 1:
		y *= mod // x + 1
		y %= mod
		x = x - mod % x
	return y
others = 1
for i in range(k - 2):
	others *= h * w - 2 - i
	others *= inv(i + 1)
	others %= mod
hh = 0
for i in range(1, h):
	hh += i * (h - i)
	hh %= mod
ww = 0
for i in range(1, w):
	ww += i * (w - i)
	ww %= mod
ans = others * (h * h * ww + w * w * hh)
print(ans % mod)


def egcd(a, b):
	(x, lastx) = (0, 1)
	(y, lasty) = (1, 0)
	while b != 0:
		q = a // b
		(a, b) = (b, a % b)
		(x, lastx) = (lastx - q * x, x)
		(y, lasty) = (lasty - q * y, y)
	return (lastx, lasty, a)

def modinv(a):
	m = MOD
	(inv, q, gcd_val) = egcd(a, m)
	return inv % m
MOD = 10 ** 9 + 7
(N, M, K) = map(int, input().split())
b = 1
for i in range(N * M - 2, N * M - K, -1):
	b = b * i % MOD
for i in range(1, K - 1):
	b = b * modinv(i) % MOD
x = 0
for i in range(1, N):
	x = x + i * (N - i)
x = x * M * M % MOD
y = 0
for j in range(1, M):
	y = y + j * (M - j)
y = y * N * N % MOD
ans = b * (x + y) % 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)


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


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

def f(N):
	a = 1
	for i in range(1, N + 1):
		a *= i
		a %= mod
	return a % mod

def Comb(N, K):
	a = f(N)
	b = f(K)
	c = f(N - K)
	return a * pow(b, mod - 2, mod) * pow(c, mod - 2, mod) % mod
c = Comb(N * M - 2, K - 2)
ans = 0
for d in range(M):
	ans += d * (M - d) * N ** 2
	ans %= mod
for d in range(N):
	ans += d * (N - d) * M ** 2
	ans %= mod
ans *= c
ans %= mod
print(ans)


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)


import sys
sys.setrecursionlimit(2147483647)
INF = float('inf')
MOD = 10 ** 9 + 7
input = lambda : sys.stdin.readline().rstrip()

def modfact(n):
	fact = [1] * (n + 1)
	invfact = [1] * (n + 1)
	for i in range(1, n + 1):
		fact[i] = i * fact[i - 1] % MOD
	invfact[n] = pow(fact[n], MOD - 2, MOD)
	for i in range(n - 1, -1, -1):
		invfact[i] = invfact[i + 1] * (i + 1) % MOD
	return (fact, invfact)

def resolve():
	from math import factorial as fact
	(n, m, k) = map(int, input().split())
	(fact, invfact) = modfact(n * m)
	ans = 1
	ans *= fact[m * n - 2] * invfact[k - 2] * invfact[m * n - k]
	ans *= ((n - 1) * n * (n + 1) * m * m + (m - 1) * m * (m + 1) * n * n) // 6
	print(ans % MOD)
resolve()


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


MAX = 2 * 10 ** 5 + 100
MOD = 10 ** 9 + 7
fact = [0] * MAX
inv = [0] * MAX
finv = [0] * MAX
fact[0] = 1
fact[1] = 1
finv[0] = 1
finv[1] = 1
inv[1] = 1
for i in range(2, MAX):
	fact[i] = fact[i - 1] * i % MOD
	inv[i] = MOD - inv[MOD % i] * (MOD // i) % MOD
	finv[i] = finv[i - 1] * inv[i] % MOD

def C(n, r):
	if n < r:
		return 0
	if n < 0 or r < 0:
		return 0
	return fact[n] * (finv[r] * finv[n - r] % MOD) % MOD
(N, M, K) = map(int, input().split())
ans = 0
for i in range(1, N):
	ans += i * (N - i) * M * M * C(N * M - 2, K - 2)
	ans %= MOD
for i in range(1, M):
	ans += i * (M - i) * N * N * C(N * M - 2, K - 2)
	ans %= MOD
print(ans)


(N, M, K) = map(int, input().split())
MOD = 10 ** 9 + 7
a = 0
for i in range(N):
	j = N - 1 - i
	a += i * (i + 1) // 2 + j * (j + 1) // 2
b = 0
for i in range(M):
	j = M - 1 - i
	b += i * (i + 1) // 2 + j * (j + 1) // 2
MAX = N * M + 5
fac = [1, 1] + [0] * MAX
finv = [1, 1] + [0] * MAX
inv = [0, 1] + [0] * MAX
for i in range(2, MAX + 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 = (a * inv[N] ** 2 + b * inv[M] ** 2) * inv[N * M - 1] * (N * M)
ans *= K * (K - 1) // 2
ans *= ncr(N * M, K)
print(ans % MOD)


def cmb(n, r):
	return fact[n] * inv_fact[r] * inv_fact[n - r] % MOD
MOD = 10 ** 9 + 7
(N, M, K) = map(int, input().split())
fact = [1]
for i in range(1, N * M - 2 + 1):
	fact.append(fact[-1] * i % MOD)
inv_fact = [-1] * (N * M - 1)
inv_fact[-1] = pow(fact[-1], MOD - 2, MOD)
for i in range(N * M - 3, -1, -1):
	inv_fact[i] = inv_fact[i + 1] * (i + 1) % MOD
pat = cmb(N * M - 2, K - 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
ans = ans * pat % MOD
print(ans)


import math
(N, M, K) = (int(x) for x in input().split())

def comb(n, k, p):
	if k < 0 or k > n:
		return 0
	if n == 0 or k == 0:
		return 1
	a = math.factorial(n) % p
	b = math.factorial(k) % p
	c = math.factorial(n - k) % p
	return a * mypow(b, p - 2, p) * mypow(c, p - 2, p) % p

def mypow(a, b, p):
	if b == 0:
		return 1
	if b % 2 == 0:
		d = mypow(a, b // 2, p)
		return d * d % p
	if b % 2 == 1:
		return a * mypow(a, b - 1, p) % p
c_sum = (M * (N - 1) * N * (N + 1) * M + N * (M - 1) * M * (M + 1) * N) // 3
c_sum = c_sum // 2
c_sum *= comb(N * M - 2, K - 2, 1000000007)
c_sum = c_sum % 1000000007
print(c_sum)


(n, m, k) = list(map(int, input().split()))
mod = 10 ** 9 + 7
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)
for i in range(k - 2):
	ans = ans * (n * m - 2 - i) * inv(i + 1) % mod
print(ans)


import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines
from math import factorial

def main():
	(n, m, k) = map(int, readline().split())
	P = int(1000000000.0) + 7

	def C(n, m, p):
		a = 1
		for i in range(m):
			a *= n - i
			a %= p
		b = factorial(m) % p
		b = pow(b, p - 2, p)
		return a * b % p
	c = C(n * m - 2, k - 2, P)
	s_x = 0
	for d in range(1, m):
		s_x += d * (m - d) * n * n
		s_x %= P
	s_y = 0
	for d in range(1, n):
		s_y += d * (n - d) * m * m
		s_y %= P
	ans = (s_x + s_y) * c % P
	print(ans)
main()


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)


from math import factorial
MOD = 10 ** 9 + 7
(N, M, K) = map(int, input().split(' '))
(n, r) = (N * M - 2, K - 2)
P = factorial(n) * pow(factorial(n - r) % MOD, MOD - 2, MOD) * pow(factorial(r) % MOD, MOD - 2, MOD)
ans = 0
for n in range(N):
	for m in range(M):
		if n == 0 and m == 0:
			continue
		d = n + m
		cnt = (N - n) * (M - m) * (1 if n == 0 or m == 0 else 2)
		ans += d * cnt
ans = ans * P % MOD
print(ans)


from collections import deque
from heapq import heappush, heappop
import re

def int_raw():
	return int(input())

def ss_raw():
	return input().split()

def ints_raw():
	return list(map(int, ss_raw()))
INF = 1 << 29
DIV = 10 ** 9 + 7

def mod_inv(a, mod=DIV):
	return pow(a, mod - 2, mod)

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

def main():
	(N, M, K) = ints_raw()
	ans = 0
	for d in range(1, N):
		ans = (ans + d * (N - d) * M * M) % DIV
	for d in range(1, M):
		ans = (ans + d * (M - d) * N * N) % DIV
	return ans * ncr(N * M - 2, K - 2) % DIV
print(main())


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)


def comb(n, k, mod=10 ** 9 + 7):
	from functools import reduce
	from operator import mul
	k = min(k, n - k)

	def calc(x, y):
		return x * y % mod
	x = reduce(calc, range(n - k + 1, n + 1), 1) % mod
	y = reduce(calc, range(1, k + 1), 1) % mod
	return x * pow(y, mod - 2, mod) % mod

def solve(string):
	(n, m, k) = map(int, string.split())
	mod = 10 ** 9 + 7
	base = n * m * (n + m) * (n * m - 1) // 6 % mod
	return str(base * comb(n * m - 2, k - 2) % mod)
print(solve(input()))


(n, m, k) = map(int, input().split())
mod = 10 ** 9 + 7
max = n * m + 10
(fac, finv, inv) = ([0] * max, [0] * max, [0] * max)

def comInit(max):
	fac[0] = fac[1] = 1
	finv[0] = finv[1] = 1
	inv[1] = 1
	for i in range(2, max):
		fac[i] = fac[i - 1] * i % mod
		inv[i] = mod - inv[mod % i] * (mod // i) % mod
		finv[i] = finv[i - 1] * inv[i] % mod
comInit(max)

def com(n, k):
	if n < k:
		return 0
	if n < 0 | k < 0:
		return 0
	return fac[n] * (finv[k] * finv[n - k] % mod) % mod
ans = n * n * (m - 1) * m * (m + 1) // 6 + m * m * (n - 1) * n * (n + 1) // 6
ans = ans % mod
ans = com(n * m - 2, k - 2) * ans % mod
print(ans)


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


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


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


(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
fact = [1, 1]
factinv = [1, 1]
inv = [0, 1]
for i in range(2, n + 1):
	fact.append(fact[-1] * i % p)
	inv.append(-inv[p % i] * (p // i) % p)
	factinv.append(factinv[-1] * inv[-1] % p)
ans1 = 0
for i in range(1, M):
	ans1 += i * (M - i) % mod
ans1 = ans1 * N * N % mod
ans2 = 0
for i in range(1, N):
	ans2 += i * (N - i) % mod
ans2 = ans2 * M * M % mod
ans = ans1 + ans2
print(ans * cmb(M * N - 2, K - 2, p) % mod)


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

def modinv(a):
	b = MOD
	(x, lastx) = (0, 1)
	(y, lasty) = (1, 0)
	while b != 0:
		q = a // b
		(a, b) = (b, a % b)
		(x, lastx) = (lastx - q * x, x)
		(y, lasty) = (lasty - q * y, y)
	return lastx % MOD
MOD = 10 ** 9 + 7
(N, M, K) = map(int, input().split())
x = 0
for i in range(1, N):
	x = x + i * (N - i)
x = x * M * M % MOD
y = 0
for j in range(1, M):
	y = y + j * (M - j)
y = y * N * N % MOD
ans = comb(N * M - 2, K - 2) * (x + y) % MOD
print(ans)


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


mod = 10 ** 9 + 7

def inv(x):
	if x == 0:
		raise ValueError('Div 0')
	(nowx, lastx) = (0, 1)
	(nowy, lasty) = (1, 0)
	a = x
	b = mod
	while b != 0:
		q = a // b
		(a, b) = (b, a % b)
		(nowx, lastx) = (lastx - q * nowx, nowx)
		(nowy, lasty) = (lasty - q * nowy, nowy)
	return a * lastx % mod
(N, M, K) = map(int, input().split())
X = sum([d * (N - d) * M * M for d in range(N)])
Y = sum([d * (M - d) * N * N for d in range(M)])
ans = X + Y
for i in range(1, K - 1):
	ans *= N * M - 1 - i
	ans %= mod
	ans *= inv(i)
	ans %= mod
print(ans)


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

def nCr(n, r):
	return math.factorial(n) // (math.factorial(r) * math.factorial(n - r)) % MOD
ans = 0
t = nCr(N * M - 2, K - 2)
for d in range(1, M):
	ans += d * (M - d) * (N * N) * t
	ans %= MOD
for d in range(1, N):
	ans += d * (N - d) * (M * M) * t
	ans %= MOD
print(ans)


MOD = 10 ** 9 + 7
MAX = 2 * 10 ** 5
(n, m, k) = map(int, input().split())
fact = [1] * (MAX + 1)
for i in range(1, MAX + 1):
	fact[i] = fact[i - 1] * i % MOD
finv = [1] * MAX + [pow(fact[-1], MOD - 2, MOD)]
for i in range(MAX, 0, -1):
	finv[i - 1] = finv[i] * i % MOD

def comb(n, r):
	return fact[n] * finv[n - r] % MOD * finv[r] % MOD
cumsum = [0] * (MAX + 1)
for i in range(1, max(n, m) + 1):
	cumsum[i] = cumsum[i - 1] + i
ans = 0
for i in range(n):
	ans += (cumsum[i] + cumsum[n - i - 1]) * m ** 2
for i in range(m):
	ans += (cumsum[i] + cumsum[m - i - 1]) * n ** 2
ans //= 2
ans = ans * comb(n * m - 2, k - 2) % MOD
print(ans)


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

def inv(x):
	y = 1
	while x != 1:
		y *= mod // x + 1
		y %= mod
		x -= mod % x
	return y
ans = 1
for i in range(k - 2):
	ans *= n * m - 2 - i
	ans %= mod
	ans *= inv(k - 2 - i)
	ans %= mod
zsum = [0]
for i in range(1, max(n, m)):
	zsum.append(zsum[-1] + i)
ans1 = 0
for i in range(m):
	ans1 += zsum[i]
	ans1 += zsum[m - 1 - i]
	ans1 %= mod
ans1 *= n ** 2
ans1 %= mod
ans2 = 0
for j in range(n):
	ans2 += zsum[j]
	ans2 += zsum[n - 1 - j]
ans2 *= m ** 2
ans2 %= mod
ans *= ans1 + ans2
ans *= inv(2)
ans %= mod
print(ans)


from math import factorial as f
(N, M, K) = map(int, input().split())
mod = 10 ** 9 + 7

def combi(n, r):
	return f(n) // (f(n - r) * f(r)) % mod
res = 0
c = combi(N * M - 2, K - 2)
for i in range(1, N):
	res += (N - i + 1) * (N - i) // 2 * c * M * M % mod
	res %= mod
for i in range(1, M):
	res += (M - i + 1) * (M - i) // 2 * c * N * N % mod
	res %= mod
print(res % 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()


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

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

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

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

def __starting_point():
	main()
__starting_point()


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

def mod_inv(a, m=MOD):

	def mod_inv_sub(a, m):
		if a == 1:
			return (1, 0)
		else:
			(d, r) = (m // a, m % a)
			(x, y) = mod_inv_sub(r, a)
			return (y - d * x, x)
	x = mod_inv_sub(a, m)[0]
	if x < 0:
		x += (-x // m + 1) * m
	return x
factor = 1
L = min(K - 2, M * N - K)
for i in range(M * N - L - 1, M * N - 1):
	factor = factor * i % MOD
for i in range(1, L + 1):
	factor = factor * mod_inv(i) % MOD
ans = (M * (N - 1) * (N + 1) + N * (M - 1) * (M + 1)) * M * N // 6 * factor % MOD
print(ans)


def inverse(a, mod):
	while a < 0:
		a += mod
	a %= mod
	p = mod
	ya = 0
	yp = 1
	while a > 0:
		n = p // a
		r = p % a
		yr = ya - n * yp
		ya = yp
		yp = yr
		p = a
		a = r
	return ya % mod if ya >= 0 else ya + mod
mod = 1000000007
(n, m, k) = map(int, input().split())
inv = 1
for i in range(1, k - 1):
	inv = inv * i % mod
inv = inverse(inv, mod)
combination = inv
for i in range(n * m - k + 1, n * m - 1):
	combination = combination * i % mod
xn = n * (n + 1) * (n - 1) // 6 % mod * m % mod * m % mod
yn = m * (m + 1) * (m - 1) // 6 % mod * n % mod * n % mod
ans = (xn + yn) * combination % mod
print(ans)


import sys
from math import factorial
input = sys.stdin.readline
Mod = 10 ** 9 + 7

def Cmb(n, k):
	return factorial(n) // (factorial(n - k) * factorial(k)) % Mod
(n, m, k) = list(map(int, input().split()))
ans = 0
for i in range(1, n + 1):
	ans += i * (n - i) * m ** 2 % Mod
for i in range(1, m + 1):
	ans += i * (m - i) * n ** 2 % Mod
ans *= Cmb(n * m - 2, k - 2)
ans %= Mod
print(ans)


(N, M, K) = map(int, input().split())
mod = 10 ** 9 + 7
fact = [1] * (2 * N * M + 1)
for n in range(1, 2 * N * M + 1):
	fact[n] = n * fact[n - 1] % mod

def inv(n):
	n %= mod
	return pow(n, mod - 2, mod)
ans = M * (M ** 2 - 1) * N ** 2 % mod + N * (N ** 2 - 1) * M ** 2 % mod
ans *= inv(6) * fact[N * M - 2] * inv(fact[K - 2] * fact[N * M - K])
ans %= mod
print(ans)


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


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


(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) = [int(_) for _ in 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)

def count(n, m):
	ret = 0
	for d in range(1, n):
		ret += d * (n - d) * m ** 2
	return ret
ans = cmb(N * M - 2, K - 2, mod)
ans *= count(N, M) + count(M, N)
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
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)
(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)


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

def makefactlist(n, mod=MOD):
	ret = [1]
	for i in range(n):
		ret += [ret[-1] * (i + 1) % MOD]
	return ret
factlist = makefactlist(N * M, mod=MOD)
comb = factlist[N * M - 2] * pow(factlist[N * M - K], MOD - 2, MOD) * pow(factlist[K - 2], MOD - 2, MOD)
comb = comb % MOD
ret = 0
for x in range(1, N):
	ret += (N - x) * M * M * comb * x % MOD
	ret = ret % MOD
for y in range(1, M):
	ret += (M - y) * N * N * comb * y % MOD
	ret = ret % MOD
print(ret)


from math import factorial as f
(n, m, k) = map(int, input().split())
print(f(n * m - 2) // (f(k - 2) * f(n * m - k)) * (m * m * (n ** 3 - n) + n * n * (m ** 3 - m)) // 6 % (10 ** 9 + 7))


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


(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 = int(1000000000.0) + 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)


MOD = 10 ** 9 + 7
UPPERLIMIT = 3 * 10 ** 5
MODMUL = [1, 1] + [0] * (UPPERLIMIT - 1)
for i in range(2, UPPERLIMIT + 1):
	MODMUL[i] = MODMUL[i - 1] * i % MOD
MODDIV = [1] * UPPERLIMIT + [pow(MODMUL[-1], MOD - 2, MOD)]
for i in range(UPPERLIMIT, 0, -1):
	MODDIV[i - 1] = MODDIV[i] * i % MOD

def MODCOMB(n, r):
	if n < r or n < 1:
		return 0
	else:
		return MODMUL[n] * MODDIV[r] % MOD * MODDIV[n - r] % MOD
(N, M, K) = map(int, input().split())
ans = 0
k = MODCOMB(N * M - 2, K - 2)
for sa in range(1, N):
	ans = (ans + sa * (N - sa) * M * M * k) % MOD
for sa in range(1, M):
	ans = (ans + sa * (M - sa) * N * N * k) % MOD
print(ans)


import sys
from collections import deque
from collections import defaultdict
import math
sys.setrecursionlimit(20000000)
input = sys.stdin.readline
(n, m, k) = map(int, input().split())
mod = 10 ** 9 + 7
fact = [1] * (n * m + 10000 + 1)
rfact = [1] * (n * m + 10000 + 1)
for i in range(n * m + 10000):
	fact[i + 1] = rr = (i + 1) * fact[i] % mod

def comb(N, K, mod):
	return fact[N] * pow(fact[K], mod - 2, mod) * pow(fact[N - K], mod - 2, mod) % mod
ans = 0
ni = comb(n * m - 2, k - 2, mod)
for i in range(n):
	for j in range(m):
		ko = (1 + j) * j * n // 2
		ko += (m - j) * (m - 1 - j) * n // 2
		ko += (1 + i) * i * m // 2
		ko += (n - i) * (n - 1 - i) * m // 2
		ans += ni * ko
		ans %= mod
print(ans * pow(2, mod - 2, mod) % mod)


mod = 10 ** 9 + 7
MAX = 2 * 10 ** 5 + 100
g1 = [1, 1]
g2 = [1, 1]
for i in range(2, MAX + 1):
	num_1 = g1[-1] * i % mod
	g1.append(num_1)
	g2.append(pow(num_1, mod - 2, mod))

def cmb(n, r):
	return g1[n] * g2[r] * g2[n - r] % mod
(N, M, K) = map(int, input().split())
ans = 0
for t in range(1, M):
	ans = (ans + t * N ** 2 * (M - t)) % mod
for t in range(1, N):
	ans = (ans + t * M ** 2 * (N - t)) % mod
print(ans * cmb(N * M - 2, K - 2) % mod)


(N, M, K) = map(int, input().split())
(O, S, T) = (10 ** 9 + 7, N + M, 6)
for i in range(K):
	S = S * (N * M - i) % O
	T = T * -~i % O
print(S * K * ~-K * pow(T, O - 2, O) % O)


(N, M, K) = map(int, input().split())
MOD = 10 ** 9 + 7
L = N * M
fact = [1] * (L + 1)
rfact = [1] * (L + 1)
r = 1
for i in range(1, L + 1):
	fact[i] = r = r * i % MOD
rfact[L] = r = pow(fact[L], MOD - 2, MOD)
for i in range(L, 0, -1):
	rfact[i - 1] = r = r * i % MOD
L = max(N, M)
D = [0] * (L + 1)
for i in range(L):
	D[i + 1] = D[i] + (i + 1)
ans = 0
for i in range(N):
	for j in range(M):
		ans += D[j] * (i + 1) + D[M - 1 - j] * i + D[i] * M
	ans %= MOD
print(ans * fact[N * M - 2] * rfact[N * M - K] * rfact[K - 2] % MOD)


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


def factorials(M):
	INF = int(1000000000.0 + 7)
	fac = [0] * (M + 1)
	finv = [0] * (M + 1)
	inv = [0] * (M + 1)
	fac[0] = fac[1] = 1
	finv[0] = finv[1] = 1
	inv[1] = 1
	for i in range(2, M + 1):
		fac[i] = fac[i - 1] * i % INF
		inv[i] = INF - inv[INF % i] * (INF // i) % INF
		finv[i] = finv[i - 1] * inv[i] % INF
	return (fac, finv)

def nCr(n, r):
	INF = int(1000000000.0 + 7)
	if n < r or r < 0 or n < 0:
		return 0
	return fac[n] * (finv[r] * finv[n - r] % INF) % INF
(N, M, K) = map(int, input().split())
MOD = 10 ** 9 + 7
(fac, finv) = factorials(N * M)
ans = 0
for d in range(1, N):
	ans += d * nCr(N * M - 2, K - 2) * ((N - d) * M ** 2 % MOD)
	ans %= MOD
for d in range(1, M):
	ans += d * nCr(N * M - 2, K - 2) * ((M - d) * N ** 2 % MOD)
print(ans % MOD)


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

def combi(x, y):
	c = 1
	m = min(x - y, y)
	for iii in range(y):
		c = c * (x - iii) * pow(y - iii, M - 2, M) % M
	return c
rem_combi = combi(n * m - 2, k - 2)

def calc(num):
	out = 0
	for iii in range(num):
		out += iii * (num - iii) % M
	return out
x_cost = m * m * calc(n) % M
y_cost = n * n * calc(m) % M
ans = rem_combi * (x_cost + y_cost) % M
print(int(ans))


import math
(n, m, k) = map(int, input().split())
mod = 10 ** 9 + 7
s = math.factorial(n * m - 2) // (math.factorial(n * m - k) * math.factorial(k - 2))
cnt = 0
s %= mod
a = m * m * s % mod
for d in range(1, n):
	cnt += d * (n - d) * a
	cnt %= mod
a = n * n * s % mod
for d in range(1, m):
	cnt += d * (m - d) * a
	cnt %= mod
print(cnt)


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

def comb(n, k):
	ret = 1
	for i in range(1, k + 1):
		ret = ret * (n - i + 1) % MOD
		ret = ret * pow(i, MOD - 2, MOD) % MOD
	return ret
ans = 0
for d in range(N):
	ans += d * (N - d) * M ** 2
	ans %= MOD
for d in range(M):
	ans += d * (M - d) * N ** 2
	ans %= MOD
print(ans * comb(N * M - 2, K - 2) % MOD)


def cmb(n, r, mod):
	if r < 0 or r > n:
		return 0
	r = min(r, n - r)
	return g1[n] * g2[r] * g2[n - r] % mod
mod = 10 ** 9 + 7
N = 2 * 10 ** 5
g1 = [1, 1]
g2 = [1, 1]
inverse = [0, 1]
for i in range(2, N + 1):
	g1.append(g1[-1] * i % mod)
	inverse.append(-inverse[mod % i] * (mod // i) % mod)
	g2.append(g2[-1] * inverse[-1] % mod)
(N, M, K) = map(int, input().split())
a = cmb(N * M - 2, K - 2, mod)
b = N ** 2 * M * (M + 1) * (M - 1) // 6 % (10 ** 9 + 7)
c = M ** 2 * N * (N + 1) * (N - 1) // 6 % (10 ** 9 + 7)
print(a * (b + c) % (10 ** 9 + 7))


mod = 10 ** 9 + 7

class Combination:

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


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)


(N, M, K) = map(int, input().split())
MOD = 10 ** 9 + 7
P = 200005
fact = [0] * P
fact[0] = 1
for i in range(1, P):
	fact[i] = fact[i - 1] * i
	fact[i] %= MOD
rfact = [0] * P
rfact[P - 1] = pow(fact[P - 1], MOD - 2, MOD)
for i in range(P - 2, -1, -1):
	rfact[i] = rfact[i + 1] * (i + 1)
	rfact[i] %= MOD

def comb(n, k):
	return fact[n] * rfact[k] * rfact[n - k] % MOD
num = 0
for d in range(1, N):
	num += d * (N - d) * M ** 2
	num %= MOD
for d in range(1, M):
	num += d * (M - d) * N ** 2
	num %= MOD
num *= comb(N * M - 2, K - 2)
print(num % MOD)


def cmb(n, r, mod):
	if r < 0 or r > n:
		return 0
	r = min(r, n - r)
	return g1[n] * g2[r] * g2[n - r] % mod
(N, M, k) = map(int, input().split())
m = 10 ** 9 + 7
g1 = [1, 1]
g2 = [1, 1]
inverse = [0, 1]
for i in range(2, N * M):
	g1.append(g1[-1] * i % m)
	inverse.append(-inverse[m % i] * (m // i) % m)
	g2.append(g2[-1] * inverse[-1] % m)
c = cmb(N * M - 2, k - 2, m)
ans = 0
for i in range(N):
	ans += (N - i) * M * M * i * c % m
	ans %= m
for i in range(M):
	ans += (M - i) * N * N * i * c % m
	ans %= m
print(ans % m)


MOD = 10 ** 9 + 7

def nCr(n, r, mod=10 ** 9 + 7):
	r = min(r, n - r)
	numer = denom = 1
	for i in range(1, r + 1):
		numer = numer * (n + 1 - i) % mod
		denom = denom * i % mod
	return numer * pow(denom, mod - 2, mod) % mod
(N, M, K) = map(int, input().split())
ans = 0
for x in range(M):
	ans += x * (M - x) * N * N
for y in range(N):
	ans += y * (N - y) * M * M
ans *= nCr(N * M - 2, K - 2, MOD)
ans %= MOD
print(ans)


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)


(N, M, K) = map(int, input().split())
P = 10 ** 9 + 7
import math
print(M * N * (M + N) * (M * N - 1) // 6 * math.factorial(N * M - 2) * pow(math.factorial(K - 2), P - 2, P) * pow(math.factorial(N * M - K), P - 2, P) % P)


from functools import reduce

def egcd(a, b):
	if a == 0:
		return (b, 0, 1)
	else:
		(g, y, x) = egcd(b % a, a)
		return (g, x - b // a * y, y)

def modinv(a, m):
	(g, x, y) = egcd(a, m)
	if g != 1:
		raise Exception('modular inverse does not exist')
	else:
		return x % m
MOD = 10 ** 9 + 7
(N, M, K) = map(int, input().split())
patterns = reduce(lambda a, b: a * b % MOD, range(N * M - K + 1, N * M + 1))
patterns *= modinv(reduce(lambda a, b: a * b % MOD, range(1, K + 1)), MOD)
patterns *= (N + M) * K * (K - 1)
patterns *= modinv(6, MOD)
patterns %= MOD
print(patterns)


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

def __starting_point():
	solve()
__starting_point()


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


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


(N, M, K) = map(int, input().split())
P = 10 ** 9 + 7
fact = [1]
for i in range(N * M):
	fact.append(fact[-1] * (i + 1) % P)

def inv(n):
	return pow(n, P - 2, P)

def C(a, b):
	return fact[a] * inv(fact[a - b] * fact[b]) % P
print((M ** 2 * (N ** 3 - N) + N ** 2 * (M ** 3 - M)) // 6 * C(M * N - 2, K - 2) % P)


(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) = map(int, input().split())
mod = 10 ** 9 + 7
table = [1]
for i in range(N * M):
	table.append(table[-1] * (i + 1) % mod)

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

def cmb(n, r):
	return table[n] * inv(table[n - r] * table[r]) % mod
print((M ** 2 * (N ** 3 - N) + N ** 2 * (M ** 3 - M)) // 6 * cmb(M * N - 2, K - 2) % mod)


(N, M, K) = map(int, input().split())
mod = 10 ** 9 + 7
X = N * M % mod
C = [1]
for i in range(K + 1):
	C.append(C[i] * (X - 2 - i) * pow(i + 1, mod - 2, mod) % mod)
P = 0
for i in range(N):
	P = (P + i * (i + 1) // 2 + (N - i) * (N - i - 1) // 2) % mod
P = P * M * M % mod
Q = 0
for i in range(M):
	Q = (Q + i * (i + 1) // 2 + (M - i) * (M - i - 1) // 2) % mod
Q = Q * N * N % mod
print((P + Q) * C[K - 2] * ((mod + 1) // 2) % mod)


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

def combination(n, m):
	c = 1
	m = min(n - m, m)
	for i in range(m):
		c = c * (n - i) * pow(m - i, MOD - 2, MOD) % MOD
	return c
pattern = combination(N * M - 2, K - 2)

def diff_x_pattern(val):
	ret = 0
	for d in range(1, val):
		ret += (val - d) * d % MOD
	return ret
ans = M * M * diff_x_pattern(N) % MOD
ans += N * N * diff_x_pattern(M) % MOD
print(ans * pattern % 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)


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

def comb(n, r):
	if n - r > 0 and n > 0 and (r > 0):
		return kaidan[n] * gyaku[n - r] * gyaku[r] % mod
	elif n == 0 or r == 0 or n == r:
		return 1
	else:
		0
ans = 0
for d in range(1, n):
	ans += d * (n - d) * m * m * comb(n * m - 2, k - 2)
for d in range(1, m):
	ans += d * (m - d) * n * n * comb(n * m - 2, k - 2)
print(ans % mod)


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


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


(N, M, K) = map(int, input().split())
mod = 10 ** 9 + 7
n = 2 * 10 ** 5
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 comb(n, k):
	return fact[n] * rfact[k] * rfact[n - k] % mod
base = comb(N * M - 2, K - 2)
first = 0
for i in range(N):
	for j in range(M):
		first += i + j
dp = [0] * M
dp[0] = first
for i in range(1, M):
	dp[i] = dp[i - 1] - (M - 2 * i) * N
SUM = sum(dp)
dp2 = [0] * N
dp2[0] = SUM
for i in range(1, N):
	dp2[i] = dp2[i - 1] - (N - 2 * i) * M ** 2
SUM = sum(dp2)
ans = comb(N * M - 2, K - 2) * (SUM % mod) * pow(2, mod - 2, mod)
ans %= mod
print(ans)


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)


def solve():
	(N, M, K) = map(int, input().split())
	mod = 10 ** 9 + 7
	com = 1
	for i in range(1, K - 1):
		com *= N * M - 1 - i
		com *= pow(i, mod - 2, mod)
		com %= mod
	ans = com * M * N * (M * N - 1) * (M + N) // 6
	return ans % mod
print(solve())


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

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


import math

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


(N, m, k) = (int(i) for i in input().split())
(n, mod, x, y) = (N * m, 10 ** 9 + 7, 0, 0)
(fn, fk) = ([1] * n, [1] * n)
for i in range(n - 1):
	fn[i + 1] = fn[i] * (i + 2) % mod

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

def comb(n, k):
	if n == 0 or k == 0 or n == k:
		return 1
	else:
		return fn[n - 1] * fk[n - k - 1] % mod * fk[k - 1] % mod
fk[-1] = power(fn[-1], mod - 2)
for i in range(2, n + 1):
	fk[-i] = fk[-i + 1] * (n + 2 - i) % mod
for i in range(N):
	x += i * (N - i) * m ** 2
for i in range(m):
	y += i * (m - i) * N ** 2
print(comb(n - 2, k - 2) * (x + y) % mod)


(n, m, k) = map(int, input().split())
U = 2 * 10 ** 5
MOD = 10 ** 9 + 7
fact = [1] * (U + 1)
fact_inv = [1] * (U + 1)
for i in range(1, U + 1):
	fact[i] = fact[i - 1] * i % MOD
fact_inv[U] = pow(fact[U], MOD - 2, MOD)
for i in range(U, 0, -1):
	fact_inv[i - 1] = fact_inv[i] * i % MOD

def comb(n, k):
	if k < 0 or k > n:
		return 0
	x = fact[n]
	x *= fact_inv[k]
	x %= MOD
	x *= fact_inv[n - k]
	x %= MOD
	return x
ans = 0
for i in range(n):
	for j in range(m):
		if i != 0 and j != 0:
			ans += 2 * (n - i) * (m - j) * (i + j)
		else:
			ans += (n - i) * (m - j) * (i + j)
		ans %= MOD
print(ans * comb(n * m - 2, k - 2) % MOD)


(n, m, k) = map(int, input().split())
mod = 10 ** 9 + 7
A = n * m - 2
from math import factorial

def Comb(N, K):
	a = factorial(N) % mod
	b = factorial(K) % mod
	c = factorial(N - K) % mod
	return a * pow(b, mod - 2, mod) * pow(c, mod - 2, mod) % mod
v = Comb(A, k - 2)
t = 0
y = 0
for i in range(1, n):
	t += i * (n - i) * m * m % mod
for i in range(1, m):
	y += i * (m - i) * n * n % mod
ans = v * (t + y)
print(ans % mod)


mod = 10 ** 9 + 7
(n, m, k) = map(int, input().split())
nm = n * m
f = [1] * (nm + 1)
for i in range(1, nm + 1):
	f[i] = f[i - 1] * i % mod

def C(n, k):
	x = f[n]
	x *= pow(f[n - k], mod - 2, mod)
	x %= mod
	x *= pow(f[k], mod - 2, mod)
	x %= mod
	return x

def F(n):
	x = n * n * (n + 1) // 2
	y = n * (n + 1) * (2 * n + 1) // 6
	return (x - y) % mod
res = (m * m % mod * F(n) % mod + n * n % mod * F(m) % mod) % mod
res *= C(n * m - 2, k - 2)
print(res % mod)


(N, M, K) = map(int, input().split())
P = 1000000007
import math
ans = 0
for A in [N, M]:
	for i in range(1, A):
		ans = (ans + (N ^ M ^ A) ** 2 * i * (i + 1) // 2) % P
print(ans * math.factorial(N * M - 2) * pow(math.factorial(K - 2), P - 2, P) * pow(math.factorial(N * M - K), P - 2, P) % P)


MOD = 10 ** 9 + 7
list_size = 200001
f_list = [1] * list_size
f_r_list = [1] * list_size
for i in range(list_size - 1):
	f_list[i + 1] = int(f_list[i] * (i + 2) % MOD)
f_r_list[-1] = pow(f_list[-1], MOD - 2, MOD)
for i in range(2, list_size + 1):
	f_r_list[-i] = int(f_r_list[-i + 1] * (list_size + 2 - i) % MOD)

def comb(n, r):
	if n < r:
		return 0
	elif n == 0 or r == 0 or n == r:
		return 1
	else:
		return f_list[n - 1] * f_r_list[n - r - 1] % MOD * f_r_list[r - 1] % MOD
(n, m, k) = map(int, input().split())
print((n + m) * n * m * (m * n - 1) // 6 * comb(n * m - 2, k - 2) % MOD)


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)


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)


mdl = 10 ** 9 + 7
(N, M, K) = map(int, input().split())
nm = N * M
mm = M * M % mdl
nn = N * N % mdl
U = 2 * 10 ** 5
fact = [1] * (U + 1)
for i in range(1, U + 1):
	fact[i] = fact[i - 1] * i % mdl

def comb(n, k):
	x = fact[n]
	x *= pow(fact[n - k], mdl - 2, mdl)
	x %= mdl
	x *= pow(fact[k], mdl - 2, mdl)
	x %= mdl
	return x

def F(N):
	x = N * (N + 1) // 2 * N
	y = N * (N + 1) * (2 * N + 1) // 6
	return (x - y) % mdl
ans = mm * F(N) % mdl + nn * F(M) % mdl
ans %= mdl
ans *= comb(nm - 2, K - 2)
ans %= mdl
print(ans)


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)


from functools import reduce

def combination2(n, r, MOD=10 ** 9 + 7):
	if not 0 <= r <= n:
		return 0
	r = min(r, n - r)
	numerator = reduce(lambda x, y: x * y % MOD, range(n, n - r, -1), 1)
	denominator = reduce(lambda x, y: x * y % MOD, range(1, r + 1), 1)
	return numerator * pow(denominator, MOD - 2, MOD) % MOD

def q(x, y):
	return cn * x ** 2 * (y ** 3 - y) // 6
(N, M, K) = map(int, input().split())
MOD = 10 ** 9 + 7
cn = combination2(N * M - 2, K - 2, MOD)
print((q(N, M) + q(M, N)) % MOD)


from collections import defaultdict
mod = 10 ** 9 + 7
memo = defaultdict(int)

def cmb(n, r):
	if (n, r) in memo:
		return memo[n, r]
	if n < r or r < 0:
		return 0
	if r > n - r:
		return cmb(n, n - r)
	res = f[n] * pow(f[n - r], mod - 2, mod) * pow(f[r], mod - 2, mod) % mod
	memo[n, r] = res
	return res
(n, m, k) = map(int, input().split())
f = [1]
for i in range(1, 2 * 10 ** 5 + 1):
	f.append(f[-1] * i % mod)
ans = 0
for d in range(1, m):
	ans += d * (m - d) * n ** 2
	ans %= mod
for d in range(1, n):
	ans += d * (n - d) * m ** 2
	ans %= mod
ans *= cmb(n * m - 2, k - 2)
ans %= mod
print(ans)


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 ** 0.5 + 1)
	arr = np.resize(arr, Lsq ** 2).reshape(Lsq, Lsq)
	for n in range(1, Lsq):
		arr[:, n] *= arr[:, n - 1]
		arr[:, n] %= MOD
	for n in range(1, Lsq):
		arr[n] *= arr[n - 1, -1]
		arr[n] %= MOD
	return arr.ravel()[:L]

def make_fact(U, MOD):
	x = np.arange(U, dtype=np.int64)
	x[0] = 1
	fact = cumprod(x, MOD)
	x = np.arange(U, 0, -1, dtype=np.int64)
	x[0] = pow(int(fact[-1]), MOD - 2, MOD)
	fact_inv = cumprod(x, MOD)[::-1]
	return (fact, fact_inv)
U = N * M + 10
(fact, fact_inv) = make_fact(U, MOD)
(fact, fact_inv) = (fact.tolist(), fact_inv.tolist())

def mod_comb_k(n, k, mod):
	return fact[n] * fact_inv[k] % mod * fact_inv[n - k] % mod
for dx in range(1, N):
	pat = N - dx
	ans += dx * pat * M * M
	ans %= MOD
for dy in range(1, M):
	pat = M - dy
	ans += dy * pat * N * N
	ans %= MOD
ans *= mod_comb_k(N * M - 2, K - 2, MOD)
ans %= MOD
print(ans)


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

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)


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


def com(n, r):
	re = fct[n] * pow(fct[r], md - 2, md) * pow(fct[n - r], md - 2, md)
	re %= md
	return re
(n, m, k) = map(int, input().split())
md = 10 ** 9 + 7
ans = 0
fct = [1] * (n * m - 1)
f = 1
for i in range(1, n * m - 1):
	f = f * i % md
	fct[i] = f
for d in range(1, m):
	ans = (ans + d * (m - d) * n * n) % md
for d in range(1, n):
	ans = (ans + d * (n - d) * m * m) % md
ans = ans * com(n * m - 2, k - 2) % md
print(ans)


import math
import sys
sys.setrecursionlimit(10000)
INF = float('inf')

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

def comb(n, r, mod):
	return math.factorial(n) * inv(math.factorial(n - r)) * inv(math.factorial(r)) % mod
(N, M, K) = list(map(int, input().split()))

def pyramid(n):
	return n * (n + 1) * (n + 2) // 6
MOD = 10 ** 9 + 7
k2 = (pyramid(N - 1) * M ** 2 + pyramid(M - 1) * N ** 2) % MOD
c = comb(N * M - 2, K - 2, mod=MOD)
print(k2 * c % MOD)


import sys
input = sys.stdin.readline
import math
(N, M, K) = [int(item) for item in input().split()]
mod = 10 ** 9 + 7
score_c = 1
for i in range(2, K):
	score_c *= (N * M - i) % mod * pow(K - i, mod - 2, mod)
	score_c %= mod
total = 0
for d in range(1, N):
	score = (N - d) * M ** 2 * d * score_c
	score %= mod
	total += score
	total %= mod
for d in range(1, M):
	score = (M - d) * N ** 2 * d * score_c
	score %= mod
	total += score
	total %= mod
print(total)


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


MOD = 10 ** 9 + 7
(N, M, K) = map(int, input().split())
invs = [1] * (K - 1)
for x in range(2, K - 1):
	invs[x] = -(MOD // x) * invs[MOD % x] % MOD
C = 1
for x in range(N * M - K + 1, N * M - 1):
	C *= x
	C %= MOD
for x in range(1, K - 1):
	C *= invs[x]
	C %= MOD
ans1 = 0
for d in range(1, N):
	ans1 += (N - d) * d
	ans1 %= MOD
ans1 *= M ** 2
ans1 %= MOD
ans2 = 0
for d in range(1, M):
	ans2 += (M - d) * d
	ans2 %= MOD
ans2 *= N ** 2
ans2 %= MOD
print((ans1 + ans2) * C % MOD)


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

def Combi2(a, b):
	if b > a:
		return 0
	ANS = 1
	for i in range(min(b, a - b)):
		ANS = ANS * (a - i) * pow(min(b, a - b) - i, mod - 2, mod) % mod
	return ANS % mod
SUM_N = SUM_M = 0
for i in range(N):
	SUM_N += i * (N - i)
for i in range(M):
	SUM_M += i * (M - i)
print((SUM_N * M ** 2 + SUM_M * N ** 2) * pow(Combi2(N * M, 2), mod - 2, mod) * Combi2(N * M, K) * K * (K - 1) * pow(2, mod - 2, mod) % mod)


(N, M, K) = map(int, input().split())
mod = 10 ** 9 + 7
frac = [1] * 200100
for i in range(2, 200100):
	frac[i] = i * frac[i - 1] % mod
fraci = [None] * 200100
fraci[-1] = pow(frac[-1], mod - 2, mod)
for i in range(-2, -200100 - 1, -1):
	fraci[i] = fraci[i + 1] * (200100 + i + 1) % mod

def comb(a, b):
	if not a >= b >= 0:
		return 0
	return frac[a] * fraci[b] * fraci[a - b] % mod
mod = 10 ** 9 + 7
t = N ** 2 * (M + 1) * M * (M - 1) // 6 + M ** 2 * (N + 1) * N * (N - 1) // 6
print(t * comb(N * M - 2, K - 2) % mod)


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 comb(self, n, k):
		return self.f[n] * self.i[n - k] * self.i[k]
(n, m, k) = map(int, input().split())
M = 10 ** 9 + 7
print((sum((m * m * i * (n - i) for i in range(n))) + sum((n * n * i * (m - i) for i in range(m)))) * Factorial(n * m, M).comb(n * m - 2, k - 2) % M)


MOD = 10 ** 9 + 7

def com(n, k, mod):
	res = 1
	tmp = 1
	for i in range(1, k + 1):
		res = res * (n - i + 1) % mod
		tmp = tmp * i % mod
	a = pow(tmp, mod - 2, mod)
	return res * a % mod
(N, M, K) = map(int, input().split())
C = com(N * M - 2, K - 2, MOD)
(X, Y) = (0, 0)
for d in range(1, N):
	X = (X + (N - d) * M ** 2 * d) % MOD
for d in range(1, M):
	Y = (Y + (M - d) * N ** 2 * d) % MOD
print((X + Y) * C % MOD)


def LI():
	return list(map(int, input().split()))
(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()


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

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

def ncr(n, r):
	res = 1
	for i in range(1, r + 1):
		res = res * n * modinv(i, div) % div
		n = n - 1
	return res % div
result = 0
comb = ncr(N * M - 2, K - 2)
for d in range(1, N):
	result += d * comb * M * M * (N - d)
for d in range(1, M):
	result += d * comb * (M - d) * N * N
print(result % div)


import sys

def MI():
	return map(int, sys.stdin.readline().rstrip().split())
(N, M, K) = MI()
mod = 10 ** 9 + 7

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

def nCr(n, r):
	return kaijou(n) * pow(kaijou(r), mod - 2, mod) * pow(kaijou(n - r), mod - 2, mod) % mod
a = N * (N + 1) * (N - 1) // 6 * M ** 2 + M * (M + 1) * (M - 1) // 6 * N ** 2
a %= mod
print(nCr(N * M - 2, K - 2) * a % mod)


import math

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


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)


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

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


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


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)


class Calc:

	def __init__(self, max_value, mod):
		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()
__starting_point()


def cmb(n, r, mod):
	if r < 0 or r > n:
		return 0
	r = min(r, n - r)
	return g1[n] * g2[r] * g2[n - r] % mod
(n, m, k) = [int(_) for _ in input().split()]
mod = 10 ** 9 + 7
N = n * m
g1 = [1, 1]
g2 = [1, 1]
inverse = [0, 1]
for i in range(2, N + 1):
	g1.append(g1[-1] * i % mod)
	inverse.append(-inverse[mod % i] * (mod // i) % mod)
	g2.append(g2[-1] * inverse[-1] % mod)
ans = 0
C = cmb(n * m - 2, k - 2, mod)
for d in range(1, n):
	ans += (n - d) * m ** 2 * d
	ans %= mod
for d in range(1, m):
	ans += (m - d) * n ** 2 * d
	ans %= mod
ans = ans * C % mod
print(ans)


(n, m, k) = list(map(int, input().split()))
MAX = 2 * 10 ** 5 + 1
MOD = 10 ** 9 + 7
fact = [0 for i in range(MAX)]
inv = [0 for i in range(MAX)]
invfact = [0 for i in range(MAX)]

def comb_build(n):
	fact[0] = inv[0] = invfact[0] = 1
	fact[1] = inv[1] = invfact[1] = 1
	for i in range(2, n):
		fact[i] = fact[i - 1] * i % MOD
		inv[i] = inv[MOD % i] * (MOD - MOD // i) % MOD
		invfact[i] = invfact[i - 1] * inv[i] % MOD

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