Buggy python submission (verdict: Runtime Error):
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)*L%mod

print(cou)

Fixed submission:
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 = 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)

Buggy python submission (verdict: Time Limit Exceeded):
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, range(n, n - r, -1))
    under = reduce(mul, 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)

Fixed submission:
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 = 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)

Buggy python submission (verdict: Time Limit Exceeded):
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, range(n, n - r, -1))
    under = reduce(mul, 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)

Fixed submission:
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 = 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)

Buggy python submission (verdict: Time Limit Exceeded):
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, range(n, n - r, -1))
    under = reduce(mul, range(1, r + 1))
    return over // under

N,M,K = 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)


Fixed submission:
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 = 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)

Buggy python submission (verdict: Time Limit Exceeded):
from math import factorial

def comb(n, r):
    return factorial(n) // (factorial(n - r) * factorial(r))

N,M,K = 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)


Fixed submission:
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 = 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)

Buggy python submission (verdict: Runtime Error):
from scipy.special import comb

N,M,K = 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)
print(ans)

Fixed submission:
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 = 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)