CodeEff/ECCO huggingface dataset:
# Original solution
def combination(n, r):
    r = min(n - r, r)
    result = 1
    for i in range(n, n - r, -1):
        result *= i
    for i in range(1, r + 1):
        result //= i
    return result


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

NM = N * M
fa = [1] * (NM + 1)

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

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

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

ans *= combination(NM-2, K-2)
ans %= mod
print(ans)


# Optimized solution
class FactMod:

    def __init__(self, n, mod):
        self.mod = mod
        self.f = [1] * (n + 1)
        for i in range(1, n + 1):
            self.f[i] = self.f[i - 1] * i % mod

        self.inv = [pow(self.f[-1], mod - 2, mod)]
        for i in reversed(list(range(1, n + 1))):
            self.inv.append(self.inv[-1] * i % mod)
        self.inv.reverse()

    def fact(self, n):
        """
        :return: n!
        """
        return self.f[n]

    def comb(self, n, r):
        """
        :return: nCr
        """
        return self.f[n] * self.inv[n - r] * self.inv[r] % self.mod

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

NM = N * M
fa = [1] * (NM + 1)

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

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

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

factmod = FactMod(NM-2, mod)
ans *= factmod.comb(NM-2, K-2)
ans %= mod
print(ans)
# Original solution
class FactMod:

    def __init__(self, n, mod):
        self.mod = mod
        self.f = [1] * (n + 1)
        for i in range(1, n + 1):
            self.f[i] = self.f[i - 1] * i % mod

        self.inv = [pow(self.f[-1], mod - 2, mod)]
        for i in reversed(list(range(1, n + 1))):
            self.inv.append(self.inv[-1] * i % mod)
        self.inv.reverse()

    def fact(self, n):
        """
        :return: n!
        """
        return self.f[n]

    def comb(self, n, r):
        """
        :return: nCr
        """
        return self.f[n] * self.inv[n - r] * self.inv[r] % self.mod

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

NM = N * M
fa = [1] * (NM + 1)

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

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

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

factmod = FactMod(NM-2, mod)
ans *= factmod.comb(NM-2, K-2)
ans %= mod
print(ans)


# Optimized solution
class FactMod:

    def __init__(self, n, mod):
        self.mod = mod
        self.f = [1] * (n + 1)
        for i in range(1, n + 1):
            self.f[i] = self.f[i - 1] * i % mod

        self.inv = [pow(self.f[-1], mod - 2, mod)]
        for i in reversed(list(range(1, n + 1))):
            self.inv.append(self.inv[-1] * i % mod)
        self.inv.reverse()

    def fact(self, n):
        """
        :return: n!
        """
        return self.f[n]

    def comb(self, n, r):
        """
        :return: nCr
        """
        return self.f[n] * self.inv[n - r] * self.inv[r] % self.mod


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

NM = N * M

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

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

factmod = FactMod(NM - 2, mod)
ans *= factmod.comb(NM - 2, K - 2)
ans %= mod
print(ans)
# Original solution
import sys
sys.setrecursionlimit(2147483647)
INF=float("inf")
MOD=10**9+7
input=lambda:sys.stdin.readline().rstrip()

class modfact(object):
    def __init__(self,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
        self.__fact=fact; self.__invfact=invfact

    def inv(self,n):
        assert(n>0)
        return self.__fact[n-1]*self.__invfact[n]%MOD

    def fact(self,n):
        return self.__fact[n]

    def invfact(self,n):
        return self.__invfact[n]

    def comb(self,n,k):
        if(k<0 or n<k): return 0
        return self.__fact[n]*self.__invfact[k]*self.__invfact[n-k]%MOD

    def perm(self,n,k):
        if(k<0 or n<k): return 0
        return self.__fact[n]*self.__invfact[n-k]%MOD

def resolve():
    n,m,k=list(map(int,input().split()))
    N=n*m
    ans=0
    mf=modfact(N)
    for d in range(1,n):
        ans+=d*(m**2)*(n-d)*mf.comb(N-2,k-2)
        ans%=MOD
    for d in range(1,m):
        ans+=d*(n**2)*(m-d)*mf.comb(N-2,k-2)
        ans%=MOD
    print(ans)
resolve()

# Optimized solution
import sys
sys.setrecursionlimit(2147483647)
INF = 1 << 60
MOD = 10**9 + 7 # 998244353
input = lambda:sys.stdin.readline().rstrip()

class modfact(object):
    def __init__(self, n):
        fact, invfact = [1] * (n + 1), [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
        self._fact, self._invfact = fact, invfact

    def inv(self, n):
        return self._fact[n - 1] * self._invfact[n] % MOD

    def fact(self, n):
        return self._fact[n]

    def invfact(self, n):
        return self._invfact[n]

    def comb(self, n, k):
        if k < 0 or n < k: return 0
        return self._fact[n] * self._invfact[k] % MOD * self._invfact[n - k] % MOD

    def perm(self, n, k):
        if k < 0 or n < k: return 0
        return self._fact[n] * self._invfact[n - k] % MOD

def resolve():
    m, n, k = list(map(int, input().split()))
    mf = modfact(m * n - 2)
    res = 0
    for s in range(1, m):
        res += n * n * (m - s) * s
        res %= MOD
    for t in range(1, n):
        res += m * m * (n - t) * t
        res %= MOD
    res *= mf.comb(m * n - 2, k - 2)
    res %= MOD
    print(res)
resolve()
# Original solution
from math import factorial

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

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


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)

# Optimized solution
from math import factorial

N, M, K = list(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)
# Original solution
N,M,K=list(map(int,input().split()))
P=1000000007
import math
ans=M**2*(N-1)*N*(N+1)//6
ans+=N**2*(M-1)*M*(M+1)//6
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))


# Optimized solution
N,M,K=list(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))
# Original solution
import sys
input = sys.stdin.readline

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

def main():
    mod=10**9+7
    N,M,K=MI()
    #2この組み合わせを全部足して，平均ぽく考える
    
    def sum2(x):
        return ((x*(x+1))//2)%mod
    
    #i,jますにいるときの相手．（i,j）よりも奥にあるやつの全組み合わせの総和
    def calc(i,j):
        temp=0
        temp=sum2(M-j-1)#該当する行
        temp=(temp+((sum2(M-j-1)+sum2(j))*(N-i-1)))%mod#下の行の階差部分，左右方向にあるから
        temp=(temp+sum2(N-i-1)*M)%mod
        return temp
        
    S=0
    for i in range(N):
        for j in range(M):
            S=(S+calc(i,j))%mod
            
    rem=(N*M)-2
    
    def cmb(n, r, mod):
        if (r < 0) or (n < r):
            return 0
        r = min(r, n - r)
        return (fact[n] * factinv[r] * factinv[n-r])%mod

    fact=[1,1]
    factinv=[1,1]
    inv=[0,1]
    
    for i in range(2, rem + 1):
        fact.append((fact[-1] * i) % mod)
        inv.append((-inv[mod % i] * (mod // i)) % mod)
        factinv.append((factinv[-1] * inv[-1]) % mod)  
        
    #print(rem,K-2)
    S=(S*cmb(rem,K-2,mod))%mod
            
    print(S)

main()


# Optimized solution
import sys
input = sys.stdin.readline

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

def main():
    mod=10**9+7
    """
    まず，xとyを分離する=>xだけを考えるので1行とみなす．
    いっぺんにいくつも考えられないので，コマを2つだけ置いた時を考える．
    2この組み合わせで距離が1=>M-1通り，2=>M-2通り,...,M-1=>1通り
    (以下のように考えてΣしても良いがやや面倒
        (j列目にコマを固定すると，相手との距離はj,j-1,...,1,0,1,2,...,M-j))
    実際にはN列あるので，選び方がN^2
        同じ行・列を選んだ場合は不可能だけど，どうせ距離も0なので関係ない
        
    ある2つの組みに注目した時，それを何回数えるかは，外野の並べ方に依存
    (N*M-2)C(K-2)通り
    """
    N,M,K=MI()
    
    
    def cmb(n, r, mod):
        if (r < 0) or (n < r):
            return 0
        r = min(r, n - r)
        return (fact[n] * factinv[r] * factinv[n-r])%mod

    fact=[1,1]
    factinv=[1,1]
    inv=[0,1]
    
    for i in range(2, N*M+5+ 1):
        fact.append((fact[-1] * i) % mod)
        inv.append((-inv[mod % i] * (mod // i)) % mod)
        factinv.append((factinv[-1] * inv[-1]) % mod)
    
    tempj=0
    for j in range(M):
        tempj+=(M-j)*j
        
    tempi=0
    for i in range(N):
        tempi+=(N-i)*i

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

main()
# Original solution
N,M,K = list(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)
    b *= i+1

comb = a//b
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)%mod
ans *= comb
print((ans%mod))

# Optimized solution
N,M,K = list(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)
# Original solution
N, M, K = list(map(int, input().split()))
MOD = 10 ** 9 + 7

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

    a, b = 1, 1
    for i in range(r) :
        a *= n - i
        b *= i + 1

    return a // b

res = 0

for i in range(1, N + 1) :
    for j in range(1, M + 1) :
        if i > 1 and j > 1 :
            res += 2 * (N + 1 - i) * (M + 1 - j) * (i + j - 2)
        else :
            res += (N + 1 - i) * (M + 1 - j) * (i + j - 2)

res *= comb(N * M - 2, K - 2)

print((res % MOD))

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

def pow(x, n) :
    y = 1
    while n > 0 :
        if n & 1 :
            y = (y * x) % MOD
        x = (x * x) % MOD
        n >>= 1

    return y

def comb(n, r) :
    if n < r or n < 0 :
        return 0

    factorial = [0] * (n + 1)
    factorial[0] = 1

    for i in range(1, n + 1) :
        factorial[i] = (factorial[i-1] * i) % MOD

    return (factorial[n] * pow(factorial[n-r], MOD - 2) * pow(factorial[r], MOD - 2)) % MOD

res = 0

for i in range(N) :
    for j in range(M) :
        if i > 0 and j > 0 :
            res += 2 * (N - i) * (M - j) * (i + j)
        else :
            res += (N - i) * (M - j) * (i + j)

res *= comb(N * M - 2, K - 2)

print((res % MOD))
# Original solution
n, m, k = list(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

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

        

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


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

A=2*10**5
mod = 10**9+7
g1 = [1, 1]
g2 = [1, 1]
inverse = [0, 1]

for i in range( 2, A + 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=list(map(int,input().split()))
keisuu=cmb(N*M-2,K-2,mod)

sum_=0
for i in range(N):
    a=min(abs(i),abs(N-i-1))
    b=max(abs(i),abs(N-i-1))
    sum_+=(M**2)*((a*(a+1)//2)+(b*(b+1)//2))
for i in range(M):
    a=min(abs(i),abs(M-i-1))
    b=max(abs(i),abs(M-i-1))
    sum_+=(N**2)*((a*(a+1)//2)+(b*(b+1)//2))
print(((keisuu * (sum_//2))%mod))

# Optimized solution
def power_mod(a,b,mod=10**9+7):
    i,temp,box=0,b,[]
    while(2**i<=b):
        i+=1
    for j in range(i-1,-1,-1):
        box=[[j,temp//2**j]]+box
        temp-=2**j*(temp//2**j)
    box[0].append(a)
    ans=box[0][1]*a%mod
    for j in range(1,i):
        box[j].append(box[j-1][2]**2%mod)
        if box[j][1]==1:
            ans=(ans*box[j][2])%mod
    return ans
def n_func(n,mod=10**9+7):
    ans=1
    for i in range(1,n+1):
        ans=(ans*i)%mod
    return ans
def nPr(n,r,mod=10**9+7):
    ans=n_func(n-r,mod)
    ans=power_mod(ans,mod-2,mod)
    return ans*n_func(n,mod)%mod
def nCr(n,r,mod=10**9+7):
    ans=n_func(n-r,mod)*n_func(r,mod)%mod
    ans=power_mod(ans,mod-2,mod)
    return ans*n_func(n,mod)%mod

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

sum_=0
for i in range(N):
    a=min(abs(i),abs(N-i-1))
    b=max(abs(i),abs(N-i-1))
    sum_+=(M**2)*((a*(a+1)//2)+(b*(b+1)//2))
for i in range(M):
    a=min(abs(i),abs(M-i-1))
    b=max(abs(i),abs(M-i-1))
    sum_+=(N**2)*((a*(a+1)//2)+(b*(b+1)//2))
print(((keisuu * (sum_//2))%mod))
# Original solution
#E問題
N,M,K = list(map(int,input().split()))
mod = 1000000007

#factorial(階乗)
def factorial(n):
    factorial_list=[]
    for i in range(n+1):
        if i == 0:
            factorial_list.append(1)
        elif i == 1:
            factorial_list.append(1)
        else:
            factorial_list.append(factorial_list[-1]*i)
    return factorial_list[-1]
#combination(組み合わせC)
def comb(comb1,comb2):
    x=min(comb1,comb2)
    y=max(comb1,comb2)
    bunshi=factorial(y)
    bunbo=factorial(x)*factorial(y-x)
    return bunshi//bunbo

C = comb(N*M-2,K-2)


X = 0
for i in range(N-1):
    d = i+1
    X+=(N-d)*M*M*d
    X%=mod
Y = 0
for i in range(M-1):
    d = i+1
    Y+=(M-d)*N*N*d
    Y%=mod
print((C*(X+Y)%mod))

# Optimized solution
#E問題
N,M,K = list(map(int,input().split()))
mod = 1000000007

zen = 1
for i in range(N*M-2):
    zen*=(i+1)
kk = 1
for i in range(K-2):
    kk*=(i+1)
zk = 1
for i in range(N*M-K):
    zk*=(i+1)
    
C = zen//(kk*zk)


X = 0
for i in range(N-1):
    d = i+1
    X+=(N-d)*M*M*d
    X%=mod
Y = 0
for i in range(M-1):
    d = i+1
    Y+=(M-d)*N*N*d
    Y%=mod
print((C*(X+Y)%mod))
# Original solution
mod = 1000000007
import sys
input = sys.stdin.readline
N, M, K = list(map(int, input().split()))

N2 = N**2
M2 = M**2
Xpattern = [(x * (N-x)) %mod for x in range(N)]
Ypattern = [(y * (M-y)) %mod for y in range(M)]
coef = 1
num = min(K-2, N*M-K)
for i in range(N*M-2, N*M-2-num, -1):
    coef *= i
for j in range(1, num+1):
    coef //= j
coef %= mod

print((((sum(Xpattern)*M2+sum(Ypattern)*N2)*coef)%mod))

# Optimized solution
mod = 1000000007
import sys
input = sys.stdin.readline
N, M, K = list(map(int, input().split()))


def cmb(n, r):
    if n<0 or r<0 or r>n:
        return None
    r = min(r, n-r)
    if r == 0:
        return 1
    if r == 1:
        return n
    numerator = [0]*r
    denominator = [0]*r
    for k in range(r):
        numerator[k] = n-r+k+1
        denominator[k] = k+1
    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
    res = 1
    for n in numerator:
        if n > 1:
            res *= n
    return res

N2 = N**2
M2 = M**2
Xpattern = [(x * (N-x)) % mod for x in range(N)]
Ypattern = [(y * (M-y)) % mod for y in range(M)]
coef = cmb(N*M-2, K-2)

print((((sum(Xpattern)*M2+sum(Ypattern)*N2)*coef)%mod))
# Original solution

"""
xi,xjのペアにおいて、このペアは何通り出てくるか
sum(全通り)sum(i)sum(j){|xi-xj|+|yi-yj|}
和なので、どんな順序でもいい

sum(i)sum(j)sum(全通り){|xi-xj|}
sum(i)sum(j)sum(全通り){|yi-yj|}

あるi,jで固定して、sum(全通り){|xi-xj|}を求めるヤツ
固定したi,jが現れる回数は、

pat = NM-2_C_K-2 こいつはどこでも一緒

pat*sum(i)sum(j){|xi-xj|}を計算する
すべてのi,j(i>j)ペアにおいて、xi,xjの和は
sum(i=1,N){i*(N-i)}
"""
from math import factorial

N,M,K = list(map(int,input().split()))
mod = 10**9+7
fa = 1
num = N*M - 2
# 固定したときの通り数
cut = min(K-2, N*M-K)
# このfor文は少ないほうが良い
for i in range(num,1,-1):
    if num-i == cut:
        break
    fa *= i
pat = fa//factorial(cut)

# print(pat)
# xi-xj
n_x = M**2*((N-1)*N*(N+1))//6
n_y = N**2*((M-1)*M*(M+1))//6
print((pat*(n_x+n_y)%mod))

# Optimized solution

"""
xi,xjのペアにおいて、このペアは何通り出てくるか
sum(全通り)sum(i)sum(j){|xi-xj|+|yi-yj|}
和なので、どんな順序でもいい

sum(i)sum(j)sum(全通り){|xi-xj|}
sum(i)sum(j)sum(全通り){|yi-yj|}

あるi,jで固定して、sum(全通り){|xi-xj|}を求めるヤツ
固定したi,jが現れる回数は、

pat = NM-2_C_K-2 こいつはどこでも一緒

pat*sum(i)sum(j){|xi-xj|}を計算する
すべてのi,j(i>j)ペアにおいて、xi,xjの和は
sum(i=1,N){i*(N-i)}
"""
from math import factorial

N,M,K = list(map(int,input().split()))
mod = 10**9+7
fa = 1
num = N*M - 2
# 固定したときの通り数
# フェルマーの小定理をつかって、mod(M)でのcombinationが高速で求まる
# http://keita-matsushita.hatenablog.com/entry/2016/12/05/184011
def mod_combination(n,r,mod):
    if n<0 or r<0 or n<r:
        return 0
    if n==0 or r == 0:
        return 1
    a = factorial(n) % mod
    b = factorial(r) % mod
    c = factorial(n-r) % mod
    return (a*pow(b,mod-2,mod)*pow(c,mod-2,mod)) % mod 

pat = mod_combination(N*M-2,K-2,mod)

n_x = M**2*((N-1)*N*(N+1))//6
n_y = N**2*((M-1)*M*(M+1))//6
print((pat*(n_x+n_y)%mod))
# Original solution
"""
組み合わせ(nCr)計算
https://qiita.com/derodero24/items/91b6468e66923a87f39f#%E7%95%AA%E5%A4%96%E7%B7%A8
逆元の求め方
http://drken1215.hatenablog.com/entry/2018/06/08/210000
解法@公式放送
https://www.youtube.com/watch?v=SS6kW-d-rJ0&feature=youtu.be&t=7376
"""


def cmb(n, r):
    return (fact[n] * finv_t[r] * finv_t[n - r]) % MOD


MOD = 10 ** 9 + 7

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

fact = [1, 1]
inv_t = [0, 1]
finv_t = [1, 1]
for i in range(2, N * M - 2 + 1):
    fact.append((fact[-1] * i) % MOD)
    inv_t.append((-inv_t[MOD % i] * (MOD // i)) % MOD)
    finv_t.append((finv_t[-1] * inv_t[-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)


# Optimized solution
"""
https://atcoder.jp/contests/abc127/submissions/5606400
"""


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
        # n!/(n-r)!

    d = 1
    for k in range(2, r + 1):
        d = (d * k) % MOD
        # r!

    return (f * inv(d)) % MOD


MOD = 10 ** 9 + 7

N, M, K = list(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)
# Original solution
# https://atcoder.jp/contests/abc127/submissions/5606400

from functools import reduce


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


def mod_add(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, list(range(n - r + 1, n + 1)))  # n!/(n-r)!
    d = reduce(mod_mul, list(range(2, r + 1)))  # r!
    return (f * pow(d, MOD - 2, MOD)) % MOD


MOD = 10 ** 9 + 7

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

ans_x = reduce(mod_add, (i * (M - i) for i in range(M))) * N ** 2
ans_y = reduce(mod_add, (i * (N - i) for i in range(N))) * M ** 2
ans = ans_x + ans_y
pat = cmb(N * M - 2, K - 2)
ans = (ans * pat) % MOD

print(ans)


# Optimized solution
# https://atcoder.jp/contests/abc127/submissions/5606400

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, list(range(n - r + 1, n + 1)))  # n!/(n-r)!
    d = reduce(mod_mul, list(range(2, r + 1)))  # r!
    return (f * pow(d, MOD - 2, MOD)) % MOD


MOD = 10 ** 9 + 7

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


# Optimized solution
def solve():
    mod = 10 ** 9 + 7

    def cmb(n, r):
        def build():
            fs = [1]
            t = 1
            for x in range(1, h * w):
                t = (t * x) % mod
                fs.append(t)

            invs = [1] * (h * w)
            t = pow(t, mod - 2, mod)
            invs[-1] = t
            for x in range(h * w - 2, 0, -1):
                t = (t * (x + 1)) % mod
                invs[x] = t
            return fs, invs

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

        fs, invs = build()
        return fs[n] * invs[n - r] * invs[r] % mod

    h, w, k = list(map(int, input().split()))
    c = cmb(h * w - 2, k - 2)
    ret = (pow(w, 2) * sum((h - d) * d for d in range(h))
           + pow(h, 2) * sum((w - d) * d for d in range(w))) * c % mod

    return ret


print((solve()))