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):
class FactMod:

    def __init__(self, n, mod):
        self.mod = mod
        self.f = [1] * (n + 1)
        for i in range(1, n + 1):
            self.f[i] = self.f[i - 1] * i % mod

        self.inv = [pow(self.f[-1], mod - 2, mod)]
        for i in reversed(range(1, n + 1)):
            self.inv.append(self.inv[-1] * i % mod)
        self.inv.reverse()

    def fact(self, n):
        """
        :return: n!
        """
        return self.f[n]

    def comb(self, n, r):
        """
        :return: nCr
        """
        return self.f[n] * self.inv[n - r] * self.inv[r] % self.mod


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

NM = N * M

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

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

factmod = FactMod(NM - 2, mod)
ans *= factmod.comb(NM - 2, K - 2)
ans %= mod
print(ans)

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

def mod_inv(a, m=MOD):
    def mod_inv_sub(a, m):
        if a == 1:
            return 1, 0
        else:
            d, r = m//a, m % a
            x, y = mod_inv_sub(r, a)
            return y - d*x, x

    x = mod_inv_sub(a,m)[0]
    if x < 0:
        x += (-x//m + 1) * m
    return x

factor = 1

L = min(K-2,M*N-K)
for i in range(M*N-L-1,M*N-1):
    factor = factor * i % MOD
for i in range(1,L+1):
    factor = factor * mod_inv(i) % MOD
ans = ((M*(N-1)*(N+1) + N*(M-1)*(M+1))*M*N //6 * factor) % MOD
print(ans)

Solution (Python):

"""
xi,xjのペアにおいて、このペアは何通り出てくるか
sum(全通り)sum(i)sum(j){|xi-xj|+|yi-yj|}
和なので、どんな順序でもいい

sum(i)sum(j)sum(全通り){|xi-xj|}
sum(i)sum(j)sum(全通り){|yi-yj|}

あるi,jで固定して、sum(全通り){|xi-xj|}を求めるヤツ
固定したi,jが現れる回数は、

pat = NM-2_C_K-2 こいつはどこでも一緒

pat*sum(i)sum(j){|xi-xj|}を計算する
すべてのi,j(i>j)ペアにおいて、xi,xjの和は
sum(i=1,N){i*(N-i)}
"""
from math import factorial

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



def mod_combination(n,r,mod):
    if n<0 or r<0 or n<r:
        return 0
    if n==0 or r == 0:
        return 1
    a = factorial(n) % mod
    b = factorial(r) % mod
    c = factorial(n-r) % mod
    return (a*pow(b,mod-2,mod)*pow(c,mod-2,mod)) % mod 

pat = mod_combination(N*M-2,K-2,mod)

n_x = M**2*((N-1)*N*(N+1))//6
n_y = N**2*((M-1)*M*(M+1))//6
print(pat*(n_x+n_y)%mod)

Solution (Python):

import copy
import random
import bisect 
import fractions 
import math
import sys

mod = 10**9+7
sys.setrecursionlimit(mod) 

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

def mCn(m,n):
    return math.factorial(m)//(math.factorial(n)*math.factorial(m-n)) % mod

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

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

for i in range( 2, N1 + 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 i in range(1,M):
    ans += (M-i)*N**2*i*(cmb(M*N-2,K-2,mod))
    if ans > mod:
        ans = ans % mod

for i in range(1,N):
    ans += (N-i)*M**2*i*(cmb(M*N-2,K-2,mod))
    if ans > mod:
        ans = ans % mod

print(ans)

Solution (Python):
import sys, re
from collections import deque, defaultdict, Counter
from math import ceil, sqrt, hypot, factorial, pi, sin, cos, radians, log2
from itertools import permutations, combinations, product
from operator import itemgetter, mul
from copy import deepcopy
from string import ascii_lowercase, ascii_uppercase, digits

def input(): return sys.stdin.readline().strip()
def INT(): return int(input())
def MAP(): return map(int, input().split())
def LIST(): return list(map(int, input().split()))
sys.setrecursionlimit(10 ** 9)
INF = float('inf')
MOD = 10 ** 9 + 7

N, M, K = MAP()

def C(n, r):
	return (((factorial[n] * power(factorial[r], MOD-2) % MOD) * power(factorial[n-r], MOD-2)) % MOD)

def power(x, y):
	if y == 0:
		return 1
	elif y == 1:
		return x % MOD
	elif y % 2 == 0:
		return power(x, int(y/2)) ** 2 % MOD
	else:
		return power(x, int((y-1)/2)) ** 2 * x % MOD

factorial = [1]
for n in range(1, 2*10**5+1):
    factorial.append(factorial[n-1]*n%MOD)

pat = C(N*M-2, K-2)

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

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

print(pat*(X+Y)%MOD)

Solution (Python):
import math
import sys

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


def ncr(n, r, mod=None):
    """
    scipy.misc.comb または scipy.special.comb と同じ
    組み合わせの数 nCr
    :param int n:
    :param int r:
    :param int mod: 3 以上の素数であること
    :rtype: int
    """
    assert n >= r

    def inv(a):
        """
        a の逆元
        :param a:
        :return:
        """
        return pow(a, mod - 2, mod)

    if mod:
        return math.factorial(n) * inv(math.factorial(r)) * inv(math.factorial(n - r)) % mod
    else:
        return math.factorial(n) // math.factorial(r) // math.factorial(n - r)


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


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


MOD = 10 ** 9 + 7


k2 = (pyramid(N - 1) * M ** 2 + pyramid(M - 1) * N ** 2) % MOD

c = ncr(N * M - 2, K - 2, mod=MOD)
print(k2 * c % MOD)

Solution (Python):
def powerf(a, b, md):
    """a^b mod md"""
    if b == 0:
        return 1

    if b % 2 == 0:
        pf2 = powerf(a, b // 2, md)
        return pf2**2 % md
    else:
        return (a * powerf(a, b-1, md)) % md


def comb(n, k, md):
    """nCk mod md"""
    from math import factorial

    if n < 0 or k < 0 or n < k:
        return 0
    if n == 0 and k == 0:
        return 1

    a = factorial(n) % md
    b = factorial(k) % md
    c = factorial(n-k) % md
    return a * powerf(b, md-2, md) * powerf(c, md-2, md) % md


M = 10**9 + 7

n, m, k = list(map(int, input().split()))
ans = 0
for dx in range(1, m):
    ans += dx * (m - dx) * n**2

for dy in range(1, n):
    ans += dy * (n - dy) * m**2

print(ans * comb(n*m-2, k-2, M) % M)

Solution (Python):




MOD = 10**9+7
Facts = [1]*(2*10**5+1)
for i in range(2,2*10**5+1):
  Facts[i] = (Facts[i-1]*i) % MOD

def mul(a,b):
  return (a*b) % MOD

def power(a,b):
  if   b==0   : return 1
  elif b==1   : return a % MOD
  elif b%2==0 : return power(a,b//2)**2 % MOD
  else        : return power(a,b//2)**2 * a % MOD

def div(a,b):
  return mul(a,power(b,MOD-2))


def combi(n,m,k):
  if k <= 2 or k >= n*m:
    return 1
  else:
    return div(Facts[n*m-2],(Facts[k-2]*Facts[n*m-k]))


def main():
  N,M,K = map(int, input().split())
  com = combi(N,M,K)
  ans = 0
  for d in range(1,N):
    ans = (ans + d * (N-d)*M**2 * com) % MOD
  for d in range(1,M):
    ans = (ans + d * (M-d)*N**2 * com) % MOD
  return ans

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

Solution (Python):
def combnk_mod(n, k, p):
    
    
    
    

    if  2 * k > n:
        k = n - k
    x = 1 
    y = 1 
    for i in range(1, k + 1):
        x *= pow(i, p - 2, p)
        y *= n - i + 1
        x %= p
        y %= p
    return (x * y) % p

n, m, k = map(int, input().split())
MOD = 10 ** 9 + 7
p = combnk_mod(n * m - 2, k - 2, MOD)
x, y = 0, 0
for i in range(2):
    for d in range(1, n):
        x += d * (n - d) * m ** 2
        x %= MOD
    x, y = y, x
    n, m = m, n
print(p * (x + y) % MOD)

Solution (Python):
import sys
from math import factorial

input = sys.stdin.readline

Mod = 10**9 + 7


def Cmb(n, k):
    return factorial(n) // (factorial(n - k) * factorial(k)) % Mod


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

ans = 0

for i in range(1, n + 1):
    ans += i * (n - i) * m**2 % Mod

for i in range(1, m + 1):
    ans += i * (m - i) * n**2 % Mod

ans *= Cmb(n * m - 2, k - 2)
ans %= Mod

print(ans)

Solution (Python):
class FactMod:

    def __init__(self, n, mod):
        self.mod = mod
        self.f = [1] * (n + 1)
        for i in range(1, n + 1):
            self.f[i] = self.f[i - 1] * i % mod

        self.inv = [pow(self.f[-1], mod - 2, mod)]
        for i in reversed(range(1, n + 1)):
            self.inv.append(self.inv[-1] * i % mod)
        self.inv.reverse()

    def fact(self, n):
        """
        :return: n!
        """
        return self.f[n]

    def comb(self, n, r):
        """
        :return: nCr
        """
        return self.f[n] * self.inv[n - r] * self.inv[r] % self.mod

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

NM = N * M
fa = [1] * (NM + 1)

for i in range(1, NM + 1):
    fa[i] = (fa[i - 1] * i) % mod

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

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

factmod = FactMod(NM-2, mod)
ans *= factmod.comb(NM-2, K-2)
ans %= mod
print(ans)

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

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(a,mod = DIV):
    return pow(a,mod-2,mod)

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 * pow(i, mod-2, 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
print(main())

Solution (Python):
from math import factorial

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

n, r = N * M - 2, K - 2

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)

Solution (Python):
mod = 10**9 + 7
n,m,k = 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)

Solution (Python):
import math
N,M,k = map(int,input().split())
p = N*M
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

P = cmb(p-2,k-2)
Q = N*(N-1)*(N+1)*M*M//6
R = M*(M-1)*(M+1)*N*N//6
print((P*(Q+R))%(10**9+7))

Solution (Python):

import sys
INF = float("inf")

MOD = 1000000007  


def cmb(n, r):
    r = min(r, n-r)
    if r == 0:
        return 1
    if r == 1:
        return n

    numer = [n - r + k + 1 for k in range(r)]
    denom = [k + 1 for k in range(r)]

    for p in range(2, r+1):
        pivot = denom[p - 1]
        if pivot > 1:
            offset = (n - r) % p
            for k in range(p-1, r, p):
                numer[k - offset] /= pivot
                denom[k] /= pivot

    result = 1
    for k in range(r):
        if numer[k] > 1:
            result *= int(numer[k])
            result %= MOD

    return result


def solve(N: int, M: int, K: int):

    sx = (M*(M+1)*(M-1)//6) % MOD
    sx = (sx * N*N) % MOD
    sy = (N*(N+1)*(N-1)//6) % MOD
    sy = (sy * M*M) % MOD
    s = (sx + sy) % MOD
    s = (s*cmb(N*M-2, K-2)) % MOD
    print(s)

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

Solution (Python):

import sys
INF = float("inf")

MOD = 1000000007  


def cmb(n, r):
    r = min(r, n-r)
    if r == 0:
        return 1
    if r == 1:
        return n

    numer = [n - r + k + 1 for k in range(r)]
    denom = [k + 1 for k in range(r)]

    for p in range(2, r+1):
        pivot = denom[p - 1]
        if pivot > 1:
            offset = (n - r) % p
            for k in range(p-1, r, p):
                numer[k - offset] /= pivot
                denom[k] /= pivot

    result = 1
    for k in range(r):
        if numer[k] > 1:
            result *= int(numer[k])

    return result


def solve(N: int, M: int, K: int):

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

    sy = 0
    for y in range(1, N+1):
        sy += y*(y-1)//2
        sy %= MOD
    sy *= M*M
    sy %= MOD

    s = (sx+sy) % MOD
    s = (s*cmb(N*M-2, K-2)) % MOD
    print(s)
    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 math
N,M,k = map(int,input().split())
p = N*M
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

P = cmb(p-2,k-2)
Q = N*(N-1)*(N+1)*M*M//6
R = M*(M-1)*(M+1)*N*N//6
print((P*(Q+R))%(10**9+7))

Solution (Python):
import numpy as np
MOD = 1000000007

def mod_add(a, b):
  return (a+b)%MOD

def mod_sub(a, b):
  return (a+MOD-b)%MOD

def mod_mul(a, b):
  return a*b%MOD


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):
    x,  = extgcd(a, MOD)
    return (MOD+x%MOD)%MOD

def mod_div(a, b):
    return a*mod_inv(b)%MOD

def mod_pow(a, b):
    res = 1
    while b > 0:
        if b & 1:
            res = res * a * MOD
        a = a * a % MOD
        b >>= 1
    return res

MAX = 510000

fac  = [1,1]
finv = [1,1]
inv  = [1,1]

def COMinit():
    for i in range(2, MAX):
        fac.append(fac[i - 1] * i % MOD)
        inv.append(MOD - inv[MOD%i] * (MOD // i) % MOD)
        finv.append(finv[i - 1] * inv[i] % MOD)

def mod_COM(n, k):
    if n < k:
        return 0
    if n < 0 or k < 0:
        return 0
    return fac[n] * (finv[k] * finv[n - k] % MOD) % MOD

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

res = 0
COMinit()
for d in range(1, N):
    tmp = mod_COM(M*N-2, K-2)
    tmp = mod_mul(tmp, d)
    tmp = mod_mul(tmp, (N-d)*M*M)
    res = mod_add(res, tmp)

for d in range(1, M):
    tmp = mod_COM(M*N-2, K-2)

    tmp = mod_mul(tmp, d)
    tmp = mod_mul(tmp, (M-d)*N*N)
    res = mod_add(res, tmp)


print(res)

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

P = 200005
fact = [0] * P
fact[0] = 1
for i in range(1, P):
    fact[i] = fact[i - 1] * i
    fact[i] %= MOD

rfact = [0] * P
rfact[P - 1] = pow(fact[P - 1], MOD - 2, MOD)
for i in range(P - 2, -1, -1):
    rfact[i] = rfact[i + 1] * (i + 1)
    rfact[i] %= MOD


def comb(n, k):
    return fact[n] * rfact[k] * rfact[n - k] % MOD


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

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

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

Solution (Python):
N, M, K = [int(_) for _ in input().split()]


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  
g1 = [1, 1]  
g2 = [1, 1]  
inverse = [0, 1]  

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


def count(n, m):
    ret = 0
    for d in range(1, n):
        ret += d * (n - d) * m**2
    return ret


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

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

MOD = 10**9+7


class ModComb:
    
    def __init__(self, n_max, m):
        self._n_max = n_max
        self._m = m
        self._mod_facts = [1] * (n_max+1)       
        self._mod_inv_facts = [1] * (n_max+1)   

        if not m >= 3:
            raise ValueError('Require mod number m >= 3')

        
        for i in range(2, n_max+1):
            self._mod_facts[i] = self._mod_facts[i-1] * i % m

        
        self._mod_inv_facts[n_max] = pow(self._mod_facts[n_max], m-2, m)
        for i in range(n_max-1, 1, -1):
            self._mod_inv_facts[i] = self._mod_inv_facts[i+1] * (i+1) % m

    def calc(self, n, r):
        
        if n > self._n_max:
            raise Exception('n is larger than n_max')

        if r == 0:
            return 1

        m = self._m
        return self._mod_facts[n] * self._mod_inv_facts[r] * self._mod_inv_facts[n-r] % m

modcomb = ModComb(N*M-2, MOD)
p = modcomb.calc(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)

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

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

ans = 0
t = nCr(N*M-2,K-2)
for d in range(1,M):
    ans += d*(M-d)*(N*N)*t
    ans %= MOD

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

print(ans)

Solution (Python):

import bisect
import heapq
import math
import random
import sys
from collections import Counter, defaultdict, deque
from decimal import ROUND_CEILING, ROUND_HALF_UP, Decimal
from functools import lru_cache, reduce
from itertools import combinations, combinations_with_replacement, product, permutations
from operator import add, mul, sub

sys.setrecursionlimit(100000)


def read_int():
    return int(input())


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


def read_float():
    return float(input())


def read_float_n():
    return list(map(float, input().split()))


def read_str():
    return input().strip()


def read_str_n():
    return list(map(str, input().split()))


def error_print(*args):
    print(*args, file=sys.stderr)


def mt(f):
    import time

    def wrap(*args, **kwargs):
        s = time.time()
        ret = f(*args, **kwargs)
        e = time.time()

        error_print(e - s, 'sec')
        return ret

    return wrap


class Mod:
    def __init__(self, m):
        self.m = m

    def add(self, a, b):
        return (a + b) % self.m

    def sub(self, a, b):
        return (a - b) % self.m

    def mul(self, a, b):
        return ((a % self.m) * (b % self.m)) % self.m

    def div(self, a, b):
        return self.mul(a, self.pow(b, self.m-2))

    def pow(self, x, y):
        if y == 0:
            return 1
        elif y == 1:
            return x % self.m
        elif y % 2 == 0:
            return self.pow(x, y//2)**2 % self.m
        else:
            return self.pow(x, y//2)**2 * x % self.m


class Combination:
    def __init__(self, n, 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)
        self.MOD = mod
        self.N = n
        self.g1 = g1
        self.g2 = g2
        self.inverse = inverse

    def __call__(self, n, r):
        if (r < 0 or r > n):
            return 0
        r = min(r, n-r)
        return self.g1[n] * self.g2[r] * self.g2[n-r] % self.MOD


@mt
def slv(N, M, K):
    m = Mod(10**9+7)
    b = Combination(N*M, 10**9+7)(N*M-2, K-2)

    MM = m.mul(M, M)
    NN = m.mul(N, N)
    ans = 0
    for d in range(1, N):
        ans = m.add(ans, m.mul(m.mul(m.mul(d, N-d), MM), b))

    for d in range(1, M):
        ans = m.add(ans, m.mul(m.mul(m.mul(d, M-d), NN), b))

    return ans


def main():
    N, M, K = read_int_n()
    print(slv(N, M, K))


if __name__ == '__main__':
    main()

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

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

def calc_rev(n, law):
	if n == 1:
		return 1
	
	pre_r = law
	pre_x = 0
	pre_y = 1
	cur_r = n
	cur_x = 1
	cur_y = 0
	
	while True:
		d = pre_r // cur_r
		r = pre_r % cur_r
		
		nxt_r = r
		nxt_x = pre_x - d * cur_x
		nxt_y = pre_y - d * cur_y
		
		if nxt_r == 1:
			ret = nxt_x
			while ret < 0:
				ret += law;
			return ret;
		elif nxt_r == 0:
			return None
		
		pre_r = cur_r;
		pre_x = cur_x;
		pre_y = cur_y;
		cur_r = nxt_r;
		cur_x = nxt_x;
		cur_y = nxt_y;

def calc_dsum(n, m):
	dx = 0
	for i in range(1, n):
		dx += i * (n - i)
	
	dy = 0
	for i in range(1, m):
		dy += i * (m - i)
	dx *= m * m
	dy *= n * n
	
	return dx + dy

def calc_csum(n, k):
	csum = 1
	for i in range(k):
		csum = (csum * (n - i)) % LAW
	for i in range(1, k + 1):
		csum = (csum * calc_rev(i, LAW)) % LAW
	

	
	csum = (csum * k * (k - 1)) % LAW
	csum = (csum * calc_rev(n * (n - 1), LAW)) % LAW
	
	return csum

def c(n, m, k):
	
	return (calc_csum(n * m, k) * calc_dsum(n, m)) % LAW

print(c(N, M, K))

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

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

Solution (Python):
from itertools import*
import math
from collections import*
from heapq import*
from bisect import bisect_left,bisect_right
from copy import deepcopy
inf = float("inf")
mod = 10**9+7
from functools import reduce
import sys
sys.setrecursionlimit(10**7)

n,m,k = map(int,input().split())
Max = n*m

fac = [0]*(Max)
finv = [0]*(Max)
inv = [0]*(Max)
fac[0]=fac[1]=1
finv[0]=finv[1]=1
inv[1]=1
for i in range(2,Max):
        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 C(n,r):
    if n < r:
        return 0
    if n < 0 or r < 0 :
        return 0
    return fac[n]*(finv[r]*finv[n-r]%mod)%mod

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


c = C(n*m-2,k-2)
res = 0

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

for i in range(1,m):
    res += i*(m-i)*(n**2)*c
    res %= mod
print(res)

Solution (Python):
import math
n,m,k=map(int,input().split())
mod=10**9+7
s=math.factorial(n*m-2)//(math.factorial(n*m-k)*math.factorial(k-2))
cnt=0
s%=mod
a=(m*m*s)%mod
for d in range(1,n):
        cnt+=d*(n-d)*a
        cnt%=mod
a=(n*n*s)%mod
for d in range(1,m):
        cnt+=d*(m-d)*a
        cnt%=mod
print(cnt)

Solution (Python):
def main():
     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 = 10**6
     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 = map(int,input().split())
     nm = n*m%(10**9+7)
     a = (cmb(nm-2,k-2,mod)%(10**9+7))*(nm*(m*(n**2-1)+n*(m**2-1)))//6
     print(a%(10**9+7))
     
if __name__=='__main__':
     main()

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










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


ans = 0
C = combi(N*M-2, K-2, ppp)

for i in range(1,N):
    ans += C * i * (N-i) * M**2
for j in range(1,M):
    ans += C * j * (M-j) * N**2


print(ans % ppp)

Solution (Python):
N,M,K=map(int,input().split())
O,S,T=10**9+7,N+M,6
for i in range(K):S=S*(N*M-i)%O;T=T*-~i%O
print(S*K*~-K*pow(T,O-2,O)%O)

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

mod = 10**9+7 
N = 10**6
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 = map(int,input().split())
nm = n*m%(10**9+7)
k = k%(10**9+7)
a = (cmb(nm-2,k-2,mod)%(10**9+7))*(nm*(m*(n**2-1)+n*(m**2-1)))//6
print(a%(10**9+7))

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



mod = int(1e9+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 )

def combination(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

  


var_1 = M**2 * (N-1) * N * (N+1) // 3
var_2 = N**2 * (M-1) * M * (M+1) // 3
cost_allpairs = int((var_1 + var_2) // 2)
cost_allpairs = cost_allpairs % mod 


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


ans = (cost_allpairs * cmb) % mod
print(int(ans))

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

def modinv(a, m):
    b = m
    u = 1
    v = 0
    while (b):
        t = a // b
        a -= t * b 
        a = a^b
        b = a^b
        a = a^b
        u -= t * v
        u = u^v
        v = u^v
        u = u^v
    u %= m 
    if (u < 0): u += m
    return u

def ncr(n, r):
    res = 1
    for i in range(1, r+1):
        res = res*n*modinv(i, div) % div
        n = n-1
    return res % div

result = 0
comb = ncr(N*M-2, K-2)
for d in range(1, N):
    result += d*comb*M*M*(N-d)

for d in range(1, M):
    result += d*comb*(M-d)*N*N

print(result % div)

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

def modinv(a, m):
    b = m
    u = 1
    v = 0
    while (b):
        t = a // b
        a -= t * b 
        a = a^b
        b = a^b
        a = a^b
        u -= t * v
        u = u^v
        v = u^v
        u = u^v
    u %= m 
    if (u < 0): u += m
    return u

def ncr(n, r):
    res = 1
    for i in range(1, r+1):
        res = res*n*modinv(i, div) % div
        n = n-1
    return res % div

result = 0
comb = ncr(N*M-2, K-2)
for d in range(1, N):
    result += d*comb*M*M*(N-d)

for d in range(1, M):
    result += d*comb*(M-d)*N*N

print(result % div)

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


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

for i in range(1, k - 2 + 1):
    ans = (ans * pow(i, MOD - 2, MOD)) % MOD

print((ans * (m ** 2 * n * (n ** 2 - 1) // 6 + n * (n - 1) * m
              * (m ** 2 - 1) // 6 + n * m * (m ** 2 - 1) // 6)) % MOD)

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

fact = [1]*(MAX+1)
for i in range(1, MAX+1):
    fact[i] = (fact[i-1]*i)%MOD
finv = [1]*MAX + [pow(fact[-1], MOD-2, MOD)]
for i in range(MAX, 0, -1):
    finv[i-1] = (finv[i]*i)%MOD

def comb(n, r):
    return ((fact[n]*finv[n-r])%MOD*finv[r])%MOD

cumsum = [0]*(MAX+1)
for i in range(1, max(n, m)+1):
    cumsum[i] = cumsum[i-1]+i

ans = 0
for i in range(n):
    ans += (cumsum[i] + cumsum[n-i-1])*(m**2)
for i in range(m):
    ans += (cumsum[i] + cumsum[m-i-1])*(n**2)

ans //= 2
ans = (ans*comb(n*m-2, k-2))%MOD
print(ans)