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

def nCk(n, k, mod):
    nu = 1
    de = 1 
    for i in range(1, k+1):
        nu *= (n-i+1)
        de *= i
        nu = nu % mod
        de = de % mod
    ans = nu * pow(de, mod-2, mod) % mod
    return ans

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

Solution (Python):
N, M, K = map(int, input().split())
P = int(1e9+7)
N -= 1
M -= 1

def comb(n, r):
    if n < r or r < 0:
        return 0
    nume = 1
    deno = 1
    for i in range(1, r + 1):
        nume = (nume * (n - i + 1)) % P
        deno = (deno * i) % P
    deno_inv = pow(deno, P-2, P)
    return (nume*deno_inv)%P

ans = 0

def sum_man_dist_ur(u, r):
    ret = 0
    u, r = min(u, r), max(u, r)
    if u == 0:
        ret =  r * (r + 1) // 2
    else:
        ret += (u * (u + 1) * (u - 1)) // 3
        ret += (u+1)*(r-u+1)*(r+u)//2
        ret += u * (u + 1) * (u + 3 * r + 2) // 6
    return ret%P

def sum_man_dist(x, y):
    ret = 0
    ret += sum_man_dist_ur(N - x, M - y)
    ret = (ret + sum_man_dist_ur(N - x, y)) % P
    ret = (ret + sum_man_dist_ur(x, M - y)) % P
    ret = (ret + sum_man_dist_ur(x, y)) % P
    ret -= (N-x)*(N-x+1)//2
    ret -= (x)*(x+1)//2
    ret -= (M-y)*(M-y+1)//2
    ret -= (y) * (y + 1) // 2
    return ret%P

ans = 0
c = comb((N+1)*(M+1)-2, K-2)
for x in range(N + 1):
    for y in range(M + 1):
        ans = (ans + c * sum_man_dist(x, y)) % P
print(ans*pow(2, P-2, P)%P)

Solution (Python):
ma = lambda :map(int,input().split())
lma = lambda :list(map(int,input().split()))
tma = lambda :tuple(map(int,input().split()))
ni = lambda:int(input())
yn = lambda fl:print("Yes") if fl else print("No")
import collections
import math
import itertools
import heapq as hq
import sys
import numpy as np
mod=10**9+7
def ncr_1(n,r,p):
    x = 1
    y = 1
    for i in range(r):
        x*=(n-i)
        y*=(i+1)
        x%=p
        y%=p
    return x*modinv(y,p) %p
def modinv(a,p):
    b,u,v = p,1,0
    while b:
        t = a//b
        a -= t*b
        a,b = b,a
        u -=t*v
        u,v = v,u
    u %=p
    return u
n,m,k = ma()
ans=0
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
print((ans*ncr_1(n*m-2,k-2,mod))%mod)

Solution (Python):
N=10**6
p=10**9+7
f,finv,inv=[0]*(N*5),[0]*(N*5),[0]*(N*5)

def nCr_init(L,p):
    for i in range(L+1):
        if i<=1:
            f[i],finv[i],inv[i]=1,1,1
        else:
            f[i]=(f[i-1]*i)%p
            inv[i]=(p-inv[p%i]*(p//i))%p
            finv[i]=(finv[i-1]*inv[i])%p

def nPr(n,r,p):
    if 0<=r<=n and 0<=n:
        return (f[n]*finv[n-r])%p
    else:
        return 0
    
nCr_init(N,p)

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

print((nPr(N*M,K,p)*((M+N)*inv[6]*finv[K-2])%p)%p)

Solution (Python):
"""
基本の方針は「各マンハッタン距離が何回足されるか」という考え方で距離の値を積み上げて行けばよい

まず、XとYは独立に考えることにする。
X軸の距離の方から積み上げる。
X1とX2の距離が足される回数は、
A.X1の列とX2の列の座標の組み合わせの数　×　B.その組み合わせが含まれるような座標の組み合わせの数
である。
A　＝　行数**２
B　＝　CMB(N*M-2 , K-2)である

同様の計算をY軸に関しても行えばよい。
"""
def cmb(n, r, p):
    if (r < 0) or (n < r):
        return 0
    r = min(r, n - r)
    return fact[n] * factinv[r] * factinv[n-r] % p

p = 10 ** 9 + 7
N = 10 ** 6  
fact = [1, 1]  
factinv = [1, 1]  
inv = [0, 1]  
 
for i in range(2, N + 1):
    fact.append((fact[-1] * i) % p)
    inv.append((-inv[p % i] * (p // i)) % p)
    factinv.append((factinv[-1] * inv[-1]) % p)

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

ans = 0


A = M*M
B = cmb(N*M-2,K-2,p)
for i in range(1,N):
    ans += ((i*i)-(i-1+0)*i//2)*A*B%p
    ans %= p


A = N*N
B = cmb(N*M-2,K-2,p)
for i in range(1,M):
    ans += ((i*i)-(i-1+0)*i//2)*A*B%p
    ans %= p

print(ans)

Solution (Python):
mod = 10 ** 9 + 7
class nCr():
    def __init__(self, n):
        self.n = n
        self.fa = [1] * (self.n + 1)
        self.fi = [1] * (self.n + 1)
        for i in range(1, self.n + 1):
            self.fa[i] = self.fa[i - 1] * i % mod
            self.fi[i] = pow(self.fa[i], mod - 2, mod)
    def comb(self, n, r):
        if n < r:return 0
        if n < 0 or r < 0:return 0
        return self.fa[n] * self.fi[r] % mod * self.fi[n - r] % mod
n, m, k = map(int, input().split())
ans = 0
c = nCr(n * m)
for i in range(n):
    ans += i * (n - i) * m * m
    ans %= mod
for i in range(m):
    ans += i * (m - i) * n * n
    ans %= mod
print(ans * c.comb(n * m - 2, k - 2) % mod)

Solution (Python):




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


MOD = 10 ** 9 + 7
table_len = 2 * 10 ** 5 + 10

fac = [1, 1]
for i in range(2, table_len):
    fac.append(fac[-1] * i % MOD)

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








distsum1 = 0
for l in range(1, N):
    distsum1 += l * (N - l)
    distsum1 %= MOD
distsum1 *= pow(M, 2, MOD)
distsum1 %= MOD


distsum2 = 0
for l in range(1, M):
    distsum2 += l * (M - l)
    distsum2 %= MOD
distsum2 *= pow(N, 2, MOD)
distsum2 %= MOD

print((distsum1 + distsum2) % MOD * fac[N * M - 2] % MOD * finv[K - 2] % MOD * finv[N * M - K] % MOD)

Solution (Python):
from itertools import combinations
import math
def cmb(n, r, p):
    if (r < 0) or (n < r):
        return 0
    r = min(r, n - r)
    return fact[n] * factinv[r] * factinv[n-r] % p

p = 10 ** 9 + 7
N = 10 ** 6  
fact = [1, 1]  
factinv = [1, 1]  
inv = [0, 1]  
 
for i in range(2, N + 1):
    fact.append((fact[-1] * i) % p)
    inv.append((-inv[p % i] * (p // i)) % p)
    factinv.append((factinv[-1] * inv[-1]) % p)


N,M,K = map(int,input().split())
ans = 0
mod = 10**9 +7
by = cmb(N*M-2,K-2,mod)


for i in range(N):
    for j in range(M):
        cnt = M-j-1
        add = (((j+1)+(M-1))*cnt//2 - j*cnt)*N%mod
        ans += add*by%mod
        ans %= mod

for i in range(M):
    for j in range(N):
        cnt = N-j-1
        add = (((j+1)+(N-1))*cnt//2 - j*cnt)*M%mod
        ans += add*by%mod
        ans %= mod


print(ans)

Solution (Python):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools
from collections import deque

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

DR = [1, -1, 0, 0]
DC = [0, 0, 1, -1]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
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 LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
     

fac = [-1] * (10**6+1)
inv = [-1] * (10**6+1)
finv = [-1] * (10**6+1)

fac[0] = fac[1] = 1
inv[1] = 1
finv[0] = finv[1] = 1

def initNCMMod(limit):
    for i in range(2, limit):
        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 NCMMod(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
initNCMMod(10**6)

def checkf(R, C):
    ans = 0
    for r1 in range(R):
        for c1 in range(C):
            for r2 in range(R):
                for c2 in range(C):
                    ans += abs(r1-r2) + abs(c1-c2)
    return ans


def main():
    R, C, K = LI()
    def f(i):
        return (i * (i+1) // 2) % mod

    ans = 0
    for i in range(C):
        ans += f(C-i-1) * R * R % mod
        ans += f(i) * R * R % mod

    for i in range(R):
        ans += f(R-i-1) * C * C % mod
        ans += f(i) * C * C % mod
    ans *= inv[2]
    ans %= mod
    ans *= NCMMod(R*C-2, K-2)
    ans %= mod

    print(ans)

main()

Solution (Python):
class ModCalc:
    mod = 10**9 + 7

    def __init__(self, n):
        self.n = n
        self.inv = [1] * (self.n + 1)
        self.fac = [1] * (self.n + 1)
        self.finv = [1] * (self.n + 1)
        self.inv[0] = 0
        self.inv_table()

    def inv_table(self):
        for i in range(2, self.n+1):
            self.inv[i] = self.inv[self.mod % i] * (self.mod - self.mod // i) % self.mod
            self.fac[i] = self.fac[i-1] * i % self.mod
            self.finv[i] = self.finv[i-1] * self.inv[i] % self.mod

    def comb(self, n, r):
        if n < r:
            return 0
        if n < 0 or r < 0:
            return 0

        return self.fac[n] * self.finv[r] % self.mod * self.finv[n-r] % self.mod

    def perm(self, n, r):
        if n < r:
            return 0
        if n < 0 or r < 0:
            return 0

        return self.fac[n] * self.finv[n-r] % self.mod

    def fact(self, n):
        return self.fac[n]


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

mc = ModCalc(n*m)
mod = 10**9 + 7

ret = 0
for d in range(1, n):
    ret += d * (n - d) % mod * m % mod * m % mod * mc.comb(n * m - 2, k-2) % mod
    ret %= mod

for d in range(1, m):
    ret += d * (m - d) % mod * n % mod * n % mod * mc.comb(n * m - 2, k - 2) % mod
    ret %= mod
print(ret)

Solution (Python):
import sys

read = sys.stdin.read
readline = sys.stdin.readline
readlines = sys.stdin.readlines
sys.setrecursionlimit(10 ** 9)
INF = 1 << 60
MOD = 1000000007


def main():
    N, M, K = map(int, readline().split())

    COM_MAX = N * M

    fac, finv, inv = [0] * (COM_MAX + 1), [0] * (COM_MAX + 1), [0] * (COM_MAX + 1)
    fac[0] = fac[1] = finv[0] = finv[1] = inv[1] = 1

    for i in range(2, COM_MAX + 1):
        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 com(n, r):
        if n < 0 or r < 0 or n < r:
            return 0
        return fac[n] * (finv[r] * finv[n - r] % MOD) % MOD

    def solve(n, m, k):
        ans = 0
        for d in range(1, n):
            ans = (ans + com(n * m - 2, k - 2) * (m ** 2) * (n - d) * d) % MOD

        return ans

    ans = (solve(N, M, K) + solve(M, N, K)) % MOD
    print(ans)
    return


if __name__ == '__main__':
    main()

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










factorial = [1 for i in range(N*M + 1)]
for i in range(1, N * M + 1):
	if i==1:
		factorial[i] = 1
	else:
		factorial[i] = factorial[i-1]*i % mod

def comb(n, k):
	return (factorial[n]*pow(factorial[n-k]*factorial[k], -1, mod)) % mod

xscore = 0
for i in range(1, N):
	xscore += i * (N-i)
xscore %= mod
xscore *= M**2 * comb(N*M-2, K-2)
xscore %= mod

yscore = 0
for i in range(1, M):
	yscore += i * (M-i)
yscore %= mod
yscore *= N**2 * comb(N*M-2, K-2)
yscore %= mod

print((xscore + yscore) % mod)

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

factorial=[1 for i in range(N*M+1)]
for i in range(1,N*M+1):
    if i==1:factorial[i]=1
    else:factorial[i] = factorial[i-1]*i % mod

def comb(n,k):
    return (factorial[n]*pow(factorial[n-k]*factorial[k], -1, mod))%mod
  
xscore = 0
for i in range(1,N):
  xscore += i * (N-i)
xscore %= mod
xscore *= M**2 * comb(N*M-2,K-2)
xscore %= mod

yscore = 0
for i in range(1,M):
  yscore += i * (M-i)
yscore %= mod
yscore *= N**2 * comb(N*M-2,K-2)
yscore %= mod

print((xscore + yscore) %mod)

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

M = [1]  
mul = 1
for i in range(1, l):
    mul = (mul * i) % mod
    M.append(mul)

def pow(x, y, mod):  
    ans = 1
    while y > 0:
        if y % 2 == 1:
            ans = (ans * x) % mod
        x = (x**2) % mod
        y //= 2
    return ans

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



MI = [0] * (l-1) +[inv(M[l-1], mod)]  
for i in range(l-2, -1, -1):
    MI[i] = MI[i+1] * (i+1) % mod

def C(x, y):  
    if y < 0 or y > x:
        return 0
    elif x > l:  
        y = min(y, x-y)
        ans = 1
        for i in range(x, x-y, -1):
            ans = (ans * i) % mod
        return (ans * MI[y]) % mod
    else:  
        ans = M[x]
        ans = (ans * MI[y]) % mod
        return (ans * MI[x-y]) % mod

def P(x, y):  
    if y < 0 or y > x:
        return 0
    elif x > l:  
        y = min(y, x-y)
        ans = 1
        for i in range(x, x-y, -1):
            ans = (ans * i) % mod
        return ans % mod
    else:  
        ans = M[x]
        return (ans * MI[x-y]) % mod

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

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

Solution (Python):






import sys
input = sys.stdin.readline

import math

import bisect
import heapq
from collections import defaultdict
from collections import deque
from collections import Counter
from functools import lru_cache





MOD = 10**9+7
INF = float('inf')
AZ = "abcdefghijklmnopqrstuvwxyz"






def I(): return int(input().strip())
def S(): return input().strip()
def IL(): return list(map(int,input().split()))
def SL(): return list(map(str,input().split()))
def ILs(n): return list(int(input()) for _ in range(n))
def SLs(n): return list(input().strip() for _ in range(n))
def ILL(n): return [list(map(int, input().split())) for _ in range(n)]
def SLL(n): return [list(map(str, input().split())) for _ in range(n)]



def DD(arg): return defaultdict(arg)


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


kaijo_memo = []
def kaijo(n):
  if(len(kaijo_memo) > n): return kaijo_memo[n]
  if(len(kaijo_memo) == 0): kaijo_memo.append(1)
  while(len(kaijo_memo) <= n): kaijo_memo.append(kaijo_memo[-1] * len(kaijo_memo) % MOD)
  return kaijo_memo[n]

gyaku_kaijo_memo = []
def gyaku_kaijo(n):
  if(len(gyaku_kaijo_memo) > n): return gyaku_kaijo_memo[n]
  if(len(gyaku_kaijo_memo) == 0): gyaku_kaijo_memo.append(1)
  while(len(gyaku_kaijo_memo) <= n): gyaku_kaijo_memo.append(gyaku_kaijo_memo[-1] * pow(len(gyaku_kaijo_memo),MOD-2,MOD) % MOD)
  return gyaku_kaijo_memo[n]

def nCr(n,r):
  if n == r: return 1
  if n < r or r < 0: return 0
  ret = 1
  ret = ret * kaijo(n) % MOD
  ret = ret * gyaku_kaijo(r) % MOD
  ret = ret * gyaku_kaijo(n-r) % MOD
  return ret


def factorization(n):
  arr = []
  temp = n
  for i in range(2, int(-(-n**0.5//1))+1):
    if temp%i==0:
      cnt=0
      while temp%i==0: 
        cnt+=1 
        temp //= i
      arr.append([i, cnt])
  if temp!=1: arr.append([temp, 1])
  if arr==[]: arr.append([n, 1])
  return arr


def make_divisors(n):
  divisors = []
  for i in range(1, int(n**0.5)+1):
    if n % i == 0:
      divisors.append(i)
      if i != n // i: 
        divisors.append(n//i)
  return divisors


def make_primes(N):
  max = int(math.sqrt(N))
  seachList = [i for i in range(2,N+1)]
  primeNum = []
  while seachList[0] <= max:
    primeNum.append(seachList[0])
    tmp = seachList[0]
    seachList = [i for i in seachList if i % tmp != 0]
  primeNum.extend(seachList)
  return primeNum


def gcd(a, b):
    while b: a, b = b, a % b
    return a


def lcm(a, b):
    return a * b // gcd (a, b)


def count_bit(n):
  count = 0
  while n:
    n &= n-1
    count += 1
  return count


def base_10_to_n(X, n):
  if X//n: return base_10_to_n(X//n, n)+[X%n]
  return [X%n]

def base_n_to_10(X, n):
  return sum(int(str(X)[-i-1])*n**i for i in range(len(str(X))))

def base_10_to_n_without_0(X, n):
  X -= 1
  if X//n: return base_10_to_n_without_0(X//n, n)+[X%n]
  return [X%n]


def int_log(n, a):
  count = 0
  while n>=a:
    n //= a
    count += 1
  return count




N,M,K = IL()

print((((N*(N+1)*(N-1)//6)*M**2+(M*(M+1)*(M-1)//6)*N**2) * nCr(K,2)*nCr(N*M,K)*inv(nCr(N*M,2)))%MOD)

Solution (Python):







class ModInt:
    def __init__(self, x, MOD=10 ** 9 + 7):
        '''
        pypyじゃないと地獄みたいな遅さなので注意
        '''
        self.mod = MOD
        self.x = x % MOD

    def __str__(self):
        return str(self.x)

    __repr__ = __str__

    def __add__(self, other):
        if isinstance(other, ModInt):
            return ModInt(self.x + other.x, self.mod)
        else:
            return ModInt(self.x + other, self.mod)

    def __sub__(self, other):
        if isinstance(other, ModInt):
            return ModInt(self.x - other.x, self.mod)
        else:
            return ModInt(self.x - other, self.mod)

    def __mul__(self, other):
        if isinstance(other, ModInt):
            return ModInt(self.x * other.x, self.mod)
        else:
            return ModInt(self.x * other, self.mod)

    def __truediv__(self, other):
        if isinstance(other, ModInt):
            return ModInt(self.x * pow(other.x, self.mod - 2, self.mod), self.mod)
        else:
            return ModInt(self.x * pow(other, self.mod - 2, self.mod), self.mod)

    def __pow__(self, other):
        if isinstance(other, ModInt):
            return ModInt(pow(self.x, other.x, self.mod), self.mod)
        else:
            return ModInt(pow(self.x, other, self.mod), self.mod)

    __radd__ = __add__

    def __rsub__(self, other):  
        if isinstance(other, ModInt):
            return ModInt(other.x - self.x, self.mod)
        else:
            return ModInt(other - self.x, self.mod)

    __rmul__ = __mul__

    def __rtruediv__(self, other):
        if isinstance(other, ModInt):
            return ModInt(other.x * pow(self.x, self.mod - 2, self.mod), self.mod)
        else:
            return ModInt(other * pow(self.x, self.mod - 2, self.mod), self.mod)

    def __rpow__(self, other):
        if isinstance(other, ModInt):
            return ModInt(pow(other.x, self.x, self.mod), self.mod)
        else:
            return ModInt(pow(other, self.x, self.mod), self.mod)


def combination(n, r):
    
    if r > n:
        return 0  
    r = min(r, n - r)
    nf = rf = ModInt(1)
    for i in range(r):
        nf = nf * (n - i)
        rf = rf * (i + 1)
    return nf / rf


N, M, K = map(int, input().split())
C = combination(N * M - 2, K - 2)
x_term = ModInt(0)
for dx in range(1, N):
    x_term += (N - dx) * dx
x_term *= M**2

y_term = ModInt(0)
for dy in range(1, M):
    y_term += (M - dy) * dy
y_term *= N**2

print(C * (x_term + y_term))

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

class Comb:
    def __init__(self, n, mod=pow(10, 9)+7, f=True, i=False):
        self.n = n
        self.mod = mod
        self.factorial_table = None
        self.inverse_table = None
        if f: self.get_factorial()
        if i: self.get_inverse()
        
    def get_factorial(self):
        self.factorial_table = [1]*(self.n+1)
        for i in range(1, self.n+1):
            self.factorial_table[i] = self.factorial_table[i-1]*i%self.mod
        self.inverse_table = None
    
    def get_inverse(self):
        mod = self.mod
        self.inverse_table = [1]*(self.n+1)
        for i in range(self.n+1):
            self.inverse_table[i] = pow(self.factorial_table[i], -1, mod)
    
    def factorial(self, n):
        if self.factorial_table == None:
            return None
        return self.factorial_table[n]
    
    def inverse(self, n):
        if self.inverse_table == None:
            return None
        return self.inverse_table[n]   
    
    def permutation(self, n, k):
        if self.factorial_table == None:
            return None
        mod = self.mod
        if self.inverse_table == None:
            return self.factorial_table[n]*pow(self.factorial_table[n-k], -1, mod)%mod
        return self.factorial_table[n]*self.inverse_table[n-k]%mod
            
    def comb(self, n, k):
        mod = self.mod
        if self.factorial_table == None:
            res = 1
            for i in range(k):
                res *= n-i
                res *= i+1
                res %= mod
            return res
        if self.inverse_table == None:
            res = self.factorial_table[n]*pow(self.factorial_table[n-k], -1, mod)%mod\
                                    *pow(self.factorial_table[k], -1, mod)%mod
        else:
            res = self.factorial_table[n]*self.inverse_table[n-k]%mod\
                                    *self.inverse_table[k]%mod
        return res
    
    def recomb(self, n, k):
        if n+k-1 > self.n:
            return None
        return self.comb(n+k-1, k)
        


def main():
    n, m, k = map(int, input().split())
    mod = pow(10, 9) + 7
    
    comb = Comb(n*m, i=True)
    
    ans = 0
    
    for i in range(n):
        ans += (i) * (n-i) * m**2
        ans %= mod
        
    
    for i in range(m):
        ans += (i) * (m-i) * n**2
        ans %= mod
    
    
    
    ans *= comb.comb(n*m-2, k-2)
    ans %= mod
    
    print(ans)
    
    
if __name__ == "__main__":
    main()

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


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

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




def RD(): return input().rstrip().decode()
def II(): return int(input())
def FI(): return int(input())
def MI(): return map(int,input().split())
def MF(): return map(float,input().split())
def LI(): return list(map(int,input().split()))
def LF(): return list(map(float,input().split()))
def TI(): return tuple(map(int,input().split()))



def main():
	n,m,k=MI()
	mod=10**9+7

	ans=0

	cnt=0
	for i in range(1,m):
		cnt+=i*(m-i)
		cnt%=mod

	ans+=cnt*n*n
	ans%=mod

	cnt=0
	for i in range(1,n):
		cnt+=i*(n-i)
		cnt%=mod

	ans+=cnt*m*m
	ans%=mod

	f=[1]
	for i in range(1,n*m+5):
		f.append((f[-1]*i)%mod)

	finv=[]
	for i in f:
		finv.append(pow(i,mod-2,mod))
	

	

	ans*=k*(k-1)

	ans*=f[n*m-2]
	ans%=mod
	


	ans*=finv[n*m-2-(k-2)]
	ans%=mod
	


	ans*=finv[k]
	ans%=mod



	print(ans)






























if __name__ == "__main__":
	main()

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

mod = 10**9 + 7





ans = 0
for y in range(1, n):
    ans += y * (n - y) * m * m % mod
for x in range(1, m):
    ans += x * (m - x) * n * n % mod


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



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

print(ans)

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

def com(n,k,mod):
    res=1
    tmp=1
    for i in range(1,k+1):
        res=res*(n-i+1)%mod
        tmp=tmp*i%mod
    a=pow(tmp,mod-2,mod)
    return res*a%mod

N,M,K=map(int,input().split())
C=com(N*M-2,K-2,MOD)

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

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


def nCr(n, r, MOD=10**9 + 7):
    from functools import reduce
    if r == 0:
        return 1
    num = reduce(lambda x, y: x * y % MOD, range(n, n - r, -1))
    den = reduce(lambda x, y: x * y % MOD, range(1, r + 1))
    return num * pow(den, MOD - 2, MOD) % MOD


ncr = nCr(n * m - 2, k - 2)
ans = 0
for d in range(1, m):
    temp = d * (m - d) * n**2
    temp %= mod
    ans += temp
    ans %= mod
for d in range(1, n):
    temp = d * (n - d) * m**2
    temp %= mod
    ans += temp
    ans %= mod
ans *= ncr
ans %= mod
print(ans)

Solution (Python):

def xgcd(a, b):
    x0, y0, x1, y1 = 1, 0, 0, 1
    while b != 0:
        q, a, b = a // b, b, a % b
        x0, x1 = x1, x0 - q * x1
        y0, y1 = y1, y0 - q * y1
    return x0

def modinv(a, m):
    x = xgcd(a, m)
    return x % m


max_num = 10 ** 9 + 7
res = 0
n, m, k = map(int, input().strip().split())

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

r = min(n*m-2, k-2)
for i in range(1, r+1):
    res *= n * m - 1 - i
    res %= max_num
for i in range(1, r+1):
    res *= modinv(i, max_num)
    res %= max_num

print(res)

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

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

Solution (Python):



def combination(n, r, mod=10**9+7):
    
    
    n1, r = n+1, min(r, n-r)
    numer = denom = 1
    for i in range(1, r+1):
        numer = numer * (n1-i) % mod
        denom = denom * i % mod
    return numer * pow(denom, mod-2, mod) % mod

def main():
    N, M, K = map(int, input().split())
    mod = 10**9 + 7
    ans1 = 0
    for h in range(1, N):
        ans1 += (N - h) * h
    ans1 = ans1 * M * M % mod
    ans2 = 0
    for w in range(1, M):
        ans2 += (M - w) * w
    ans2 = ans2 * N * N % mod
    ans = (ans1 + ans2) * combination(N*M-2, K-2) % mod
    print(ans)

main()

Solution (Python):
N,M,K=map(int,input().split())
P=10**9+7
import math
print(M*N*(M+N)*(M*N-1)//6*math.factorial(N*M-2)*pow(math.factorial(K-2),P-2,P)*pow(math.factorial(N*M-K),P-2,P)%P)

Solution (Python):
N,M,K=map(int,input().split())
P=1000000007
import math
print(M*N*(M+N)*(M*N-1)//6*math.factorial(N*M-2)*pow(math.factorial(K-2),P-2,P)*pow(math.factorial(N*M-K),P-2,P)%P)

Solution (Python):
N,M,K=map(int,input().split())
P=1000000007
import math
ans=M**2*(N-1)*N*(N+1)//6
ans+=N**2*(M-1)*M*(M+1)//6
print(ans*math.factorial(N*M-2)*pow(math.factorial(K-2),P-2,P)*pow(math.factorial(N*M-K),P-2,P)%P)

Solution (Python):
N,M,K=map(int,input().split())
P=1000000007
import math
ans=0
for A in [N,M]:
    for i in range(1,A):
        ans = (ans+(N^M^A)**2*i*(i+1)//2)%P
print(ans*math.factorial(N*M-2)*pow(math.factorial(K-2),P-2,P)*pow(math.factorial(N*M-K),P-2,P)%P)

Solution (Python):
N,M,K=map(int,input().split())
P=1000000007
f=[0]*(N*M+1)
f[0]=1
for i in range(N*M):
    f[i+1]=(i+1)*f[i]%P
ans=0
for A in [N,M]:
    for i in range(1,A):
        ans = (ans+(N^M^A)**2*i*(i+1)//2)%P
print(ans*f[N*M-2]*pow(f[K-2],P-2,P)*pow(f[N*M-K],P-2,P)%P)

Solution (Python):
N,M,K = map(int, input().split())
P=1000000007
MPL = 1
f = 1
for i in range(1, N*M-1):
    MPL = MPL * i % P
    MPL %= P
    if i == K - 2:
        f *= pow(MPL, P - 2, P)
    if i == N*M - K:
        f *= pow(MPL, P - 2, P)
MPL = MPL * f % P
ans = 0
for A in [N, M]:
    width = N ^ M ^ A
    for i in range(1, A):
        ans = (ans + width * width * i * (i + 1) // 2) % P
print(ans * MPL % P)

Solution (Python):
from functools import reduce

def combination2(n, r, MOD=10**9+7):
    if not 0 <= r <= n: return 0
    r = min(r, n - r)
    numerator = reduce(lambda x, y: x * y % MOD, range(n, n - r, -1), 1)
    denominator = reduce(lambda x, y: x * y % MOD, range(1, r + 1), 1)
    return numerator * pow(denominator, MOD - 2, MOD) % MOD

def q(x,y):
    return cn*(x**2)*(y**3-y)//6

N, M, K = map(int, input().split())
MOD=10**9+7
cn = combination2(N*M - 2, K-2, MOD)
print((q(N,M) + q(M,N))%MOD)

Solution (Python):
N,M,K=map(int,input().split())
c=[0 for i in range(N*M)]
c[1]=1
c[0]=1
mod=10**9+7
for i in range(1,N*M-1):
  c[i+1]=(i+1)*c[i]%mod
def comb(a,b):
  return c[a]*pow(c[b],mod-2,mod)*pow(c[a-b],mod-2,mod)%mod
ans=0
for i in range(N):
  ans+=i*(N-i)*pow(M,2)
  ans%=mod
for i in range(M):
  ans+=i*(M-i)*pow(N,2)%mod
  ans%=mod
ans*=comb(N*M-2,K-2)

ans%=mod
print(ans)

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




def RD(): return input().rstrip().decode()
def II(): return int(input())
def FI(): return int(input())
def MI(): return map(int,input().split())
def MF(): return map(float,input().split())
def LI(): return list(map(int,input().split()))
def LF(): return list(map(float,input().split()))
def TI(): return tuple(map(int,input().split()))


def main():
	n,m,k=MI()
	mod=10**9+7

	f=[1]
	for i in range(1,n*m+10):
		f.append((f[-1]*i)%mod)
	

	finv=[]
	for i in f:
		finv.append(pow(i,mod-2,mod))
	


	ans=0

	cnt=0
	for i in range(1,m):
		cnt+=i*(m-i)
		cnt%=mod
	ans+=n*n*cnt

	cnt=0
	for i in range(1,n):
		cnt+=i*(n-i)
		cnt%=mod
	ans+=m*m*cnt

	ans%=mod
	

	ans*=k*(k-1)
	ans%=mod

	ans*=f[n*m-2]*finv[n*m-2-(k-2)]
	ans*=finv[k]

	print(ans%mod)


















if __name__ == "__main__":
	main()

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




def RD(): return input().rstrip().decode()
def II(): return int(input())
def FI(): return int(input())
def MI(): return map(int,input().split())
def MF(): return map(float,input().split())
def LI(): return list(map(int,input().split()))
def LF(): return list(map(float,input().split()))
def TI(): return tuple(map(int,input().split()))


def main():
	n,m,k=MI()
	mod=10**9+7

	f=[1]
	for i in range(1,n*m+10):
		f.append((f[-1]*i)%mod)
	

	finv=[]
	for i in f:
		finv.append(pow(i,mod-2,mod))
	


	ans=0

	cnt=0
	for i in range(1,m):
		cnt+=i*(m-i)
		cnt%=mod
	ans+=n*n*cnt

	cnt=0
	for i in range(1,n):
		cnt+=i*(n-i)
		cnt%=mod
	ans+=m*m*cnt

	ans%=mod
	

	ans*=k*(k-1)
	ans%=mod

	ans*=f[n*m-2]*finv[n*m-2-(k-2)]
	ans*=finv[k]

	print(ans%mod)


















if __name__ == "__main__":
	main()

Solution (Python):
def solve():
  N, M, K = map(int, input().split())
  mod = 10**9+7
  com = 1
  for i in range(1,K-1):
    com *= N*M-1-i
    com *= pow(i,mod-2,mod)
    com %= mod
  ans = com*M*N*(M*N-1)*(M+N)//6
  return ans%mod
print(solve())

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

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

MOD = 10 ** 9 + 7
N = n * m + 10
fact = [0 for _ in range(N)]
invfact = [0 for _ in range(N)]
fact[0] = 1
for i in range(1, N):
    fact[i] = fact[i - 1] * i % MOD

invfact[N - 1] = pow(fact[N - 1], MOD - 2, MOD)

for i in range(N - 2, -1, -1):
    invfact[i] = invfact[i + 1] * (i + 1) % MOD
def nCk(n, k):
    if k < 0 or n < k:
        return 0
    else:
        return fact[n] * invfact[k] * invfact[n - k] % MOD
        
ans = 0
for i in range(n-1, 0, -1):
    ans += m * m * (i + 1) * i // 2 
    ans %= MOD
for i in range(m - 1, 0, -1):
    ans += n * n * (i + 1) * i // 2
    ans %= MOD

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

Solution (Python):
N,M,K=map(int,input().split())
mod=10**9+7
X=N*M%mod
C=[1]
for i in range(K+1):
  C.append(C[i]*(X-2-i)*pow(i+1,mod-2,mod)%mod)
P=0
for i in range(N):
  P=(P+(i*(i+1)//2)+((N-i)*(N-i-1)//2))%mod
P=P*M*M%mod
Q=0
for i in range(M):
  Q=(Q+(i*(i+1)//2)+((M-i)*(M-i-1)//2))%mod
Q=Q*N*N%mod
print((P+Q)*C[K-2]*((mod+1)//2)%mod)