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.



Sample input:
2 2 2


Sample output:
8

Submission status: Accepted

Solution code:
MAX = 2*10**5+100
MOD = 10**9+7
fact = [0]*MAX #fact[i]: i!
inv = [0]*MAX #inv[i]: iの逆元
finv = [0]*MAX #finv[i]: i!の逆元
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 i in range(1, N):
    ans += i*(N-i)*M*M*C(N*M-2, K-2)
    ans %= MOD

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

print(ans)

Submission status: Accepted

Solution code:
import sys
readline = sys.stdin.readline
sys.setrecursionlimit(10**8)
mod = 10**9+7
#mod = 998244353
INF = 10**18
eps = 10**-7

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

def comb(n, r, mod):
    r = min(r, n-r)
    mol = 1
    deno = 1
    for i in range(1, r+1):
        mol = mol * (n-r+i) % mod
        deno = deno * i % mod
    ret = mol * pow(deno, mod-2, mod) % mod
    return ret

def f1(x,y,a,b):
    return (a+b-x-y)*(a-x+1)*(b-y+1)//2

def f2(x,a):
    return (a-x)*(a-x+1)//2

ans = 0

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

ans = (ans * comb(m*n-2,k-2,mod))%mod
ans = ans*pow(2,mod-2,mod)%mod

print(ans)

Submission status: Accepted

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

# cmbが10**10くらいだけど求められるか？って感じ
# 問題読み違えていた。。。N*M <= 2*10**5だ。。。まあ普通だ。
# もうライブラリ使おう。昔nCrのrが小さい時の工夫とかあったけど今回は大丈夫だ。

# https://ikatakos.com/pot/programming_algorithm/number_theory/mod_combination
import numpy as np
 
def prepare(n, MOD):
    nrt = int(n ** 0.5) + 1
    nsq = nrt * nrt
    facts = np.arange(nsq, dtype=np.int64).reshape(nrt, nrt)
    facts[0, 0] = 1
    for i in range(1, nrt):
        facts[:, i] = facts[:, i] * facts[:, i - 1] % MOD
    for i in range(1, nrt):
        facts[i] = facts[i] * facts[i - 1, -1] % MOD
    facts = facts.ravel().tolist()
 
    invs = np.arange(1, nsq + 1, dtype=np.int64).reshape(nrt, nrt)
    invs[-1, -1] = pow(facts[-1], MOD - 2, MOD)
    for i in range(nrt - 2, -1, -1):
        invs[:, i] = invs[:, i] * invs[:, i + 1] % MOD
    for i in range(nrt - 2, -1, -1):
        invs[i] = invs[i] * invs[i + 1, 0] % MOD
    invs = invs.ravel().tolist()
 
    return facts, invs

facts, invs = prepare(N*M+10,mod)
def cmb(n,r,MOD):
  return (((facts[n] * invs[n-r]) % mod) * invs[r]) % mod

## シグマを分離すること。２点i,jの選び方数、かけることのi,jの距離とする。あとはこれに２点以外の選び方のパターン数をかければ良い
# X,Yを独立して考える
# 例えばXについて考えると、使う列の選び方がN*2, 使う行の選び方はM-dis(1ならM-1だし、M-1離れてるのは１個しか選べない）
# disは距離が答えに寄与するから
# cmb()はその２点以外の選び方の個数の話。
#  もちろん他のi,jについて見たら同じ盤面何度も見てるけど、今のi,jに関する値しか計算してないからOK
ans = 0

for dis in range(M):
  ans += N**2 * (M-dis) * dis * cmb(N*M-2, K-2, mod)
  ans %= mod
for dis in range(N):
  ans += M**2 * (N-dis) * dis * cmb(N*M-2, K-2, mod)
  ans %= mod

print(ans)

Submission status: Accepted

Solution code:
def sum_distance(a, b):
    x = 0
    for i in range(a):
        x += ((i + 1) * (a - i))
        x %= mod
    x *= pow(b, 2)
    return x % mod

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

n, m, k = map(int, input().split())
mod = pow(10, 9) + 7
ans = 0
ans += sum_distance(n - 1, m)
ans += sum_distance(m - 1, n)
ans %= mod
ans *= comb(n * m - 2, k - 2)
ans %= mod
print(ans)

Submission status: Accepted

Solution code:
class Combination:
    def __init__(self, mod, max_n):
        self.MOD = mod
        self.MAX_N = max_n

        self.f = self.factorial(self.MAX_N)
        self.f_inv = [self.inv(x) for x in self.f]

    def inv(self,x):
        return pow(x, self.MOD-2, self.MOD)
    
    def factorial(self, n):
        res = [1]
        for i in range(1,n+1):
            res.append(res[-1] * i % self.MOD)
        return res
    
    def comb(self, n, r):
        return (self.f[n] * self.f_inv[r] % self.MOD) * self.f_inv[n-r] % self.MOD

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

ans = 0

for d in range(N):
    ans += (((CB.comb(N*M-2, K-2) * (N-d) % CB.MOD ) * M % CB.MOD ) * M % CB.MOD) * d % CB.MOD
    ans %= CB.MOD

for d in range(M):
    ans += (((CB.comb(N*M-2, K-2) * (M-d) % CB.MOD ) * N % CB.MOD ) * N % CB.MOD) * d % CB.MOD
    ans %= CB.MOD

print(ans)

Submission status: Accepted

Solution code:
from collections import defaultdict,deque
import math
def gcd(a,b):
    return b if not a%b else gcd(b,a%b)
def lcm(a,b):
    return (a*b)//gcd(a,b)
def prepare(n, MOD):
    # 1! - n! の計算
    f = 1
    factorials = [1]  # 0!の分
    for m in range(1, n + 1):
        f *= m
        f %= MOD
        factorials.append(f)
    # n!^-1 の計算
    inv = pow(f, MOD - 2, MOD)
    # n!^-1 - 1!^-1 の計算
    invs = [1] * (n + 1)
    invs[n] = inv
    for m in range(n, 1, -1):
        inv *= m
        inv %= MOD
        invs[m - 1] = inv

    return factorials, invs
    #nCr = (factorials[n]%MOD * invs[r]%MOD * invs[n-r]%MOD)%MOD
def main():
    n,m,k=map(int, input().split())
    MOD = 10**9 + 7
    #k=2の場合nについて最大の差n-1は1通り、n-2は2通り・・・n-iはi通り(i<n)
    facs,invs = prepare(max(k,n*m),MOD)
    res = 0
    tmp = min(k,n)
    for i in range(1,n):
        res = res + ((n-i)*i*m*m)*(facs[n*m-2]%MOD * invs[k-2]%MOD * invs[n*m-2-(k-2)]%MOD)%MOD
        res%=MOD
    tmp = min(k,m)
    for i in range(1,m):
        res = res + ((m-i)*i*n*n)*(facs[n*m-2]%MOD * invs[k-2]%MOD * invs[n*m-2-(k-2)]%MOD)%MOD
        res%=MOD
    print(res)


if __name__ == '__main__':
    main()

Submission status: Accepted

Solution code:
import sys
def MI(): return map(int,sys.stdin.readline().rstrip().split())

N,M,K = MI()
mod = 10**9+7


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

def nCr(n,r):
    return (kaijou(n)*pow(kaijou(r),mod-2,mod)*pow(kaijou(n-r),mod-2,mod)) % mod


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

Submission status: Accepted

Solution code:
def inverse(a, mod):
    while(a < 0):
        a += mod
    a %= mod
    p = mod
    ya = 0
    yp = 1
    while (a > 0):
        n = p // a
        r = p % a
        yr = ya - n * yp
        ya = yp
        yp = yr
        p = a
        a = r
    return ya % mod if ya >= 0 else ya + mod

mod = 1000000007

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

inv = 1
for i in range(1, k-1):
    inv = inv * i % mod
inv = inverse(inv, mod)

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


xn = n * (n + 1) * (n - 1) // 6 % mod * m % mod * m % mod
yn = m * (m + 1) * (m - 1) // 6 % mod * n % mod * n % mod
ans = (xn + yn) * combination % mod
print(ans)

Submission status: Accepted

Solution code:
p,q,r=map(int,input().split())

    
def comb_mod(N,K,modp):

    K_fact=1
    NK_fact=1
    N_fact=1

    for i in range(1,K+1):
        K_fact*=i
        K_fact=K_fact%modp
    for i in range(1,N-K+1):
        NK_fact*=i
        NK_fact=NK_fact%modp
    for i in range(1,N+1):
        N_fact*=i
        N_fact=N_fact%modp

    K_inv=pow(K_fact, -1, modp)
    NK_inv=pow(NK_fact, -1, modp)
    pat=(N_fact*K_inv)%modp*NK_inv%modp
    return pat

pat=comb_mod(p*q-2,r-2,10**9+7)

ans=0
for i in range(0,q+1):
    ans+=i*(q-i)*p*p
    ans=ans%(10**9+7)
    
for k in range(0,p+1):
    ans+=k*(p-k)*q*q%(10**9+7)
    ans=ans%(10**9+7)

print(ans*pat%(10**9+7))

Submission status: Accepted

Solution code:
# nCk ( O(N) )
N_MAX = 10**6
MOD = 10**9 + 7
fac, finv, inv = [0]*N_MAX ,[0]*N_MAX, [0]*N_MAX
def com_init():
    fac[0], fac[1] = 1, 1
    finv[0], finv[1] = 1, 1
    inv[1] = 1
    for i in range(2, N_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 com(n, k):
    if n < k: return 0
    if n < 0 or k < 0: return 0
    return fac[n] * (finv[k] * finv[n - k] % MOD) % MOD


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

ans = 0
for dist in range(1,m):
    comb = (m-dist)*n*n
    other_comb = com(n*m-2,k-2)
    ans += comb*other_comb*dist
    ans %= MOD

for dist in range(1,n):
    comb = (n-dist)*m*m
    other_comb = com(n*m-2,k-2)
    ans += comb*other_comb*dist
    ans %= MOD

print(ans)

Submission status: Accepted

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

MOD = 10**9 + 7
num = 2000020
fact = [1] * (num+1)
ifact = [1] * (num+1)

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

memo = {}
def nCr2(n,r):
    if (n,r) in memo:
        return memo[(n,r)]
    else:
        memo[(n,r)] = (fact[n] * ifact[r] * ifact[n-r]) % MOD
        return memo[(n,r)]

def nCr(n,r):
    if r > n:
        return 0
    return (fact[n] * ifact[r] * ifact[n-r]) % MOD

ans = 0
ans1 = 0
ans2 = 0

for j in range(2,K+1):
    ans1 += M**2 * nCr(2*(M-1),j-2)* nCr((N-2)*M, K-j) 
    ans1 %= MOD
ans1 *= sum((i*(N-i) for i in range(1,N)))

for j in range(2,K+1):
    ans2 += N**2 * nCr(2*(N-1),j-2) * nCr((M-2)*N, K-j)
    ans2 %= MOD
ans2 *= sum((i*(M-i) for i in range(1,M)))

ans = (ans1 + ans2) % MOD
'''

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

for i in range(1,M):
    print(i)
    for j in range(2,min(K,2*N)+1):
        for k in range(1,j):
            ans += i * (M-i) * k* (j-k) * nCr(N,k) * nCr(N,j-k) * nCr((M-2)*N, K-j)
            ans %= MOD
'''

print(ans)

Submission status: Accepted

Solution code:
import sys
import heapq
import math
import fractions
import bisect
import itertools
from collections import Counter
from collections import deque
from operator import itemgetter
def input(): return sys.stdin.readline().strip()
def mp(): return map(int,input().split())
def lmp(): return list(map(int,input().split()))

mod = 10**9 + 7

# x ** a をmodで割った余りを、O(log(a))時間で求める。
def power(x, a):
    if a == 0:
        return 1
    elif a == 1:
        return x
    elif a % 2 == 0:
        return power(x, a//2) **2 % mod
    else:
        return power(x, a//2) **2 * x % mod

# xの逆元を求める。フェルマーの小定理より、 x の逆元は x ^ (mod - 2) に等しい。計算時間はO(log(mod))程度。
# https://qiita.com/Yaruki00/items/fd1fc269ff7fe40d09a6
def modinv(x):
    return power(x, mod-2)

def binomial_coefficients(n, k):
    numera = 1  # 分子
    denomi = 1  # 分母

    for i in range(k):
        numera *= n-i
        numera %= mod
        denomi *= i+1
        denomi %= mod
    return numera * modinv(denomi) % mod

ans=0
n,m,k=mp()
for i in range(n):
    for j in range(m):
        if i!=0 or j!=0:
            if i==0 or j==0:
                ans+=(n-i)*(m-j)*(i+j)%mod
            else:
                ans+=2*(n-i)*(m-j)*(i+j)%mod
            ans%=mod
print(ans*binomial_coefficients(n*m-2, k-2)%mod)

Submission status: Accepted

Solution code:
def ext_euc(a, b):
    # return x, y, gcd(a, b) s.t. ax + by = gcd(a, b)
    if b == 0:
        return 1, 0, a
    y, x, v = ext_euc(b, a % b)
    y -= (a // b) * x
    return x, y, v


def mod_inv(a, mod):
    x, _, _ = ext_euc(a, mod)
    return x % mod


def comb(n, k, mod):
    if k >= n // 2:
        return comb(n, n - k, mod)
    ret = 1
    for i in range(k):
        ret *= n - i
        ret *= mod_inv(k - i, mod)
        ret %= mod
    return ret


def main():
    MOD = 10 ** 9 + 7
    N, M, K = list(map(int, input().split(' ')))
    ans = comb(N * M - 2, K - 2, MOD)
    ans *= ((N ** 2) * (M - 1) * M * (M + 1) + (M ** 2) * (N - 1) * N * (N + 1)) // 6
    ans %= MOD
    print(ans)


if __name__ == '__main__':
    main()

Submission status: Accepted

Solution code:
import sys
input = sys.stdin.readline

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

def main():
    mod=10**9+7
    N,M,K=MI()
    #2この組み合わせを全部足して，平均ぽく考える
    
    def sum2(x):
        return ((x*(x+1))//2)%mod
    
    #i,jますにいるときの相手．（i,j）よりも奥にあるやつの全組み合わせの総和
    def calc(i,j):
        temp=0
        temp=sum2(M-j-1)#該当する行
        temp=(temp+((sum2(M-j-1)+sum2(j))*(N-i-1)))%mod#下の行の階差部分，左右方向にあるから
        temp=(temp+sum2(N-i-1)*M)%mod
        return temp
        
    S=0
    for i in range(N):
        for j in range(M):
            S=(S+calc(i,j))%mod
            
    rem=(N*M)-2
    
    def cmb(n, r, mod):
        if (r < 0) or (n < r):
            return 0
        r = min(r, n - r)
        return (fact[n] * factinv[r] * factinv[n-r])%mod

    fact=[1,1]
    factinv=[1,1]
    inv=[0,1]
    
    for i in range(2, rem + 1):
        fact.append((fact[-1] * i) % mod)
        inv.append((-inv[mod % i] * (mod // i)) % mod)
        factinv.append((factinv[-1] * inv[-1]) % mod)  
        
    #print(rem,K-2)
    S=(S*cmb(rem,K-2,mod))%mod
            
    print(S)

main()

Submission status: Accepted

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

mod = 10**9+7 #出力の制限
N = 2*(10**5)
g1 = [1, 1] # 元テーブル
g2 = [1, 1] #逆元テーブル
inverse = [0, 1] #逆元テーブル計算用テーブル

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


N,M,K = map(int,input().split())
a=cmb(N*M-2, K-2, mod)
b=((N**2)*M*(M+1)*(M-1)//6)%(10**9+7)
c=((M**2)*N*(N+1)*(N-1)//6)%(10**9+7)
print((a*(b+c))%(10**9+7))

Submission status: Accepted

Solution code:
def cmb(a,b,c):
    b = min(b,a-b)
    num = 1
    for i in range(b):
        num = num*(a-i) % c
    den = 1
    for i in range(b):
        den = den*(i+1) % c
    return num * pow(den,c-2,c) % c
n,m,k = map(int,input().split())
mod = 10**9+7
ans = 0
C = cmb(n*m-2,k-2,mod)
for i in range(1,n):
    ans += i*(n-i)*m**2*C
    ans %= mod
for i in range(1,m):
    ans += i*(m-i)*n**2*C
    ans %= mod
print(ans)

Submission status: Accepted

Solution code:
mod = 10 ** 9 + 7
MAX = 2 * 10 ** 5

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

inv = [1] * (MAX + 1)
for i in range(2, MAX + 1):
    inv[i] = inv[mod % i] * (mod - mod // i) % mod

fact_inv = [1] * (MAX + 1)
for i in range(1, MAX + 1):
    fact_inv[i] = fact_inv[i-1] * inv[i] % mod


def comb(n, k):
    if n < k:
        return 0
    return fact[n] * fact_inv[n-k] * fact_inv[k] % mod


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

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

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

print(ans)

Submission status: Accepted

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

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

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

res = res_n + res_m
#print(res)

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

res = res * f[n*m-2] * pow(f[k-2], mod-2, mod) * pow(f[n*m-k], mod-2, mod) % mod

print(res)

Submission status: Accepted

Solution code:
import sys
sys.setrecursionlimit(10000000)
MOD = 10 ** 9 + 7
INF = 10 ** 15

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

    MAXN = N*M
    fact = [1]
    for i in range(1,MAXN + 1):
        fact.append(fact[-1]*i%MOD)
    inv_fact = [-1] * (MAXN + 1)
    inv_fact[-1] = pow(fact[-1],MOD - 2,MOD)
    for i in range(MAXN - 1,-1,-1):
        inv_fact[i] = inv_fact[i + 1]*(i + 1)%MOD
    
    nck = lambda N,K: 0 if K > N or K < 0 else fact[N]*inv_fact[N - K]*inv_fact[K]%MOD

    ans = 0
    const = nck(N*M - 2,K - 2)
    for i in range(1,N):
        ans += (N - i) * M*M * const * i % MOD
        ans %= MOD
    
    for i in range(1,M):
        ans += (M - i) * N*N * const * i % MOD
        ans %= MOD
    print(ans)
if __name__ == '__main__':
    main()

Submission status: Accepted

Solution code:
class Calc:
    def __init__(self, max_value, mod):
        """combination(max_value, all)"""
        fact = [-1] * (max_value + 1)
        fact[0] = 1
        fact[1] = 1
        for x in range(2, max_value + 1):
            fact[x] = x * fact[x - 1] % mod

        invs = [1] * (max_value + 1)
        invs[max_value] = pow(fact[max_value], mod - 2, mod)
        for x in range(max_value - 1, 0, -1):
            invs[x] = invs[x + 1] * (x + 1) % mod

        self.fact = fact
        self.invs = invs
        self.mod = mod

    def nCr(self, n, r):
        r = min(n - r, r)
        if r < 0: return 0
        if r == 0: return 1
        if r == 1: return n
        return self.fact[n] * self.invs[r] * self.invs[n - r] % self.mod

    def nHr(self, n, r):
        return self.nCr(n - 1 + r, r)


def main():
    MOD = 10 ** 9 + 7

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

    calc = Calc(max_value=N * M - 2, mod=MOD)

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

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

    ans = (t + u) % MOD

    print(ans)


if __name__ == '__main__':
    main()

Submission status: Accepted

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

Submission status: Accepted

Solution code:
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():
        # 素数のmod取るときのみ　速い
        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)
    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
"""

Submission status: Accepted

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

Submission status: Accepted

Solution code:
# E - Cell Distance
# https://atcoder.jp/contests/abc127/tasks/abc127_e

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)

Submission status: Accepted

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

Submission status: Accepted

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

Submission status: Accepted

Solution code:
# -*- coding: utf-8 -*-
import sys
# sys.setrecursionlimit(10**6)
# buff_readline = sys.stdin.buffer.readline
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()

Submission status: Accepted

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

Submission status: Accepted

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

Submission status: Accepted

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

Submission status: Accepted

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