Problem:
 Score : 500 points    Problem Statement We have a grid of squares with N rows and M columns. Let (i, j) denote the square at the i -th row from the top and j -th column from the left. We will choose K of the squares and put a piece on each of them.  If we place the K pieces on squares (x_1, y_1) , (x_2, y_2) , ..., and (x_K, y_K) , the cost of this arrangement is computed as:  \sum_{i=1}^{K-1} \sum_{j=i+1}^K (|x_i - x_j| + |y_i - y_j|)  Find the sum of the costs of all possible arrangements of the pieces. Since this value can be tremendous, print it modulo 10^9+7 .  We consider two arrangements of the pieces different if and only if there is a square that contains a piece in one of the arrangements but not in the other.      Constraints  2 \leq N \times M \leq 2 \times 10^5  2 \leq K \leq N \times M  All values in input are integers.         Input Input is given from Standard Input in the following format:  N  M  K       Output Print the sum of the costs of all possible arrangements of the pieces, modulo 10^9+7 .        Sample Input 1 2 2 2      Sample Output 1 8  There are six possible arrangements of the pieces, as follows:   ((1,1),(1,2)) , with the cost |1-1|+|1-2| = 1  ((1,1),(2,1)) , with the cost |1-2|+|1-1| = 1  ((1,1),(2,2)) , with the cost |1-2|+|1-2| = 2  ((1,2),(2,1)) , with the cost |1-2|+|2-1| = 2  ((1,2),(2,2)) , with the cost |1-2|+|2-2| = 1  ((2,1),(2,2)) , with the cost |2-2|+|1-2| = 1   The sum of these costs is 8 .       Sample Input 2 4 5 4      Sample Output 2 87210       Sample Input 3 100 100 5000      Sample Output 3 817260251  Be sure to print the sum modulo 10^9+7 .

Solution (Python):
MAX = 10**9 + 7

def mod_inv(i, mod):
    return pow(i, mod-2, mod)

def combination(_N, _K):
    if _K > (_N/2):
        _K = _N - _K
    up = 1
    low = 1
    for i in range (_K):
        up = up * (_N - i) % MAX
        low = low * (i + 1) % MAX
    ret = up * mod_inv(low, MAX) % MAX
    return ret

def solve():
    N, M, K = [int(_) for _ in raw_input().split()]
    disSum = M*N* (M*(N+1)*(N-1) + N*(M+1)*(M-1)) // 6
    ret = disSum * combination(M*N-2, K-2) % MAX
    print ret

if __name__ == '__main__':
    solve()

Solution (Python):
MAX = 10**9 + 7

def mod_inv(i, mod):
    return pow(i, mod-2, mod)

def combination(_N, _K):
    if _K > (_N/2):
        _K = _N - _K
    up = 1
    low = 1
    for i in range (_K):
        up = up * (_N - i) % MAX
        low = low * (i + 1) % MAX
    ret = up * mod_inv(low, MAX) % MAX
    return ret

def solve():
    N, M, K = [int(_) for _ in raw_input().split()]
    disSum = (M*N* ((M*(N+1)*(N-1))%MAX + (N*(M+1)*(M-1)) % MAX)) % MAX
    disSum *= mod_inv(6, MAX) % MAX
    ret = disSum * combination(M*N-2, K-2) % MAX
    print ret

if __name__ == '__main__':
    solve()

Solution (Python):
import sys
def LI(): return [int(x) for x in sys.stdin.readline().split()]
def II(): return int(sys.stdin.readline())
def LS(): return sys.stdin.readline().split()
sys.setrecursionlimit(10**7)
INF = 10 ** 18
MOD = 10 ** 9 + 7
def LI_(): return [int(x) - 1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def SI(): return input()
YN = lambda b: print('YES') if b else print('NO')
yn = lambda b: print('Yes') if b else print('No')

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 = 2*10**5
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 = LI()

    rows = M*M*sum([d * (N-d) for d in range(1, N)])
    cols = N*N*sum([d * (M-d) for d in range(1, M)])
    ans = (rows + cols) * cmb(M*N-2, K-2, MOD) % MOD

    print(ans)

main()

Solution (Python):
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
    

N, M, K = map(lambda x: int(x), 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)

Solution (Python):
def mfactorial(n: int, mod: int) -> int:
  cache = mfactorial.__cache if hasattr(mfactorial, '__cache') else [1, 1]
  n_len = len(cache)
  n_new = n - n_len + 1
  if n_new > 0:
    cache = cache + ([0] * n_new)
    for i in range(n_len, n + 1):
      cache[i] = cache[i - 1] * i % mod
    mfactorial.__cache = cache

  return mfactorial.__cache[n]


def cmb(n, r, mod):
  result = mfactorial(n, mod)
  result *= pow(mfactorial(r, mod), mod - 2, mod) % mod
  result %= mod
  result *= pow(mfactorial(n - r, mod), mod - 2, mod) % mod
  result %= mod

  return result


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

  N, M, K = [int(i) for i in input().strip().split(' ')]
  c = cmb(N * M - 2, K - 2, MOD)

  X = 0
  Y = 0
  for i in range(M):
    X += (M - i) * i
  for i in range(N):
    Y += (N - i) * i

  total = (X * (N ** 2) + Y * (M ** 2)) % MOD * c % MOD
  print(total)


main()

Solution (Python):
from itertools import product, combinations
N,M,K=map(int, raw_input().split())


mod=10**9+7
fact, inv_fact = [1], [1]
fact_tmp = 1
for i in range(1, N*M+1):
    fact_tmp *= i
    fact_tmp %= mod
    fact.append(fact_tmp)
    inv_fact.append(pow(fact_tmp, mod-2, mod))
 
 
def ncr(n,r):
	if n < 0 or r < 0 or n < r:	return 0
	else:	return (fact[n] * inv_fact[r] * inv_fact[n-r]) %mod


torisuu=ncr(N*M-2,K-2)

ans=0
for i in range(N):
	for j in range(M):
		l_sum=(1+j)*j/2*N		
		r_sum=(1+(M-j-1))*(M-j-1)/2*N	
		u_sum=(1+i)*i/2*M		
		d_sum=(1+(N-i-1))*(N-i-1)/2*M	
		ans+=(l_sum+r_sum+u_sum+d_sum)*torisuu
		ans%=mod

print ( ans*pow(2,mod-2,mod) ) %mod

Solution (Python):
import sys
stdin = sys.stdin

sys.setrecursionlimit(10 ** 7)

def li(): return map(int, stdin.readline().split())
def li_(): return map(lambda x: int(x) - 1, stdin.readline().split())
def lf(): return map(float, stdin.readline().split())
def ls(): return stdin.readline().split()
def ns(): return stdin.readline().rstrip()
def lc(): return list(ns())
def ni(): return int(stdin.readline())
def nf(): return float(stdin.readline())


n, m, k = li()
MOD = 10**9 + 7

ans = 0



def inv_mod(n: int, mod: int) -> list:
    inv = [0, 1]
    for i in range(2, n + 1):
        inv.append(mod - ((mod // i) * inv[mod % i]) % mod)
    return inv



def fact(n: int, mod: int) -> list:
    fac = [1, 1]
    res = 1
    for i in range(2, n + 1):
        res = res * i % mod
        fac.append(res)
    return fac



def fact_inv(n: int, inv: list, mod: int) -> list:
    facInv = [1, 1]
    for i in range(2, n + 1):
        facInv.append(facInv[i - 1] * inv[i] % mod)
    return facInv



def nCr(n: int, r: int, mod: int, fac: list, facInv: list) -> int:
    if not (0 <= r and r <= n):
        return 0

    return ((fac[n] * facInv[r]) % mod) * facInv[n - r] % mod


def onedim(width: int, depth: int, mod: int):
    ret = 0
    for d in range(width):
        ret += depth**2 * (width - d) * d * nCr(m*n-2, k-2, mod, fac, facInv)
        ret %= mod

    return ret

fac = fact(n*m+1, MOD)
facInv = fact_inv(n*m, inv_mod(m*n+1, MOD), MOD)

print((onedim(n, m, MOD) + onedim(m, n, MOD)) % MOD)

Solution (Python):
N,M,K = map(int,input().split())

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

def inv(n):
    n %= mod
    return pow(n,mod-2,mod)

ans = M*(M**2-1)*N**2%mod+N*(N**2-1)*M**2%mod
ans *= inv(6)*fact[N*M-2]*inv(fact[K-2]*fact[N*M-K])
ans %= mod

print(ans)

Solution (Python):
N,M,K=map(int,input().split())
L=N*M+10
mod=10**9+7
table=[1]*(L+3)
t=1
for i in range(1,L+3):
    t*=i
    t%=mod
    table[i]=t
rtable=[1]*(L+3)
t=1
for i in range(1,L+3):
    t*=pow(i,mod-2,mod)
    t%=mod
    rtable[i]=t
ans=0
t=table[N*M-2]*rtable[K-2]*rtable[N*M-K]
for i in range(1,M):
    ans+=N*N*(M-i)*i
    ans%=mod
for i in range(1,N):
    ans+=M*M*(N-i)*i
    ans%=mod
print((ans*t)%mod)

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


MOD = 10 ** 9 + 7

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

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

inv_fact = [-1] * (N * M - 1)
inv_fact[-1] = pow(fact[-1], MOD - 2, MOD)
for i in range(N * M - 3, -1, -1):
    inv_fact[i] = inv_fact[i + 1] * (i + 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)

Solution (Python):
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)


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)

Solution (Python):
class Factorial:
    def __init__(self,n,mod):
        self.f=[1]
        for i in range(1,n+1):
            self.f.append(self.f[-1]*i%mod)
        self.i=[pow(self.f[-1],mod-2,mod)]
        for i in range(1,n+1)[::-1]:
            self.i.append(self.i[-1]*i%mod)
        self.i.reverse()
    def factorial(self,i):
        return self.f[i]
    def ifactorial(self,i):
        return self.i[i]
    def comb(self,n,k):
        return self.f[n]*self.i[n-k]*self.i[k]
n,m,k=map(int,input().split())
M=10**9+7
print((sum(m*m*i*(n-i)for i in range(n))+sum(n*n*i*(m-i)for i in range(m)))*Factorial(n*m,M).comb(n*m-2,k-2)%M)

Solution (Python):
class Factorial:
    def __init__(self,n,mod):
        self.f=[1]
        for i in range(1,n+1):
            self.f.append(self.f[-1]*i%mod)
        self.i=[pow(self.f[-1],mod-2,mod)]
        for i in range(1,n+1)[::-1]:
            self.i.append(self.i[-1]*i%mod)
        self.i.reverse()
    def factorial(self,i):
        return self.f[i]
    def ifactorial(self,i):
        return self.i[i]
    def comb(self,n,k):
        return self.f[n]*self.i[n-k]*self.i[k]
n,m,k=map(int,input().split())
M=10**9+7
print((sum(m*m*i*(n-i)for i in range(n))+sum(n*n*i*(m-i)for i in range(m)))*Factorial(n*m,M).comb(n*m-2,k-2)%M)

Solution (Python):
from collections import deque
from heapq import heappush,heappop
import re

def int_raw():
    return int(input())

def ss_raw():
    return input().split()

def ints_raw():
    return list(map(int, ss_raw()))

INF = 1<<29

DIV = 10**9+7

def mod_inv_prime(a,mod = DIV):
    return pow(a,mod-2,mod)

def mod_inv(a,b):
    r = a
    w = b
    u = 1
    v = 0
    while w!=0:
        t = r//w
        r -=t*w
        r,w = w,r
        u -= t*v
        u,v = v,u
    return (u%b+b)%b

def ncr(n,r, mod =DIV):
	r = min(r,n-r)
	ret = 1
	for i in range(1,r+1):
		ret = ret * (n-i+1) % mod
		ret = ret * mod_inv(i, mod) % mod
	return ret


def main():
    N, M,K = ints_raw()
    ans = 0
    for d in range(1,N):
        ans =(ans+d*(N-d)*M*M)%DIV
    for d in range(1,M):
        ans =(ans+d*(M-d)*N*N)%DIV
    return (ans*ncr(N*M-2,K-2))%DIV

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

Solution (Python):
n, m, k = [int(i) for i in input().split()]
p = 10 ** 9 + 7

S = 0

for d in range(n):
    S += m ** 2 * d * (n - d)
    S %= p

for d in range(m):
    S += n ** 2 * d * (m - d)
    S %= p

def fact(n, p=10**9 + 7):
    f = [1]
    for i in range(1, n+1):
        f.append(f[-1]*i%p)
    return f

def invfact(n, f, p=10**9 + 7):
    inv = [pow(f[n], p-2, p)]
    for i in range(n, 0, -1):
        inv.append(inv[-1]*i%p)
    return inv[::-1]

f = fact(n * m - 2)
invf = invfact(n * m - 2, f)

def comb(a, b):
    if a < b:
        return 0
    if a < 0 or b < 0:
        return 0
    return f[a] * invf[b] * invf[a-b] % p

S *= comb(n * m - 2, k - 2)

print(S % p)

Solution (Python):
H, W, K = map(int, input().split())
MOD = 10 ** 9 + 7

Hd = sum([d * (H - d) * (W ** 2) for d in range(H)])
Wd = sum([d * (W - d) * (H ** 2) for d in range(W)])



factorial = [1]
inverse = [1]
for i in range(1, H * W + 2):
    factorial.append(factorial[-1] * i % MOD)
    inverse.append(pow(factorial[-1], MOD - 2, MOD))


ans = (Hd + Wd) * factorial[H * W - 2] * inverse[K - 2] * inverse[H * W - K] % MOD
print(ans)

Solution (Python):

import sys

MOD = 1000000007  

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

N = 2*10**5
g1 = [1, 1]
g2 = [1, 1]
inv = [0, 1]

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

def solve(N: int, M: int, K: int):
    ans = 0
    for i in range(1, M):
        ans += i * N**2 * (M-i)
        ans %= MOD
    for i in range(1, N):
        ans += i * M**2 * (N-i)
        ans %= MOD
    ans *= cmb(N*M-2, K-2, MOD)
    ans %= MOD
    print(ans)
    return



def main():
    def iterate_tokens():
        for line in sys.stdin:
            for word in line.split():
                yield word
    tokens = iterate_tokens()
    N = int(next(tokens))  
    M = int(next(tokens))  
    K = int(next(tokens))  
    solve(N, M, K)

if __name__ == '__main__':
    main()

Solution (Python):
LARGE = 10**9 + 7


def p_inv(k, p):
    return pow(k, p-2, p)


def nck_mod(n, k, m):
    res = 1
    for i in range(n, n-k, -1):
        res = res * i % m
    for i in range(1, k+1):
        res = res * p_inv(i, m) % m
    return res


def solve(n, m, k):
    res_1 = 0
    r = n * (n - 1) // 2
    res_1 += r
    for i in range(1, n):
        r += i
        r -= n - i
        res_1 += r
    res_2 = 0
    r = m * (m - 1) // 2
    res_2 += r
    for i in range(1, m):
        r += i
        r -= m - i
        res_2 += r
    res = (res_1 * m * m + res_2 * n * n) % LARGE

    weight = nck_mod(n * m - 1, k - 1, LARGE) * p_inv(n * m - 1, LARGE) * (k - 1) * p_inv(2, LARGE) % LARGE

    return res * weight % LARGE


if __name__ == "__main__":
    N, M, K = map(int, input().split())
    print(solve(N, M, K))

Solution (Python):

import math


def calc_C(n, r):
    return (math.factorial(n) // (math.factorial(n - r) * math.factorial(r))) % (10**9 + 7)


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

ans = 0


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


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

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

Solution (Python):
n, m, k = map(int, input().split())
MOD = (10**9) + 7


S = 0

def stick(n):
    return n * (n+1) // 2

def calc(w, h):
    L = w * (w-1) // 2
    R = L + w * (h-1)
    return (L + R) * h // 2

for i in range(n):
    for j in range(m):
        S += calc(i+1, j+1)
        S += calc(i+1, m-j)
        S += calc(n-i, j+1)
        S += calc(n-i, m-j)
        S -= stick(i)
        S -= stick(j)
        S -= stick(n-i-1)
        S -= stick(m-j-1)
        S %= MOD
        
if S % 2 == 0:S //= 2
else: S = (S + MOD) // 2
S %= MOD
if S < 0:S += MOD

def fact(x):
    res = 1
    for i in range(1, x+1):
        res = res * i % MOD
    return res


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

def nCk(n, k):
    return fact(n) * rev(fact(k)) % MOD * rev(fact(n-k)) % MOD

print(S * nCk(n*m-2, k-2) % MOD)

Solution (Python):
class BigCombination(object):
    __slots__ = ["mod", "factorial", "inverse"]

    def __init__(self, mod: int = 10**9+7, max_n: int = 10**6):
        fac, inv = [1], []
        fac_append, inv_append = fac.append, inv.append

        for i in range(1, max_n+1):
            fac_append(fac[-1] * i % mod)

        inv_append(pow(fac[-1], mod-2, mod))

        for i in range(max_n, 0, -1):
            inv_append(inv[-1] * i % mod)

        self.mod, self.factorial, self.inverse = mod, fac, inv[::-1]

    def get_combination(self, n, r):
        return self.factorial[n] * self.inverse[r] * self.inverse[n-r] % self.mod

    def get_permutation(self, n, r):
        return self.factorial[n] * self.inverse[n-r] % self.mod



import math
Big = BigCombination()
s = 0
m = 10**9+7


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

res =0
for i in range(N-1):
  res += (i+1)*(N-i-1)*(M**2)
for i in range(M-1):
  res += (i+1)*(M-i-1)*(N**2)
if K>2:
  anp = Big.get_combination(N*M-2,K-2)
else:anp = 1
print((res*anp)%m)

Solution (Python):
class BigCombination(object):
    __slots__ = ["mod", "factorial", "inverse"]

    def __init__(self, mod: int = 10**9+7, max_n: int = 10**6):
        fac, inv = [1], []
        fac_append, inv_append = fac.append, inv.append

        for i in range(1, max_n+1):
            fac_append(fac[-1] * i % mod)

        inv_append(pow(fac[-1], mod-2, mod))

        for i in range(max_n, 0, -1):
            inv_append(inv[-1] * i % mod)

        self.mod, self.factorial, self.inverse = mod, fac, inv[::-1]

    def get_combination(self, n, r):
        return self.factorial[n] * self.inverse[r] * self.inverse[n-r] % self.mod

    def get_permutation(self, n, r):
        return self.factorial[n] * self.inverse[n-r] % self.mod



import math
Big = BigCombination()
s = 0
m = 10**9+7


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

res =0
for i in range(N-1):
  res += (i+1)*(N-i-1)*(M**2)
for i in range(M-1):
  res += (i+1)*(M-i-1)*(N**2)
if K>2:
  anp = Big.get_combination(N*M-2,K-2)
else:anp = 1
print((res*anp)%m)

Solution (Python):

"""
Created on Wed May 29 22:49:19 2019

@author: Yamazaki Kenichi
"""


"""
Created on Sat May 25 21:24:29 2019

@author: Yamazaki Kenichi
"""

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

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

g1 = [1, 1] 
g2 = [1, 1] 
inv = [0, 1] 

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

ans1 = 0
def c(N,M):
    res = 0
    for i in range(1,N):
        res += i*(N-i)*M*M
    return res
ans1 += c(N,M) + c(M,N)
ans1 %= mod

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

print(ans1)

Solution (Python):

"""
Created on Sat May 25 21:24:29 2019

@author: Yamazaki Kenichi
"""

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

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

g1 = [1, 1] 
g2 = [1, 1] 
inv = [0, 1] 

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

def a(N,M):
    return N*(N+1)//2*(M+1) + M*(M+1)//2*(N+1)
def b(N):
    return N*(N+1)//2

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

ans *= g2[2]

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

print(ans)

Solution (Python):
mod=10**9+7
def inv(x):
    if x==0:
        raise ValueError("Div 0")
    (nowx,lastx) = (0,1)
    (nowy,lasty) = (1,0)
    a=x;b=mod
    while b!=0:
        q=a//b
        (a,b)=(b,a%b)
        (nowx,lastx)=(lastx-q*nowx,nowx)
        (nowy,lasty)=(lasty-q*nowy,nowy)
    return (a*lastx)%mod
N,M,K=map(int,input().split())
X=sum([d*(N-d)*M*M for d in range(N)])
Y=sum([d*(M-d)*N*N for d in range(M)])
ans=X+Y
for i in range(1,K-1):
    ans*=(N*M-1-i)
    ans%=mod
    ans*=inv(i)
    ans%=mod
print(ans)

Solution (Python):
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=map(int,input().split())


x=0
for i in range(1,N//2+1):
   x+=2*i*(N-i)
if N%2==0: x-=N*N//4

y=0
for j in range(1,M//2+1):
    y+=2*j*(M-j)
if M%2==0: y-=M*M//4 

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

Solution (Python):
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=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)

Solution (Python):
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=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)

Solution (Python):
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=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)

Solution (Python):
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)


def modinv(a):
    m=MOD
    (inv, q, gcd_val) = egcd(a, m)
    return inv % m

MOD=10**9+7
N,M,K=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)

Solution (Python):
N, M, K = map(int, input().split())

mod = 10**9+7
frac = [1]*200100
for i in range(2,200100):
    frac[i] = i * frac[i-1]%mod
fraci = [None]*200100
fraci[-1] = pow(frac[-1], mod -2, mod)
for i in range(-2, -200100 - 1, -1):
    fraci[i] = fraci[i+1] * (200100 + i + 1) % mod
def comb(a, b):
    if not a >= b >= 0:
        return 0
    return frac[a]*fraci[b]*fraci[a-b]%mod


mod = 10**9+7

t = N**2*(M+1)*M*(M-1)//6 + M**2*(N+1)*N*(N-1)//6
print(t*comb(N*M-2,K-2)%mod)

Solution (Python):
n,m,k=map(int,input().split())
M=10**9+7
def combi(x,y):
  c=1
  m=min(x-y,y)
  for iii in range(y):
    c=c*(x-iii)*pow(y-iii,M-2,M)%M
  return c

rem_combi=combi(n*m-2,k-2)

def calc(num):
  out=0
  for iii in range(num):
    out+=iii*(num-iii)%M
  return out

x_cost=m*m*calc(n)%M
y_cost=n*n*calc(m)%M
ans=rem_combi*(x_cost+y_cost)%M
print(int(ans))

Solution (Python):
mod = 10**9+7
n,m,k = map(int,input().split())



x = k
fact = [1] * (x + 1)
for i in range(1, x + 1):
    fact[i] = fact[i - 1] * i % mod

def mod_comb_k(n,k,mod):
    ret = 1
    if n - k < k:
        k = n - k
    for i in range(n, n - k, -1):
        ret = ret * i % mod
    return ret % mod * pow(fact[k], mod - 2, mod) % mod

res = 0
for d in range(1,n):
    res += m*m*(n-d)*d
    res = res%mod

n,m = m,n
for d in range(1,n):
    res += m*m*(n-d)*d
    res = res%mod

res *= mod_comb_k(n*m-2,k-2,mod)
'''
for d in range(1,m):
    res += n*n*(m-d)*mod_comb_k(n*m-2,k-2,mod)%mod
    res = res%mod
'''
print(res%mod)

Solution (Python):
import sys, math, collections, heapq, itertools
F = sys.stdin
def single_input(): return F.readline().strip("\n")
def line_input(): return F.readline().strip("\n").split()
mod = 7 + 10 ** 9
factorial = [1] * (2 * 10 ** 5 + 1)
for i in range(1, 2 * 10 ** 5 + 1):
    factorial[i] = (factorial[i-1] * i) % mod
revfactorial = [1] * (2 * 10 ** 5 + 1)
revfactorial[-1] = pow(factorial[-1], mod - 2, mod)
for i in reversed(range(1, 2 * 10**5)):
    revfactorial[i] = (revfactorial[i+1] * (i+1)) % mod

def comb(n, r):
    return (factorial[n] * revfactorial[n-r] * revfactorial[r]) % mod

def solve():
    N, M, K = map(int, line_input())
    xsum, ysum = 0, 0
    for i in range(1, N): 
        xsum += ((N - i) * M * M * i) % mod
        xsum %= mod
    for j in range(1, M):
        ysum += ((M - j) * N * N * j) % mod
        ysum %= mod
    print(((xsum + ysum) * comb(N*M-2, K-2)) % mod)
    return 0
  
if __name__ == "__main__":
    solve()

Solution (Python):
n, m, k = map(int, input().split())

MOD = 10 ** 9 + 7

class ModCmb:
    def __init__(self, size):
        self.inv = [1] * (size + 1)
        self.fact = [1] * (size + 1)
        temp_inv = [1] * (size + 1)
        for i in range(2, size + 1):
            temp_inv[i] = ( -(MOD // i) * temp_inv[MOD%i] ) % MOD
        for i in range(2, size + 1):
            self.fact[i] = self.fact[i-1] * i % MOD
            self.inv[i] = self.inv[i-1] * temp_inv[i] % MOD

    def cmb(self, n, r):
        if n <= 1 or n == r:
            return 1
        elif r == 1:
            return n
        else:
            return  (self.fact[n] * self.inv[r] % MOD ) * self.inv[n-r] % MOD

mc = ModCmb(n*m)

temp = mc.cmb(n*m-2, k-2)

ans = 0

for h in range(n):
    for w in range(m):
        if h != 0 and w != 0:
            ans += (h + w) * temp * (n - h) * (m - w) * 2
        else:
            ans += (h + w) * temp * (n - h) * (m - w)
        ans %= MOD

print(ans)

Solution (Python):
n, m, k = list(map(int, input().split()))

MAX = 2*10**5+1
MOD = 10**9+7
fact = [0 for i in range(MAX)]
inv = [0 for i in range(MAX)]
invfact = [0 for i in range(MAX)]

def comb_build(n):
  fact[0] = inv[0] = invfact[0] = 1
  fact[1] = inv[1] = invfact[1] = 1
  for i in range(2, n):
    fact[i] = fact[i-1]*i%MOD
    inv[i] = inv[MOD%i]*(MOD-MOD//i)%MOD
    invfact[i] = invfact[i-1]*inv[i]%MOD

def nCk(n, k):
  if n<k or n<0 or k<0:
    return 0
  return (((fact[n]*invfact[k])%MOD)*invfact[n-k])%MOD

comb_build(MAX)
ans = 0
t = nCk(n*m-2, k-2)
for i in range(1, n+1):
  ans += (((i*(n-i)*m*m)%MOD)*t)%MOD
for i in range(1, m+1):
  ans += (((i*(m-i)*n*n)%MOD)*t)%MOD
print(ans%MOD)

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

tn = 0
tm = 0

for n in range(N):
    tn = (tn + n*(n+1)//2) % mod
tn = tn*M*M % mod

for m in range(M):
    tm = (tm + m*(m+1)//2) % mod
tm = tm*N*N % mod

class combinationmod():
    def __init__(self, MAX, mod):
        self.mod = mod
        self.fac = [0]*(MAX+1)
        self.finv = [0]*(MAX+1)
        self.inv = [0]*(MAX+1)
        self.fac[0] = self.fac[1] = 1
        self.finv[0] = self.finv[1] = 1
        self.inv[1] = 1
        for i in range(2, MAX+1):
            self.fac[i] = self.fac[i-1]*i % mod
            self.inv[i] = (-self.inv[mod % i]*(mod//i)) % mod
            self.finv[i] = self.finv[i-1] * self.inv[i] % mod
    def nCk(self,n,k):
        if n<k or n<0 or k<0:
            return 0
        return self.fac[n]*(self.finv[k]*self.finv[n-k]%self.mod)%self.mod
 
c = combinationmod(200000,int(1e9+7))

ans = (tn + tm)*c.nCk(N*M-2, K-2) %mod
print(ans)

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

def makefactlist(n, mod=MOD):
    ret = [1]
    for i in range(n):
        ret += [ (ret[-1]*(i+1))%MOD ]
    return ret

factlist = makefactlist(N*M, mod=MOD)
comb = factlist[N*M-2]*pow(factlist[N*M-K], MOD-2, MOD)*pow(factlist[K-2], MOD-2, MOD)
comb = comb % MOD

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

Solution (Python):
N, M, K = 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)

Solution (Python):
import sys
input = sys.stdin.readline
inf = float('inf')
mod = 10**9+7


def INT_(n): return int(n)-1


def MI(): return map(int, input().split())


def MF(): return map(float, input().split())


def MI_(): return map(INT_, input().split())


def LI(): return list(MI())


def LI_(): return [int(x) - 1 for x in input().split()]


def LF(): return list(MF())


def LIN(n: int): return [input() for _ in range(n)]


def LLIN(n: int): return [LI() for _ in range(n)]


def LLIN_(n: int): return [LI_() for _ in range(n)]


def LLI(): return [list(map(int, l.split())) for l in input()]


def I(): return int(input())


def F(): return float(input())


def ST(): return input().replace('\n', '')


def main():

    class BigCombination(object):

        def __init__(self, mod: int = 10**9+7, max_n: int = 10**6):
            fac, inv = [1], []
            fac_append, inv_append = fac.append, inv.append

            for i in range(1, max_n+1):
                fac_append(fac[-1] * i % mod)

                inv_append(pow(fac[-1], mod-2, mod))

            for i in range(max_n, 0, -1):
                inv_append(inv[-1] * i % mod)

            self.mod, self.factorial, self.inverse = mod, fac, inv[::-1]

        def get_combination(self, n, r):
            return self.factorial[n] * self.inverse[r] * self.inverse[n-r] % self.mod

        def get_permutation(self, n, r):
            return self.factorial[n] * self.inverse[n-r] % self.mod
    N, M, K = MI()
    comb = BigCombination(max_n=N*M)
    get_comb = comb.get_combination

    ans = 0
    for distansssss in range(N):
        ans += distansssss * (M ** 2) * (N - distansssss)
        ans %= mod
    for distansssss in range(M):
        ans += distansssss * (N ** 2) * (M - distansssss)
        ans %= mod
    ans *= get_comb(N*M-2, K-2)
    print(ans % mod)


if __name__ == '__main__':
    main()

Solution (Python):
N, M, K = map(int, input().split())

mod = 10**9 + 7
table = [1]
for i in range(N*M):
    table.append(table[-1] * (i+1) % mod)


def inv(N):
    return pow(N, mod-2, mod)  


def cmb(n, r):
    return table[n] * inv(table[n-r]*table[r]) % mod


print((M ** 2 * (N ** 3 - N) + N ** 2 * (M ** 3 - M)) // 6 * cmb(M * N - 2, K - 2) % mod)