CodeEff/ECCO huggingface dataset:
# Original solution
def egcd(a, b):
    (x, lastx) = (0, 1)
    (y, lasty) = (1, 0)
    while b != 0:
        q = a // b
        (a, b) = (b, a % b)
        (x, lastx) = (lastx - q * x, x)
        (y, lasty) = (lasty - q * y, y)
    return (lastx, lasty, a)

# ax ≡ 1 (mod m)なるxを返す
def modinv(a):
    m=MOD
    (inv, q, gcd_val) = egcd(a, m)
    return inv % m

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

b=1
for i in range(N*M-2,N*M-K,-1):
    b=(b*i)%MOD

for i in range(1,K-1):
    b=(b*modinv(i))%MOD

x=0
for i in range(1,N):
   x=x+i*(N-i)

x=(x*M*M)%MOD

y=0
for j in range(1,M):
    y=y+j*(M-j)
y=(y*N*N)%MOD
    
ans=b*(x+y)%MOD
print(ans)

# Optimized solution
def modinv(a):
    b=MOD
    (x, lastx) = (0, 1)
    (y, lasty) = (1, 0)
    while b != 0:
        q = a // b
        (a, b) = (b, a % b)
        (x, lastx) = (lastx - q * x, x)
        (y, lasty) = (lasty - q * y, y)
    return lastx%MOD

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

b=1
for i in range(1,K-1):
    b=(b*(N*M-1-i)*modinv(i))%MOD

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

y=0
for j in range(1,M):
    y=y+j*(M-j)
y=(y*N*N)%MOD
    
ans=b*(x+y)%MOD
print(ans)
# Original solution
def modinv(a):
    b=MOD
    (x, lastx) = (0, 1)
    (y, lasty) = (1, 0)
    while b != 0:
        q = a // b
        (a, b) = (b, a % b)
        (x, lastx) = (lastx - q * x, x)
        (y, lasty) = (lasty - q * y, y)
    return lastx%MOD

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

b=1
for i in range(1,K-1):
    b=(b*(N*M-1-i)*modinv(i))%MOD

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

y=0
for j in range(1,M):
    y=y+j*(M-j)
y=(y*N*N)%MOD
    
ans=b*(x+y)%MOD
print(ans)

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

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

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


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

y=0
for j in range(1,M):
    y=y+j*(M-j)
y=(y*N*N)%MOD
    
ans=comb(N*M-2,K-2)*(x+y)%MOD
print(ans)
# Original solution
def comb(n, k):
    x, y = 1, 1
    k = min(k, n-k)
    for i in range(k):
        x *= n - i
        x %= MOD
        y *= i + 1
        y %= MOD
    return x*pow(y, MOD-2, MOD) % MOD

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

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


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

y=0
for j in range(1,M):
    y=y+j*(M-j)
y=(y*N*N)%MOD
    
ans=comb(N*M-2,K-2)*(x+y)%MOD
print(ans)

# Optimized solution
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*modinv(y) % MOD

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

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


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

y=0
for j in range(1,M):
    y=y+j*(M-j)
y=(y*N*N)%MOD
    
ans=comb(N*M-2,K-2)*(x+y)%MOD
print(ans)
# Original solution
def cmb(a,b,c):
    b = min(b,a-b)
    num = 1
    for i in range(b):
        num = num*(a-i) % c
    den = 1
    for i in range(b):
        den = den*(i+1) % c
    return num * pow(den,c-2,c) % c
n,m,k = list(map(int,input().split()))
mod = 10**9+7
ans = 0
for i in range(1,n):
    ans += i*(n-i)*m**2*cmb(n*m-2,k-2,mod)
    ans %= mod
for i in range(1,m):
    ans += i*(m-i)*n**2*cmb(n*m-2,k-2,mod)
    ans %= mod
print(ans)


# Optimized solution
def cmb(a,b,c):
    b = min(b,a-b)
    num = 1
    for i in range(b):
        num = num*(a-i) % c
    den = 1
    for i in range(b):
        den = den*(i+1) % c
    return num * pow(den,c-2,c) % c
n,m,k = list(map(int,input().split()))
mod = 10**9+7
ans = 0
C = cmb(n*m-2,k-2,mod)
for i in range(1,n):
    ans += i*(n-i)*m**2*C
    ans %= mod
for i in range(1,m):
    ans += i*(m-i)*n**2*C
    ans %= mod
print(ans)
# Original solution
N,M,K=list(map(int,input().split()))
import math
def com(n, r):
  return math.factorial(n) // (math.factorial(n - r) * math.factorial(r))
A=com(M*N-2, K-2)
MOD=10**9+7
ans=0
for d in range(1,N):
  ans+=d*(N-d)*M*M*A%MOD
for d in range(1,M):
  ans+=d*(M-d)*N*N*A%MOD
print((ans%MOD))

# Optimized solution
N,M,K=list(map(int,input().split()))
def cmb(x,y):
  X=1
  Y=1
  for i in range(1,y+1):
    X= X*(x+1-i)%MOD
    Y= Y*i%MOD
  return (X*pow(Y,MOD-2,MOD))%MOD
MOD=10**9+7
A=cmb(M*N-2, K-2)
ans=0
for d in range(1,N):
  ans+=d*(N-d)*M*M*A%MOD
for d in range(1,M):
  ans+=d*(M-d)*N*N*A%MOD
print((ans%MOD))
# Original solution
n,m,k=list(map(int,input().split()))
mod=10**9+7
f=[1]
for i in range(2*10**5+7):f+=[f[-1]*(i+1)%mod]
def comb(a,b):return f[a]*pow(f[b],mod-2,mod)*pow(f[a-b],mod-2,mod)%mod
ans=0
for i in range(1,n):
  ans+=(i*(n-i)*m**2)*(comb(n*m-2,k-2))
  ans%=mod
for i in range(1,m):
  ans+=(i*(m-i)*n**2)*(comb(n*m-2,k-2))
  ans%=mod
print(ans)

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

mod = 10**9+7
f = [1]
for i in range(n*m):
  f += [f[-1]*(i+1) % mod]
def comb(a, b):
  return f[a]*pow(f[b], mod-2, mod)*pow(f[a-b], mod-2, mod) % mod

ans=0
for dist in range(1,n):
  c=n-dist
  p=(c*m*m)%mod
  ans+=p*dist
for dist in range(1,m):
  c=m-dist
  p=(c*n*n)%mod
  ans+=p*dist
ans%=mod
print((ans*comb(n*m-2,k-2)%mod))
# Original solution
import math

MOD = int(10**9 + 7)

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

def modinv(a, m):
    g, x, y = egcd(a, m)
    if g != 1:
        raise Exception("No inverse exists")
    else:
        return x % m

def combination(n, r):
    return math.factorial(n) * modinv(math.factorial(r), MOD) * modinv(math.factorial(n-r), MOD)
    # return math.factorial(n) // math.factorial(r) // math.factorial(n-r)

N, M, K = [int(x) for x in input().split()]

ans = 0

k = combination(N*M-2, K-2)

for d in range(1, N):
    ans += d * k * M * M * (N - d)

for d in range(1, M):
    ans += d * k * N * N * (M - d)

print((ans % MOD))


# Optimized solution
import math

MOD = int(10**9 + 7)

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

def modinv(a, m):
    g, x, y = egcd(a, m)
    if g != 1:
        raise Exception("No inverse exists")
    else:
        return x % m

def combination(n, r):
    return math.factorial(n)%MOD * modinv(math.factorial(r), MOD) * modinv(math.factorial(n-r), MOD) % MOD
    # return math.factorial(n) // math.factorial(r) // math.factorial(n-r)

N, M, K = [int(x) for x in input().split()]

ans = 0

k = combination(N*M-2, K-2)

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

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

print((ans % MOD))
# Original solution
import sys
input = sys.stdin.readline
import math

N, M, K = [int(item) for item in input().split()]
mod = 10**9 + 7

total = 0
for d in range(1, N):
    score = (N - d) * M**2 * d
    score %= mod
    for i in range(2, K):
        score *= ((N * M - i) % mod * pow(K - i, mod - 2, mod))
        score %= mod
    total += score
    total %= mod

for d in range(1, M):
    score = (M - d) * N**2 * d
    score %= mod
    for i in range(2, K):
        score *= ((N * M - i) % mod * pow(K - i, mod - 2, mod))
        score %= mod
    total += score
    total %= mod


print(total)


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

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

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

res = 0

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

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

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

# Optimized solution
mod = 10**9 + 7
n,m,k = list(map(int,input().split()))
nm = n*m
f = [1]*(nm+1)
for i in range(1,nm+1):
    f[i] = (f[i-1]*i)%mod
def C(n,k):
    x = f[n]
    x *= pow(f[n-k],mod-2,mod)
    x %= mod
    x *= pow(f[k],mod-2,mod)
    x %= mod
    return x
def F(n):
    x = n*n*(n+1)//2
    y = n*(n+1)*(2*n+1)//6
    return (x-y)%mod
res = (m*m%mod*F(n)%mod + n*n%mod*F(m)%mod)%mod
res *= C(n*m-2,k-2)
print((res%mod))
# 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()))

ans = 0

for i in range(1, M):
    ans += (i * (M-i) * N**2) % (10**9 + 7) 
for i in range(1, N):
    ans += (i * (N-i) * M**2) % (10**9 + 7)

ans *= cmb((N*M-2)%(10**9+7) , K-2%(10**9+7))

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


# Optimized solution
def framod(n, mod, a=1):
    for i in range(1,n+1):
        a=a * i % mod
    return a

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

def comb(n, k, mod):
    a=framod(n, mod)
    b=framod(k, mod)
    c=framod(n-k, mod)
    return (a * power(b, mod-2, mod) * power(c, mod-2, mod)) % mod

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

ans = 0

for i in range(1, M):
    ans += (i * (M-i) * N**2) % (10**9 + 7) 
for i in range(1, N):
    ans += (i * (N-i) * M**2) % (10**9 + 7)

ans *= comb(N*M-2 , K-2, 10**9+7)

print((int(ans % (10**9+7))))
# Original solution
MOD = 10**9 + 7

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

invs = [1] * (K-1)
for x in range(2, K-1):
    invs[x] = (-(MOD//x) * invs[MOD%x]) % MOD

C = 1
for x in range(N*M-K+1, N*M-1):
    C *= x
    C %= MOD
for x in range(1, K-1):
    C *= invs[x]
    C %= MOD

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

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

print(((ans1+ans2) * C % MOD))


# Optimized solution
MOD = 10**9 + 7

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

def getComb(n, k, MOD):
    if n < k:
        return 0
    if n-k < k:
        k = n-k
    comb = 1
    for x in range(n-k+1, n+1):
        comb = (comb * x) % MOD
    d = 1
    for x in range(1, k+1):
        d = (d * x) % MOD
    comb *= pow(d, MOD-2, MOD)
    return comb % MOD

comb = getComb(N*M-2, K-2, MOD)

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

for j in range(M):
    ans += j*(M-j)*N*N
    ans %= MOD

ans *= comb
ans %= MOD

print(ans)
# Original solution
import math

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

def nCm(n,m):
    if m < 0 or m > n:
        return 0
    nf = math.factorial(n)
    nm = math.factorial(m)
    nm2 = math.factorial(n-m)
    return (nf//nm//nm2) % 1000000007

c_sum = (M*(N-1)*N*(N+1)*M+N*(M-1)*M*(M+1)*N)//3
c_sum = c_sum//2
c_sum *= nCm(N*M-2, K-2)
c_sum = c_sum % 1000000007

print(c_sum)


# Optimized solution
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)
# Original solution
def combi(n, r) :
    ret = 1
    for i in range(n, n-r, -1) :
        ret *= i
        while r > 1 and ret % r == 0 :
            ret //= r
            r -= 1

    return ret

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

n, r = N * M - 2, K - 2
P = combi(n, r) % MOD

ans = 0
for n in range(N) :
    for m in range(M) :
        if n == 0 and m == 0 :
            continue

        d = n + m
        cnt = (N - n) * (M - m) * (1 if n == 0 or m == 0 else 2)

        ans += d * cnt * P
        ans %= MOD

print(ans)


# Optimized solution
from math import factorial

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

n, r = N * M - 2, K - 2
#P = factorial(n) // factorial(n-r) // factorial(r) % MOD
P = factorial(n) * pow(factorial(n-r) % MOD, MOD - 2, MOD) * pow(factorial(r) % MOD, MOD - 2, MOD)

ans = 0
for n in range(N) :
    for m in range(M) :
        if n == 0 and m == 0 :
            continue

        d = n + m
        cnt = (N - n) * (M - m) * (1 if n == 0 or m == 0 else 2)

        ans += d * cnt

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

MOD = 10**9+7
fact = [1 for i in range(N*M+1)]
for i in range(2, N*M+1):
    fact[i] = fact[i-1] * i

def comb(n, r):
    if n == 0 or r == 0:
        return 1
    return fact[n] // fact[n-r] // fact[r]

p = comb(N*M-2, K-2)
ans = 0
for i in range(N):
    for j in range(M):
        if i == 0 and j == 0:
            continue
        d = i+j
        #print(i, j, calc_rec_num(N, M, 0, 0, i, j))
        cnt = (N-i) * (M-j)
        if i != 0 and j != 0:
            cnt *= 2
        ans += d * cnt
        #if ans > DIV:
        #    ans %= DIV
        #print(ans)

print(((ans*p)%MOD))


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

MOD = 10**9+7

def comb(n, r):
    r = n-r if n-r < r else r
    if r == 0:
        return 1
    ndfact = 1
    for i in range(n, n-r, -1):
        ndfact *= i
    return ndfact // factorial(r)

p = comb(N*M-2, K-2)
ans = 0
for i in range(N):
    for j in range(M):
        if i == 0 and j == 0:
            continue
        d = i+j
        #print(i, j, calc_rec_num(N, M, 0, 0, i, j))
        cnt = (N-i) * (M-j)
        if i != 0 and j != 0:
            cnt *= 2
        ans += d * cnt
        #if ans > DIV:
        #    ans %= DIV
        #print(ans)

print(((ans*p)%MOD))
# Original solution
from math import factorial
N, M, K = list(map(int, input().split()))

MOD = 10**9+7

def comb(n, r):
    r = n-r if n-r < r else r
    if r == 0:
        return 1
    ndfact = 1
    for i in range(n, n-r, -1):
        ndfact *= i
    return ndfact // factorial(r)

p = comb(N*M-2, K-2) % MOD
ans = 0
for i in range(N):
    for j in range(M):
        if i == 0 and j == 0:
            continue
        d = i+j
        cnt = (N-i) * (M-j)
        if i != 0 and j != 0:
            cnt *= 2
        ans += d * cnt
        ans %= MOD

print(((ans*p)%MOD))


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

MOD = 10**9+7

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

def comb(n, r):
    r = n-r if n-r < r else r
    return fact[n] * pow(fact[r], MOD-2, MOD) * pow(fact[n-r], MOD-2, MOD) % MOD

p = comb(N*M-2, K-2)
ans = 0
for i in range(N):
    for j in range(M):
        if i == 0 and j == 0:
            continue
        d = i+j
        cnt = (N-i) * (M-j)
        if i != 0 and j != 0:
            cnt *= 2
        ans += d * cnt % MOD
ans = ans * p % MOD

print(ans)
# Original solution
def com(n,r):
    if n-r<r:
        r=n-r
    re=1
    for i in range(n,n-r,-1):
        re*=i
    for i in range(1,r+1):
        re//=i
    return re

def f(n,m,k):
    md = 10 ** 9 + 7
    ans = 0
    for d in range(1, m):
        ans = (ans + d * (m - d) * n * n) % md
    for d in range(1, n):
        ans = (ans + d * (n - d) * m * m) % md
    ans = (ans * com(n * m - 2, k - 2)) % md
    print(ans)

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


# Optimized solution
def com(n,r):
    re=fct[n]*pow(fct[r],md-2,md)*pow(fct[n-r],md-2,md)
    re%=md
    return re

n,m,k=list(map(int,input().split()))
md=10**9+7
ans=0
fct=[1]*(n*m-1)
f=1
for i in range(1,n*m-1):
    f=f*i%md
    fct[i]=f
for d in range(1,m):
    ans=(ans+d*(m-d)*n*n)%md
for d in range(1,n):
    ans=(ans+d*(n-d)*m*m)%md
ans=(ans*com(n*m-2,k-2))%md
print(ans)
# Original solution
from math import factorial
mod=10**9+7
def cmb(n, r):
    return factorial(n) // factorial(r) // factorial(n - r)%mod

n,m,k=list(map(int,input().split()))
print((((pow(m,2,mod)*cmb(n*m-2,k-2)*sum((n-i)*i%mod for i in range(1,n))%mod)+pow(n,2,mod)*cmb(n*m-2,k-2)*sum((m-i)*i%mod for i in range(1,m))%mod)%mod))

# Optimized solution
mod=10**9+7
def mod_combination(n, k, mod):
    def ext_gcd(a, b):
        if b == 0:return a, 1, 0
        else:d,x,y = ext_gcd(b,a%b);x-=(a//b)*y 
        return d,y,x
    p,q=1,1
    for i in range(n-k+1,n+1):
        p=(p*i)%mod
    for i in range(2,k+1):
        q=(q*i)%mod
    return int(p*(ext_gcd(q, mod)[1]%mod)%mod)


n,m,k=list(map(int,input().split()))
print((((pow(m,2,mod)*mod_combination(n*m-2,k-2,mod)*sum((n-i)*i%mod for i in range(1,n))%mod)+pow(n,2,mod)*mod_combination(n*m-2,k-2,mod)*sum((m-i)*i%mod for i in range(1,m))%mod)%mod))
# Original solution
# https://atcoder.jp/contests/abc127/tasks/abc127_e

import itertools
from collections import Counter
from collections import defaultdict
import bisect
import math

# Calculate count of combination
def combination(n, r):
    if r == 0:
        return 1
    a = 1
    b = 1
    for i in range(r):
        a *= (n - i)
        if a % (i + 1):
            a //= (i + 1)
        else:
            b *= (i + 1)
    return a // b


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

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

    ans = 0

    for dx in range(1, M):
        ans += dx * (M - dx) * pattern * (N * N)
        ans = ans % MOD

    # y
    for dy in range(1, N):
        ans += dy * (N - dy) * pattern * (M * M)
        ans = ans % MOD

    ans = ans % MOD
    print(ans)



if __name__ == '__main__':
    main()


# Optimized solution
# https://atcoder.jp/contests/abc127/tasks/abc127_e

import itertools
from collections import Counter
from collections import defaultdict
import bisect
import math

# Calculate count of combination
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


def combination(n, r):
    if r == 0:
        return 1
    a = 1
    b = 1
    for i in range(r):
        a *= (n - i)
        if a % (i + 1):
            a //= (i + 1)
        else:
            b *= (i + 1)
    return a // b


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

    pattern = cmb(N * M - 2, K - 2)
    # print(com)
    pattern = pattern % MOD
    # print(com)

    ans = 0

    for dx in range(1, M):
        ans += dx * (M - dx) * pattern * (N * N)
        ans = ans % MOD

    # y
    for dy in range(1, N):
        ans += dy * (N - dy) * pattern * (M * M)
        ans = ans % MOD

    ans = ans % MOD
    print(ans)


if __name__ == '__main__':
    main()
# Original solution
import sys
input = sys.stdin.readline
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()))
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 = cmb(N*M-2, K-2) * (x + y) % mod
print(ans)

# Optimized solution
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 = list(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)
# Original solution
import itertools
import sys
from functools import reduce
from operator import mul

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


def comb(n, r):
    """
    scipy.misc.comb または scipy.special.comb と同じ
    組み合わせの数 nCr
    :param int n:
    :param int r:
    :rtype: int
    """
    assert n >= r
    r = min(n - r, r)
    if r == 0:
        return 1
    return reduce(mul, list(range(n, n - r, -1))) // reduce(mul, list(range(r, 0, -1)))


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

def pyramid(n):
    return n * (n + 1) * (n + 2) // 6


# K ==2 の場合
# print(pyramid(N - 1) * M ** 2 + pyramid(M - 1) * N ** 2)
a = comb(N * M-2, K-2)
b = pyramid(N - 1) * M ** 2 + pyramid(M - 1) * N ** 2
print((a*b% (10**9+7)))


# Optimized solution
import math
import sys

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


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


def comb(n, r, mod):
    # https://fushime2.hatenablog.com/entry/2016/11/12/155543
    return (math.factorial(n) * inv(math.factorial(n - r)) * inv(math.factorial(r))) % mod


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


def pyramid(n):
    return n * (n + 1) * (n + 2) // 6


MOD = 10 ** 9 + 7

# K ==2 の場合
# print(pyramid(N - 1) * M ** 2 + pyramid(M - 1) * N ** 2)
k2 = (pyramid(N - 1) * M ** 2 + pyramid(M - 1) * N ** 2) % MOD
# k2 を何回使うか
c = comb(N * M - 2, K - 2, mod=MOD)
print((k2 * c % MOD))
# Original 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
        y *= i + 1
    return x*pow(y, mod-2, mod) % mod

N, M, K = list(map(int, input().split()))
print(((sum(i*(N-i)*M**2 for i in range(N))+sum(i*(M-i)*N**2 for i in range(M)))*comb(N*M-2, K-2) % mod))

# 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()))
print(((sum(i*(N-i)*M**2 for i in range(N))+sum(i*(M-i)*N**2 for i in range(M)))*comb(N*M-2, K-2) % mod))
# Original 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()))
print(((sum(i*(N-i)*M**2 for i in range(N))+sum(i*(M-i)*N**2 for i in range(M)))*comb(N*M-2, K-2) % mod))

# 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()))
print(((M*M%mod*((N**3-N)//6%mod)%mod+N*N%mod*((M**3-M)//6%mod)%mod)%mod*comb(N*M-2,K-2)%mod))