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):
mod = 10 ** 9 + 7


class Combination:
    """
    SIZEが10^6程度以下の二項係数を何回も呼び出したいときに使う
    使い方:
    comb = Combination(SIZE, MOD)
    comb(10, 3) => 120
    """

    def __init__(self, N, MOD=10 ** 9 + 7):
        self.MOD = MOD
        self.fact, self.inv = self._make_factorial_list(N)

    def __call__(self, n, k):
        if k < 0 or k > n:
            return 0
        res = self.fact[n] * self.inv[k] % self.MOD
        res = res * self.inv[n - k] % self.MOD
        return res

    def _make_factorial_list(self, N):
        fact = [1] * (N + 1)
        inv = [1] * (N + 1)
        MOD = self.MOD
        for i in range(1, N + 1):
            fact[i] = (fact[i - 1] * i) % MOD
        inv[N] = pow(fact[N], MOD - 2, MOD)
        for i in range(N, 0, -1):
            inv[i - 1] = (inv[i] * i) % MOD
        return fact, inv


if __name__ == "__main__":
    N, M, K = map(int, input().split())
    comb = Combination(N * M + 10, mod)

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

    for d in range(1, M):
        tmp = (M - d) * N * N % mod
        ans += d * tmp * comb(N * M - 2, K - 2)
        ans %= mod

    print(ans)

Solution (Python):
def examA():
    C1 = SI()
    C2 = SI()
    ans = "YES"
    if C1[0]!=C2[2] or C1[1]!=C2[1] or C1[2]!=C2[0]:
        ans = "NO"
    print(ans)
    return

def examB():
    A, B, C, X, Y = LI()
    loop = max(X,Y)
    ans = inf
    for i in range(loop+1):
        cost = i*2*C
        if X>i:
            cost += A*(X-i)
        if Y>i:
            cost += B*(Y-i)
        ans = min(ans,cost)
    print(ans)
    return

def examC():
    def gcd(x, y):
        if y == 0:
            return x
        while (y != 0):
            x, y = y, x % y
        return x
    N, x = LI()
    X = LI()
    if N==1:
        print(abs(X[0]-x))
        return
    for i in range(N):
        X[i] = abs(X[i]-x)
    now = gcd(X[0],X[1])
    for i in range(2,N):
        now = gcd(now,X[i])
    ans = now
    print(ans)
    return

def examD():
    N = I()
    C = [LI()for _ in range(N-1)]
    ans = [inf]*N
    for i in range(N):
        cur = 0
        for j in range(i,N-1):
            c, s, f = C[j]
            next = 0
            if cur<=s:
                cur = s+c
            else:
                cur = (1+(cur-1)//f)*f + c
        ans[i] = cur
    for v in ans:
        print(v)
    return

def examE():
    N, C = LI()
    S = defaultdict(list)
    S = [LI()for _ in range(N)]
    ans = C


    print(ans)
    return

def examF():
    class combination():
        
        def __init__(self, n, mod):
            self.n = n
            self.mod = mod
            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):
            if r > n or n < 0 or r < 0:
                return 0
            return self.fac[n] * self.inv[n - r] * self.inv[r] % self.mod
    N, M, K = LI()
    B = N*M
    C = combination(B,mod)
    ans = 0
    for i in range(N):
        for j in range(M):
            if i>0 and j>0:
                ans += (i+j)*(N - i) * (M - j) * C.comb(B - 2, K-2) *2
            else:
                ans += (i+j)*(N - i) * (M - j) * C.comb(B - 2, K-2)
            ans %= mod
            
    print(ans)
    return

from decimal import Decimal as dec
import sys,bisect,itertools,heapq,math,random
from copy import deepcopy
from heapq import heappop,heappush,heapify
from collections import Counter,defaultdict,deque
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines
def I(): return int(input())
def LI(): return list(map(int,sys.stdin.readline().split()))
def DI(): return dec(input())
def LDI(): return list(map(dec,sys.stdin.readline().split()))
def LSI(): return list(map(str,sys.stdin.readline().split()))
def LS(): return sys.stdin.readline().split()
def SI(): return sys.stdin.readline().strip()
global mod,mod2,inf,alphabet,_ep
mod = 10**9 + 7
mod2 = 998244353
inf = 10**18
_ep = 10**(-12)
alphabet = [chr(ord('a') + i) for i in range(26)]

sys.setrecursionlimit(10**7)

if __name__ == '__main__':
    examF()

"""
142
12 9 1445 0 1
asd dfg hj o o
aidn
"""

Solution (Python):
mod = 10 ** 9 + 7
kkai = [1]
for i in range(1, 3 * 10 ** 5):
    kkai.append(kkai[-1] * i % mod)


def kai(x, p=mod):
    return kkai[x]


def comb(a, b, p=mod):
    if a < 0 or b < 0:
        return 0
    elif a < b:
        return 0
    c = 1
    c *= kai(a, p)
    c *= pow(kai(b, p), p - 2, p)
    c *= pow(kai(a - b, p), p - 2, p)
    return c % p


n, m, k = [int(i) for i in input().split()]

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

Solution (Python):



def comb(n,k,mod):
    if n<k:
        return 0
    if n<0 or k<0:
        return 0
    k=min(n-k,k)
    ans=1
    inv=[1]*(k+1)
    if k>=1:
        ans*=(n-k+1)%mod
    for i in range(2,k+1):
        inv[i]=mod-inv[mod%i]*(mod//i)%mod
        ans=ans*(n-k+i)*inv[i]%mod
    return ans

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

MOD = 1000000007

distX = 0
for d in range(n):
    distX += d * (n - d) * m ** 2

distY = 0
for d in range(m):
    distY += d * (m - d) * n ** 2

ans = ((distX + distY) * comb(n * m - 2, k - 2, MOD)) % MOD

print(ans)

Solution (Python):
MOD = 10 ** 9 + 7
n, m, k = map(int, input().split())
def f(i, j):
    return i*(i-1)*(i+1)//6 * j ** 2
def c(a, b):
    r = 1
    for i in range(a, a - b, -1):
        r = r * i % MOD
    t = 1
    for i in range(1, b + 1):
        t = t * i % MOD
    return r * pow(t, MOD - 2, MOD)
print((f(n, m) + f(m, n)) * c(n*m - 2, k - 2) % MOD)

Solution (Python):
import sys
input = lambda : sys.stdin.readline().rstrip()
mod = 10**9 + 7

def make_tables(n, mod = 10 ** 9 + 7):
    fac = [1, 1] 
    finv = [1, 1] 
    inv = [0, 1] 

    for i in range(2, n + 1):
        fac.append((fac[-1] * i) % mod)
        inv.append((mod -inv[mod % i] * (mod // i)) % mod)
        finv.append((finv[-1] * inv[-1]) % mod)
    return fac, finv

def nCk(n, k, mod = 10 ** 9 + 7):
    k = min(k, n-k)
    return fac[n] * finv[k] * finv[n-k] % mod



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

fac, finv = make_tables(n*m - 2)

ans = 0

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

print(ans)

Solution (Python):

import sys


buff_readline = sys.stdin.readline
readline = sys.stdin.readline

INF = 2**62-1


def read_int():
    return int(buff_readline())


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


def read_float():
    return float(buff_readline())


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


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


def read_str_n():
    return readline().strip().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 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):
    mod = 10**9+7
    C = Combination(N*M, mod)
    b = C(N*M-2, K-2)
    ans = 0
    m = M**2
    for d in range(1, N):
        ans += (N-d) * b * d * m
        ans %= mod
    m = N**2
    for d in range(1, M):
        ans += (M-d) * b * d * m
        ans %= mod



    return ans



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


if __name__ == '__main__':
    main()

Solution (Python):
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines
from math import factorial 
def main():
    n, m, k = map(int, readline().split())
    P = int(1e9)+7

    def C(n, m, p):
        a = 1
        for i in range(m):
            a *= (n-i)
            a %= p
        b = factorial(m)%p 
        b = pow(b, p-2, p)
        return a*b%p

    c = C(n*m-2, k-2, P)
    s_x = 0
    for d in range(1, m):
        s_x += d*(m-d)*n*n
        s_x %= P
    s_y = 0
    for d in range(1, n):
        s_y += d*(n-d)*m*m
        s_y %= P

    ans = (s_x+s_y)*c%P
    print(ans)
if __name__ == '__main__':
    main()

Solution (Python):
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines
from math import factorial 
def main():
    n, m, k = map(int, readline().split())
    P = int(1e9)+7

    def C(n, m, p):
        a = 1
        for i in range(m):
            a *= (n-i)
            a %= p
        b = factorial(m)%p 
        b = pow(b, p-2, p)
        return a*b%p

    c = C(n*m-2, k-2, P)
    s_x = 0
    for d in range(1, m):
        s_x += d*(m-d)*n*n
        s_x %= P
    s_y = 0
    for d in range(1, n):
        s_y += d*(n-d)*m*m
        s_y %= P

    ans = (s_x+s_y)*c%P
    print(ans)
if __name__ == '__main__':
    main()

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


def main() -> None:
    N, M, K = map(int, input().split())
    MOD = 10 ** 9 + 7
    ans = 0
    for d in range(1, N):
        ans += M * M * (N - d) * d
        ans %= MOD
    for d in range(1, M):
        ans += N * N * (M - d) * d
        ans %= MOD

    for i in range(K - 2):
        ans = (ans * (N * M - 2 - i)) % MOD
    for i in range(K - 2):
        ans = (ans * inv(i + 1, MOD)) % MOD
    print(ans)


if __name__ == '__main__':
    main()

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


mod = 10 ** 9 + 7
N = 10 ** 6


inv_t = [0]+[1]
for i in range(2, N):
  inv_t += [inv_t[mod % i] * (mod - int(mod / i)) % mod]


kai = [1, 1]
rev_kai = [1, inv_t[1]]
for i in range(2, N):
	kai.append(kai[-1] * i % mod)
	rev_kai.append(rev_kai[-1] * inv_t[i] % mod)


def cmb(n,r):
    res = 1
    for i in range(r):
        res *= n * inv_t[r]
        res %= mod
        n -= 1
        r -= 1
    return res

def calc(n):
    res = 0
    for i in range(1, n):
        res += (n - i) * i
        res %= mod
    return res

ans = 0
ans += m * m * cmb(n*m -2, k-2) * calc(n)
ans %= mod

ans += n * n * cmb(n*m -2, k-2) * calc(m)
ans %= mod

print(ans)

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

def MOD_perm(n, r):
    rtn = 1
    for _ in range(r):
        rtn *= n
        rtn %= MOD
        n -= 1
    return rtn

def MOD_inv(a):
    b = MOD
    u = 1
    v = 0
    while b > 0:
        t = a // b
        a -= t * b
        a, b = b, a
        u -= t * v
        u, v = v, u
    u = u % MOD
    if u < 0:
        u += MOD
    return u

def MOD_COM2(n, r):
    return MOD_perm(n, r) * MOD_inv(MOD_perm(r, r)) % MOD

def main():
    n, m, k = map(int, input().split())
    comb = MOD_COM2(n*m-2, k-2)
    ans = 0
    for i in range(1, n+1):
        U = (i-1)*i//2
        U *= m*m
        U *= comb
        U %= MOD
        ans += U
        ans %= MOD
    for j in range(1, m+1):
        L = (j-1)*j//2
        L *= n*n
        L *= comb
        L %= MOD
        ans += L
        ans %= MOD
    print(ans%MOD)

if __name__ == "__main__":
    main()

Solution (Python):


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


N, M, K = LI()

MOD = 10 ** 9 + 7
MAX = 10 ** 5 * 2
fac, finv, inv = [None] * MAX, [None] * MAX, [None] * MAX


def comb_init():
    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] = -inv[MOD%i] * int(MOD / i) % MOD
        finv[i] = finv[i - 1] * inv[i] % MOD


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


def s(x):
    return x * (x + 1) // 2


def cal(a, b):
    ret = a * s(b - 1) + b * s(a - 1)
    return ret % MOD


def main():
    comb_init()
    ans = 0
    for i in range(N):
        for j in range(M):
            v = cal(N - i, M - j) + cal(N - i, j + 1) - (s(j) + s(N - i - 1))
            ans = (ans + v) % MOD

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


if __name__ == "__main__":
    main()

Solution (Python):
N,M,K=map(int,input().split())
MOD=10**9+7
ans=0
x,y=1,1
k=min(K-2,N*M-K)
for i in range(k):
    x*=N*M-2-i
    y*=k-i
    x%=MOD
    y%=MOD
comb=x*pow(y,MOD-2,MOD)
comb%=MOD
for i in range(1,N):
    a=i*comb*(N-i)*M**2
    ans+=a
    ans%=MOD
for i in range(1,M):
    a=i*comb*(M-i)*N**2
    ans+=a
    ans%=MOD
print(ans%MOD)

Solution (Python):
MOD=10**9+7
UPPERLIMIT=3*10**5
MODMUL=[1, 1]+[0]*(UPPERLIMIT-1)
for i in range(2, UPPERLIMIT+1):
  MODMUL[i]=MODMUL[i-1]*i%MOD
MODDIV=[1]*UPPERLIMIT+[pow(MODMUL[-1], MOD-2, MOD)]
for i in range(UPPERLIMIT, 0, -1):
  MODDIV[i-1]=MODDIV[i]*i%MOD
def MODCOMB(n, r):
  if n<r or n<1:
    return 0
  else:
    return (((MODMUL[n]*MODDIV[r])%MOD)*MODDIV[n-r])%MOD

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

ans=0
k=MODCOMB(N*M-2, K-2)
for sa in range(1, N):
  ans=(ans+sa*(N-sa)*M*M*k)%MOD
    
for sa in range(1, M):
  ans=(ans+sa*(M-sa)*N*N*k)%MOD
    
print(ans)

Solution (Python):
def main():
    import sys
    input = sys.stdin.readline
    sys.setrecursionlimit(10**7)
    from collections import Counter, deque
    from itertools import combinations, permutations, accumulate, groupby, product
    from bisect import bisect_left,bisect_right
    from heapq import heapify, heappop, heappush
    import math
    

    
    mod = 10**9 + 7

    
    max_n = 4*10**5
    fac, inv = [1]*(max_n+1), [0]*(max_n+1)
    for i in range(2, max_n+1):
        fac[i] = fac[i-1] * i % mod
    inv[-1] = pow(fac[-1], mod-2, mod)
    for i in range(max_n, 0, -1):
        inv[i-1] = inv[i] * i % mod
    
    
    def ncr(n, r):
        return fac[n]*inv[r]*inv[n-r]%mod

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

    base = ncr(n*m-2, k-2)

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

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

    print(res*base%mod)

if __name__ == '__main__':
    main()

Solution (Python):
def main():
    import sys
    input = sys.stdin.readline
    sys.setrecursionlimit(10**7)
    from collections import Counter, deque
    from itertools import combinations, permutations, accumulate, groupby, product
    from bisect import bisect_left,bisect_right
    from heapq import heapify, heappop, heappush
    import math
    

    
    mod = 10**9 + 7

    
    max_n = 4*10**5
    fac, inv = [1]*(max_n+1), [0]*(max_n+1)
    for i in range(2, max_n+1):
        fac[i] = fac[i-1] * i % mod
    inv[-1] = pow(fac[-1], mod-2, mod)
    for i in range(max_n, 0, -1):
        inv[i-1] = inv[i] * i % mod
    
    
    def ncr(n, r):
        return fac[n]*inv[r]*inv[n-r]%mod

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

    base = ncr(n*m-2, k-2)

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

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

    print(res*base%mod)

if __name__ == '__main__':
    main()

Solution (Python):
def combination_mod(n,r,MOD=10**9+7):
    n_=1
    for i in range(1,n+1):
        n_=(n_*i)%MOD
    r_=1
    for i in range(1,r+1):
        r_=(r_*i)%MOD
    nr_=1
    for i in range(1,n-r+1):
        nr_=(nr_*i)%MOD

    power_r=pow(r_,MOD-2,MOD)
    power_nr=pow(nr_,MOD-2,MOD)

    return (n_*power_r*power_nr)%MOD


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

    MOD=10**9+7
    num=combination_mod(N*M-2,K-2)

    col_ans=0
    row_ans=0

    for i in range(1,M):
        col_ans+=pow(N,2)*(M-i)*num*i

    for i in range(1,N):
        row_ans+=pow(M,2)*(N-i)*num*i

    ans=(col_ans+row_ans)%MOD

    print(ans)

if __name__=="__main__":
    main()

Solution (Python):
from operator import mul
from functools import reduce
from fractions import gcd
import math
import bisect
import itertools
import sys
sys.setrecursionlimit(10**7)
input = sys.stdin.readline
INF = float("inf")

MOD = 1000000007
class Combination(object):
    def __init__(self, N, MOD):
        self.fac = [0] * (N + 1)
        self.inv = [0] * (N + 1)
        self.finv = [0] * (N + 1)
        self.fac[0] = 1
        self.finv[0] = 1
        if N > 0:
            self.fac[1] = 1
            self.inv[1] = 1
            self.finv[1] = 1
 
        
        for i in range(2, N + 1):
            self.fac[i] = self.fac[i - 1] * i % MOD
            self.inv[i] = self.inv[MOD % i] * (MOD - (MOD // i)) % MOD
            self.finv[i] = self.finv[i - 1] * self.inv[i] % MOD
 
    def com(self, N, k):
        return (self.fac[N] * self.finv[k] * self.finv[N - k]) % MOD


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

    ans = 0
    Com = Combination(N*M-2, MOD)
    a = Com.com(N*M-2, K-2)
    for dx in range(1, M):
        b = (M-dx) * N * N
        ans += a * b * dx
        ans %= MOD
    for dy in range(1, N):
        b = (N-dy) * M * M
        ans += a * b * dy
        ans %= MOD
    
    print(ans)


if __name__ == '__main__':
    main()

Solution (Python):
def comb(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())
MOD = 10**9+7
sumx = 0
sumy = 0

for d in range(1, N):
    cnt = (N-d)*M**2
    cnt %= MOD
    cnt *= comb(N*M-2, K-2, MOD)
    sumx += d*cnt
    sumx %= MOD

for d in range(1, M):
    cnt = (M-d)*N**2
    cnt %= MOD
    cnt *= comb(N*M-2, K-2, MOD)
    sumy += d*cnt
    sumy %= MOD
print((sumx + sumy)%MOD)

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

fac = [1] * n*m
finv = [1] * n*m
inv = [1] * n*m

ans = 0
def COMinit():
    for i in range(2, n*m):
        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, k):
    if n < 0 or k < 0 or n < k:
        return 0
    return fac[n] * (finv[k] * finv[n - k] % mod) % mod
COMinit()
for i in range(1, n):
    ans += (COM(n*m - 2, k-2)) * ((n-i)*i * m**2)%mod
for j in range(1, m):
    ans += (COM(n*m - 2, k-2)) * ((m-j)*j * n**2)%mod
print(ans%mod)

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

def make_array_for_comb(N, mod=10**9+7):
    fact = [1,1]
    fact_inv = [1,1]
    inv = [0,1]
    for i in range(2, N+1):
        fact.append((fact[-1]*i) % mod)
        inv.append((-inv[mod%i] * (mod//i)) % mod)
        fact_inv.append((fact_inv[-1]*inv[i]) % mod)
    return fact, fact_inv

def comb(n, r, mod=10**9+7):
  if (r < 0) or (n < r):
      return 0
  r = min(r, n - r)
  return fact[n] * fact_inv[r] * fact_inv[n-r] % mod

fact, fact_inv = make_array_for_comb(n*m,mod)

l = comb(n*m-2,k-2,mod)
ans = 0
for a in range(1,n+1):
  ans += a*m*(a-1)*m
  ans -= a*m*(n-a)*m
ans %= mod

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

print(ans*l%mod)

Solution (Python):
import sys

readline = sys.stdin.readline
MOD = 10 ** 9 + 7
INF = float('INF')
sys.setrecursionlimit(10 ** 5)


def comb_mod(n, r):
    res = 1
    r = min(n - r, r)

    for i in range(r):
        res *= (n - i)
        res %= MOD
        res *= pow((r - i), MOD - 2, MOD)

    return res


def main():
    n, m, k = list(map(int, readline().split()))

    com = comb_mod(n * m - 2, k - 2)
    x = 0
    y = 0

    for i in range(n):
        right = n - 1 - i
        x += (right + 1) * right // 2
        x %= MOD

    for i in range(m):
        down = m - 1 - i
        y += (down + 1) * down // 2
        y %= MOD

    ans = 0
    ans += m ** 2 * x * com
    ans += n ** 2 * y * com
    ans %= MOD

    print(ans)


if __name__ == '__main__':
    main()

Solution (Python):
import sys
sys.setrecursionlimit(2147483647)
INF = 1 << 60
MOD = 10**9 + 7 
input = lambda:sys.stdin.readline().rstrip()

class modfact(object):
    def __init__(self, n):
        fact, invfact = [1] * (n + 1), [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, self._invfact = fact, invfact

    def inv(self, n):
        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] % MOD * 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():
    m, n, k = map(int, input().split())
    mf = modfact(m * n - 2)
    res = 0
    for s in range(1, m):
        res += n * n * (m - s) * s
        res %= MOD
    for t in range(1, n):
        res += m * m * (n - t) * t
        res %= MOD
    res *= mf.comb(m * n - 2, k - 2)
    res %= MOD
    print(res)
resolve()

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

mod = 10**9 + 7
n = 2*10**5

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

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

base = comb(N*M-2,K-2)
first = 0
for i in range(N):
  for j in range(M):
    first += (i+j)
dp = [0]*M
dp[0] = first
for i in range(1,M):
  dp[i] = dp[i-1] - (M-2*i) * N
SUM = sum(dp)
dp2 = [0]*N
dp2[0] = SUM
for i in range(1,N):
  dp2[i] = dp2[i-1] - (N-2*i) * (M ** 2)
SUM = sum(dp2)
ans = comb(N*M-2,K-2) * (SUM % mod) * pow(2,mod-2,mod)
ans %= mod
print(ans)

Solution (Python):
import sys
def input(): return sys.stdin.readline().strip()
mod = 10**9+7

class Combination:
    """
    O(n)の前計算を1回行うことで，O(1)でnCr mod mを求められる
    n_max = 10**6のとき前処理は約950ms (PyPyなら約340ms, 10**7で約1800ms)
    使用例：
    comb = Combination(1000000)
    print(comb(5, 3))  
    """
    def __init__(self, n_max, mod=10**9+7):
        self.mod = mod
        self.modinv = self.make_modinv_list(n_max)
        self.fac, self.facinv = self.make_factorial_list(n_max)

    def __call__(self, n, r):
        return self.fac[n] * self.facinv[r] % self.mod * self.facinv[n-r] % self.mod

    def make_factorial_list(self, n):
        
        
        fac = [1]
        facinv = [1]
        for i in range(1, n+1):
            fac.append(fac[i-1] * i % self.mod)
            facinv.append(facinv[i-1] * self.modinv[i] % self.mod)
        return fac, facinv

    def make_modinv_list(self, n):
        
        modinv = [0] * (n+1)
        modinv[1] = 1
        for i in range(2, n+1):
            modinv[i] = self.mod - self.mod//i * modinv[self.mod%i] % self.mod
        return modinv

comb = Combination(2 * 100000)

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

    """
    xについてだけ考えればいいのはOK
    この後だが、xi, xjを固定したときに\Sum_{置き方} |xi - xj|を求めればよい。
    （シグマの順番を入れ替えた。つまりxiとxjに駒が置かれるような配置が何通りあるかを求めたい。）
    （iとjに騙されないこと！xi, xjと書くんじゃなくてx, x'と書いた方が分かりやすいかもしれない。）

    これは単純にcomb(NM-2, K-2)なので一応これで解けた。すなわち
        ans = \Sum_{xi, xj} |xi - xj| * comb(NM-2, K-2)

    しかしi, jに関してそれぞれforループを回しているのでこれだと間に合わない。
    まずcomb(NM-2, K-2)はi, jに依存しないので二重forループの外に出せる。
    そして\Sum_{i<=j}|xi-xj|だが、|xi-xj|=dとなる(i, j)の組は(N-d)*M^2個になる。（解説放送参照）
    """

    ans_x = 0
    for d in range(1, N):
        ans_x += d * (N - d) * M ** 2
        ans_x %= mod
    ans_y = 0
    for d in range(1, M):
        ans_y += d * (M - d) * N ** 2
        ans_y %= mod
    ans = (ans_x + ans_y) * comb(N * M - 2, K - 2)
    print(ans % mod)




if __name__ == "__main__":
    main()

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


def kai(x, p=mod):
    a = 1
    for i in range(1, x + 1):
        a *= i
        a %= p
    return a


def comb(a, b, p=mod):
    if a < 0 or b < 0:
        return 0
    elif a < b:
        return 0
    c = 1
    c *= kai(a, p)
    c *= pow(kai(b, p), p - 2, p)
    c *= pow(kai(a - b, p), p - 2, p)
    return c % p


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
ans *= comb(n * m - 2, k - 2)
ans %= mod
print(ans)

Solution (Python):
class Combination:
    """
    O(n)の前計算を1回行うことで，O(1)でnCr mod mを求められる
    n_max = 10**6のとき前処理は約950ms (PyPyなら約340ms, 10**7で約1800ms)
    使用例：
    comb = Combination(1000000)
    print(comb(5, 3))  
    """
    def __init__(self, n_max, mod=10**9+7):
        self.mod = mod
        self.modinv = self.make_modinv_list(n_max)
        self.fac, self.facinv = self.make_factorial_list(n_max)

    def __call__(self, n, r):
        if n < r:
            return 0
        return self.fac[n] * self.facinv[r] % self.mod * self.facinv[n-r] % self.mod

    def make_factorial_list(self, n):
        
        
        fac = [1]
        facinv = [1]
        for i in range(1, n+1):
            fac.append(fac[i-1] * i % self.mod)
            facinv.append(facinv[i-1] * self.modinv[i] % self.mod)
        return fac, facinv

    def make_modinv_list(self, n):
        
        modinv = [0] * (n+1)
        modinv[1] = 1
        for i in range(2, n+1):
            modinv[i] = self.mod - self.mod//i * modinv[self.mod%i] % self.mod
        return modinv

def main():
    N, M, K = map(int, input().split())
    MOD = 10 ** 9 + 7
    comb = Combination(N*M)
    ans = 0
    for i in range(M):
        ans += i * (M - i) * N * N
        ans %= MOD
    
    for i in range(N):
        ans += i * (N - i) * M * M
        ans %= MOD

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



if __name__ == "__main__":
    main()

Solution (Python):
import bisect
import copy
import heapq
import math
import sys
from collections import *
from itertools import accumulate, combinations, permutations, product

def input():
    return sys.stdin.readline()[:-1]
def ruiseki(lst):
    return [0]+list(accumulate(lst))
mod=pow(10,9)+7
al=[chr(ord('a') + i) for i in range(26)]
direction=[[1,0],[0,1],[-1,0],[0,-1]]

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


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

N = 2*10**5+10
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 )

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

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

print(ans)

Solution (Python):

import sys
sys.setrecursionlimit(10**9)
INF=10**18
MOD=10**9+7
input=lambda: sys.stdin.readline().rstrip()
YesNo=lambda b: bool([print('Yes')] if b else print('No'))
YESNO=lambda b: bool([print('YES')] if b else print('NO'))
int1=lambda x:int(x)-1

def main():
    N,M,K=map(int,input().split())
    
    def COMinit(n,MOD):
        fac,finv,inv=[0]*max(2,n+1),[0]*max(2,n+1),[0]*max(2,n+1)
        fac[0]=fac[1]=1
        finv[0]=finv[1]=1
        inv[1]=1
        for i in range(2,(n+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
        return fac,finv,inv
    
    fac,finv,inv=COMinit(N*M,MOD)
    
    def COM(n, k, MOD):
        if n<k or n<0 or k<0:
            return 0
        return fac[n]*(finv[k]*finv[n-k]%MOD)%MOD
    
    ans=0
    for a,b in [(N,M),(M,N)]:
        for i in range(1,a):
            ans+=(i*(a-i)*b**2)
            ans%=MOD
    ans*=COM(N*M-2,K-2,MOD)
    ans%=MOD
    print(ans)

if __name__ == '__main__':
    main()

Solution (Python):
MAX = 2*10**5+100
MOD = 10**9+7
fact = [0]*MAX 
inv = [0]*MAX 
finv = [0]*MAX 
fact[0] = 1
fact[1] = 1
finv[0] = 1
finv[1] = 1
inv[1] = 1
    
for i in range(2, MAX):
    fact[i] = fact[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 fact[n]*(finv[r]*finv[n-r]%MOD)%MOD

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

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

ans *= C(N*M-2, K-2)
ans %= MOD

print(ans)

Solution (Python):
N,M,K=map(int,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)
ans=0
for i in range(1,N):
  ans+=i*(N-i)*M*M
for j in range(1,M):
  ans+=j*(M-j)*N*N
ans*=cmb(N*M-2,K-2,mod)
print(ans%mod)