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:
import sys
s2nn = lambda s: [int(c) for c in s.split(' ')]
ss2nn = lambda ss: [int(s) for s in ss]
ss2nnn = lambda ss: [s2nn(s) for s in ss]
i2s = lambda: sys.stdin.readline().rstrip()
i2n = lambda: int(i2s())
i2nn = lambda: s2nn(i2s())
ii2ss = lambda n: [sys.stdin.readline().rstrip() for _ in range(n)]
ii2sss = lambda n: [list(sys.stdin.readline().rstrip()) for _ in range(n)]
ii2nn = lambda n: ss2nn(ii2ss(n))
ii2nnn = lambda n: ss2nnn(ii2ss(n))

MOD = int(1e+9) + 7

def comb(n, k):
    r = 1
    for i in range(n, max(n-k, k), -1):
        r *= i
    for i in range(1, min(n-k, k) + 1, 1):
        r //= i
    return r

class Comb:
    def __init__(self, MAX, MOD):
        MAX += 1
        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
        self.fac = fac
        self.finv = finv
        self.MOD = MOD
    
    def comb(self, n, k):
        fac = self.fac
        finv = self.finv
        MOD = self.MOD
        if (n < k):
            return 0
        if (n < 0 or k < 0):
            return 0
        return fac[n] * (finv[k] * finv[n - k] % MOD) % MOD

def main():
    # 行数, 列数, 駒数
    N, M, K = i2nn()
    # 配置コストの総和 % MOD

    # 例）N, M, K = 2, 2, 2
    # 駒が置ける場所は4箇所: N * M
    # 駒の配置は6通り: comb(4, 2) = comb(N*M, K)
    # K=2 の場合のコスト
    xcost = 0
    for i in range(1, M):
        # iをx座標の差とする
        # x座標のみに注目した場合、差がiとなるのは(M-i)通り
        #   M=4, i=1, => (1,2), (2,3), (3,4) = 3通り
        #   M=4, i=2, => (1,3), (2,4) = 2通り
        #   M=4, i=3, => (1,4) = 1通り
        # さらにy座標の組み合わせは, 左側がN通り, 右側がN通り
        #   (1,2) について M=2 の場合 ((1,1), (2,1)), ((1,1), (2,2)), ((1,2), (2,1)), ((1,2), (2,2))
        xcost = (xcost + i * (M-i) * (N*N)) % MOD
    ycost = 0
    for i in range(1, N):
        ycost = (ycost + i * (N-i) * (M*M)) % MOD
    cost = (xcost + ycost) % MOD
    # K=3 の場合、K=2 の場合のコストを NM-2 倍すればよい
    # 汎用的に考えると、 K=2 の場合のコストを comb(MN-2, K-2) 倍すればよい
    #cost = (cost * comb(M*N-2, K-2)) % MOD
    cost = (cost * Comb(N*M-2, MOD).comb(N*M-2, K-2)) % MOD
    print(cost)

main()

Submission status: Accepted

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

# 二項係数を 10^9+7 で割った余りを求める
# https://qiita.com/drken/items/3b4fdf0a78e7a138cd9a#5-%E4%BA%8C%E9%A0%85%E4%BF%82%E6%95%B0-ncr
# https://qiita.com/Yaruki00/items/fd1fc269ff7fe40d09a6
# https://www.hamayanhamayan.com/entry/2018/06/06/210256

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)

# 二項係数の左側の数字の最大値を max_len　とする。nとかだと他の変数と被りそうなので。
# factori_table = [1, 1, 2, 6, 24, 120, ...] 要は factori_table[n] = n!
# 計算時間はO(max_len * log(mod))

max_len = n * m #適宜変更する

factori_table = [1] * (max_len + 1)
factori_inv_table = [1] * (max_len + 1)
for i in range(1, max_len + 1):
	factori_table[i] = factori_table[i-1] * (i) % mod

# 	factori_inv_table[i] = modinv(factori_table[i])
#	↑が遅いので、使うところだけ逆元を計算する↓
for i in [k-2, n*m-k]:
	factori_inv_table[i] = modinv(factori_table[i])

def binomial_coefficients(n, k):
	# n! / (k! * (n-k)! )
	return factori_table[n] * factori_inv_table[k] * factori_inv_table[n-k]

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

print(binomial_coefficients(n*m-2, k-2) * ans_1 % mod)

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 = 3 * 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())

mod = 10**9 + 7

base_n = 0
base_m = 0

for i in range(M):
    if i == 0:
        s = (M * (M - 1) // 2) % mod
    else:
        s += (i - (M - i))
    base_m = (base_m + s + mod) % mod
for i in range(N):
    if i == 0:
        s = (N * (N - 1) // 2) % mod
    else:
        s += (i - (N - i))
    base_n = (base_n + s + mod) % mod
base = (base_m * N * N + base_n * M * M) * g2[2] % mod

ans = base * cmb((N * M) - 2, K - 2, mod) % mod

print(ans)

Submission status: Accepted

Solution code:
# using main() makes code faster from the point of view of "access to variables in global name-space"
# def main():

""" convenient functions
# for i, a in enumerate(iterable)
# q, mod = divmod(a, b)
# divmod(x, y) returns the tuple (x//y, x%y)
# manage median(s) using two heapq https://atcoder.jp/contests/abc127/tasks/abc127_f
"""

import sys
sys.setrecursionlimit(10**7)
from itertools import accumulate, combinations, permutations, product # https://docs.python.org/ja/3/library/itertools.html
# accumulate() returns iterator! to get list: list(accumulate())
from math import factorial, ceil, floor
def factorize(n):
    """return the factors of the Arg and count of each factor
    
    Args:
        n (long): number to be resolved into factors
    
    Returns:
        list of tuples: factorize(220) returns [(2, 2), (5, 1), (11, 1)]
    """
    fct = []  # prime factor
    b, e = 2, 0  # base, exponent
    while b * b <= n:
        while n % b == 0:
            n = n // b
            e = e + 1
        if e > 0:
            fct.append((b, e))
        b, e = b + 1, 0
    if n > 1:
        fct.append((n, 1))
    return fct
def combinations_count(n, r):
    """Return the number of selecting r pieces of items from n kinds of items.
    
    Args:
        n (long): number
        r (long): number
    
    Raises:
        Exception: not defined when n or r is negative
    
    Returns:
        long: number
    """
    # TODO: How should I do when n - r is negative?
    if n < 0 or r < 0:
        raise Exception('combinations_count(n, r) not defined when n or r is negative')
    if n - r < r: r = n - r
    if r < 0: return 0
    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 combinations_with_replacement_count(n, r):
    """Return the number of selecting r pieces of items from n kinds of items allowing individual elements to be repeated more than once.
    
    Args:
        n (long): number
        r (long): number
    
    Raises:
        Exception: not defined when n or r is negative
    
    Returns:
        long: number
    """
    if n < 0 or r < 0:
        raise Exception('combinations_with_replacement_count(n, r) not defined when n or r is negative')
    elif n == 0:
        return 1
    else:
        return combinations_count(n + r - 1, r)
from bisect import bisect_left, bisect_right
from collections import deque, Counter, defaultdict # https://docs.python.org/ja/3/library/collections.html#collections.deque
from heapq import heapify, heappop, heappush, heappushpop, heapreplace,nlargest,nsmallest # https://docs.python.org/ja/3/library/heapq.html
from copy import deepcopy, copy # https://docs.python.org/ja/3/library/copy.html
from operator import itemgetter
# ex1: List.sort(key=itemgetter(1))
# ex2: sorted(tuples, key=itemgetter(1,2))
from functools import reduce
def chmin(x, y):
    """change minimum
    if x > y, x = y and return (x, True).
    convenient when solving problems of dp[i]
    
    Args:
        x (long): current minimum value
        y (long): potential minimum value
    
    Returns:
        (x, bool): (x, True) when updated, else (x, False)
    """
    if x > y:
        x = y
        return (x, True)
    else:
        return (x, False)
def chmax(x, y):
    """change maximum
    if x < y, x = y and return (x, True).
    convenient when solving problems of dp[i]
    
    Args:
        x (long): current maximum value
        y (long): potential maximum value
    
    Returns:
        (x, bool): (x, True) when updated, else (x, False)
    """
    if x < y:
        x = y
        return (x, True)
    else:
        return (x, False)

from fractions import gcd # Deprecated since version 3.5: Use math.gcd() instead.
def gcds(numbers):
    return reduce(gcd, numbers)
def lcm(x, y):
    return (x * y) // gcd(x, y)
def lcms(numbers):
    return reduce(lcm, numbers, 1)

# first create factorial_list
# fac_list = mod_factorial_list(n)
INF = 10 ** 18
MOD = 10 ** 9 + 7
modpow = lambda a, n, p = MOD: pow(a, n, p) # Recursive function in python is slow!
def modinv(a, p = MOD):
    # evaluate reciprocal using Fermat's little theorem:
    # a**(p-1) is identical to 1 (mod p) when a and p is coprime
    return modpow(a, p-2, p)
def modinv_list(n, p = MOD):
    if n <= 1:
        return [0,1][:n+1]
    else:
        inv_t = [0,1]
        for i in range(2, n+1):
            inv_t += [inv_t[p % i] * (p - int(p / i)) % p]
        return inv_t
def modfactorial_list(n, p = MOD):
    if n == 0:
        return [1]
    else:
        l = [0] * (n+1)
        tmp = 1
        for i in range(1, n+1):
            tmp = tmp * i % p
            l[i] = tmp
        return l
def modcomb(n, k, fac_list = [], p = MOD):
    # fac_list = modfactorial_list(100)
    # print(modcomb(100, 5, modfactorial_list(100)))
    from math import factorial
    if n < 0 or k < 0 or n < k: return 0
    if n == 0 or k == 0: return 1
    if len(fac_list) <= n:
        a = factorial(n) % p
        b = factorial(k) % p
        c = factorial(n-k) % p
    else:
        a = fac_list[n]
        b = fac_list[k]
        c = fac_list[n-k]
    return (a * modpow(b, p-2, p) * modpow(c, p-2, p)) % p
def modadd(a, b, p = MOD):
    return (a + b) % MOD
def modsub(a, b, p = MOD):
    return (a - b) % p
def modmul(a, b, p = MOD):
    return ((a % p) * (b % p)) % p
def moddiv(a, b, p = MOD):
    return modmul(a, modpow(b, p-2, p))

""" initialize variables and set inputs
# initialize variables
    # to initialize list, use [0] * n
    # to initialize two dimentional array, use [[0] * N for _ in range(N)]
# set inputs
    # open(0).read() is a convenient method:
    # ex) n, m, *x = map(int, open(0).read().split())
    #     min(x[::2]) - max(x[1::2])
    # ex2) *x, = map(int, open(0).read().split())
    #     don't forget to add comma after *x if only one variable is used
# preprocessing
    # transpose = [x for x in zip(*data)]
    # ex) [[1, 2, 3], [4, 5, 6], [7, 8, 9]] => [(1, 4, 7), (2, 5, 8), (3, 6, 9)]
    # flat = [flatten for inner in data for flatten in inner]
    # ex) [[1, 2, 3], [4, 5, 6], [7, 8, 9]] => [1, 2, 3, 4, 5, 6, 7, 8, 9]
# calculate and output
    # output pattern
    # ex1) print(*l) => when l = [2, 5, 6], printed 2 5 6
"""

# functions used
r = lambda: sys.stdin.readline().strip()
r_int = lambda: int(r())
R = lambda: list(map(int, r().split()))
Rfloat = lambda: list(map(float, r().split()))
Rtuple = lambda: tuple(map(int, r().split()))
Rmap = lambda: map(int, r().split())

""" how to treat input
# single int: int(r())
# single string: r()
# single float: float(r())
# line int: R()
# line string: r().split()
# line (str, int, int): [j if i == 0 else int(j) for i, j in enumerate(r().split())]
# lines int: [R() for _ in range(n)]
"""

# for test
# sys.stdin = open('sample.txt') 

# ----------------------------------
# main
N, M, K = R()

s = modcomb(N*M-2,K-2)
ans = 0
for i in range(1,M):
    ans += modmul(modmul(modmul(i,M-i),modmul(N,N)),s)
N, M = M, N
for i in range(1,M):
    ans += modmul(modmul(modmul(i,M-i),modmul(N,N)),s)

print(ans%MOD)
# end of main
# ----------------------------------

"""memo: how to use defaultdict of list
# initialize
Dic = defaultdict(list)
# append / extend
Dic[x].append(y)
# for
for k, v in Dic.items():
"""

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

Submission status: Accepted

Solution code:
#!usr/bin/env python3
from collections import defaultdict
from collections import deque
from heapq import heappush, heappop
import sys
import math
import bisect
import random
def LI(): return list(map(int, sys.stdin.readline().split()))
def I(): return int(sys.stdin.readline())
def LS():return list(map(list, sys.stdin.readline().split()))
def S(): return list(sys.stdin.readline())[:-1]
def IR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = I()
    return l
def LIR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = LI()
    return l
def SR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = S()
    return l
def LSR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = LS()
    return l
sys.setrecursionlimit(1000000)
mod = 1000000007

#A
def A():
    return

#B
def B():
    return

#C
def C():
    n,m = LI()
    t = [0 for i in range(n+1)]
    for i in range(m):
        l,r = LI()
        t[l-1] += 1
        t[r] -= 1
    for i in range(n):
        t[i+1] += t[i]
        if t[i] == m:
            t[i] = 1
        else:
            t[i] = 0
    print(sum(t[:n]))
    return

#D
def D():
    n,m = LI()
    a = LI()
    q = []
    for i in range(n):
        heappush(q,a[i])
    l = LIR(m)
    l.sort(key = lambda x:-x[1])
    for b,c in l:
        for i in range(b):
            x = heappop(q)
            if x < c:
                x = c
                heappush(q,x)
            else:
                heappush(q,x)
                break
    print(sum(q))
    return

#E
def E():
    n,m,k = LI()
    M = n*m
    fact = [1]*(M+1)
    for i in range(M):
        fact[i+1] = fact[i]*(i+1)%mod
    inv = [1]*(M+1)
    inv[M] = pow(fact[M],mod-2,mod)
    for i in range(M)[::-1]:
        inv[i] = inv[i+1]*(i+1)%mod
    ans = fact[M-2]*inv[k-2]*inv[M-k]%mod
    s = 0
    for i in range(1,M):
        s += i*(max(0,n-i)*m**2+max(0,m-i)*n**2)%mod
        s %= mod
    ans = ans*s%mod
    print(ans)
    return

#F
def F():
    q = I()
    m = 0
    l = [-float("inf"),float("inf")]
    x = 0
    num = 0
    for _ in range(q):
        Q = LI()
        if len(Q) == 1:
            print(x,m)
        else:
            z,a,b = Q
            m += b
            if num == 0:
                num += 1
                l = [a,float("inf")]
                x = a
            elif num == 1:
                if a != l[0]:
                    num += 1
                    l = [min(a,l[0]),max(a,l[0])]
                    x = l[0]
                    m += l[1]-l[0]
            else:
                num += 1
                if num % 2:
                    if l[0] <= a < l[1]:
                        m += min(a-l[0],l[1]-a)
                        x = a
                    elif a < l[0]:
                        m += l[0]-a
                        x = l[0]
                        l = [a,l[1]]
                    else:
                        m += a-l[1]
                        x = l[1]
                        l = [l[0],a]
                else:
                    if l[0] <= a < l[1]:
                        if a < x:
                            m += x-a
                            l = [a,x]
                            x = a
                        else:
                            m += a-x
                            l = [x,a]
                    elif a < l[0]:
                        m += x-a
                        x = l[0]
                        l = [l[0],x]
                    else:
                        m += a-x
                        l = [x,l[1]]
    return

#G
def G():
    return

#H
def H():
    return

#I
def I_():
    return

#J
def J():
    return

#Solve
if __name__ == "__main__":
    E()

Submission status: Accepted

Solution code:
MOD = 10**9 + 7
MAX = int(2e5)


def div(a, b):
    return a * pow(b, MOD-2, MOD) % MOD


FACT = [1] * (MAX+1)
for i in range(1, MAX+1):
    FACT[i] = (i * FACT[i-1]) % MOD
INV = [1] * (MAX+1)
INV[MAX] = div(1, FACT[MAX])
for i in range(MAX, 0, -1):
    INV[i-1] = (INV[i] * i) % MOD


def combination(r, c):
    return FACT[r] * INV[r-c] * INV[c]


def main():
    N, M, K = map(int, input().split())
    ans1 = 0
    for d in range(1, N):
        ans1 = (ans1 + d * (N-d) * (M**2)) % MOD
    ans2 = 0
    for d in range(1, M):
        ans2 = (ans2 + d * (M-d) * (N**2)) % MOD
    ans = ((ans1 + ans2) * combination(M*N-2, K-2)) % MOD
    print(ans)


if __name__ == "__main__":
    main()

Submission status: Accepted

Solution code:
facs = []
M = 1000000007

def modulo(n):
    global facs
    facs = [1]*(n+1)
    last = 1
    for i in range(1,n+1):
        last = facs[i] = mulmod(last,i,M)
        
def mulmod(x,y,p):
    return x*y % p
def divmod(x,y,p):
    return mulmod(x,powmod(y,p-2,p),p)

def ncr(n,r):
    if n<r:
        return 0
    if n==r:
        return 1
    res = facs[n]
    res = divmod(res,facs[r],M)
    res = divmod(res,facs[n-r],M)
    return res

def powmod(x,y,p):
    if y==0:
        return 1
    elif y==1:
        return x % p
    elif (y%2)==0:
        return powmod(x,y//2,p)**2 % p
    else:
        return powmod(x,y//2,p)**2 * x % p
def resolve():
    n,m,k = map(int,input().split())
    modulo(2*10**5+5)
    kumi=ncr(n*m-2,k-2)
    ans = 0
    for i in range(1,n):
        ans = ans + ((i*(n-i)*m*m)%M)
    for i in range(1,m):
        ans = ans + ((i*(m-i)*n*n)%M)
    print(ans*kumi%M)

if __name__ == "__main__":
    resolve()

Submission status: Accepted

Solution code:
import math

N, M, K = (int(x) for x in input().split())

def comb(n,k,p):
    if k < 0 or k > n:
        return 0
    if n == 0 or k == 0:
        return 1
    a=math.factorial(n) % p
    b=math.factorial(k) % p
    c=math.factorial(n-k) % p
    return (a*mypow(b,p-2,p)*mypow(c,p-2,p)) % p

def mypow(a,b,p):
    if b == 0: return 1
    if b % 2 == 0:
        d=mypow(a,b//2,p)
        return d*d % p
    if b % 2 == 1:
        return (a*mypow(a,b-1,p )) % p

c_sum = (M*(N-1)*N*(N+1)*M+N*(M-1)*M*(M+1)*N)//3
c_sum = c_sum//2
c_sum *= comb(N*M-2, K-2, 1000000007)
c_sum = c_sum % 1000000007

print(c_sum)

Submission status: Accepted

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

a = 0
for i in range(N):
    j = N-1-i
    a += i*(i+1)//2 + j*(j+1)//2

b = 0
for i in range(M):
    j = M-1-i
    b += i*(i+1)//2 + j*(j+1)//2

MAX = N*M+5
fac = [1,1] + [0]*MAX
finv = [1,1] + [0]*MAX
inv = [0,1] + [0]*MAX
for i in range(2,MAX+2):
    fac[i] = fac[i-1] * i % MOD
    inv[i] = -inv[MOD%i] * (MOD // i) % MOD
    finv[i] = finv[i-1] * inv[i] % MOD

def ncr(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

ans = (a * inv[N]**2 + b * inv[M]**2) * inv[N*M-1] * (N*M)
ans *= K*(K-1)//2
ans *= ncr(N*M, K)
print(ans % MOD)

Submission status: Accepted

Solution code:
def pmod(a,b,m):
  ret=1
  while b>0:
    if b&1:
      ret=(ret*a)%m
    a=(a*a)%m
    b=b>>1
  return ret
  
def fmod(n,m):
  ret=1
  for i in range(n,0,-1):
    ret=(ret*i)%m
  return ret

def cmod(n,r,m):
  ret=1
  for i in range(r):
    ret=(ret*(n-i))%m
  return (ret*pmod(fmod(r,m),m-2,m))%m

mod=10**9+7
n,m,k=map(int,input().split())
ans=0
c=cmod(n*m-2,k-2,mod)
for i in range(1,n):
  ans+=i*c*(m**2)*(n-i)
  ans=ans%mod
for i in range(1,m):
  ans+=i*c*(n**2)*(m-i)
  ans=ans%mod
print(ans%mod)

Submission status: Accepted

Solution code:

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

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

comb = cmb(N*M-2, K-2) % 1000000007

cost = 0
for d in range(1, N):
    cost += d * M**2 * (N - d) * comb
    cost %= 1000000007
for d in range(1, M):
    cost += d * N**2 * (M - d) * comb
    cost %= 1000000007
    
print(cost)

Submission status: Accepted

Solution code:
def factorials(M):#M!までのリストｘと其の逆元のリストyを返すmod INF
    INF=int(1e9+7)
    fac=[0]*(M+1)
    finv=[0]*(M+1)
    inv=[0]*(M+1)
    fac[0]=fac[1]=1
    finv[0]=finv[1]=1
    inv[1]=1
    for i in range(2,M+1):
        fac[i]=(fac[i-1]*i)%INF
        inv[i]=INF-(inv[INF%i]*(INF//i))%INF
        finv[i]=finv[i-1]*inv[i]%INF
    return fac,finv
def nCr(n,r):
    INF=int(1e9+7)
    if n<r or r<0 or n<0 :
        return 0
    return fac[n]*(finv[r]*finv[n-r]%INF)%INF

N,M,K=map(int,input().split())
MOD=10**9+7
fac,finv=factorials(N*M)
ans=0
for d in range(1,N):
    ans+=d*nCr(N*M-2,K-2)*(((N-d)*M**2)%MOD)
    ans%=MOD
for d in range(1,M):
    ans+=d*nCr(N*M-2,K-2)*(((M-d)*N**2)%MOD)
print(ans%MOD)

Submission status: Accepted

Solution code:
N,M,K=map(int,input().split())
mod=10**9+7
NM=N*M
kaijo=[0]*(NM+1)
gyaku=[0]*(NM+1)
kaijo[0]=1
for i in range(1,NM+1):
    kaijo[i]=kaijo[i-1]*i%mod
gyaku[NM]=pow(kaijo[NM],mod-2,mod)
for i in range(NM,0,-1):
    gyaku[i-1]=gyaku[i]*i%mod
C=kaijo[NM-2]*gyaku[K-2]*gyaku[NM-K]%mod
ans=0
for k in range(1,M):
    x=k*(M-k)*C*N*N%mod
    ans=(ans+x)%mod
for k in range(1,N):
    x=k*(N-k)*C*M*M%mod
    ans=(ans+x)%mod
print(ans)

Submission status: Accepted

Solution code:
def mod_combination(mod: int, n: int, r: int) -> int:
  a, b, c = mod_factorials(mod, [n, r, n - r])

  return a * pow(b, mod - 2, mod) % mod * pow(c, mod - 2, mod) % mod

def mod_factorials(mod: int, N: [int]) -> [int]:
  result = [1] * len(N)
  tuples = sorted([(i, n) for i, n in enumerate(N)], key=lambda t: t[1])

  fact = 1
  i = 1
  for t in tuples:
    for j in range(i, t[1] + 1):
      fact = fact * j % mod
      i = i + 1

    result[t[0]] = fact

  return result


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

  N, M, K = [int(i) for i in input().strip().split(' ')]
  c = mod_combination(MOD, N * M - 2, K - 2)

  X = 0
  Y = 0
  for i in range(M):
    X += (M - i) * i
  for i in range(N):
    Y += (N - i) * i

  total = (X * (N ** 2) + Y * (M ** 2)) % MOD * c % MOD
  print(total)


main()

Submission status: Accepted

Solution code:
from math import factorial


def comb(n, k, mod):
  if n < 0 or k < 0 or n < k:
      return 0
  if n == 0 or k == 0:
      return 1
  
  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


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


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

Submission status: Accepted

Solution code:
from math import factorial


def comb(n, k, mod):
  if n < 0 or k < 0 or n < k:
      return 0
  if n == 0 or k == 0:
      return 1
  
  a=factorial(n) % mod
  b=factorial(k) % mod
  c=factorial(n-k) % mod
  return (a*power_func(b, mod-2, mod)*power_func(c, mod-2, mod)) % mod


def power_func(a,b,mod):
  if b == 0:
      return 1
  if b % 2 == 0:
    d = power_func(a, b//2, mod)
    return d*d % mod
  if b % 2 == 1:
    return (a*power_func(a, b-1, mod)) % mod


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


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

Submission status: Accepted

Solution code:
MAX = 10**9 + 7

def mod_inv(i, mod):
    return pow(i, mod-2, mod)

def combination(_N, _K):
    if _K > (_N/2):
        _K = _N - _K
    up = 1
    low = 1
    for i in range (_K):
        up = up * (_N - i) % MAX
        low = low * (i + 1) % MAX
    ret = up * mod_inv(low, MAX) % MAX
    return ret

def cum_sum(N):
    return N * (N+1) // 2

def solve2():
    N, M, K = [int(_) for _ in raw_input().split()]
    sumx = 0
    for xi in range(1, N+1):
        sumx = (sumx + cum_sum(N) - (N + 1 - xi) * xi)
    sumx = (sumx * M * M) // 2 %MAX
    sumy = 0
    for yi in range(1, M+1):
        sumy = (sumy + cum_sum(M) - (M + 1 - yi) * yi)
    sumy = (sumy * N * N) // 2 %MAX
    ret = (sumx + sumy) * combination(M*N-2, K-2) % MAX
    print ret

if __name__ == '__main__':
    solve2()

Submission status: Accepted

Solution code:
MAX = 10**9 + 7

def mod_inv(i, mod):
    return pow(i, mod-2, mod)

def combination(_N, _K):
    if _K > (_N/2):
        _K = _N - _K
    up = 1
    low = 1
    for i in range (_K):
        up = up * (_N - i) % MAX
        low = low * (i + 1) % MAX
    ret = up * mod_inv(low, MAX) % MAX
    return ret

def cum_sum(N):
    return N * (N+1) // 2

def solve2():
    N, M, K = [int(_) for _ in raw_input().split()]
    sumx = 0
    for xi in range(1, N+1):
        sumx = (sumx + cum_sum(N) - (N + 1 - xi) * xi) % MAX
    sumx = (sumx * M * M) // 2 %MAX
    sumy = 0
    for yi in range(1, M+1):
        sumy = (sumy + cum_sum(M) - (M + 1 - yi) * yi) % MAX
    sumy = (sumy * N * N) // 2 %MAX
    ret = (sumx + sumy) * combination(M*N-2, K-2) % MAX
    print ret

if __name__ == '__main__':
    solve2()

Submission status: Accepted

Solution code:

MAX = 10**9 + 7

def mod_inv(i, mod):
    return pow(i, mod-2, mod)

def combination(_N, _K):
    if _K > (_N/2):
        _K = _N - _K
    up = 1
    low = 1
    for i in range (_K):
        up = up * (_N - i) % MAX
        low = low * (i + 1) % MAX
    ret = up * mod_inv(low, MAX) % MAX
    return ret

def solve():
    N, M, K = [int(_) for _ in raw_input().split()]
    disSum = M*N* (M*(N+1)*(N-1) + N*(M+1)*(M-1)) / 6
    ret = disSum * combination(M*N-2, K-2) % MAX
    print ret

if __name__ == '__main__':
    solve()

Submission status: Accepted

Solution code:
MAX = 10**9 + 7

def mod_inv(i, mod):
    return pow(i, mod-2, mod)

def combination(_N, _K):
    if _K > (_N/2):
        _K = _N - _K
    up = 1
    low = 1
    for i in range (_K):
        up = up * (_N - i) % MAX
        low = low * (i + 1) % MAX
    ret = up * mod_inv(low, MAX) % MAX
    return ret

def solve():
    N, M, K = [int(_) for _ in raw_input().split()]
    disSum = M*N* (M*(N+1)*(N-1) + N*(M+1)*(M-1)) // 6
    ret = disSum * combination(M*N-2, K-2) % MAX
    print ret

if __name__ == '__main__':
    solve()

Submission status: Accepted

Solution code:
MAX = 10**9 + 7

def mod_inv(i, mod):
    return pow(i, mod-2, mod)

def combination(_N, _K):
    if _K > (_N/2):
        _K = _N - _K
    up = 1
    low = 1
    for i in range (_K):
        up = up * (_N - i) % MAX
        low = low * (i + 1) % MAX
    ret = up * mod_inv(low, MAX) % MAX
    return ret

def solve():
    N, M, K = [int(_) for _ in raw_input().split()]
    disSum = (M*N* ((M*(N+1)*(N-1))%MAX + (N*(M+1)*(M-1)) % MAX)) % MAX
    disSum *= mod_inv(6, MAX) % MAX
    ret = disSum * combination(M*N-2, K-2) % MAX
    print ret

if __name__ == '__main__':
    solve()

Submission status: Accepted

Solution code:
import sys
def LI(): return [int(x) for x in sys.stdin.readline().split()]
def II(): return int(sys.stdin.readline())
def LS(): return sys.stdin.readline().split()
sys.setrecursionlimit(10**7)
INF = 10 ** 18
MOD = 10 ** 9 + 7
def LI_(): return [int(x) - 1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def SI(): return input()
YN = lambda b: print('YES') if b else print('NO')
yn = lambda b: print('Yes') if b else print('No')

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 )

def main():
    N, M, K = LI()

    rows = M*M*sum([d * (N-d) for d in range(1, N)])
    cols = N*N*sum([d * (M-d) for d in range(1, M)])
    ans = (rows + cols) * cmb(M*N-2, K-2, MOD) % MOD

    print(ans)

main()

Submission status: Accepted

Solution code:
import math

MOD = int(10**9 + 7)

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

def modinv(a, m):
    g, x, y = egcd(a, m)
    if g != 1:
        raise Exception("No inverse exists")
    else:
        return x % m

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

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

ans = 0

k = combination(N*M-2, K-2)

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

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

print(ans % MOD)

Submission status: Accepted

Solution code:
def mfactorial(n: int, mod: int) -> int:
  cache = mfactorial.__cache if hasattr(mfactorial, '__cache') else [1, 1]
  n_len = len(cache)
  n_new = n - n_len + 1
  if n_new > 0:
    cache = cache + ([0] * n_new)
    for i in range(n_len, n + 1):
      cache[i] = cache[i - 1] * i % mod
    mfactorial.__cache = cache

  return mfactorial.__cache[n]


def cmb(n, r, mod):
  result = mfactorial(n, mod)
  result *= pow(mfactorial(r, mod), mod - 2, mod) % mod
  result %= mod
  result *= pow(mfactorial(n - r, mod), mod - 2, mod) % mod
  result %= mod

  return result


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

  N, M, K = [int(i) for i in input().strip().split(' ')]
  c = cmb(N * M - 2, K - 2, MOD)

  X = 0
  Y = 0
  for i in range(M):
    X += (M - i) * i
  for i in range(N):
    Y += (N - i) * i

  total = (X * (N ** 2) + Y * (M ** 2)) % MOD * c % MOD
  print(total)


main()