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

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


def comb(n, r):
    if n - r < r: r = n - r
    if r == 0: return 1
    if r == 1: return n

    numerator = [n - r + k + 1 for k in range(r)]
    denominator = [k + 1 for k in range(r)]

    for p in range(2,r+1):
        pivot = denominator[p - 1]
        if pivot > 1:
            offset = (n - r) % p
            for k in range(p-1,r,p):
                numerator[k - offset] /= pivot
                denominator[k] /= pivot

    result = 1
    for k in range(r):
        if numerator[k] > 1:
            result *= int(numerator[k])
            result%=MOD

    return result



ans_x=0
ans_y=0

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

Solution (Python):
MOD = 10 ** 9 + 7
MAX = 200000 + 5
fact = [1 for _ in range(MAX)]
finv = [1 for _ in range(MAX)]
for i in range(2, MAX):
    fact[i] = fact[i - 1] * i % MOD
    finv[i] = pow(fact[i], MOD-2, MOD)

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

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

ans = 0
for i in range(n):
    for j in range(m):
        a = j * (j+1) // 2
        b = (m-1-j) * (m-1-j+1) // 2
        ans += (a + b) * i + m * (i*(i+1)//2) + a
        ans %= MOD

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

print(ans)

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

import math

from operator import mul
from functools import reduce
def combinations_count(n, r):
    if n - r < r: r = n - r
    if r == 0: return 1
    if r == 1: return n

    numerator = [n - r + k + 1 for k in range(r)]
    denominator = [k + 1 for k in range(r)]

    for p in range(2,r+1):
        pivot = denominator[p - 1]
        if pivot > 1:
            offset = (n - r) % p
            for k in range(p-1,r,p):
                numerator[k - offset] /= pivot
                denominator[k] /= pivot

    result = 1
    for k in range(r):
        if numerator[k] > 1:
            result *= int(numerator[k])

    return result

ans_x = 0
c = combinations_count(N*M-2, K-2)


    
    
ans_x = (N*N*(N+1) + N*(N+1)*(2*N+1)//6 - (2*N+1)*N*(N+1)//2)//2
ans_x = M*M*ans_x*c

ans_y = 0
c = combinations_count(N*M-2, K-2)


    
    
ans_y = (M*M*(M+1) + M*(M+1)*(2*M+1)//6 - (2*M+1)*M*(M+1)//2)//2
ans_y = N*N*ans_y*c
ans = (ans_x+ans_y)%(10**9+7)
print(ans)

Solution (Python):
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)

Solution (Python):
import sys
sys.setrecursionlimit(100000000)
def input():
    return sys.stdin.readline()[:-1]
from bisect import *
from collections import *
from heapq import *
from math import *
N, M, K = map(int, input().split())
L = N*M
mod = 10**9+7
def factrial(N):
    fc = [1, 1] + [None] * N
    inv = [0, 1] + [None] * N
    fcInv = [1, 1] + [None] * N
    for i in range(2, N+1):
        fc[i] = i*fc[i-1]%mod
        inv[i] = mod-(inv[mod%i]*(mod//i))%mod
        fcInv[i] = fcInv[i-1]*inv[i]%mod
    return fc, fcInv
fc, fcInv = factrial(L)
def modCom(n, k):
    return 0 if n < k else fc[n] * fcInv[k] * fcInv[n-k] % mod
ans = 0
for d in range(1, N):
    ans += d*(N-d)*M*M*modCom(L-2, K-2)
    ans %= mod
for d in range(1, M):
    ans += d*(M-d)*N*N*modCom(L-2, K-2)
    ans %= mod
print(ans)

Solution (Python):
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)

Solution (Python):
import math, string, itertools, fractions, heapq, collections, re, array, bisect, copy, functools, random
import sys
from collections import deque, defaultdict, Counter; from heapq import heappush, heappop
from itertools import permutations, combinations, product, accumulate, groupby
from bisect import bisect_left, bisect_right, insort_left, insort_right
from operator import itemgetter as ig
sys.setrecursionlimit(10 ** 7)
inf = 10 ** 20; INF = float("INF"); ans = 0; tmp = 0; ansli = []; tmpli = []; candili = []
eps = 1.0 / 10 ** 10; mod = 10 ** 9 + 7
dd = [(-1, 0), (0, 1), (1, 0), (0, -1)]; ddn = dd + [(-1, 1), (1, 1), (1, -1), (-1, -1)]; ddn9 = ddn + [(0, 0)]
"""for dx, dy in dd:
        nx = j + dx; ny = i + dy
            if 0 <= nx < w and 0 <= ny < h:"""
def wi(): return list(map(int, sys.stdin.readline().split()))
def wip(): return [int(x) - 1 for x in sys.stdin.readline().split()]
def ws(): return sys.stdin.readline().split()
def i(): return int(sys.stdin.readline())
def s(): return input()
def hi(n): return [i() for _ in range(n)]
def hs(n): return [s() for _ in range(n)]
def mi(n): return [wi() for _ in range(n)]
def mip(n): return [wip() for _ in range(n)]
def ms(n): return [ws() for _ in range(n)]

n, m, k = wi()

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

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

MAX_N = 3 * 10 ** 5 + 100
factorial = [1] * MAX_N

def calc_factorial():
    for i in range(1, MAX_N):
        factorial[i] = i * factorial[i - 1] % mod

def comb(n, k):
    a = factorial[n] % mod
    b = (factorial[k] * factorial[n - k]) % mod
    b_ = pow(b, mod - 2, mod)
    return (a * b_) % mod



calc_factorial()

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

Solution (Python):
n,m,k=map(int,input().split())
mod=10**9+7
import math
def f(x):
  return math.factorial(x)
def c(x,y):
  return f(x)//(f(x-y)*f(y))

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

Solution (Python):
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))

Solution (Python):
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)

Solution (Python):
import sys
s2nn = lambda s: [int(c) for c in s.split(' ')]
ss2nn = lambda ss: [int(s) for s in ss]
ss2nnn = lambda ss: [s2nn(s) for s in ss]
i2s = lambda: sys.stdin.readline().rstrip()
i2n = lambda: int(i2s())
i2nn = lambda: s2nn(i2s())
ii2ss = lambda n: [sys.stdin.readline().rstrip() for _ in range(n)]
ii2sss = lambda n: [list(sys.stdin.readline().rstrip()) for _ in range(n)]
ii2nn = lambda n: ss2nn(ii2ss(n))
ii2nnn = lambda n: ss2nnn(ii2ss(n))

MOD = int(1e+9) + 7

def comb(n, k):
    r = 1
    for i in range(n, max(n-k, k), -1):
        r *= i
    for i in range(1, min(n-k, k) + 1, 1):
        r //= i
    return r

class Comb:
    def __init__(self, MAX, MOD):
        MAX += 1
        fac = [0] * MAX
        finv = [0] * MAX
        inv = [0] * 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
        self.fac = fac
        self.finv = finv
        self.MOD = MOD
    
    def comb(self, n, k):
        fac = self.fac
        finv = self.finv
        MOD = self.MOD
        if (n < k):
            return 0
        if (n < 0 or k < 0):
            return 0
        return fac[n] * (finv[k] * finv[n - k] % MOD) % MOD

def main():
    
    N, M, K = i2nn()
    

    
    
    
    
    xcost = 0
    for i in range(1, M):
        
        
        
        
        
        
        
        xcost = (xcost + i * (M-i) * (N*N)) % MOD
    ycost = 0
    for i in range(1, N):
        ycost = (ycost + i * (N-i) * (M*M)) % MOD
    cost = (xcost + ycost) % MOD
    
    
    
    cost = (cost * Comb(N*M-2, MOD).comb(N*M-2, K-2)) % MOD
    print(cost)

main()

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






mod = 10**9 + 7


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



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





max_len = n * m 

factori_table = [1] * (max_len + 1)
factori_inv_table = [1] * (max_len + 1)
for i in range(1, max_len + 1):
	factori_table[i] = factori_table[i-1] * (i) % mod



for i in [k-2, n*m-k]:
	factori_inv_table[i] = modinv(factori_table[i])

def binomial_coefficients(n, k):
	
	return factori_table[n] * factori_inv_table[k] * factori_inv_table[n-k]

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

print(binomial_coefficients(n*m-2, k-2) * ans_1 % mod)

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


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

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

mod = 10**9 + 7

base_n = 0
base_m = 0

for i in range(M):
    if i == 0:
        s = (M * (M - 1) // 2) % mod
    else:
        s += (i - (M - i))
    base_m = (base_m + s + mod) % mod
for i in range(N):
    if i == 0:
        s = (N * (N - 1) // 2) % mod
    else:
        s += (i - (N - i))
    base_n = (base_n + s + mod) % mod
base = (base_m * N * N + base_n * M * M) * g2[2] % mod

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

print(ans)

Solution (Python):



""" convenient functions




"""

import sys
sys.setrecursionlimit(10**7)
from itertools import accumulate, combinations, permutations, product 

from math import factorial, ceil, floor
def factorize(n):
    """return the factors of the Arg and count of each factor
    
    Args:
        n (long): number to be resolved into factors
    
    Returns:
        list of tuples: factorize(220) returns [(2, 2), (5, 1), (11, 1)]
    """
    fct = []  
    b, e = 2, 0  
    while b * b <= n:
        while n % b == 0:
            n = n // b
            e = e + 1
        if e > 0:
            fct.append((b, e))
        b, e = b + 1, 0
    if n > 1:
        fct.append((n, 1))
    return fct
def combinations_count(n, r):
    """Return the number of selecting r pieces of items from n kinds of items.
    
    Args:
        n (long): number
        r (long): number
    
    Raises:
        Exception: not defined when n or r is negative
    
    Returns:
        long: number
    """
    
    if n < 0 or r < 0:
        raise Exception('combinations_count(n, r) not defined when n or r is negative')
    if n - r < r: r = n - r
    if r < 0: return 0
    if r == 0: return 1
    if r == 1: return n
    numerator = [n - r + k + 1 for k in range(r)]
    denominator = [k + 1 for k in range(r)]
    for p in range(2,r+1):
        pivot = denominator[p - 1]
        if pivot > 1:
            offset = (n - r) % p
            for k in range(p-1,r,p):
                numerator[k - offset] /= pivot
                denominator[k] /= pivot
    result = 1
    for k in range(r):
        if numerator[k] > 1:
            result *= int(numerator[k])
    return result
def combinations_with_replacement_count(n, r):
    """Return the number of selecting r pieces of items from n kinds of items allowing individual elements to be repeated more than once.
    
    Args:
        n (long): number
        r (long): number
    
    Raises:
        Exception: not defined when n or r is negative
    
    Returns:
        long: number
    """
    if n < 0 or r < 0:
        raise Exception('combinations_with_replacement_count(n, r) not defined when n or r is negative')
    elif n == 0:
        return 1
    else:
        return combinations_count(n + r - 1, r)
from bisect import bisect_left, bisect_right
from collections import deque, Counter, defaultdict 
from heapq import heapify, heappop, heappush, heappushpop, heapreplace,nlargest,nsmallest 
from copy import deepcopy, copy 
from operator import itemgetter


from functools import reduce
def chmin(x, y):
    """change minimum
    if x > y, x = y and return (x, True).
    convenient when solving problems of dp[i]
    
    Args:
        x (long): current minimum value
        y (long): potential minimum value
    
    Returns:
        (x, bool): (x, True) when updated, else (x, False)
    """
    if x > y:
        x = y
        return (x, True)
    else:
        return (x, False)
def chmax(x, y):
    """change maximum
    if x < y, x = y and return (x, True).
    convenient when solving problems of dp[i]
    
    Args:
        x (long): current maximum value
        y (long): potential maximum value
    
    Returns:
        (x, bool): (x, True) when updated, else (x, False)
    """
    if x < y:
        x = y
        return (x, True)
    else:
        return (x, False)

from fractions import gcd 
def gcds(numbers):
    return reduce(gcd, numbers)
def lcm(x, y):
    return (x * y) // gcd(x, y)
def lcms(numbers):
    return reduce(lcm, numbers, 1)



INF = 10 ** 18
MOD = 10 ** 9 + 7
modpow = lambda a, n, p = MOD: pow(a, n, p) 
def modinv(a, p = MOD):
    
    
    return modpow(a, p-2, p)
def modinv_list(n, p = MOD):
    if n <= 1:
        return [0,1][:n+1]
    else:
        inv_t = [0,1]
        for i in range(2, n+1):
            inv_t += [inv_t[p % i] * (p - int(p / i)) % p]
        return inv_t
def modfactorial_list(n, p = MOD):
    if n == 0:
        return [1]
    else:
        l = [0] * (n+1)
        tmp = 1
        for i in range(1, n+1):
            tmp = tmp * i % p
            l[i] = tmp
        return l
def modcomb(n, k, fac_list = [], p = MOD):
    
    
    from math import factorial
    if n < 0 or k < 0 or n < k: return 0
    if n == 0 or k == 0: return 1
    if len(fac_list) <= n:
        a = factorial(n) % p
        b = factorial(k) % p
        c = factorial(n-k) % p
    else:
        a = fac_list[n]
        b = fac_list[k]
        c = fac_list[n-k]
    return (a * modpow(b, p-2, p) * modpow(c, p-2, p)) % p
def modadd(a, b, p = MOD):
    return (a + b) % MOD
def modsub(a, b, p = MOD):
    return (a - b) % p
def modmul(a, b, p = MOD):
    return ((a % p) * (b % p)) % p
def moddiv(a, b, p = MOD):
    return modmul(a, modpow(b, p-2, p))

""" initialize variables and set inputs

    
    

    
    
    
    
    

    
    
    
    

    
    
"""


r = lambda: sys.stdin.readline().strip()
r_int = lambda: int(r())
R = lambda: list(map(int, r().split()))
Rfloat = lambda: list(map(float, r().split()))
Rtuple = lambda: tuple(map(int, r().split()))
Rmap = lambda: map(int, r().split())

""" how to treat input







"""






N, M, K = R()

s = modcomb(N*M-2,K-2)
ans = 0
for i in range(1,M):
    ans += modmul(modmul(modmul(i,M-i),modmul(N,N)),s)
N, M = M, N
for i in range(1,M):
    ans += modmul(modmul(modmul(i,M-i),modmul(N,N)),s)

print(ans%MOD)



"""memo: how to use defaultdict of list

Dic = defaultdict(list)

Dic[x].append(y)

for k, v in Dic.items():
"""

Solution (Python):

from collections import defaultdict
from collections import deque
from heapq import heappush, heappop
import sys
import math
import bisect
import random
def LI(): return list(map(int, sys.stdin.readline().split()))
def I(): return int(sys.stdin.readline())
def LS():return list(map(list, sys.stdin.readline().split()))
def S(): return list(sys.stdin.readline())[:-1]
def IR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = I()
    return l
def LIR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = LI()
    return l
def SR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = S()
    return l
def LSR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = LS()
    return l
sys.setrecursionlimit(1000000)
mod = 1000000007


def A():
    return


def B():
    return


def C():
    n,m = LI()
    t = [0 for i in range(n+1)]
    for i in range(m):
        l,r = LI()
        t[l-1] += 1
        t[r] -= 1
    for i in range(n):
        t[i+1] += t[i]
        if t[i] == m:
            t[i] = 1
        else:
            t[i] = 0
    print(sum(t[:n]))
    return


def D():
    n,m = LI()
    a = LI()
    q = []
    for i in range(n):
        heappush(q,a[i])
    l = LIR(m)
    l.sort(key = lambda x:-x[1])
    for b,c in l:
        for i in range(b):
            x = heappop(q)
            if x < c:
                x = c
                heappush(q,x)
            else:
                heappush(q,x)
                break
    print(sum(q))
    return


def E():
    n,m,k = LI()
    M = n*m
    fact = [1]*(M+1)
    for i in range(M):
        fact[i+1] = fact[i]*(i+1)%mod
    inv = [1]*(M+1)
    inv[M] = pow(fact[M],mod-2,mod)
    for i in range(M)[::-1]:
        inv[i] = inv[i+1]*(i+1)%mod
    ans = fact[M-2]*inv[k-2]*inv[M-k]%mod
    s = 0
    for i in range(1,M):
        s += i*(max(0,n-i)*m**2+max(0,m-i)*n**2)%mod
        s %= mod
    ans = ans*s%mod
    print(ans)
    return


def F():
    q = I()
    m = 0
    l = [-float("inf"),float("inf")]
    x = 0
    num = 0
    for _ in range(q):
        Q = LI()
        if len(Q) == 1:
            print(x,m)
        else:
            z,a,b = Q
            m += b
            if num == 0:
                num += 1
                l = [a,float("inf")]
                x = a
            elif num == 1:
                if a != l[0]:
                    num += 1
                    l = [min(a,l[0]),max(a,l[0])]
                    x = l[0]
                    m += l[1]-l[0]
            else:
                num += 1
                if num % 2:
                    if l[0] <= a < l[1]:
                        m += min(a-l[0],l[1]-a)
                        x = a
                    elif a < l[0]:
                        m += l[0]-a
                        x = l[0]
                        l = [a,l[1]]
                    else:
                        m += a-l[1]
                        x = l[1]
                        l = [l[0],a]
                else:
                    if l[0] <= a < l[1]:
                        if a < x:
                            m += x-a
                            l = [a,x]
                            x = a
                        else:
                            m += a-x
                            l = [x,a]
                    elif a < l[0]:
                        m += x-a
                        x = l[0]
                        l = [l[0],x]
                    else:
                        m += a-x
                        l = [x,l[1]]
    return


def G():
    return


def H():
    return


def I_():
    return


def J():
    return


if __name__ == "__main__":
    E()

Solution (Python):
MOD = 10**9 + 7
MAX = int(2e5)


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)


if __name__ == "__main__":
    main()

Solution (Python):
facs = []
M = 1000000007

def modulo(n):
    global facs
    facs = [1]*(n+1)
    last = 1
    for i in range(1,n+1):
        last = facs[i] = mulmod(last,i,M)
        
def mulmod(x,y,p):
    return x*y % p
def divmod(x,y,p):
    return mulmod(x,powmod(y,p-2,p),p)

def ncr(n,r):
    if n<r:
        return 0
    if n==r:
        return 1
    res = facs[n]
    res = divmod(res,facs[r],M)
    res = divmod(res,facs[n-r],M)
    return res

def powmod(x,y,p):
    if y==0:
        return 1
    elif y==1:
        return x % p
    elif (y%2)==0:
        return powmod(x,y//2,p)**2 % p
    else:
        return powmod(x,y//2,p)**2 * x % p
def resolve():
    n,m,k = map(int,input().split())
    modulo(2*10**5+5)
    kumi=ncr(n*m-2,k-2)
    ans = 0
    for i in range(1,n):
        ans = ans + ((i*(n-i)*m*m)%M)
    for i in range(1,m):
        ans = ans + ((i*(m-i)*n*n)%M)
    print(ans*kumi%M)

if __name__ == "__main__":
    resolve()

Solution (Python):
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)

Solution (Python):
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)

Solution (Python):
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)

Solution (Python):
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)

Solution (Python):

import math

def cmb(n, r):
    if n - r < r: r = n - r
    if r == 0: return 1
    if r == 1: return n

    numerator = [n - r + k + 1 for k in range(r)]
    denominator = [k + 1 for k in range(r)]

    for p in range(2,r+1):
        pivot = denominator[p - 1]
        if pivot > 1:
            offset = (n - r) % p
            for k in range(p-1,r,p):
                numerator[k - offset] /= pivot
                denominator[k] /= pivot

    result = 1
    for k in range(r):
        if numerator[k] > 1:
            result *= int(numerator[k])

    return result

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

comb = cmb(N*M-2, K-2) % 1000000007

cost = 0
for d in range(1, N):
    cost += d * M**2 * (N - d) * comb
    cost %= 1000000007
for d in range(1, M):
    cost += d * N**2 * (M - d) * comb
    cost %= 1000000007
    
print(cost)

Solution (Python):
def factorials(M):
    INF=int(1e9+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(1e9+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)

Solution (Python):
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)

Solution (Python):
def mod_combination(mod: int, n: int, r: int) -> int:
  a, b, c = mod_factorials(mod, [n, r, n - r])

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

def mod_factorials(mod: int, N: [int]) -> [int]:
  result = [1] * len(N)
  tuples = sorted([(i, n) for i, n in enumerate(N)], key=lambda t: t[1])

  fact = 1
  i = 1
  for t in tuples:
    for j in range(i, t[1] + 1):
      fact = fact * j % mod
      i = i + 1

    result[t[0]] = fact

  return result


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

  N, M, K = [int(i) for i in input().strip().split(' ')]
  c = mod_combination(MOD, N * M - 2, K - 2)

  X = 0
  Y = 0
  for i in range(M):
    X += (M - i) * i
  for i in range(N):
    Y += (N - i) * i

  total = (X * (N ** 2) + Y * (M ** 2)) % MOD * c % MOD
  print(total)


main()

Solution (Python):
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)

Solution (Python):
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*power_func(b, mod-2, mod)*power_func(c, mod-2, mod)) % mod


def power_func(a,b,mod):
  if b == 0:
      return 1
  if b % 2 == 0:
    d = power_func(a, b//2, mod)
    return d*d % mod
  if b % 2 == 1:
    return (a*power_func(a, b-1, 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)

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

def mod_inv(i, mod):
    return pow(i, mod-2, mod)

def combination(_N, _K):
    if _K > (_N/2):
        _K = _N - _K
    up = 1
    low = 1
    for i in range (_K):
        up = up * (_N - i) % MAX
        low = low * (i + 1) % MAX
    ret = up * mod_inv(low, MAX) % MAX
    return ret

def cum_sum(N):
    return N * (N+1) // 2

def solve2():
    N, M, K = [int(_) for _ in raw_input().split()]
    sumx = 0
    for xi in range(1, N+1):
        sumx = (sumx + cum_sum(N) - (N + 1 - xi) * xi)
    sumx = (sumx * M * M) // 2 %MAX
    sumy = 0
    for yi in range(1, M+1):
        sumy = (sumy + cum_sum(M) - (M + 1 - yi) * yi)
    sumy = (sumy * N * N) // 2 %MAX
    ret = (sumx + sumy) * combination(M*N-2, K-2) % MAX
    print ret

if __name__ == '__main__':
    solve2()

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

def mod_inv(i, mod):
    return pow(i, mod-2, mod)

def combination(_N, _K):
    if _K > (_N/2):
        _K = _N - _K
    up = 1
    low = 1
    for i in range (_K):
        up = up * (_N - i) % MAX
        low = low * (i + 1) % MAX
    ret = up * mod_inv(low, MAX) % MAX
    return ret

def cum_sum(N):
    return N * (N+1) // 2

def solve2():
    N, M, K = [int(_) for _ in raw_input().split()]
    sumx = 0
    for xi in range(1, N+1):
        sumx = (sumx + cum_sum(N) - (N + 1 - xi) * xi) % MAX
    sumx = (sumx * M * M) // 2 %MAX
    sumy = 0
    for yi in range(1, M+1):
        sumy = (sumy + cum_sum(M) - (M + 1 - yi) * yi) % MAX
    sumy = (sumy * N * N) // 2 %MAX
    ret = (sumx + sumy) * combination(M*N-2, K-2) % MAX
    print ret

if __name__ == '__main__':
    solve2()

Solution (Python):

MAX = 10**9 + 7

def mod_inv(i, mod):
    return pow(i, mod-2, mod)

def combination(_N, _K):
    if _K > (_N/2):
        _K = _N - _K
    up = 1
    low = 1
    for i in range (_K):
        up = up * (_N - i) % MAX
        low = low * (i + 1) % MAX
    ret = up * mod_inv(low, MAX) % MAX
    return ret

def solve():
    N, M, K = [int(_) for _ in raw_input().split()]
    disSum = M*N* (M*(N+1)*(N-1) + N*(M+1)*(M-1)) / 6
    ret = disSum * combination(M*N-2, K-2) % MAX
    print ret

if __name__ == '__main__':
    solve()