CodeEff/ECCO huggingface dataset:
# Original solution
  
    
N,M,K = list(map(int,input().split()))


mod = 10 ** 9 + 7
MAX = N*M
fac = [0] * MAX
finv = [0] * MAX
inv = [0] * MAX

def cominit():
    fac[0] = 1
    fac[1] = 1
    finv[0] = 1
    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



def cmb(n,r):
    if n < 0 or r < 0 or r > n:return 0
    if r > n/2: r = n-r        
    return fac[n] * (finv[r] * finv[n-r] % mod) % mod
    
cominit()
ans= 0
NCK = cmb(N*M-2,K-2) 
M2 = M**2 * NCK
N2 = N**2 * NCK
for i in range(N-1):
    for j in range(i,N):
        ans += M2 * (j-i) 
        ans %= mod
        
for i in range(M-1):
    for j in range(i,M):
        ans += N2 * (j-i) 
        ans %= mod
        
print(ans)

# Optimized solution
  
    
N,M,K = list(map(int,input().split()))


mod = 10 ** 9 + 7
MAX = N*M
fac = [0] * MAX
finv = [0] * MAX
inv = [0] * MAX

def cominit():
    fac[0] = 1
    fac[1] = 1
    finv[0] = 1
    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



def cmb(n,r):
    if n < 0 or r < 0 or r > n:return 0
    if r > n/2: r = n-r        
    return fac[n] * (finv[r] * finv[n-r] % mod) % mod
    
cominit()
ans= 0
NCK = cmb(N*M-2,K-2) 
M2 = M**2 * NCK
N2 = N**2 * NCK
for i in range(1,N):
    ans += M2 * i * (N-i)
    ans %= mod
    
for i in range(1,M):
    ans += N2 * i * (M-i)
    ans %= mod
        
print(ans)
# Original solution
#!/usr/bin/python3
# -*- coding: utf-8 -*-


MOD = 10**9+7
Facts = [1]*(2*10**5+1)
for i in range(2,2*10**5+1):
  Facts[i] = (Facts[i-1]*i) % MOD

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

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

def div(a,b):
  return mul(a,power(b,MOD-2))


def combi(n,m,k):
  if k <= 2 or k >= n*m:
    return 1
  else:
    return div(Facts[n*m-2],(Facts[k-2]*Facts[n*m-k]))


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

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

# Optimized solution
#!/usr/bin/python3
# -*- coding: utf-8 -*-


MOD = 10**9+7
Facts = [1]*(2*10**5+1)
for i in range(2,2*10**5+1):
  Facts[i] = (Facts[i-1]*i) % MOD

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

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

def div(a,b):
  return mul(a,power(b,MOD-2))


def combi(n,m,k):
  if k <= 2 or k >= n*m:
    return 1
  else:
    return div(Facts[n*m-2],(Facts[k-2]*Facts[n*m-k]))


def main():
  N,M,K = list(map(int, input().split()))
  com = combi(N,M,K)
  ans = 0
  for d in range(1,N):
    ans = (ans + d * (N-d)*M**2 * com) % MOD
  for d in range(1,M):
    ans = (ans + d * (M-d)*N**2 * com) % MOD
  return ans

if __name__ == "__main__":
  print((main()))
# Original solution
import sys
readline = sys.stdin.readline
sys.setrecursionlimit(10**8)
mod = 10**9+7
#mod = 998244353
INF = 10**18
eps = 10**-7

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

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

def f1(x,y,a,b):
    return (a+b-x-y)*(a-x+1)*(b-y+1)//2

def f2(x,a):
    return (a-x)*(a-x+1)//2

ans = 0

for x in range(1,n+1):
    for y in range(1,m+1):
        ans += f1(x,y,n,m)+f1(1,1,x,y)+f1(1,y,x,m)+f1(x,1,n,y)
        ans -= f2(x,n) + f2(1,x) + f2(1,y) + f2(y,m)
        ans %= mod

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

print(ans)


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

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

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

def f1(x,y,a,b):
    return (a+b-x-y)*(a-x+1)*(b-y+1)//2

def f2(x,a):
    return (a-x)*(a-x+1)//2

ans = 0

for x in range(1,n+1):
    for y in range(1,m+1):
        a = n-x+1
        b = m-y+1
        c = m*(a*(a-1)+x*(x-1))
        d = n*(b*(b-1)+y*(y-1))
        ans += (c+d)//2
        ans %= mod

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

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

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

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

def f1(x,y,a,b):
    return (a+b-x-y)*(a-x+1)*(b-y+1)//2

def f2(x,a):
    return (a-x)*(a-x+1)//2

ans = 0

for x in range(1,n+1):
    for y in range(1,m+1):
        a = n-x+1
        b = m-y+1
        c = m*(a*(a-1)+x*(x-1))
        d = n*(b*(b-1)+y*(y-1))
        ans += (c+d)//2
        ans %= mod

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

print(ans)



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

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

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

ans = m*n*(m+n)*(m*n-1)//3

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

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

mod = 10**9+7 #出力の制限
N = 10**6
g1 = [1, 1] # 元テーブル
g2 = [1, 1] #逆元テーブル
inverse = [0, 1] #逆元テーブル計算用テーブル

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


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

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

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



if __name__ == "__main__":
    main()

# Optimized solution
class Combination:
    """
    O(n)の前計算を1回行うことで，O(1)でnCr mod mを求められる
    n_max = 10**6のとき前処理は約950ms (PyPyなら約340ms, 10**7で約1800ms)
    使用例：
    comb = Combination(1000000)
    print(comb(5, 3))  # 10
    """
    def __init__(self, n_max, mod=10**9+7):
        self.mod = mod
        self.modinv = self.make_modinv_list(n_max)
        self.fac, self.facinv = self.make_factorial_list(n_max)

    def __call__(self, n, r):
        if n < r:
            return 0
        return self.fac[n] * self.facinv[r] % self.mod * self.facinv[n-r] % self.mod

    def make_factorial_list(self, n):
        # 階乗のリストと階乗のmod逆元のリストを返す O(n)
        # self.make_modinv_list()が先に実行されている必要がある
        fac = [1]
        facinv = [1]
        for i in range(1, n+1):
            fac.append(fac[i-1] * i % self.mod)
            facinv.append(facinv[i-1] * self.modinv[i] % self.mod)
        return fac, facinv

    def make_modinv_list(self, n):
        # 0からnまでのmod逆元のリストを返す O(n)
        modinv = [0] * (n+1)
        modinv[1] = 1
        for i in range(2, n+1):
            modinv[i] = self.mod - self.mod//i * modinv[self.mod%i] % self.mod
        return modinv

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

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



if __name__ == "__main__":
    main()
# Original solution
class Combination:

    def __init__(self, n, MOD):
        self.fact = [1]
        for i in range(1, n + 1):
            self.fact.append(self.fact[-1] * i % MOD)
        self.inv_fact = [pow(self.fact[i] ,MOD - 2 ,MOD) for i in range(n + 1)]
        self.MOD = MOD

    def factorial(self, k):
        '''k!を求める'''
        return self.fact[k]

    def inverse_factorial(self, k):
        '''k!の逆元を求める'''
        return self.inv_fact[k]

    def combination(self, k, r):
        '''kCrを求める'''
        return (self.fact[k] * self.inv_fact[k - r] * self.inv_fact[r]) % self.MOD
      
h, w, k = list(map(int, input().split()))
MOD = 10**9 + 7
comb = Combination(h*w, MOD)

kumi = comb.combination(h*w - 2, k - 2)
ans = 0

for i in range(h):
    for j in range(i, h):
        diff_h = abs(i - j)
        ans += diff_h * kumi * w * w
        ans %= MOD

for i in range(w):
    for j in range(i, w):
        diff_w = abs(i - j)
        ans += diff_w * kumi * h * h
        ans %= MOD

print(ans)

# Optimized solution
class Combination:

    def __init__(self, n, MOD):
        self.fact = [1]
        for i in range(1, n + 1):
            self.fact.append(self.fact[-1] * i % MOD)
        self.inv_fact = [pow(self.fact[i] ,MOD - 2 ,MOD) for i in range(n + 1)]
        self.MOD = MOD

    def factorial(self, k):
        '''k!を求める'''
        return self.fact[k]

    def inverse_factorial(self, k):
        '''k!の逆元を求める'''
        return self.inv_fact[k]

    def combination(self, k, r):
        '''kCrを求める'''
        return (self.fact[k] * self.inv_fact[k - r] * self.inv_fact[r]) % self.MOD
      
h, w, k = list(map(int, input().split()))
MOD = 10**9 + 7
comb = Combination(h*w, MOD)

kumi = comb.combination(h*w - 2, k - 2)
ans = 0

for d in range(0, h):
    ans += d * (h - d) * kumi * w * w
    ans %= MOD

for d in range(0, w):
    ans += d * (w - d) * kumi * h * h
    ans %= MOD

print(ans)
# Original solution
n, m, k = list(map(int, input().split()))

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

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

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

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

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

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

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

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

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

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

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

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

print(((x + y) % mod))
# Original solution
N,M,K=list(map(int,input().split()))
MOD=10**9+7

from fractions import Fraction

# 組合せの総数
def comb(n, k):
    
    if n < 0 or k < 0:
        pdt = 0
    
    else:
        pdt = 1
        for s in range(1, k + 1):
            pdt *= Fraction(n - s + 1, s)
    
    # 戻り値をint型にするために分子を取り出す
    return pdt.numerator


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

# Optimized solution
N,M,K=list(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)))
# Original solution
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 = list(map(int, input().split()))
s = 0
mod = 10 ** 9 + 7
for x in range(m):
    for i in range(x, m):
        s += abs(x-i) * n **2
for y in range(n):
    for j in range(y, n):
        s += abs(y-j) * m ** 2


print((s * comb_mod(n * m - 2, k - 2, mod) % mod))


# Optimized solution
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 = list(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))
# Original solution
N, M, K = list(map(int, input().split()))
P = 10**9+7
fact = [1]
for i in range(N*M):
    fact.append(fact[-1]*(i+1))
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))

# Optimized solution
N, M, K = list(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))
# Original solution
import math
n,m,k = list(map(int,input().split()))
mod = 10**9+7
memo = {}
def factorial(num):
    if num == 0:
        return 1
    if not num in memo:
        memo[num] = num*factorial(num-1)
    return memo[num]

temp = math.factorial(n*m-2)//(math.factorial(k-2)*math.factorial(n*m-k))
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))

# Optimized solution
n,m,k = list(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))
# Original solution
n,m,k=list(map(int,input().split()))
ans=0
INF=10**9+7
if n<=m:
    M=(m-1)*m*(2*m-1)//3
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
            else:
                d=j-i
                ans=ans+d*m+2*d*m*(m-1)+(m-d)*m*(m-1)-M
            ans=ans%INF
else:
    (n,m)=(m,n)
    M=(m-1)*m*(2*m-1)//3
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
            else:
                d=j-i
                ans=ans+d*m+2*d*m*(m-1)+(m-d)*m*(m-1)-M
            ans=ans%INF
import math
ans=(math.factorial(n*m-2)//(math.factorial(n*m-k)*math.factorial(k-2)))*ans
ans=ans%INF
print(ans)


# Optimized solution
n,m,k=list(map(int,input().split()))
ans=0
INF=10**9+7
if n<=m:
    M=(m-1)*m*(2*m-1)//3
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
            else:
                d=j-i
                ans=ans+d*m+(m+d)*m*(m-1)-M
            ans=ans%INF
else:
    (n,m)=(m,n)
    M=(m-1)*m*(2*m-1)//3
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
            else:
                d=j-i
                ans=ans+d*m+(m+d)*m*(m-1)-M
            ans=ans%INF
import math
ans=(math.factorial(n*m-2)//(math.factorial(n*m-k)*math.factorial(k-2)))*ans
ans=ans%INF
print(ans)
# Original solution
n,m,k=list(map(int,input().split()))
ans=0
INF=10**9+7
if n<=m:
    M=-(m-1)*m*(2*m-1)//3+m*m*(m-1)
    Ms=m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
    M=M%INF
    Ms=Ms%INF
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+Ms
            else:
                d=j-i
                ans=ans+d*m**2+M
            ans=ans%INF
else:
    (n,m)=(m,n)
    M=-(m-1)*m*(2*m-1)//3+m*m*(m-1)
    Ms=m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
    M=M%INF
    Ms=Ms%INF
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+Ms
            else:
                d=j-i
                ans=ans+d*m**2+M
            ans=ans%INF
a=1
b=1
x=min(n*m-2,k-2)
for i in range(x):
    a=a*(n*m-2-i)
    b=b*(i+1)
ans=(ans*a//b)%INF
print(ans)


# Optimized solution
n,m,k=list(map(int,input().split()))
ans=0
INF=10**9+7
if n<=m:
    M=-(m-1)*m*(2*m-1)//3+m*m*(m-1)
    Ms=m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
    M=M%INF
    Ms=Ms%INF
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+Ms
            else:
                d=j-i
                ans=ans+d*m**2+M
            ans=ans%INF
else:
    (n,m)=(m,n)
    M=-(m-1)*m*(2*m-1)//3+m*m*(m-1)
    Ms=m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
    M=M%INF
    Ms=Ms%INF
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+Ms
            else:
                d=j-i
                ans=ans+d*m**2+M
            ans=ans%INF
import math
ans=(math.factorial(n*m-2)//(math.factorial(n*m-k)*math.factorial(k-2)))*ans
ans=ans%INF
print(ans)
# Original solution
n,m,k=list(map(int,input().split()))
ans=0
INF=10**9+7
if n<=m:
    M=-(m-1)*m*(2*m-1)//3+m*m*(m-1)
    Ms=m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
    M=M%INF
    Ms=Ms%INF
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+Ms
            else:
                d=j-i
                ans=ans+d*m**2+M
            ans=ans%INF
else:
    (n,m)=(m,n)
    M=-(m-1)*m*(2*m-1)//3+m*m*(m-1)
    Ms=m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
    M=M%INF
    Ms=Ms%INF
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+Ms
            else:
                d=j-i
                ans=ans+d*m**2+M
            ans=ans%INF
a=1
b=1
c=n*m-2
x=min(c,k-2)
for i in range(x):
    a=a*(c-i)
    b=b*(i+1)
ans=(ans*a//b)%INF
print(ans)


# Optimized solution
n,m,k=list(map(int,input().split()))
ans=0
INF=10**9+7
if n<=m:
    M=-(m-1)*m*(2*m-1)//3+m*m*(m-1)
    Ms=m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
    M=M%INF
    Ms=Ms%INF
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+Ms
            else:
                d=j-i
                ans=ans+d*m**2+M
            ans=ans%INF
else:
    (n,m)=(m,n)
    M=-(m-1)*m*(2*m-1)//3+m*m*(m-1)
    Ms=m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
    M=M%INF
    Ms=Ms%INF
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+Ms
            else:
                d=j-i
                ans=ans+d*m**2+M
            ans=ans%INF
def modcomb(n,k,m):
    fac=[0]*(n+1)
    finv=[0]*(n+1)
    inv=[0]*(n+1)
    fac[0]=fac[1]=1
    finv[0]=finv[1]=1
    inv[1]=1
    for i in range(2,n+1):
        fac[i]=fac[i-1]*i%m
        inv[i]=m-inv[m%i]*(m//i)%m
        finv[i]=finv[i-1]*inv[i]%m
    return fac[n]*(finv[k]*finv[n-k]%m)%m
ans=ans*modcomb(n*m-2,k-2,INF)%INF
print(ans)
# Original solution
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)


# Optimized solution
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)
# Original solution
inpl = lambda: list(map(int,input().split()))
N, M ,K = inpl()
MOD = 10**9+7
factor = 1
# MN-2 C K-2 = MN-2 C MN-K = MN-2 C L
L = min(K-2,M*N-K)
for i in range(M*N-L-1,M*N-1):
    factor *= i
for i in range(1,L+1):
    factor //= i
factor %= MOD
ans = ((M*(N-1)*(N+1) + N*(M-1)*(M+1))*M*N //6 * factor) % MOD
print(ans)

# Optimized solution
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
# MN-2 C K-2 = MN-2 C MN-K = MN-2 C L
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)