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


import itertools
from collections import Counter
from collections import defaultdict
import bisect
import math


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


def combination(n, r):
    if r == 0:
        return 1
    a = 1
    b = 1
    for i in range(r):
        a *= (n - i)
        if a % (i + 1):
            a //= (i + 1)
        else:
            b *= (i + 1)
    return a // b


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

    pattern = cmb(N * M - 2, K - 2)
    
    pattern = pattern % MOD
    

    ans = 0

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

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

    ans = ans % MOD
    print(ans)


if __name__ == '__main__':
    main()

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

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

def Combi2(a,b):
    if b>a:
        return 0
    ANS=1
    for i in range(min(b,a-b)):
        ANS=ANS*(a-i)*pow(min(b,a-b)-i,mod-2,mod)%mod

    return ANS%mod

SUM_N=SUM_M=0
for i in range(N):
    SUM_N+=i*(N-i)

for i in range(M):
    SUM_M+=i*(M-i)


print((SUM_N*(M**2)+SUM_M*(N**2))*pow(Combi2(N*M,2),mod-2,mod)*Combi2(N*M,K)*K*(K-1)*pow(2,mod-2,mod)%mod)

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

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

def Combi2(a,b):
    if b>a:
        return 0
    ANS=1
    for i in range(min(b,a-b)):
        ANS=ANS*(a-i)*pow(min(b,a-b)-i,mod-2,mod)%mod

    return ANS%mod

SUM_N=SUM_M=0
for i in range(N):
    SUM_N+=i*(N-i)

for i in range(M):
    SUM_M+=i*(M-i)


print((SUM_N*(M**2)+SUM_M*(N**2))*pow(Combi2(N*M,2),mod-2,mod)*Combi2(N*M,K)*K*(K-1)*pow(2,mod-2,mod)%mod)

Solution (Python):

mod = 10 ** 9 + 7
n, m, k = map(int, input().split())
nm = n * m
fact = [1] * (nm + 1)

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


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


def F(x):
    s1 = x * x * (x + 1) // 2
    s2 = x * (x + 1) * (2 * x + 1) // 6
    return (s1 - s2) % mod


res = (m ** 2 % mod * F(n) % mod + n ** 2 * F(m) % mod) % mod
res *= C(n * m - 2, k - 2)
res %= mod
print(res)

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

def comb(n, k):
    x, y = 1, 1
    k = min(k, n-k)
    for i in range(k):
        x *= n - i
        x %= mod
        y *= i + 1
        y %= mod
    return x*pow(y, mod-2, mod) % mod

N, M, K = map(int, input().split())
print((M*M%mod*((N**3-N)//6%mod)%mod+N*N%mod*((M**3-M)//6%mod)%mod)%mod*comb(N*M-2,K-2)%mod)

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

def comb(n, k):
    x, y = 1, 1
    k = min(k, n-k)
    for i in range(k):
        x *= n - i
        x %= mod
        y *= i + 1
        y %= mod
    return x*pow(y, mod-2, mod) % mod

N, M, K = map(int, input().split())
print((sum(i*(N-i)*M**2 for i in range(N))+sum(i*(M-i)*N**2 for i in range(M)))*comb(N*M-2, K-2) % mod)

Solution (Python):
from math import factorial

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

def f(N):
    a = 1
    for i in range(1, N+1):
        a *= i
        a %= mod
    return a % mod

def Comb(N, K):
    a = f(N)
    b = f(K)
    c = f(N-K)
    return (a * pow(b, mod-2, mod) * pow(c, mod-2, mod)) % mod


c = Comb(N*M-2, K-2)

ans = 0

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

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

ans *= c
ans %= mod

print(ans)

Solution (Python):
from math import factorial

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

def f(N):
    a = 1
    for i in range(1, N+1):
        a *= i
        a %= mod
    return a % mod

def Comb(N, K):
    a = f(N)
    b = f(K)
    c = f(N-K)
    return (a * pow(b, mod-2, mod) * pow(c, mod-2, mod)) % mod


c = Comb(N*M-2, K-2)

ans = 0

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

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

ans *= c
ans %= mod

print(ans)

Solution (Python):
from math import factorial

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

def Comb(N, K):
    a = factorial(N) % mod
    b = factorial(K) % mod
    c = factorial(N-K) % mod
    return (a * pow(b, mod-2, mod) * pow(c, mod-2, mod)) % mod


c = Comb(N*M-2, K-2)

ans = 0

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

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

ans *= c
ans %= mod

print(ans)

Solution (Python):
from math import factorial

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

def Comb(N, K):
    a = factorial(N) % mod
    b = factorial(K) % mod
    c = factorial(N-K) % mod
    return (a * pow(b, mod-2, mod) * pow(c, mod-2, mod)) % mod


c = Comb(N*M-2, K-2)

ans = 0

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

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

ans *= c
ans %= mod

print(ans)

Solution (Python):
n,m,k=map(int,input().split())
ans=0
INF=10**9+7
if n<=m:
    M=-(m-1)*m*(2*m-1)//3+m*m*(m-1)
    Ms=m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
    M=M%INF
    Ms=Ms%INF
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+Ms
            else:
                d=j-i
                ans=ans+d*m**2+M
            ans=ans%INF
else:
    (n,m)=(m,n)
    M=-(m-1)*m*(2*m-1)//3+m*m*(m-1)
    Ms=m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
    M=M%INF
    Ms=Ms%INF
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+Ms
            else:
                d=j-i
                ans=ans+d*m**2+M
            ans=ans%INF
def modcomb(n,k,m):
    fac=[0]*(n+1)
    finv=[0]*(n+1)
    inv=[0]*(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%m
        inv[i]=m-inv[m%i]*(m//i)%m
        finv[i]=finv[i-1]*inv[i]%m
    return fac[n]*(finv[k]*finv[n-k]%m)%m
ans=ans*modcomb(n*m-2,k-2,INF)%INF
print(ans)

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


mod = 10 ** 9 + 7
MAX = N*M
fac = [0] * MAX
finv = [0] * MAX
inv = [0] * MAX

def cominit():
    fac[0] = 1
    fac[1] = 1
    finv[0] = 1
    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 cmb(n,r):
    if n < 0 or r < 0 or r > n:return 0
    if r > n/2: r = n-r        
    return fac[n] * (finv[r] * finv[n-r] % mod) % mod
    
cominit()
ans= 0
NCK = cmb(N*M-2,K-2) 
M2 = M**2 * NCK
N2 = N**2 * NCK
for i in range(1,N):
    ans += M2 * i * (N-i)
    ans %= mod
    
for i in range(1,M):
    ans += N2 * i * (M-i)
    ans %= mod
        
print(ans)

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

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

nm = N*M
mm = M*M % mdl
nn = N*N % mdl

U = 2*10**5
fact = [1] * (U+1)

for i in range(1,U+1):
    fact[i] = (fact[i-1]*i) % mdl

def comb(n,k):
    x = fact[n]
    x *= pow(fact[n-k], mdl-2, mdl)
    x %= mdl
    x *= pow(fact[k], mdl-2, mdl)
    x %= mdl
    return x

def F(N):
    x = N*(N+1) // 2 * N
    y = N*(N+1)*(2*N+1) // 6
    return (x-y) % mdl

ans = mm*F(N) % mdl + nn*F(M) % mdl
ans %= mdl
ans *= comb(nm-2,K-2)
ans %= mdl

print(ans)

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

def pow(x, n) :
    y = 1
    while n > 0 :
        if n & 1 :
            y = (y * x) % MOD
        x = (x * x) % MOD
        n >>= 1

    return y

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

    factorial = [0] * (n + 1)
    factorial[0] = 1

    for i in range(1, n + 1) :
        factorial[i] = (factorial[i-1] * i) % MOD

    return (factorial[n] * pow(factorial[n-r], MOD - 2) * pow(factorial[r], MOD - 2)) % MOD

res = 0

for i in range(N) :
    for j in range(M) :
        if i > 0 and j > 0 :
            res += 2 * (N - i) * (M - j) * (i + j)
        else :
            res += (N - i) * (M - j) * (i + j)

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

print(res % MOD)

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

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

invs = [1] * (K-1)
for x in range(2, K-1):
    invs[x] = (-(MOD//x) * invs[MOD%x]) % MOD

C = 1
for x in range(N*M-K+1, N*M-1):
    C *= x
    C %= MOD
for x in range(1, K-1):
    C *= invs[x]
    C %= MOD

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

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

print((ans1+ans2) * C % MOD)

Solution (Python):
def com(n,r):
    re=fct[n]*pow(fct[r],md-2,md)*pow(fct[n-r],md-2,md)
    re%=md
    return re

n,m,k=map(int,input().split())
md=10**9+7
ans=0
fct=[1]*(n*m-1)
f=1
for i in range(1,n*m-1):
    f=f*i%md
    fct[i]=f
for d in range(1,m):
    ans=(ans+d*(m-d)*n*n)%md
for d in range(1,n):
    ans=(ans+d*(n-d)*m*m)%md
ans=(ans*com(n*m-2,k-2))%md
print(ans)

Solution (Python):
class PERM_COMB_MOD():
    
    def __init__(self, max_n=510000, mod=10**9+7):
        self.fac = [0]*max_n
        self.finv = [0]*max_n
        self.inv = [0]*max_n
        self.fac[0] = self.fac[1] = 1
        self.finv[0] = self.finv[1] = 1
        self.inv[1] = 1
        self.max = max_n
        self.mod = mod
        self._maesyori()

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

    def perm(self, n, k):
        if n < k : return 0
        if n < 0 or k < 0:return 0
        return self.fac[n] * self.finv[n-k] % self.mod % self.mod

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

n,m,k = map(int, input().split())
mod = 10**9+7
p = PERM_COMB_MOD()
ans = 0


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

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

Solution (Python):


from functools import reduce
from operator import add


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


def cmb(n, r):
    r = min(r, n - r)
    if r == 0:
        return 1
    if r == 1:
        return n
    f = reduce(mod_mul, range(n - r + 1, n + 1))  
    d = reduce(mod_mul, range(2, r + 1))  
    return (f * pow(d, MOD - 2, MOD)) % MOD


MOD = 10 ** 9 + 7

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

ans_x = reduce(add, (i * (M - i) for i in range(M))) * N ** 2
ans_y = reduce(add, (i * (N - i) for i in range(N))) * M ** 2
ans = (ans_x + ans_y) % MOD
pat = cmb(N * M - 2, K - 2)
ans = (ans * pat) % MOD

print(ans)

Solution (Python):


from functools import reduce


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


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


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

    f = reduce(mod_mul, range(n - r + 1, n + 1))  
    d = reduce(mod_mul, range(2, r + 1))  
    return (f * pow(d, MOD - 2, MOD)) % MOD


MOD = 10 ** 9 + 7

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

ans_x = reduce(mod_add, (i * (M - i) for i in range(M))) * N ** 2
ans_y = reduce(mod_add, (i * (N - i) for i in range(N))) * M ** 2
ans = ans_x + ans_y
pat = cmb(N * M - 2, K - 2)
ans = (ans * pat) % MOD

print(ans)

Solution (Python):


from functools import reduce


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


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

    f = reduce(mod_mul, range(n - r + 1, n + 1))  
    d = reduce(mod_mul, range(2, r + 1))  
    return (f * pow(d, MOD - 2, MOD)) % MOD


MOD = 10 ** 9 + 7

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

ans_x = 0
for i in range(M):
    ans_x += i * (M - i)
ans_x *= N ** 2

ans_y = 0
for i in range(N):
    ans_y += i * (N - i)
ans_y *= M ** 2

ans = ans_x + ans_y

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

print(ans)

Solution (Python):
"""
https://atcoder.jp/contests/abc127/submissions/5606400
"""


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


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

    f = 1
    for k in range(n, n - r, -1):
        f = (f * k) % MOD
        

    d = 1
    for k in range(2, r + 1):
        d = (d * k) % MOD
        

    return (f * inv(d)) % MOD


MOD = 10 ** 9 + 7

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

ans = 0

for i in range(M):
    ans += i * (M - i) * N ** 2

for i in range(N):
    ans += i * (N - i) * M ** 2

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

print(ans)

Solution (Python):
"""
組み合わせ(nCr)計算
https://qiita.com/derodero24/items/91b6468e66923a87f39f
逆元の求め方
http://drken1215.hatenablog.com/entry/2018/06/08/210000
解法@公式放送
https://www.youtube.com/watch?v=SS6kW-d-rJ0&feature=youtu.be&t=7376
"""


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


MOD = 10 ** 9 + 7

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

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

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

ans = 0

for i in range(M):
    ans += i * (M - i) * N ** 2

for i in range(N):
    ans += i * (N - i) * M ** 2

ans = (ans * pat) % MOD

print(ans)

Solution (Python):
from functools import reduce


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

MOD = 10**9 + 7


class ModInt:
    def __init__(self, x):
        self.x = x % MOD

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

    __repr__ = __str__

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

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

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

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

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

    __radd__ = __add__

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

    __rmul__ = __mul__

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

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


ans = (N * M * (M * (N**2 - 1) + N * (M**2 - 1)) // 6) * reduce(
    lambda acc, i: acc * (N * M - 2 - i) / (K - 2 - i),
    range(K - 2),
    ModInt(1)
)
print(ans)

Solution (Python):
N,M,K = map(int, input().split())
mod = 10**9+7
nm = N*M
n2 = (N**2)%mod
m2 = (M**2)%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

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

for i in range( 2, 10**6 +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 x in range(1,M):
    ans += x*(M-x)*n2*cmb(nm-2,K-2,mod)
    ans %= mod

for y in range(1,N):
    ans += y*(N-y)*m2*cmb(nm-2,K-2,mod)
    ans %= mod

print(int(ans))

Solution (Python):
import sys
from itertools import combinations

def comb_mod(n, r, MOD=10**9+7):
    if n < 0 or r < 0 or n < r:
        return 0

    factrial = [1] * (n+1)
    for k in range(1, n+1):
        factrial[k] = (factrial[k-1] * k) % MOD

    fact_inv = [1] * (n+1)
    fact_inv[n] = pow(factrial[n], MOD - 2, MOD)
    for k in range(n-1, -1, -1):
        fact_inv[k] = (fact_inv[k+1] * (k+1)) % MOD

    return (factrial[n] * fact_inv[r] * fact_inv[n-r]) % MOD


def main():
    input = sys.stdin.readline
    MOD = 10**9 + 7
    N, M, K = map(int, input().split())

    
    num = comb_mod(N*M-2, K-2)

    ans = 0
    
    for i in range(M):
        ans += i * (M - i) * N**2 * num
        ans %= MOD

    
    for i in range(N):
        ans += i * (N - i) * M**2 * num
        ans %= MOD

    return ans


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

Solution (Python):


from sys import stdout
from bisect import bisect_left as binl
from copy import copy, deepcopy
from collections import defaultdict
import math


mod = 1


def intin():
    input_tuple = input().split()
    if len(input_tuple) <= 1:
        return int(input_tuple[0])
    return tuple(map(int, input_tuple))


def intina():
    return [int(i) for i in input().split()]


def intinl(count):
    return [intin() for _ in range(count)]


def modadd(x, y):
    global mod
    return (x + y) % mod


def modmlt(x, y):
    global mod
    return (x * y) % mod


def lcm(x, y):
    while y != 0:
        z = x % y
        x = y
        y = z
    return x


def combination(x, y):
    assert(x >= y)

    ret = math.factorial(x)
    ret = ret // (math.factorial(x - y) * math.factorial(y))

    return ret


def get_divisors(x):
    retlist = []
    for i in range(1, int(x**0.5) + 3):
        if x % i == 0:
            retlist.append(i)
            retlist.append(x // i)
    return retlist


def get_factors(x):
    retlist = []
    for i in range(2, int(x**0.5) + 3):
        while x % i == 0:
            retlist.append(i)
            x = x // i
    retlist.append(x)
    return retlist


def make_linklist(xylist):
    linklist = {}
    for a, b in xylist:
        linklist.setdefault(a, [])
        linklist.setdefault(b, [])
        linklist[a].append(b)
        linklist[b].append(a)
    return linklist


def calc_longest_distance(linklist, v=1):
    distance_list = {}
    distance_count = 0
    distance = 0
    vlist_previous = []
    vlist = [v]
    nodecount = len(linklist)

    while distance_count < nodecount:
        vlist_next = []
        for v in vlist:
            distance_list[v] = distance
            distance_count += 1
            vlist_next.extend(linklist[v])
        distance += 1
        vlist_to_del = vlist_previous
        vlist_previous = vlist
        vlist = list(set(vlist_next) - set(vlist_to_del))

    max_distance = -1
    max_v = None
    for v, distance in distance_list.items():
        if distance > max_distance:
            max_distance = distance
            max_v = v

    return (max_distance, max_v)


def calc_tree_diameter(linklist, v=1):
    _, u = calc_longest_distance(linklist, v)
    distance, _ = calc_longest_distance(linklist, u)
    return distance


class UnionFind:

    def __init__(self, n):
        self.parent = [i for i in range(n)]

    def root(self, i):
        if self.parent[i] == i:
            return i
        self.parent[i] = self.root(self.parent[i])
        return self.parent[i]

    def unite(self, i, j):
        rooti = self.root(i)
        rootj = self.root(j)
        if rooti == rootj:
            return
        if rooti < rootj:
            self.parent[rootj] = rooti
        else:
            self.parent[rooti] = rootj

    def same(self, i, j):
        return self.root(i) == self.root(j)





mod = 10**9 + 7


def main():
    n, m, k = intin()

    c = combination(n * m - 2, k - 2) % mod

    x_ans = ((n - 1) * n * (n + 1)) // 6
    x_ans = modmlt(x_ans, m)
    x_ans = modmlt(x_ans, m)
    x_ans = modmlt(x_ans, c)

    y_ans = ((m - 1) * m * (m + 1)) // 6
    y_ans = modmlt(y_ans, n)
    y_ans = modmlt(y_ans, n)
    y_ans = modmlt(y_ans, c)

    print(modadd(x_ans, y_ans))


if __name__ == '__main__':
    main()

Solution (Python):
import sys
sys.setrecursionlimit(10**8)

MOD = 10**9 + 7
def power_mod(a,b,p):
    if b==0:
        return 1
    if b%2==0:
        d=power_mod(a,b//2,p)
        return (d*d)%p
    else:
        return(a*power_mod(a,b-1,p)%p)

def inverse(a,p):
    return power_mod(a,p-2,p)

def combination(n, r, p): 
    r = min(n-r, r)
    result = 1
    for i in range(n, n-r, -1):
        result *= i
        result %=p
    for i in range(1, r+1):
        result *= inverse(i,p)
        result %= p
    return result

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

ans = 0
for i in range(n):
    ans +=(n-i)*(m**2)*i
    ans %=MOD
for i in range(m):
    ans +=(m-i)*(n**2)*i
    ans %=MOD

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

print(ans)

Solution (Python):
N,m,k = (int(i) for i in input().split())
n,mod,x,y = N*m,10**9+7,0,0
fn,fk = [1]*n,[1]*n
for i in range(n-1): fn[i+1] = (fn[i]*(i+2))%mod
def power(n,k):
	if k==1: return n
	elif k%2==0: return power((n**2)%mod,k//2)
	else: return (n*power(n,k-1))%mod
def comb(n,k):
	if n==0 or k==0 or n==k: return 1
	else: return (((fn[n-1]*fk[n-k-1])%mod)*fk[k-1])%mod
fk[-1] = power(fn[-1],mod-2)
for i in range(2,n+1): fk[-i] = (fk[-i+1]*(n+2-i))%mod
for i in range(N): x+=i*(N-i)*(m**2)
for i in range(m): y+=i*(m-i)*(N**2)
print(comb(n-2,k-2)*(x+y)%mod)

Solution (Python):
def cmb(n,r):
    if n-r<r:
        r=n-r
    if r==0:
        return 1
    if r==1:
        return n
    
    numerator=[n-r+k+1 for k in range(r)]
    denominator=[k+1 for k in range(r)]
    
    for p in range(2,r+1):
        pivot=denominator[p-1]
        if pivot>1:
            offset=(n-r)%p
            for k in range(p-1,r,p):
                numerator[k-offset]/=pivot
                denominator[k]/=pivot
    result=1
    for k in range(r):
        if numerator[k]>1:
            result*=int(numerator[k])
    return result
mod=10**9+7
N,M,K=map(int,input().split())
retu=0
gyou=0
for i in range(1,M):
    retu+=i*(M-i)*N*N
for j in range(1,N):
    gyou+=j*(N-j)*M*M
cost=retu+gyou
num=cmb(N*M-2,K-2)%mod
ans=int((num*(cost %mod))%mod)
print(ans)

Solution (Python):


import sys


class bigint:
    primeset = set()
    def __init__(self,n,mod):
        if not mod in bigint.primeset:
            assert mod**2 < sys.maxsize, "modulo is too large"
            
            bigint.primeset.add(mod)
        
        self.n = n % mod
        self.mod = mod
        return
    
    def __str__(self):
        return str(self.n)
    
    def __add__(self,other):
        assert self.mod==other.mod, "modulos of the two not match"
        return bigint(self.n+other.n,self.mod)
    
    def __sub__(self,other):
        assert self.mod==other.mod, "modulos of the two not match"
        return bigint(self.n-other.n,self.mod)
    
    def __mul__(self,other):
        assert self.mod==other.mod, "modulos of the two not match"
        return bigint(self.n*other.n,self.mod)
    
    def __pow__(self,other):
        assert type(other)==int, "power should be integer"
        
        if other<0: return self ** ( other % ( self.mod-1 ) )
        elif other==0: return 1
        elif other==1: return self
        elif other%2==0:
            tmp= self ** (other//2)
            return bigint( ( tmp.n**2 ) % self.mod,self.mod )
        else:
            tmp= self ** (other//2)
            return bigint( ( tmp.n**2 * self.n ) % self.mod,self.mod )
    
    def __truediv__(self,other):
        assert self.mod==other.mod, "modulos of the two not match"
        return self * ( other ** (other.mod-2) )


mod= 10**9 + 7

def combination(n,m):
    global mod
    assert type(n)==type(m)==int, "type error"
    
    mm= min(m,n-m)
    if mm==0: return bigint(1,mod)
    elif mm>0:
        a = bigint(n,mod)
        b = bigint(1,mod)
        for i in range(1,mm):
            a = a * bigint(n-i,mod)
            b = b * bigint(i+1,mod)
        return a / b
    
    raise ValueError

n,m,k = [ int(x) for x in sys.stdin.readline().split() ]


nn= n * n
mm= m * m
cc= combination(n*m-2,k-2)


score=bigint(0,mod)


for i in range(1,n):
    score= score + bigint(i * mm * (n-i),mod) * cc


for i in range(1,m):
    score= score + bigint(i * nn * (m-i),mod) * cc

print( score )

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

f_list = [1] * list_size
f_r_list = [1] * list_size

for i in range(list_size - 1):
	f_list[i + 1] = int((f_list[i] * (i + 2)) % MOD)

f_r_list[-1] = pow(f_list[-1], MOD - 2, MOD)

for i in range(2, list_size + 1):
	f_r_list[-i] = int((f_r_list[-i + 1] * (list_size + 2 - i)) % MOD)

def comb(n, r):
	if n < r:
		return 0
	elif n == 0 or r == 0 or n == r:
		return 1
	else:
		return (((f_list[n - 1] * f_r_list[n - r - 1]) % MOD) * f_r_list[r - 1]) % MOD 

n, m, k = map(int, input().split())
print(((n+m)*k*(k-1))*comb(n*m, k)*f_r_list[2]%MOD)

Solution (Python):
def e_cell_distance(N, M, K, MOD=10**9 + 7):
    class Combination(object):
        """
        組み合わせ
        参考: https://harigami.net/contents?id=5f169f85-5707-4137-87a5-f0068749d9bb
        """
        __slots__ = ["mod", "factorial", "inverse"]

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

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

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

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

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

        def combination(self, n, r):
            if r == n or r == 0:
                return 1
            if r > n:
                return 0
            return self.factorial[n] * self.inverse[r] * self.inverse[n - r] % self.mod

    comb = Combination(N * M).combination(N * M - 2, K - 2)
    t1 = sum([d * (M - d) * (N**2) for d in range(1, M)])
    t2 = sum([d * (N - d) * (M**2) for d in range(1, N)])
    t = (t1 + t2) % MOD
    return (comb * t) % MOD

N, M, K = [int(i) for i in input().split()]
print(e_cell_distance(N, M, K))

Solution (Python):
def power_mod(a,b,mod=10**9+7):
    i,temp,box=0,b,[]
    while(2**i<=b):
        i+=1
    for j in range(i-1,-1,-1):
        box=[[j,temp//2**j]]+box
        temp-=2**j*(temp//2**j)
    box[0].append(a)
    ans=box[0][1]*a%mod
    for j in range(1,i):
        box[j].append(box[j-1][2]**2%mod)
        if box[j][1]==1:
            ans=(ans*box[j][2])%mod
    return ans
def n_func(n,mod=10**9+7):
    ans=1
    for i in range(1,n+1):
        ans=(ans*i)%mod
    return ans
def nPr(n,r,mod=10**9+7):
    ans=n_func(n-r,mod)
    ans=power_mod(ans,mod-2,mod)
    return ans*n_func(n,mod)%mod
def nCr(n,r,mod=10**9+7):
    ans=n_func(n-r,mod)*n_func(r,mod)%mod
    ans=power_mod(ans,mod-2,mod)
    return ans*n_func(n,mod)%mod

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

sum_=0
for i in range(N):
    a=min(abs(i),abs(N-i-1))
    b=max(abs(i),abs(N-i-1))
    sum_+=(M**2)*((a*(a+1)//2)+(b*(b+1)//2))
for i in range(M):
    a=min(abs(i),abs(M-i-1))
    b=max(abs(i),abs(M-i-1))
    sum_+=(N**2)*((a*(a+1)//2)+(b*(b+1)//2))
print((keisuu * (sum_//2))%mod)

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

f_list = [1] * list_size
f_r_list = [1] * list_size

for i in range(list_size - 1):
	f_list[i + 1] = int((f_list[i] * (i + 2)) % MOD)

f_r_list[-1] = pow(f_list[-1], MOD - 2, MOD)

for i in range(2, list_size + 1):
	f_r_list[-i] = int((f_r_list[-i + 1] * (list_size + 2 - i)) % MOD)

def comb(n, r):
	if n < r:
		return 0
	elif n == 0 or r == 0 or n == r:
		return 1
	else:
		return (((f_list[n - 1] * f_r_list[n - r - 1]) % MOD) * f_r_list[r - 1]) % MOD 

n, m, k = map(int, input().split())
print((((n+m)*n*m*(m*n-1))//6)*comb(n*m-2, k-2)%MOD)

Solution (Python):
class Combination:

    def __init__(self, n, MOD):
        self.fact = [1]
        for i in range(1, n + 1):
            self.fact.append(self.fact[-1] * i % MOD)
        self.inv_fact = [pow(self.fact[i] ,MOD - 2 ,MOD) for i in range(n + 1)]
        self.MOD = MOD

    def factorial(self, k):
        '''k!を求める'''
        return self.fact[k]

    def inverse_factorial(self, k):
        '''k!の逆元を求める'''
        return self.inv_fact[k]

    def combination(self, k, r):
        '''kCrを求める'''
        return (self.fact[k] * self.inv_fact[k - r] * self.inv_fact[r]) % self.MOD
      
h, w, k = map(int, input().split())
MOD = 10**9 + 7
comb = Combination(h*w, MOD)

kumi = comb.combination(h*w - 2, k - 2)
ans = 0

for d in range(0, h):
    ans += d * (h - d) * kumi * w * w
    ans %= MOD

for d in range(0, w):
    ans += d * (w - d) * kumi * h * h
    ans %= MOD

print(ans)