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):
n,m,k=map(int,input().split())
mod = 10**9+7
fac = [1]
caf = [1]
def pow(n,p):
    res=1
    while p >0:
        if p%2==0:
            n = n**2 % mod
            p //= 2
        else:
            res = res * n % mod
            p -= 1
    return res % mod

for i in range(1,n*m):
    fac.append((fac[i-1]*i)%mod)
    caf.append(pow(fac[i],mod-2))

def ncr(n,r):
    if n < r or r < 0: return 0
    return fac[n]*caf[r]*caf[n-r]

div2 = caf[2]

sums = [0]
for i in range((max(n,m))):
    sums.append((sums[-1]+i+1) % mod)

res1 = 0
for i in range(m):
    res1 += sums[i]+sums[m-1-i]
    res1 %= mod
res1 *= (n*n)%mod
res1 %= mod

res2 = 0
for i in range(n):
    res2 += sums[i]+sums[n-1-i]
    res2 %= mod
res2 *= (m*m)%mod
res2 %= mod

res = ((res1+res2)*div2)%mod
kc2 = ncr(n*m-2,k-2)
print((res*kc2)%mod)

Solution (Python):



import sys
import math
from bisect import bisect_right as br
from bisect import bisect_left as bl
sys.setrecursionlimit(1000000)
from heapq import heappush, heappop,heappushpop
from collections import defaultdict
from itertools import accumulate
from collections import Counter
from collections import deque
from operator import itemgetter
from itertools import permutations
mod = 10**9 + 7
inf = float('inf')
def I(): return int(sys.stdin.readline())
def LI(): return list(map(int,sys.stdin.readline().split()))

def C(n,r):
    return fact[n]*pow(fact[n-r],mod-2,mod)*pow(fact[r],mod-2,mod) % mod

n,m,k = LI()

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

ans = 0
"""xの差"""
for d in range(1,m):
    ans += d*(m-d)*n*n
    ans %= mod

"""yの差"""
for d in range(1,n):
    ans += d*(n-d)*m*m
    ans %= mod

ans *= C(n*m-2,k-2)
ans %= mod
print(ans)

Solution (Python):
n,m,k = map(int, input().split())
n_ = n*m
mod = 10**9+7
fun = [1]*(n_+1)
for i in range(1,n_+1):
    fun[i] = fun[i-1]*i%mod
rev = [1]*(n_+1)
rev[n_] = pow(fun[n_],mod-2,mod)
for i in range(n_-1,0,-1):
    rev[i] = rev[i+1]*(i+1)%mod

def nCr(n,r):
    if r > n:
        return 0
    return fun[n]*rev[r]%mod*rev[n-r]%mod
  
x = sum(i*(n-i) for i in range(n))%mod
y = sum(i*(m-i) for i in range(m))%mod
ans = (x*m**2 + y*n**2)%mod*nCr(n*m-2,k-2)%mod
print(ans)

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

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

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

def inv_mod(a, p=MOD):
    def inv_mod_sub(a, p):
        if a == 1:
            return 1, 0
        else:
            d, r = p//a, p%a
            x, y = inv_mod_sub(r, a)
            return y-d*x, x
    if p < 0: p = -p
    a %= p
    return inv_mod_sub(a,p)[0] % p

def comb_mod(n, k):
    if k < 0 or k > n:
        return 0
    else:
        return f_mod[n]*f_mod_inverse[k]*f_mod_inverse[n-k] % MOD
f_mod=[1]*(num+1)
f_mod_inverse=[1]*(num+1)
for i in range(1,num+1):
  f_mod[i]=f_mod[i-1]*i%MOD
  f_mod_inverse[i]=f_mod_inverse[i-1]*inv_mod(i)%MOD
N=0;M=0
for i in range(1,n):
  N+=i*(n-i)
  N%=MOD
for i in range(1,m):
  M+=i*(m-i)
  M%=MOD
print((N*(m**2)%MOD+M*(n**2)%MOD)*comb_mod(n*m-2,k-2)%MOD)

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

print((n*m*(m-1)*(m+1)*k*(k-1)*cmb(n*m,k)//(6*m*(n*m-1))+m*n*(n-1)*(n+1)*k*(k-1)*cmb(n*m,k)//(6*n*(n*m-1)))%(10**9+7))

Solution (Python):

n,m,k = list(map(int,input().split()))
mod = 10 ** 9 + 7
def inv(x):
  y = 1
  while x != 1:
    y *= mod//x+1
    y %= mod
    x -= mod%x
  return y

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

zsum = [0]
for i in range(1,max(n,m)):
  zsum.append(zsum[-1]+i)

ans1 = 0
for i in range(m):
  ans1 += zsum[i]
  ans1 += zsum[m-1-i]
  ans1 %= mod
ans1 *= n ** 2
ans1 %= mod
ans2 = 0
for j in range(n):
  ans2 += zsum[j]
  ans2 += zsum[n-1-j]
ans2 *= m ** 2
ans2 %= mod
ans *= ans1 + ans2
ans *= inv(2)
ans %= mod


print(ans)

Solution (Python):

import sys
input = sys.stdin.buffer.readline

con = 10 ** 9 + 7

from collections import defaultdict

class Combination(object):
	def __init__(self, N, con):
		self.fac = [0] * (N + 1)
		self.inv = [0] * (N + 1)
		self.finv = [0] * (N + 1)
		self.fac[0], self.fac[1] = 1, 1
		self.inv[1] = 1
		self.finv[0], self.finv[1] = 1, 1

		
		for i in range(2, N + 1):
			self.fac[i] = self.fac[i - 1] * i % con
			self.inv[i] = self.inv[con % i] * (con - (con // i)) % con
			self.finv[i] = self.finv[i - 1] * self.inv[i] % con

	def com(self, N, k):
		return (self.fac[N] * self.finv[k] * self.finv[N - k]) % con


def getlist():
	return list(map(int, input().split()))


def main():
	N, M, K = getlist()
	C = Combination(N * M - 2, con)
	c = C.com(N * M - 2, K - 2)
	s = (N * (N + 1) * (N - 1)) // 6 * (M ** 2)
	t = (M * (M + 1) * (M - 1)) // 6 * (N ** 2)
	ans = c * (s + t)
	print(ans % con)

if __name__ == '__main__':
	main()

Solution (Python):
from math import factorial as f
n,m,k=map(int, input().split())
print((f(n*m-2)//(f(k-2)*f(n*m-k))*(m*m*(n**3-n)+n*n*(m**3-m))//6)%(10**9+7))

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

L = int(0)



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





C = combination(N*M-2,mod)
cou = C.comb(N*M-2,K-2,mod)*L%mod

print(cou)

Solution (Python):
def comb(n, k, p):
    """power_funcを用いて(nCk) mod p を求める"""
    from math import factorial
    if n < 0 or k < 0 or n < k: return 0
    if n == 0 or k == 0: return 1
    a = factorial(n) % p
    b = factorial(k) % p
    c = factorial(n - k) % p
    return (a * power_func(b, p - 2, p) * power_func(c, p - 2, p)) % p


def power_func(a, b, p):
    """a^b mod p を求める"""
    if b == 0: return 1
    if b % 2 == 0:
        d = power_func(a, b // 2, p)
        return d * d % p
    if b % 2 == 1:
        return (a * power_func(a, b - 1, p)) % p


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

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

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

ans = (((sumx + sumy) % mod) * c) % mod
print(ans)

Solution (Python):
n,m,k=map(int,input().split())
mod=10**9+7
A=n*m-2
from math import factorial
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
v=Comb(A,k-2)
t=0
y=0
for i in range(1,n):
  t+=i*(n-i)*m*m%mod
for i in range(1,m):
  y+=i*(m-i)*n*n%mod

ans=v*(t+y)
print(ans%mod)

Solution (Python):
n,m,k=map(int,input().split())
mod=10**9+7
A=n*m-2
from math import factorial
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
v=Comb(A,k-2)
t=0
y=0
for i in range(1,n):
  t+=i*(n-i)*m*m%mod
for i in range(1,m):
  y+=i*(m-i)*n*n%mod

ans=v*(t+y)
print(ans%mod)

Solution (Python):

mod=10**9+7
MAX=2*10**5+100

g1=[1,1]
g2=[1,1]
for i in range(2,MAX+1):
    num_1=g1[-1]*i%mod
    g1.append(num_1)
    g2.append(pow(num_1,mod-2,mod))
    
def cmb(n,r,MOD):
    return g1[n]*g2[r]*g2[n-r]%MOD

n,m,k=map(int,input().split())
s=0
t=0
for i in range(1,n):
    s+=i*(n-i)
for i in range(1,m):
    t+=i*(m-i)

ans=(n*n*t+m*m*s)%mod
ans=(ans*cmb(n*m-2,k-2,mod))%mod
print(ans)

Solution (Python):
import sys
sys.setrecursionlimit(2147483647)
INF=float("inf")
MOD=10**9+7
input=lambda :sys.stdin.readline().rstrip()
def modfact(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
    return fact,invfact

def resolve():
    n,m,k=map(int,input().split())
    fact,invfact=modfact(n*m-2)
    ans=1
    ans*=fact[m*n-2]*invfact[k-2]*invfact[m*n-k]
    ans*=((n-1)*n*(n+1)*m*m+(m-1)*m*(m+1)*n*n)//6
    print(ans%MOD)
resolve()

Solution (Python):
import sys
sys.setrecursionlimit(2147483647)
INF=float("inf")
MOD=10**9+7
input=lambda :sys.stdin.readline().rstrip()
def modfact(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
    return fact,invfact

def resolve():
    from math import factorial as fact
    n,m,k=map(int,input().split())
    fact,invfact=modfact(n*m)
    ans=1
    ans*=fact[m*n-2]*invfact[k-2]*invfact[m*n-k]
    ans*=((n-1)*n*(n+1)*m*m+(m-1)*m*(m+1)*n*n)//6
    print(ans%MOD)
resolve()

Solution (Python):
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,m,k=map(int,input().split())
mod=10**9+7
 
g1=[1,1]
g2=[1,1]
inverse=[0,1]
 
for i in range(2,n*m):
  g1.append((g1[-1]*i)%mod)
  inverse.append((-inverse[mod%i]*(mod//i))%mod)
  g2.append((g2[-1]*inverse[-1])%mod)

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

Solution (Python):
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,M,k=map(int,input().split())
m=10**9+7
 
g1=[1,1]
g2=[1,1]
inverse=[0,1]
 
for i in range(2,N*M):
  g1.append((g1[-1]*i)%m)
  inverse.append((-inverse[m%i]*(m//i))%m)
  g2.append((g2[-1]*inverse[-1])%m)

c=cmb(N*M-2,k-2,m)
ans=0
for i in range(N):
  ans+=((N-i)*M*M*i*c)%m
  ans%=m
for i in range(M):
  ans+=((M-i)*N*N*i*c)%m
  ans%=m
print(ans%m)

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


nmax = 2 * 10 ** 5 + 1  
fac = [0] * (nmax + 1)
finv = [0] * (nmax + 1)
inv = [0] * (nmax + 1)
fac[0] = 1
fac[1] = 1
finv[0] = 1
finv[1] = 1
inv[1] = 1
for i in range(2, nmax):
    fac[i] = fac[i - 1] * i % mod
    inv[i] = mod - inv[mod % i] * (mod // i) % mod
    finv[i] = finv[i - 1] * inv[i] % mod


def comb(n, r):
    if n < r:
        return 0
    else:
        return (fac[n] * ((finv[r] * finv[n - r]) % mod)) % mod


def s(k):
    return ((k * (k+1) * (2 * k + 1))%mod * inv[6] + (k * (k+1) * inv[2])%mod)%mod


C = comb(N*M-2, K-2)
print(((C * ((N ** 2 * s(M-1))%mod + (M ** 2 * s(N-1))%mod)%mod)%mod) * inv[2]%mod)

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

def cmb(n, r, mod=mod):
    if ( r<0 or r>n ):
        return 0
    r = min(r, n-r)
    return g1[n] * g2[r] * g2[n-r] % mod

NN = 10**6 
g1 = [1, 1] 
g2 = [1, 1] 
inverse = [0, 1] 

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

ans = 0
for n in range(1, N+1):
    le = n*(n-1)//2 + (N-n)*(N-n+1)//2
    le *= M**2
    ans += le

for m in range(1, M+1):
    le = m*(m-1)//2 + (M-m)*(M-m+1)//2
    le *= N**2
    ans += le

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

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

L = int(0)



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)






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

ans = cou *L%mod

print(ans)

Solution (Python):
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

import numpy as np



MOD = 10**9 + 7

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

def cumprod(arr,MOD):
    L = len(arr); Lsq = int(L**.5+1)
    arr = np.resize(arr,Lsq**2).reshape(Lsq,Lsq)
    for n in range(1,Lsq):
        arr[:,n] *= arr[:,n-1]; arr[:,n] %= MOD
    for n in range(1,Lsq):
        arr[n] *= arr[n-1,-1]; arr[n] %= MOD
    return arr.ravel()[:L]

def make_fact(U,MOD):
    x = np.arange(U,dtype=np.int64); x[0] = 1
    fact = cumprod(x,MOD)
    x = np.arange(U,0,-1,dtype=np.int64); x[0] = pow(int(fact[-1]),MOD-2,MOD)
    fact_inv = cumprod(x,MOD)[::-1]
    return fact,fact_inv

U = 10**6
fact, fact_inv = make_fact(U,MOD)

coef = fact[N*M-2] * fact_inv[K-2] % MOD * fact_inv[N*M-K] % MOD

x = (N-1)*N*(N+1)*M*M//6
y = (M-1)*M*(M+1)*N*N//6
answer = coef * ((x+y) % MOD) % MOD
print(answer)

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

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

def mod_combination(n, k):
    def extended_gcd(a, b):
        if b == 0:
            return (a, 1, 0)
        d, x, y = extended_gcd(b, a % b)
        return (d, y, x - (a // b) * y)
    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 p * (extended_gcd(q, MOD)[1] % MOD) % MOD

def mod_pow(a, b):
    return pow(a, b, MOD)

def mod_inv(a):
    return pow(a, MOD - 2, MOD)


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

X = mod_pow(M, 2) * (mod_pow(N, 3) - N) % MOD
Y = mod_pow(N, 2) * (mod_pow(M, 3) - M) % MOD

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

ans = (X + Y) * mod_inv(6) % MOD * comb % MOD

print(ans)

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

def mod_combination(n, k):
    def extended_gcd(a, b):
        if b == 0:
            return (a, 1, 0)
        d, x, y = extended_gcd(b, a % b)
        return (d, y, x - (a // b) * y)
    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 p * (extended_gcd(q, MOD)[1] % MOD) % MOD

def mod_pow(a, b):
    return pow(a, b, MOD)


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

X = mod_pow(M, 2) * (mod_pow(N, 3) - N) % MOD
Y = mod_pow(N, 2) * (mod_pow(M, 3) - M) % MOD

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

ans = (X + Y) * mod_pow(6, MOD - 2) * comb % MOD

print(ans)

Solution (Python):
n, m, k = 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)

Solution (Python):
import sys
import math
from collections import defaultdict

sys.setrecursionlimit(10**7)
def input():
    return sys.stdin.readline()[:-1]

mod = 10**9 + 7

def I(): return int(input())
def II(): return map(int, input().split())
def III(): return list(map(int, input().split()))
def Line(N,num):
    if N<=0:
        return [[]]*num
    elif num==1:
        return [I() for _ in range(N)]
    else:
        read_all = [tuple(II()) for _ in range(N)]
        return map(list, zip(*read_all))




def fact(a,M=mod):
    ans = 1
    for i in range(2,a+1):
        ans = ans*i
        ans = ans%M
    return ans


def extgcd(a,b):
    r = [1,0,a]
    w = [0,1,b]
    while w[2]!=1:
        q = r[2]//w[2]
        r2 = w
        w2 = [r[0]-q*w[0],r[1]-q*w[1],r[2]-q*w[2]]
        r = r2
        w = w2
    return [w[0],w[1]]



def mod_inv(a,M=mod):
    x = extgcd(a,M)[0]
    return (M+x%M)%M

N,M,K = II()
ans = 0
for i in range(N):
    x = i*(i+1)//2 + (N-1-i)*(N-i)//2
    for j in range(M):
        ans += M*x + N*(j*(j+1)+(M-1-j)*(M-j))//2
        ans %= mod

ans = ans*fact(N*M-2)*mod_inv(fact(K-2))*mod_inv(fact(N*M-K))*mod_inv(2)
ans %= mod
print(ans)

Solution (Python):
import sys
mod = pow(10, 9) + 7
def C(n, r):
    if ( r<0 or r>n ):
        return 0
    r = min(r, n-r)
    return g1[n] * g2[r] * g2[n-r] % mod
NNN = (10**6)
g1 = [1, 1]
g2 = [1, 1]
inverse = [0, 1]
for i in range(2, NNN+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=map(int,input().split())

print(C(n*m,k)*k*(k-1)%mod*(n+m)*inverse[6]%mod)

Solution (Python):
mod = 10**9+7
inv_2 = pow(2, mod-2, mod)
N, M, K = map(int, input().split())
def make_fact(n):
    fact = [1]*(n+1)
    for i in range(1, n+1):
        fact[i] = fact[i-1]*i%mod
    return fact

def make_fact_inv(n):
    fact_inv = [1]*(n+1)
    fact_inv[n] = pow(fact[n], mod-2, mod)
    for i in range(n, 0, -1):
        fact_inv[i-1] = fact_inv[i]*i%mod
    return fact_inv

def comb(n, k):
    return fact[n]*fact_inv[k]*fact_inv[n-k]%mod
fact = make_fact(N*M)
fact_inv = make_fact_inv(N*M)
Sum = [0]*(max(N, M)+1)
for i in range(1, len(Sum)):
    Sum[i] = (Sum[i-1]+i)%mod
ans = 0

for x in range(1, N+1):
    for y in range(1, M+1):
        ans+=(Sum[x-1]+Sum[N-x])*M+(Sum[y-1]+Sum[M-y])*N
        ans%=mod
print(ans*comb(N*M-2, K-2)*inv_2%mod)

Solution (Python):

h,w,k = map(int,input().split())
mod = 10 ** 9 + 7
def inv(x):
  y = 1
  while x != 1:
    y *= mod//x + 1
    y %= mod
    x = x - mod%x
  return y
others = 1
for i in range(k-2):
  others *= h*w-2-i
  others *= inv(i+1)
  others %= mod
hh = 0
for i in range(1,h):
  hh += i * (h-i)
  hh %= mod
ww = 0
for i in range(1,w):
  ww += i * (w-i)
  ww %= mod
ans = others * ( h * h * ww + w * w * hh )
print(ans % mod)

Solution (Python):
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import math
import bisect
import random
from itertools import permutations, accumulate, combinations
import sys
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor
from functools import reduce



INF = float('inf')
def LI(): return list(map(int, sys.stdin.readline().split()))
def I(): return int(sys.stdin.readline())
def LS(): return sys.stdin.readline().split()
def S(): return sys.stdin.readline().strip()
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]
mod = 1000000007


def mod_combination(n, r, mod=10 ** 9 + 7):
    fact = 1
    for i in range(n - r + 1, n + 1):
        fact = fact * i % mod

    divisor = 1
    for j in range(1, r + 1):
        divisor = divisor * j % mod

    return fact * pow(divisor, mod - 2, mod) % mod


def main():
    mod = 10**9+7
    n, m, k = LI()
    
    combination = mod_combination(n*m-2, k-2)
    x=y=0
    for i in range(1, n):
        x += i * (n-i)
    for j in range(1, m):
        y += j * (m-j)
    return (((m ** 2) % mod * x % mod) + ((n ** 2) % mod * y % mod)) * combination % mod


print(main())

Solution (Python):
from math import factorial as f
n,m,k=map(int,input().strip().split())
print((f(n*m-2)//(f(k-2)*f(n*m-k))*(m*m*(n**3-n)+n*n*(m**3-m))//6)%(10**9+7))

Solution (Python):
mod = 1000000007
import sys
input = sys.stdin.readline
N, M, K = 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)

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

Solution (Python):
def mulmod(x, y, mod=10**9 + 7):
    return x * y % mod


def comb(n, k, mod=10**9 + 7):
    from functools import reduce
    k = min(k, n - k)
    x = reduce(mulmod, range(n - k + 1, n + 1), 1)
    y = reduce(mulmod, range(1, k + 1), 1)
    return x * pow(y, mod - 2, mod) % mod


def solve(string):
    n, m, k = map(int, string.split())
    mod = 10**9 + 7
    base = n * m * (n + m) * (n * m - 1) // 6 % mod
    return str(base * comb(n * m - 2, k - 2) % mod)


if __name__ == '__main__':
    print(solve(input()))

Solution (Python):
def comb(n, k, mod=10**9 + 7):
    from functools import reduce
    from operator import mul
    k = min(k, n - k)

    def calc(x, y):
        return x * y % mod

    x = reduce(calc, range(n - k + 1, n + 1), 1) % mod
    y = reduce(calc, range(1, k + 1), 1) % mod
    return x * pow(y, mod - 2, mod) % mod


def solve(string):
    n, m, k = map(int, string.split())
    mod = 10**9 + 7
    base = n * m * (n + m) * (n * m - 1) // 6 % mod
    return str(base * comb(n * m - 2, k - 2) % mod)


if __name__ == '__main__':
    print(solve(input()))

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


def solve(string):
    n, m, k = map(int, string.split())
    base = n * m * (n + m) * (n * m - 1) // 6 % mod
    return str(base * comb(n * m - 2, k - 2) % mod)


if __name__ == '__main__':
    print(solve(input()))

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


def solve(string):
    n, m, k = map(int, string.split())
    base = n * m * ((n + 1) * (m + 1) * (n + m - 2) - n**2 - m**2 + 2) // 6 % mod
    return str(base * comb(n * m - 2, k - 2) % mod)


if __name__ == '__main__':
    print(solve(input()))

Solution (Python):

import sys

MOD = 1000000007  

def mfact():
    f = [1]
    for i in range(1, 200007):
        f.append(f[-1] * i % MOD)
    return f

def inv(x):
    res = 1
    k = MOD - 2
    y = x
    while (k):
        if(k % 2 == 1):
            res = (res * y) % MOD
        y = (y ** 2) % MOD
        k = k // 2
    return res

def comb(n, k, f):
    a = f[n]
    b = f[n-k]
    c = f[k]
    bc = (b * c) % MOD
    return (a * inv(bc)) % MOD

def solve(N: int, M: int, K: int):
    ans = 0
    for i in range(N):
        for j in range(M):
            a = j * (j+1) // 2
            b = (M-1-j) * (M-1-j+1) // 2
            ans += (a+b) * i + M * (i*(i+1)//2) + a
            ans %= MOD
    f = mfact()
    ans *= comb(N*M-2, K-2, f)
    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):
import sys,collections as cl,bisect as bs
sys.setrecursionlimit(100000)
mod = 10**9+7
Max = sys.maxsize
def l(): 
    return list(map(int,input().split()))
def m(): 
    return map(int,input().split())
def onem(): 
    return int(input())
def s(x): 
    a = []
    aa = x[0]
    su = 1
    for i in range(len(x)-1):
        if aa != x[i+1]:
            a.append([aa,su])
            aa = x[i+1]
            su = 1
        else:
            su += 1
    a.append([aa,su])
    return a
def jo(x): 
    return " ".join(map(str,x))
def max2(x): 
    return max(map(max,x))
mod = 10**9+7 
N = 2*(10**5)
def cmb(n, r):
    if ( r<0 or r>n ):
        return 0
    r = min(r, n-r)
    return g1[n] * g2[r] * g2[n-r] % mod
def p(n,r):
    if ( r<0 or r>n ):
        return 0
    return g1[n] * g2[n-r] % mod

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 )



n,m,k = m()

sx = 0
sy = 0

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

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

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

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

mod=10**9+7
a=cmb(n*m-2,k-2)%mod
res=0
for d in range(1,n):
	res+=(d*(n-d)*m**2*a)%mod
	res%=mod
for d in range(1,m):
	res+=(d*(m-d)*n**2*a)%mod
	res%=mod
print(res)

Solution (Python):
from math import factorial as f
N, M, K = map(int, input().split())
mod = 10 ** 9 + 7

def combi(n, r):
  return f(n) // (f(n - r) * f(r)) % mod

res = 0
c = combi(N * M - 2, K - 2)
for i in range(1, N):
  res += ((N - i + 1) * (N - i) // 2) * c * M * M % mod
  res %= mod

for i in range(1, M):
  res += ((M - i + 1) * (M - i) // 2) * c * N * N % mod
  res %= mod

print(res % mod)