CodeEff/ECCO huggingface dataset:
# Original 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()))

# Optimized solution
class Calc:
    def __init__(self, max_value, mod):
        """combination(max_value, all)"""
        fact = [-1] * (max_value + 1)
        fact[0] = 1
        fact[1] = 1
        for x in range(2, max_value + 1):
            fact[x] = x * fact[x - 1] % mod

        invs = [1] * (max_value + 1)
        invs[max_value] = pow(fact[max_value], mod - 2, mod)
        for x in range(max_value - 1, 0, -1):
            invs[x] = invs[x + 1] * (x + 1) % mod

        self.fact = fact
        self.invs = invs
        self.mod = mod

    def combination(self, n, r):
        if n - r < r:
            return self.combination(n, n - r)
        if r < 0:
            return 0
        if r == 0:
            return 1
        if r == 1:
            return n
        return self.fact[n] * self.invs[r] * self.invs[n - r] % self.mod


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

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

    calc = Calc(max_value=N * M - 2, mod=MOD)

    xsum = 0
    for d in range(1, N):
        xsum = (xsum + d * (N - d)) % MOD
    xsum = (xsum * M * M) % MOD

    ysum = 0
    for d in range(1, M):
        ysum = (ysum + d * (M - d)) % MOD
    ysum = (ysum * N * N) % MOD

    ans = (xsum + ysum) * calc.combination(N * M - 2, K - 2) % MOD

    print(ans)


if __name__ == '__main__':
    main()

# import sys
#
# sys.setrecursionlimit(10 ** 7)
#
# input = sys.stdin.readline
# rstrip()
# int(input())
# map(int, input().split())
# Original solution
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(N-1):
        for j in range(i+1,N):
            d = abs(i - j)
            ans += d * M**2
    for i in range(M-1):
        for j in range(i+1,M):
            d = abs(i - j)
            ans += d * N**2
    print(((ans*complement) % mod))

main()


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

MOD = 10 ** 9 + 7

NM = N * M

fact = [1 for _ in range(NM+1)]
for i in range(1, NM+1):
    fact[i] = fact[i-1] * i

def combination(n, m):
    return fact[n] // fact[n-m] //fact[m]

pattern = combination(NM-2, K-2) % MOD

def diff_x_pattern(val):
    # for d in range(1, val):
    #     ret += (val - d) * d
    return val * (val - 1) * val // 2 - (val - 1) * (2 * val - 1) * val // 6

ans = M * M * diff_x_pattern(N)
ans += N * N * diff_x_pattern(M)
print(((ans*pattern)%MOD))


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

MOD = 10 ** 9 + 7

def combination(n, m):
    m = min(m, n-m)
    c = 1
    for i in range(m):
        c *= n - i
    for i in range(2, m+1):
        c //= i
    return c % MOD

pattern = combination(N*M-2, K-2) % MOD

def diff_x_pattern(val):
    # for d in range(1, val):
    #     ret += (val - d) * d
    return (val * (val - 1) * val // 2 - (val - 1) * (2 * val - 1) * val // 6) % MOD

ans = (M * M * diff_x_pattern(N)) % MOD
ans += (N * N * diff_x_pattern(M)) % MOD
print(((ans*pattern)%MOD))
# Original solution
N, M, K = list(map(int, input().split()))

MOD = 10 ** 9 + 7

fact = [1 for _ in range(N*M+1)]
for i in range(1, N*M+1):
    fact[i] = fact[i-1] * i

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

def combination(n, m):
    return fact[n] * inv(fact[n-m]*fact[m]) % MOD

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


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


# Optimized solution
N, M, K = list(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))
# Original solution
import sys
from itertools import combinations

def comb_mod(n, r, MOD=10**9+7):
    if n < 0 or r < 0 or n < r:
        return 0

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

    fact_inv = [1] * (n+1)
    fact_inv[n] = pow(factrial[n], MOD - 2, MOD)
    for k in range(n-1, -1, -1):
        fact_inv[k] = (fact_inv[k+1] * (k+1)) % MOD

    return (factrial[n] * fact_inv[r] * fact_inv[n-r]) % MOD


def main():
    input = sys.stdin.readline
    MOD = 10**9 + 7
    N, M, K = list(map(int, input().split()))

    # The number of times that pieces are located on certain two tiles.
    num = comb_mod(N*M-2, K-2)

    ans = 0
    # Column-wise distance
    for a, b in combinations(list(range(M)), 2):
        ans += abs(a - b) * N**2 * num
        ans %= MOD

    # Row-wise distance
    for a, b in combinations(list(range(N)), 2):
        ans += abs(a - b) * M**2 * num
        ans %= MOD

    return ans


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


# Optimized solution
import sys
from itertools import combinations

def comb_mod(n, r, MOD=10**9+7):
    if n < 0 or r < 0 or n < r:
        return 0

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

    fact_inv = [1] * (n+1)
    fact_inv[n] = pow(factrial[n], MOD - 2, MOD)
    for k in range(n-1, -1, -1):
        fact_inv[k] = (fact_inv[k+1] * (k+1)) % MOD

    return (factrial[n] * fact_inv[r] * fact_inv[n-r]) % MOD


def main():
    input = sys.stdin.readline
    MOD = 10**9 + 7
    N, M, K = list(map(int, input().split()))

    # The number of times that pieces are located on certain two tiles.
    num = comb_mod(N*M-2, K-2)

    ans = 0
    # Column-wise distance
    for i in range(M):
        ans += i * (M - i) * N**2 * num
        ans %= MOD

    # Row-wise distance
    for i in range(N):
        ans += i * (N - i) * M**2 * num
        ans %= MOD

    return ans


if __name__ == '__main__':
    print((main()))
# Original solution
p,q,r=list(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
    for i in range(1,N-K+1):
        NK_fact*=i
    for i in range(1,N+1):
        N_fact*=i

    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

for k in range(0,p+1):
    ans+=k*(p-k)*q*q


print((ans*pat%(10**9+7)))

# Optimized solution
p,q,r=list(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)))
# Original solution
from operator import mul
from functools import reduce

def cmb(n,r):
    r = min(n-r,r)
    if r == 0: return 1
    over = reduce(mul, list(range(n, n - r, -1)))
    under = reduce(mul, list(range(1,r + 1)))
    return over // under

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

c = cmb(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)) % 1000000007

print(ans)

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

def expo_div_p(a: int,n: int,p: int):
    if n == 1:
        d = a
    elif n % 2 == 0:
        d = expo_div_p(a, n//2, p)**2 % p
    else:
        d = a * expo_div_p(a, n//2, p)**2 % p
    return d

def combi(x, y, p):
    c = 1
    for i in range(1,y+1):
        c *= (x-i+1) * expo_div_p(i, p-2, p)
        c %= p
    return c

ans = 0
C = 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(C)
print((ans % ppp))


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

#def expo_div_p(a: int,n: int,p: int):
#    if n == 1:
#        d = a
#    elif n % 2 == 0:
#        d = expo_div_p(a, n//2, p)**2 % p
#    else:
#        d = a * expo_div_p(a, n//2, p)**2 % p
#    return d
#
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(C)
print((ans % ppp))
# Original solution
from operator import mul
from functools import reduce

def comb(n,r):
    r = min(n - r, r)
    if r == 0: return 1
    over = reduce(mul, list(range(n, n - r, -1)))
    under = reduce(mul, list(range(1, r + 1)))
    return over // under

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

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

A = int(0)
for d in range(N):
    a = (N-d)*(M**2)

    A += int(d*a)

B = int(0)
for d in range(M):
    b = (M-d)*(N**2)
    B += int(d*b)

ans = int((A+B)*C)
ans = ans%(1000000007)
print(ans)


# Optimized solution
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])

    return result

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

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

A = int(0)
for d in range(N):
    a = (N-d)*(M**2)

    A += int(d*a)

B = int(0)
for d in range(M):
    b = (M-d)*(N**2)
    B += int(d*b)

ans = int((A+B)*C)
ans = ans%(1000000007)
print(ans)
# Original solution
from operator import mul
from functools import reduce
def comb2(n,r):
    r = min(n - r, r)
    if r == 0: return 1
    over = reduce(mul, list(range(n, n - r, -1)))
    under = reduce(mul, list(range(1, r + 1)))
    return over // under

import sys
def I(): return int(sys.stdin.readline())
def LI(): return list(map(int,sys.stdin.readline().split()))
mod = 10**9 + 7
inf = float('inf')
ans = int(0)

N, M, K = LI()

cou = int(0)
#maxL = N + M -2
L = int(0)
#距離i+1の組み合わせ

#for i in range(maxL):
for j in range(N):
    L += (N - j - 1) * (M ** 2) * (j + 1)
for j in range(M):
    L += (M - j - 1) * (N ** 2) * (j + 1)
#    for j in range(i+2):
#        if i>1 and i+1>j>0:
#            L[i] += max((N-j)*(M-(1+i)+j),0)*(i+1)*2
#        else:
#            L[i] += max(((N - j) * (M - (1 + i) + j)), 0) * (i + 1)

cou = comb2(N*M-2,K-2)%mod

ans = cou *L%mod

print(ans)

# Optimized solution
from operator import mul
from functools import reduce
def comb2(n,r):
    r = min(n - r, r)
    if r == 0: return 1
    over = reduce(mul, list(range(n, n - r, -1)))
    under = reduce(mul, list(range(1, r + 1)))
    return over // under
def comb3(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


import sys
def I(): return int(sys.stdin.readline())
def LI(): return list(map(int,sys.stdin.readline().split()))
mod = 10**9 + 7
inf = float('inf')
ans = int(0)

N, M, K = LI()

cou = int(0)
#maxL = N + M -2
L = int(0)
#距離i+1の組み合わせ

#for i in range(maxL):
for j in range(N):
    L += (N - j - 1) * (M ** 2) * (j + 1)
for j in range(M):
    L += (M - j - 1) * (N ** 2) * (j + 1)
#    for j in range(i+2):
#        if i>1 and i+1>j>0:
#            L[i] += max((N-j)*(M-(1+i)+j),0)*(i+1)*2
#        else:
#            L[i] += max(((N - j) * (M - (1 + i) + j)), 0) * (i + 1)

cou = comb3(N*M-2,K-2)%mod

ans = cou *L%mod

print(ans)
# Original solution
from operator import mul
from functools import reduce
def comb2(n,r):
    r = min(n - r, r)
    if r == 0: return 1
    over = reduce(mul, list(range(n, n - r, -1)))
    under = reduce(mul, list(range(1, r + 1)))
    return over // under
def comb3(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


import sys
def I(): return int(sys.stdin.readline())
def LI(): return list(map(int,sys.stdin.readline().split()))
mod = 10**9 + 7
inf = float('inf')
ans = int(0)

N, M, K = LI()

cou = int(0)
#maxL = N + M -2
L = int(0)
#距離i+1の組み合わせ

#for i in range(maxL):
for j in range(N):
    L += (N - j - 1) * (M ** 2) * (j + 1)
for j in range(M):
    L += (M - j - 1) * (N ** 2) * (j + 1)
#    for j in range(i+2):
#        if i>1 and i+1>j>0:
#            L[i] += max((N-j)*(M-(1+i)+j),0)*(i+1)*2
#        else:
#            L[i] += max(((N - j) * (M - (1 + i) + j)), 0) * (i + 1)

cou = comb3(N*M-2,K-2)%mod

ans = cou *L%mod

print(ans)

# Optimized solution
class combination():
    #素数のmod取るときのみ　速い
    def __init__(self, n, mod):
        self.n = n
        self.fac = [1] * (n + 1)
        self.inv = [1] * (n + 1)
        for j in range(1, n + 1):
            self.fac[j] = self.fac[j - 1] * j % mod

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

    def comb(self, n, r, mod):
        if r > n or n < 0 or r < 0:
            return 0
        return self.fac[n] * self.inv[n - r] * self.inv[r] % mod


import sys
def I(): return int(sys.stdin.readline())
def LI(): return list(map(int,sys.stdin.readline().split()))
mod = 10**9 + 7
inf = float('inf')
ans = int(0)

N, M, K = LI()

cou = int(0)
#maxL = N + M -2
L = int(0)
#距離i+1の組み合わせ

#for i in range(maxL):
for j in range(N):
    L += (N - j - 1) * (M ** 2) * (j + 1)
    L %= mod
for j in range(M):
    L += (M - j - 1) * (N ** 2) * (j + 1)
    L %= mod
#    for j in range(i+2):
#        if i>1 and i+1>j>0:
#            L[i] += max((N-j)*(M-(1+i)+j),0)*(i+1)*2
#        else:
#            L[i] += max(((N - j) * (M - (1 + i) + j)), 0) * (i + 1)
C = combination(N*M-2,mod)
cou = C.comb(N*M-2,K-2,mod)*L%mod

print(cou)
# Original solution
N,M,K=list(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*=i+1
print((S*K*~-K*pow(T,O-2,O)%O))

# Optimized solution
N,M,K=list(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))
# Original solution
import math
n,m,k = list(map(int,input().split()))
ss = (math.factorial(m*n)//math.factorial(k))//math.factorial(m*n-k)
s = (ss*k*(k-1))//2
r = 10**9 +7
xav = ((m**3-m) *n**2)//6
yav = ((n**3-n)*m**2)//6
av = xav + yav
d = (m*n*(m*n-1))//2
print((int(((av*s)//d)%r)))



# Optimized solution
n,m,k = list(map(int,input().split()))
x = m*n
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
ss = cmb(x,k)
s = (ss*k*(k-1))//2
r = 10**9 +7
av = (((m**3-m) *n**2)//6 + ((n**3-n)*m**2)//6)%r
d = (x*(x-1))//2
print((int(((av*s)//d)%r)))
# Original solution
# -*- coding: utf-8 -*-
mod = 10 ** 9 + 7
n, m, k = list(map(int, input().split()))
nm = n * m
fact = [1] * (nm + 1)

for i in range(1, nm + 1):
    fact[i] = fact[i - 1] * i


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)


# Optimized solution
# -*- coding: utf-8 -*-
mod = 10 ** 9 + 7
n, m, k = list(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)
# Original solution
mod = 10**9 + 7
n, m, k = list(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]
comb = fact[nm-2] * inv_fact[k-2] * inv_fact[nm-k] % mod
cost_x = (n-1) * n * (n+1) // 6 * m * m % mod
cost_y = (m-1) * m * (m+1) // 6 * n * n % mod
print((comb * (cost_x + cost_y) % mod))

# Optimized solution
mod = 10**9 + 7
n, m, k = list(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))
# Original solution
import sys
input = sys.stdin.readline

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

def Combi2(a,b):#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)

    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))/Combi2(N*M,2)*Combi2(N*M,K)*K*(K-1)//2%mod)
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))


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

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

def Combi2(a,b):#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))/Combi2(N*M,2)*Combi2(N*M,K)*K*(K-1)//2%mod)
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))