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:
def combnk_mod(n, k, p):
    # pが素数かつaとpが互いに素であるとき
    # フェルマーの小定理より以下が成り立つ.
    # 1 / a = a ** (p - 2) (mod p)
    # これを使って (1 / k!) mod p を計算する.

    if  2 * k > n:
        k = n - k
    x = 1 # 分母にあたる数 1 / n! mod p
    y = 1 # 分子にあたる数. n! / (n - k)! mod p
    for i in range(1, k + 1):
        x *= pow(i, p - 2, p)
        y *= n - i + 1
        x %= p
        y %= p
    return (x * y) % p

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

Submission status: Accepted

Solution code:
class FactMod:

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

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

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

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

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

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

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

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

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

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

Submission status: Accepted

Solution code:
mod=10**9+7
def mod_combination(n, k, mod):
    def ext_gcd(a, b):
        if b == 0:return a, 1, 0
        else:d,x,y = ext_gcd(b,a%b);x-=(a//b)*y 
        return d,y,x
    p,q=1,1
    for i in range(n-k+1,n+1):
        p=(p*i)%mod
    for i in range(2,k+1):
        q=(q*i)%mod
    return int(p*(ext_gcd(q, mod)[1]%mod)%mod)


n,m,k=map(int,input().split())
print(((pow(m,2,mod)*mod_combination(n*m-2,k-2,mod)*sum((n-i)*i%mod for i in range(1,n))%mod)+pow(n,2,mod)*mod_combination(n*m-2,k-2,mod)*sum((m-i)*i%mod for i in range(1,m))%mod)%mod)

Submission status: Accepted

Solution code:
from math import factorial

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

n, r = N * M - 2, K - 2
#P = factorial(n) // factorial(n-r) // factorial(r) % MOD
P = factorial(n) * pow(factorial(n-r) % MOD, MOD - 2, MOD) * pow(factorial(r) % MOD, MOD - 2, MOD)

ans = 0
for n in range(N) :
    for m in range(M) :
        if n == 0 and m == 0 :
            continue

        d = n + m
        cnt = (N - n) * (M - m) * (1 if n == 0 or m == 0 else 2)

        ans += d * cnt

ans = (ans * P) % MOD
print(ans)

Submission status: Accepted

Solution code:
mod = 10**9 + 7
n,m,k = map(int,input().split())
nm = n*m
f = [1]*(nm+1)
for i in range(1,nm+1):
    f[i] = (f[i-1]*i)%mod
def C(n,k):
    x = f[n]
    x *= pow(f[n-k],mod-2,mod)
    x %= mod
    x *= pow(f[k],mod-2,mod)
    x %= mod
    return x
def F(n):
    x = n*n*(n+1)//2
    y = n*(n+1)*(2*n+1)//6
    return (x-y)%mod
res = (m*m%mod*F(n)%mod + n*n%mod*F(m)%mod)%mod
res *= C(n*m-2,k-2)
print(res%mod)

Submission status: Accepted

Solution code:
import math
N,M,k = map(int,input().split())
p = N*M
def cmb(n, r):
    if n - r < r: r = n - r
    if r == 0: return 1
    if r == 1: return n

    numerator = [n - r + k + 1 for k in range(r)]
    denominator = [k + 1 for k in range(r)]

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

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

    return result

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

Submission status: Accepted

Solution code:
#!/usr/bin/env python3
import sys
INF = float("inf")

MOD = 1000000007  # type: int


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

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

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

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

    return result


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

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

    return


def main():

    def iterate_tokens():
        for line in sys.stdin:
            for word in line.split():
                yield word
    tokens = iterate_tokens()
    N = int(next(tokens))  # type: int
    M = int(next(tokens))  # type: int
    K = int(next(tokens))  # type: int
    solve(N, M, K)


if __name__ == '__main__':
    main()

Submission status: Accepted

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

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

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

res = 0

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

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

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

Submission status: Accepted

Solution code:
#!/usr/bin/env python3
import sys
INF = float("inf")

MOD = 1000000007  # type: int


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

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

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

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

    return result


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

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

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

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


def main():

    def iterate_tokens():
        for line in sys.stdin:
            for word in line.split():
                yield word
    tokens = iterate_tokens()
    N = int(next(tokens))  # type: int
    M = int(next(tokens))  # type: int
    K = int(next(tokens))  # type: int
    solve(N, M, K)


if __name__ == '__main__':
    main()

Submission status: Accepted

Solution code:
import math
N,M,k = map(int,input().split())
p = N*M
def cmb(n, r):
    if n - r < r: r = n - r
    if r == 0: return 1
    if r == 1: return n

    numerator = [n - r + k + 1 for k in range(r)]
    denominator = [k + 1 for k in range(r)]

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

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

    return result

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

Submission status: Accepted

Solution code:
import numpy as np
MOD = 1000000007

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

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

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

#互いに素なa,bについて、a*x+b*y=1の一つの解
def extgcd(a,b):
    r = (1,0,a)
    w = (0,1,b)
    while w[2]!=1:
        q = r[2]//w[2]
        r2 = w
        w2 = [r[0]-q*w[0],r[1]-q*w[1],r[2]-q*w[2]]
        r = r2
        w = w2
    return (w[0],w[1])

def mod_inv(a):
    x,  = extgcd(a, MOD)
    return (MOD+x%MOD)%MOD

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

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

MAX = 510000

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

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

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

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

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

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

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


print(res)

Submission status: Accepted

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

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

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


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


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

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

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

Submission status: Accepted

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


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


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

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


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


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

Submission status: Accepted

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

MOD = 10**9+7

# WARN: Require mod number m is prime number
class ModComb:
    # For nCr mod m
    def __init__(self, n_max, m):
        self._n_max = n_max
        self._m = m
        self._mod_facts = [1] * (n_max+1)       # list of n! mod m
        self._mod_inv_facts = [1] * (n_max+1)   # list of inverse r! mod m

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

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

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

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

        if r == 0:
            return 1

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

modcomb = ModComb(N*M-2, MOD)
p = modcomb.calc(N*M-2, K-2)

ans = 0
for i in range(N):
    for j in range(M):
        if i == 0 and j == 0:
            continue
        d = i+j
        cnt = (N-i) * (M-j)
        if i != 0 and j != 0:
            cnt *= 2
        ans += d * cnt % MOD
ans = ans * p % MOD

print(ans)

Submission status: Accepted

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

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

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

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

print(ans)

Submission status: Accepted

Solution code:
# -*- coding: utf-8 -*-
import bisect
import heapq
import math
import random
import sys
from collections import Counter, defaultdict, deque
from decimal import ROUND_CEILING, ROUND_HALF_UP, Decimal
from functools import lru_cache, reduce
from itertools import combinations, combinations_with_replacement, product, permutations
from operator import add, mul, sub

sys.setrecursionlimit(100000)


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


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


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


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


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


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


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


def mt(f):
    import time

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

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

    return wrap


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

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

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

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

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

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


class Combination:
    def __init__(self, n, mod):

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

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


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

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

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

    return ans


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


if __name__ == '__main__':
    main()

Submission status: Accepted

Solution code:
LAW = 10 ** 9 + 7

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

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

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

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

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

print(c(N, M, K))

Submission status: Accepted

Solution code:
n,m,k = map(int,input().split())
ans = 0
mod = 10**9+7
max = 2*10**5
g1 = [1, 1] # 元テーブル
g2 = [1, 1] #逆元テーブル
inverse = [0, 1] #逆元テーブル計算用テーブル
for i in range( 2, max + 1 ):
    g1.append( ( g1[-1] * i ) % mod )
    inverse.append( ( -inverse[mod % i] * (mod//i) ) % mod )
    g2.append( (g2[-1] * inverse[-1]) % mod )
def cmb(n, r, mod):
    if ( r<0 or r>n ):
        return 0
    r = min(r, n-r)
    return g1[n] * g2[r] * g2[n-r] % mod
C = cmb(m*n-2,k-2,mod)
for i in range(1,n):
    ans += (i*(n-i)*(m**2))*C
for i in range(1,m):
    ans += (i*(m-i)*(n**2))*C
print(ans%mod)

Submission status: Accepted

Solution code:
from math import factorial
N, M, K = map(int, input().split())

MOD = 10**9+7

fact = [1] * (N*M+1)
for i in range(2, N*M+1):
    fact[i] = fact[i-1] * i % MOD

def comb(n, r):
    r = n-r if n-r < r else r
    return fact[n] * pow(fact[r], MOD-2, MOD) * pow(fact[n-r], MOD-2, MOD) % MOD

p = comb(N*M-2, K-2)
ans = 0
for i in range(N):
    for j in range(M):
        if i == 0 and j == 0:
            continue
        d = i+j
        cnt = (N-i) * (M-j)
        if i != 0 and j != 0:
            cnt *= 2
        ans += d * cnt % MOD
ans = ans * p % MOD

print(ans)

Submission status: Accepted

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

n,m,k = map(int,input().split())
Max = n*m
#二項係数とその逆元テーブルを作る前処理
fac = [0]*(Max)
finv = [0]*(Max)
inv = [0]*(Max)
fac[0]=fac[1]=1
finv[0]=finv[1]=1
inv[1]=1
for i in range(2,Max):
        fac[i] = fac[i-1] * i % mod
        inv[i] = mod - inv[mod%i]*(mod//i)%mod
        finv[i] = finv[i-1]*inv[i]%mod
#O(1)でmod計算した組合せ数を計算
def C(n,r):
    if n < r:
        return 0
    if n < 0 or r < 0 :
        return 0
    return fac[n]*(finv[r]*finv[n-r]%mod)%mod
#mod計算した順列の計算
def P(n,r):
    if n < r:
        return 0
    if n<0 or r<0:
        return 0
    return (fac[n]*finv[n-r])%mod


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

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

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

Submission status: Accepted

Solution code:
def main():
     def cmb(n, r, mod):
          if ( r<0 or r>n ):
               return 0
          r = min(r, n-r)
          return g1[n] * g2[r] * g2[n-r] % mod
    
     mod = 10**9+7 #出力の制限
     N = 10**6
     g1 = [1, 1] # 元テーブル
     g2 = [1, 1] #逆元テーブル
     inverse = [0, 1] #逆元テーブル計算用テーブル
     
     for i in range( 2, N + 1 ):
          g1.append( ( g1[-1] * i ) % mod )
          inverse.append( ( -inverse[mod % i] * (mod//i) ) % mod )
          g2.append( (g2[-1] * inverse[-1]) % mod )

     n,m,k = map(int,input().split())
     nm = n*m%(10**9+7)
     a = (cmb(nm-2,k-2,mod)%(10**9+7))*(nm*(m*(n**2-1)+n*(m**2-1)))//6
     print(a%(10**9+7))
     
if __name__=='__main__':
     main()

Submission status: Accepted

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

#def expo_div_p(a: int,n: int,p: int):
#    if n == 1:
#        d = a
#    elif n % 2 == 0:
#        d = expo_div_p(a, n//2, p)**2 % p
#    else:
#        d = a * expo_div_p(a, n//2, p)**2 % p
#    return d
#
def combi(x, y, p):
    c = 1
    for i in range(1,y+1):
        c *= (x-i+1) * pow(i, p-2, p)
        c %= p
    return c


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

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

#print(C)
print(ans % ppp)

Submission status: Accepted

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

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 = 10**6
g1 = [1, 1] # 元テーブル
g2 = [1, 1] #逆元テーブル
inverse = [0, 1] #逆元テーブル計算用テーブル

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

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

Submission status: Accepted

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

## コンビネーション計算用 ========
## https://qiita.com/derodero24/items/91b6468e66923a87f39f
mod = int(1e9+7)
g1 = [1, 1] # 元テーブル
g2 = [1, 1] # 逆元テーブル
inverse = [0, 1] # 逆元テーブル計算用テーブル

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

def combination(n, r, mod):
    if ( r<0 or r>n ):
        return 0
    r = min(r, n-r)
    return g1[n] * g2[r] * g2[n-r] % mod
## コンビネーション計算用 --------
  

## 二つの駒の距離の総和 = cost_allpairs
var_1 = M**2 * (N-1) * N * (N+1) // 3
var_2 = N**2 * (M-1) * M * (M+1) // 3
cost_allpairs = int((var_1 + var_2) // 2)
cost_allpairs = cost_allpairs % mod ## ここでmod取っておく

## コンビネーションの計算
cmb = combination(N*M-2, K-2, mod)

## 答え
ans = (cost_allpairs * cmb) % mod
print(int(ans))

Submission status: Accepted

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

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

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

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

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

print(result % div)

Submission status: Accepted

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


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

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

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


# pairs = []
# for i in range(1, n + 1):
#     for j in range(1, m + 1):
#         pairs.append((i, j))
#
# pairs.sort()
# print(pairs)
#
# brute = 0
# for i in range(len(pairs)):
#     for j in range(i + 1, len(pairs)):
#         brute += abs(pairs[i][0] - pairs[j][0]) + \
#                      abs(pairs[i][1] - pairs[j][1])
#
# print(ans * brute)

Submission status: Accepted

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

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

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

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

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

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

Submission status: Accepted

Solution code:
N, M, K = map(int, input().split())
ans=0
mod = 10**9+7 #出力の制限
max = 2*10**5
g1 = [1, 1] # 元テーブル
g2 = [1, 1] #逆元テーブル
inverse = [0, 1] #逆元テーブル計算用テーブル
for i in range( 2, max + 1 ):
    g1.append( ( g1[-1] * i ) % mod )
    inverse.append( ( -inverse[mod % i] * (mod//i) ) % mod )
    g2.append( (g2[-1] * inverse[-1]) % mod )
def cmb(n, r, mod):
    if ( r<0 or r>n ):
        return 0
    r = min(r, n-r)
    return g1[n] * g2[r] * g2[n-r] % mod
for d in range(1,N):
    ans += d*(N-d)*M**2 % mod
for d in range(1,M):
    ans += d*(M-d)*N**2 % mod
ans = ans*cmb(N*M-2, K-2, mod) % mod
print(ans)

Submission status: Accepted

Solution code:
N, M, K = [int(item) for item in input().split()]
MOD = 10**9 + 7

n = 2 * 10**5 + 10
fac = [1] + [0] * n
fac_inv = [1] + [0] * n
for i in range(1, n+1):
    fac[i] = fac[i-1] * (i) % MOD
    # Fermat's little theorem says
    # a**(p-1) mod p == 1
    # then, a * a**(p-2) mod p == 1
    # it means a**(p-2) is inverse element
    fac_inv[i] = fac_inv[i-1] * pow(i, MOD-2, MOD) % MOD

def mod_nCr(n, r):
    if n == 0 and r == 0:
        return 1
    if n < r or n < 0:
        return 0
    tmp = fac_inv[n-r] * fac_inv[r] % MOD
    return tmp * fac[n] % MOD 

total = 0
for i in range(M):
    total += (M - i) * N * N * mod_nCr(M * N - 2, K - 2) * i
    total %= MOD  
for i in range(N):
    total += (N - i) * M * M * mod_nCr(M * N - 2, K - 2) * i
    total %= MOD  
print(total)

Submission status: Accepted

Solution code:
N,M,K=map(int,input().split())
cost=0
mod=10**9+7
NM=N*M
kaijo=[0]*(NM-1)
kaijo[0]=1
for i in range(1,NM-1):
    kaijo[i]=(kaijo[i-1]*i)%mod
gyaku=[0]*(NM-1)
gyaku[NM-2]=pow(kaijo[NM-2],mod-2,mod)
for i in range(NM-2,0,-1):
    gyaku[i-1]=(gyaku[i]*i)%mod
def conb(a,b):
    return (kaijo[a]*gyaku[b]*gyaku[a-b])%mod
x=conb(NM-2,K-2)
for d in range(1,N):
    cost=(cost+x*d*(N-d)*(M**2))%mod
for d in range(1,M):
    cost=(cost+x*d*(M-d)*(N**2))%mod
print(cost)

Submission status: Accepted

Solution code:
ri = lambda: int(input())
rl = lambda: list(map(int,input().split()))
rr = lambda N: [ri() for _ in range(N)]
YN = lambda b: print('YES') if b else print('NO')
yn = lambda b: print('Yes') if b else print('No')
OE = lambda x: print('Odd') if x%2 else print('Even')
INF = 10**18
MOD=10**9+7

N,M,K=rl()

#nCr % MOD
def cmb(n, r, MOD=0):
    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])
            if MOD != 0:
                result %= MOD
 
    return result


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

Submission status: Accepted

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

factorial=[1]
for i in range(1,s+1):
    factorial.append(factorial[i-1]*i%mod)
inverse=[0]*(s+1)
inverse[-1]=pow(factorial[-1],mod-2,mod)
for i in range(s-1,-1,-1):
    inverse[i]=inverse[i+1]*(i+1)%mod
def comb(n,r):
    if n<r:
        return 0
    return factorial[n]*inverse[r]*inverse[n-r]%mod

c=comb(n*m-2,k-2)
ans=0
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:
import sys,heapq
from collections import deque,defaultdict
printn = lambda x: sys.stdout.write(x)
inn = lambda : int(input())
inl   = lambda: list(map(int, input().split()))
inm   = lambda:      map(int, input().split())
R = 10**9+7
DBG = True  and False
def ddprint(x):
  if DBG:
    print(x)

def exgcd(x,y):
    if y == 0:
        return (x,1,0)
    else:
        g,b,a = exgcd(y, x%y)
        return (g, a, b-(x//y)*a)

# return y s.t. x*y mod p == 1
def modinv(x,p):
    g,a,b = exgcd(x,p)  # ax+bp = g, g=1 if mutually prime
    return a%p

def comb(p,q):
    if q==0 or q==p:
        return 1
    if q > p-q:
        q = p-q
    num = p
    den = q
    for i in range(1,q):
        num = (num * (p-i))%R
        den = (den * (q-i))%R
    xden = modinv(den,R)
    return (num*xden)%R

n,m,k = inm()
x = ( (n*n*(m-1)*m*(m+1)+m*m*(n-1)*n*(n+1))//6 )%R
y = comb(n*m-2,k-2)
ddprint(x)
ddprint(y)
print((x*y)%R)

Submission status: Accepted

Solution code:
from itertools import accumulate
def inpl(): return list(map(int, input().split()))
MOD = 10**9 + 7
N, M, K = inpl()

X = sum([n*(n+1)//2 for n in range(1, N)])*M*(M)%MOD
Y = sum([m*(m+1)//2 for m in range(1, M)])*N*(N)%MOD


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


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 )

print((X+Y)*cmb(N*M-2, K-2, MOD)%MOD)

Submission status: Accepted

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

mod = 10 ** 9 + 7


class comb(object):
    def __init__(self, N, mod):
        self.mod = mod
        self.fac = [1] * N
        for i in range(1, N):
            self.fac[i] = (self.fac[i - 1] * i) % mod
        self.inv = [1] * N
        self.inv[N - 1] = pow(self.fac[N - 1], mod - 2, mod)
        for i in range(N - 1, 0, -1):
            self.inv[i - 1] = (self.inv[i] * i) % mod

    def calc(self, n, k):
        return self.fac[n] * (self.inv[k] * self.inv[n - k] %
                              self.mod) % self.mod


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