Problem:
 Score : 500 points    Problem Statement We have a grid of squares with N rows and M columns. Let (i, j) denote the square at the i -th row from the top and j -th column from the left. We will choose K of the squares and put a piece on each of them.  If we place the K pieces on squares (x_1, y_1) , (x_2, y_2) , ..., and (x_K, y_K) , the cost of this arrangement is computed as:  \sum_{i=1}^{K-1} \sum_{j=i+1}^K (|x_i - x_j| + |y_i - y_j|)  Find the sum of the costs of all possible arrangements of the pieces. Since this value can be tremendous, print it modulo 10^9+7 .  We consider two arrangements of the pieces different if and only if there is a square that contains a piece in one of the arrangements but not in the other.      Constraints  2 \leq N \times M \leq 2 \times 10^5  2 \leq K \leq N \times M  All values in input are integers.         Input Input is given from Standard Input in the following format:  N  M  K       Output Print the sum of the costs of all possible arrangements of the pieces, modulo 10^9+7 .        Sample Input 1 2 2 2      Sample Output 1 8  There are six possible arrangements of the pieces, as follows:   ((1,1),(1,2)) , with the cost |1-1|+|1-2| = 1  ((1,1),(2,1)) , with the cost |1-2|+|1-1| = 1  ((1,1),(2,2)) , with the cost |1-2|+|1-2| = 2  ((1,2),(2,1)) , with the cost |1-2|+|2-1| = 2  ((1,2),(2,2)) , with the cost |1-2|+|2-2| = 1  ((2,1),(2,2)) , with the cost |2-2|+|1-2| = 1   The sum of these costs is 8 .       Sample Input 2 4 5 4      Sample Output 2 87210       Sample Input 3 100 100 5000      Sample Output 3 817260251  Be sure to print the sum modulo 10^9+7 .

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

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

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

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


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)

Solution (Python):
n,m,k=map(int,input().split())
ans=0
INF=10**9+7
if n<=m:
    M=-(m-1)*m*(2*m-1)//3+m*m*(m-1)
    Ms=m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
    M=M%INF
    Ms=Ms%INF
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+Ms
            else:
                d=j-i
                ans=ans+d*m**2+M
            ans=ans%INF
else:
    (n,m)=(m,n)
    M=-(m-1)*m*(2*m-1)//3+m*m*(m-1)
    Ms=m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
    M=M%INF
    Ms=Ms%INF
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+Ms
            else:
                d=j-i
                ans=ans+d*m**2+M
            ans=ans%INF
import math
ans=(math.factorial(n*m-2)//(math.factorial(n*m-k)*math.factorial(k-2)))*ans
ans=ans%INF
print(ans)

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

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


def modinv(x,p):
    g,a,b = exgcd(x,p)  
    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)

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

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

Solution (Python):

from collections import defaultdict
from collections import deque
from heapq import heappush, heappop
import sys
import math
import bisect
import random
import itertools
sys.setrecursionlimit(10**5)
stdin = sys.stdin
bisect_left = bisect.bisect_left
bisect_right = bisect.bisect_right
def LI(): return list(map(int, stdin.readline().split()))
def LF(): return list(map(float, stdin.readline().split()))
def LI_(): return list(map(lambda x: int(x)-1, stdin.readline().split()))
def II(): return int(stdin.readline())
def IF(): return float(stdin.readline())
def LS(): return list(map(list, stdin.readline().split()))
def S(): return list(stdin.readline().rstrip())
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = float('INF')


def A():
    a, b = LI()
    if a <= 5:
        print(0)
    elif 6 <= a <= 12:
        print(b // 2)
    else:
        print(b)
    return


def B():
    r, d, x = LI()
    for i in range(10):
        print(x * r - d)
        x = x * r - d
    return


def C():
    n, m = LI()
    imo = [0 for i in range(n+1)]
    for _ in range(m):
        l, r = LI_()
        imo[l] += 1
        imo[r + 1] -= 1
    for i in range(1,n):
        imo[i] += imo[i - 1]
    ans = 0
    
    for i in imo:
        if i == m:
            ans += 1
    print(ans)
    return


def D():
    n, m = LI()
    a = LI()
    bc = LIR(m)
    bc.sort(key=lambda x: x[1], reverse=True)
    a.sort()
    index = 0
    suma = 0
    for b, c in bc:
        if index + b >= n:
            if a[n - 1] <= c:
                suma += (n - index) * c
            else: 
                br = bisect_left(a, c)
                if br > index:
                    suma += (br - index) * c
                    index = br
                    break
                else:
                    break
            index = n
            break
        if a[index + b - 1] < c:
            suma += b * c
            index += b
        else:
            if a[index] < c:
                bl = bisect_left(a, c)
                suma += (bl - index) * c
                index = bl
            else:
                break
    suma += sum(a[index:])
    print(suma)

    return


def E():
    n, m, k = LI()
    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)
        
    ans = 0
    for d in range(1,n+1):
        ans += d*cmb(n * m - 2, k - 2, mod) * (n - d) * m ** 2 
    for d in range(1,m+1):
        ans += d*cmb(n * m - 2, k - 2, mod) * (m - d) * n ** 2
    print(ans%mod)
    

    
    return


def F():
    return


def G():
    return


def H():
    return


if __name__ == '__main__':
    E()

Solution (Python):
def factorial_mod(n, mod):
    a = 1
    for i in range(1, n + 1):
        a *= i
        a %= mod
    return a


def comb_mod(n, k, mod):
    fact_n = factorial_mod(n, mod)
    fact_k = factorial_mod(k, mod)
    fact_n_k = factorial_mod(n - k, mod)
    return (fact_n * pow(fact_k, mod - 2, mod) * pow(fact_n_k, mod - 2, mod)) % mod


n, m, k = map(int, input().split())
s = 0
mod = 10 ** 9 + 7
for x in range(m):
    s += (m - x - 1) * (m - x) // 2 * n **2
for y in range(n):
    s += (n - y - 1) * (n - y) // 2 * m ** 2


print(s * comb_mod(n * m - 2, k - 2, mod) % mod)

Solution (Python):
import math
import sys
sys.setrecursionlimit(10 ** 6 + 10)
from collections import defaultdict

from math import factorial

def combinations_count(n, r):
    return factorial(n) // (factorial(n - r) * factorial(r))

def permutations_count(n, r):
    return math.factorial(n) // math.factorial(n - r)



def main():
    mod = 10 ** 9 + 7
    high, width, num = map(int, input().split())
    toori = combinations_count(high * width - 2, num - 2)

    kyori = high + width - 2

    ans = 0
    for now_kyori in range(1, kyori + 1):
        kari = 0
        if now_kyori < high:
            kari += width * width * (high - now_kyori)
        if now_kyori < width:
            kari += high * high * (width - now_kyori)
        ans += now_kyori * kari

        if now_kyori >= high and now_kyori >= width:
            break

    ans = ans * toori % mod

    print(ans)



if __name__ == '__main__':
    main()

Solution (Python):
import sys
import bisect
from collections import defaultdict

ans=0

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



max_value = n*m+2
mod = 10 ** 9 + 7

fac = [0 for _ in range(max_value)]
finv = [0 for _ in range(max_value)]
inv = [0 for _ in range(max_value)]

def com_init():
    fac[0] = 1
    fac[1] = 1
    finv[0] = 1
    finv[1] = 1
    inv[1] = 1

    for i in range(2, max_value):
        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

def dist(a,b):
    if a<=0 or b<=0:
        return 0
    
    a,b=min(a,b),max(a,b)
    a,b=a%mod,b%mod
    '''
    rt=((a*(a+1))%mod)*(a+b)//2
    rt+=((a*(a+b))%mod)*(b-a-1) if a!=b else 0
    rt+= - a*a if a==b else 0
    '''
    rt=((b*a*(a+b))//2)%mod
    
    return rt%mod

com_init()
inv2=500000004

d=0
for i in range(n):
    for j in range(m):
        d+=dist(j,i+1)+dist(i,m-j)+dist(j+1,n-i-1)+dist(n-i,m-j-1)
        d%=mod
ans=d*com(n*m-2,k-2)
ans*=inv2

print(ans%mod)

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

MOD = 10**9 + 7
R = 1
for i in range(K - 2):
    R = (R * (N * M - 2 - i) * pow(i + 1, MOD - 2, MOD)) % MOD

ans = 0
for i in range(1, N + 1):
    for j in range(1, M + 1):
        u = (N - i + 1)
        d = i
        l = j
        r = (M - j + 1)
        D = 0
        D += u * l * (u + l - 2) // 2
        D += u * r * (u + r - 2) // 2
        D += d * l * (d + l - 2) // 2
        D += d * r * (d + r - 2) // 2
        D -= u * (u - 1) // 2
        D -= d * (d - 1) // 2
        D -= l * (l - 1) // 2
        D -= r * (r - 1) // 2
        ans += D

ans = (ans * R * pow(2, MOD - 2, MOD)) % MOD
print(ans)

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

MOD = 10**9 + 7
L = N*M

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

L = max(N, M)
D = [0]*(L+1)
for i in range(L):
    D[i+1] = D[i] + (i+1)

ans = 0
for i in range(N):
    for j in range(M):
        ans += D[j] * (i+1) + D[M-1-j] * i + D[i] * M
    ans %= MOD
print(ans * fact[N*M-2] * rfact[N*M-K] * rfact[K-2] % MOD)

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

A=2*10**5
mod = 10**9+7
g1 = [1, 1]
g2 = [1, 1]
inverse = [0, 1]

for i in range( 2, A + 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())
keisuu=cmb(N*M-2,K-2,mod)

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

Solution (Python):
import math
import sys

sys.setrecursionlimit(10000)
INF = float('inf')


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


def comb(n, r, mod):
    
    return (math.factorial(n) * inv(math.factorial(n - r)) * inv(math.factorial(r))) % mod


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


def pyramid(n):
    return n * (n + 1) * (n + 2) // 6


MOD = 10 ** 9 + 7



k2 = (pyramid(N - 1) * M ** 2 + pyramid(M - 1) * N ** 2) % MOD

c = comb(N * M - 2, K - 2, mod=MOD)
print(k2 * c % MOD)

Solution (Python):
import math
import sys

sys.setrecursionlimit(10000)
INF = float('inf')


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


def comb(n, r, mod):
    
    return (math.factorial(n) * inv(math.factorial(n - r)) * inv(math.factorial(r))) % mod


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


def pyramid(n):
    return n * (n + 1) * (n + 2) // 6


MOD = 10 ** 9 + 7



k2 = (pyramid(N - 1) * M ** 2 + pyramid(M - 1) * N ** 2) % MOD

c = comb(N * M - 2, K - 2, mod=MOD)
print(k2 * c % MOD)

Solution (Python):
def getInv(N):
    inv = [0] * (N + 1)
    inv[0] = 1
    inv[1] = 1
    for i in range(2, N + 1):
        inv[i] = (-(Q // i) * inv[Q%i]) % Q
    return inv

Q = 10**9+7
modinv = getInv(4*10**5)
N, M, K = map( int, input().split())
T = N*M
comb = 1
for i in range(K-2):
    comb *= T-2-i
    comb %= Q
    comb *= modinv[K-2-i]
    comb %= Q
ans = 0
for i in range(1, M):
    ans += N**2*(M-i)*i
    ans %= Q
for i in range(1, N):
    ans += M**2*(N-i)*i
    ans %= Q
ans %= Q
print( ans*comb%Q)

Solution (Python):
def mod_inv(n: int, mod: int)->int:
    b, u, v = mod, 1, 0
    while b > 0:
        t = n // b

        n -= t * b
        u -= t * v

        n, b = b, n
        u, v = v, u

    return (u+mod) % mod


def cell_distance(N: int, M: int, K: int)->int:
    MOD = 10**9 + 7

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

    S = 0
    for d in range(N):
        S += d * (N-d) * M * M
        S %= MOD

    for d in range(M):
        S += d * (M-d) * N * N
        S %= MOD

    S *= C
    S %= MOD

    return S


if __name__ == "__main__":
    N, M, K = map(int, input().split())

    ans = cell_distance(N, M, K)
    print(ans)

Solution (Python):
n,m,k=map(int,input().split())
ans=0
INF=10**9+7
if n<=m:
    M=-(m-1)*m*(2*m-1)//3+m*m*(m-1)
    Ms=m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
    M=M%INF
    Ms=Ms%INF
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+Ms
            else:
                d=j-i
                ans=ans+d*m**2+M
            ans=ans%INF
else:
    (n,m)=(m,n)
    M=-(m-1)*m*(2*m-1)//3+m*m*(m-1)
    Ms=m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
    M=M%INF
    Ms=Ms%INF
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+Ms
            else:
                d=j-i
                ans=ans+d*m**2+M
            ans=ans%INF
import math
ans=(math.factorial(n*m-2)//(math.factorial(n*m-k)*math.factorial(k-2)))*ans
ans=ans%INF
print(ans)

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

Solution (Python):
n,m,k=map(int,input().split())
ans=0
INF=10**9+7
if n<=m:
    M=(m-1)*m*(2*m-1)//3
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
            else:
                d=j-i
                ans=ans+d*m+(m+d)*m*(m-1)-M
            ans=ans%INF
else:
    (n,m)=(m,n)
    M=(m-1)*m*(2*m-1)//3
    for i in range(1,n+1):
        for j in range(i,n+1):
            if i==j:
                ans=ans+m*m*(m-1)//2-(m-1)*m*(2*m-1)//6
            else:
                d=j-i
                ans=ans+d*m+(m+d)*m*(m-1)-M
            ans=ans%INF
import math
ans=(math.factorial(n*m-2)//(math.factorial(n*m-k)*math.factorial(k-2)))*ans
ans=ans%INF
print(ans)

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

Solution (Python):
def dist(arr, n):
	res = 0; su = 0
	for i in range(n):
		res = (res + (arr[i] * i) - su) % mod
		su = (su + arr[i]) % mod

	return res % mod

n,m,k = map(int,input().split())
mod = int(1e9+7)
fact = [1]
for i in range(1,n*m+1):
	fact.append((fact[-1] * i)%mod)

x = []
y = []

for i in range(1,n+1):
	for j in range(1,m+1):
		x.append(i)
		y.append(j)

y.sort()

tot = (dist(x,n*m) + dist(y,n*m))%mod
N = n*m - 2
Need = k-2
num = fact[N]
den = (fact[Need] * fact[N-Need])%mod
ncr = (num * pow(den,mod-2,mod)) % mod
ans = (tot * ncr) % mod
print(ans)

Solution (Python):
n,m,k = map(int,input().split())
x = m*n
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
ss = cmb(x,k)
s = (ss*k*(k-1))//2
r = 10**9 +7
av = (((m**3-m) *n**2)//6 + ((n**3-n)*m**2)//6)%r
d = (x*(x-1))//2
print(int(((av*s)//d)%r))

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

mod = 10**9 + 7

def mpow(x, n):
  ans = 1
  while n != 0:
    if n&1:
      ans = ans * x % mod
    x = x * x % mod
    n = n // 2
  return ans

fac = [0] * (n*m+2)
finv = [0] * (n*m+2)
inv = [0] * (n*m+2)
fac[0] = 1
fac[1] = 1
finv[0] = 1
finv[1] = 1
inv[1] = 1

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

def choose(n,k):
  if n == 0 and k == 0:
    return 1
  if n < 0:
    return 0
  if n < k:
    return 0
  return fac[n] * (finv[k] * finv[n-k] % mod) % mod


sum = 0
for i in range(1, n+1):
  for j in range(1, m+1):
    score = i-1 + j-1
    cnt = (n-i+1) * (m-j+1)
    sum += (score * cnt) % mod
    

for i in range(2, n+1):
  for j in range(1, m):
    score = i-1 + m-j
    cnt = (n-i+1) * j
    sum += (score * cnt) % mod
    

sum = sum * choose(m*n-2,k-2)
sum = sum % mod
print(sum)

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

def P(n, r):
    return math.factorial(n)//math.factorial(n-r)

def C(n, r):
    return P(n, r)//math.factorial(r)


dis = M * M * N * N * (N+1) // 2 % sosu
dis = dis - M * M * N * (N+1) * ( 2* N + 1) // 6 % sosu
dis = dis + N * N * M * M * (M+1) //2 % sosu
dis = dis - N * N * M * (M+1) * (2*M + 1) // 6 % sosu

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 = sosu 
N = N * M
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 )


if K >= 3:
    dis = dis * cmb(N-2, K-2, mod) % sosu
if dis >= 0:
    print(dis)
else:
    print(dis + sosu)

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

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

def power(n, r, mod):
    if r == 0: return 1
    if r%2 == 0:
        return power(n*n % mod, r//2, mod) % mod
    if r%2 == 1:
        return n * power(n, r-1, mod) % mod

def comb(n, k, mod):
    a=framod(n, mod)
    b=framod(k, mod)
    c=framod(n-k, mod)
    return (a * power(b, mod-2, mod) * power(c, mod-2, mod)) % mod


ans = 0

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

c = comb(N*M-2,K-2,MOD)

for i in range(N*M):
    tate, yoko = i//M, i%M
    migi = f(M-1-yoko)
    hidari = f(yoko)
    sita = (hidari+migi)*(N-1-tate) + M*f(N-1-tate)
    coef = migi + sita
    ans += coef*c
    ans %= MOD
print(ans)

Solution (Python):
def init():
    global factorial
    MOD = 10 ** 9 + 7
    factorial = [0] * 200001
    factorial[0] = 1
    for i in range(1, 200001):
        factorial[i] = (factorial[i - 1] * i) % MOD

def main():
    init()
    n, m, k = map(int, input().split())
    s = 0
    MOD = 10 ** 9 + 7
    for i in range(n):
        s += (((((n - 1 - i) * (n - i)) // 2) % MOD) * ((m * m) % MOD)) % MOD
    for i in range(m):
        s += (((((m - 1 - i) * (m - i)) // 2) % MOD) * ((n * n) % MOD)) % MOD
    print((((s * ((factorial[n * m - 2] * pow(factorial[k - 2], MOD - 2, MOD)) % MOD)) % MOD) * pow(factorial[n * m - k], MOD - 2, MOD)) % MOD)

main()

Solution (Python):

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)
    d = Combination(N*M, 10**9+7)(N*M-2, K-2)

    ans = 0
    for i in range(1, N+1):
        for j in range(1, M+1):
            if i == N and j == M:
                continue
            if i == N or j == M:
                a = 1
            else:
                a = 2
            ans = m.add(ans, m.mul(i*j*a, m.mul(i + j, d)))
    ans = m.div(ans, 2)

    return ans


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


if __name__ == '__main__':
    main()

Solution (Python):
from functools import reduce

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('modular inverse does not exist')
    else:
        return x % m

MOD = 10**9+7

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

patterns = reduce(lambda a,b: (a*b)%MOD, range(N*M-K+1, N*M+1))
patterns *= modinv(reduce(lambda a,b: (a*b)%MOD, range(1, K+1)), MOD)

patterns *= (N+M)*K*(K-1)
patterns *= modinv(6, MOD)
patterns %= MOD

print(patterns)

Solution (Python):
p =10**9 +7
def comb(n,k,p):
  """power_funcを用いて(nCk) mod p を求める"""
  from math import factorial
  if n<0 or k<0 or n<k: return 0
  if n==0 or k==0: return 1
  a=factorial(n) %p
  b=factorial(k) %p
  c=factorial(n-k) %p
  return (a*power_func(b,p-2,p)*power_func(c,p-2,p))%p
def power_func(a,b,p):
  """a^b mod p を求める"""
  if b==0: return 1
  if b%2==0:
    d=power_func(a,b//2,p)
    return d*d %p
  if b%2==1:
    return (a*power_func(a,b-1,p ))%p
N,M,K =map(int,input().split())
def keisan(a,b):
    foo = a*(a**2 +3*a +2)//6
    foo *= b**2
    foo = foo%p
    return foo
ans = (keisan(N-1,M) + keisan(M-1,N))%p
foo = comb(N*M-2,K-2,p)
bar =9*17
bar =bar%p
ans = (ans*foo)%p
print(ans)

Solution (Python):


MOD=10**9+7
fact=[1]
infact=[1]
temp=1
intemp=1
for i in xrange(1,10**6):
    temp*=i
    temp%=MOD
    fact+=[temp]
    intemp*=pow(i,MOD-2,MOD)
    intemp%=MOD
    infact+=[intemp]
    

def binom(a,b):
    up=fact[a]
    down=(infact[b]*infact[a-b])%MOD
    ans=(up*down)%MOD
    return ans


N,M,k=list(map(int,raw_input().strip().split(' ')))
if N==1 or M==1:
    if M<N:
        N,M=M,N
    ans=(binom(M*N,k)*binom(k,2))%MOD
    ans*=binom(M+1,3)
    ans%=MOD
    ans*=pow(binom(M*N,2),MOD-2,MOD)
    ans%=MOD
    print(ans)

else:       


    ans=(binom(M*N,k)*binom(k,2))%MOD
    ans*=((M*M*binom(N+1,3))%MOD+(N*N*binom(M+1,3))%MOD)%MOD
    ans*=pow(binom(M*N,2),MOD-2,MOD)
    ans%=MOD
    print(ans)

Solution (Python):
def mod_exp(b,e,mod):
    r = 1
    while e > 0:
        if (e&1) == 1:
            r = (r*b)%mod
        b = (b*b)%mod
        e >>= 1

    return r


def fact_exp(n,p):
    e = 0
    u = p
    t = n
    while u <= t:
        e += t//u
        u *= p
    return e



def get_base_digits(n,b):
    d = []
    while n > 0:
        d.append(n % b)
        n  = n // b
    return d




def fermat_binom_advanced(n,k,p):
    
    num_degree = fact_exp(n,p) - fact_exp(n-k,p)
    den_degree = fact_exp(k,p)
    if num_degree > den_degree:
        return 0
    if k > n:
        return 0
    
    num = 1
    for i in range(n,n-k,-1):
        cur = i
        while cur%p == 0:
            cur //= p
        num = (num*cur)%p
    
    denom = 1
    for i in range(1,k+1):
        cur = i
        while cur%p == 0:
            cur //= p
        denom = (denom*cur)%p
    
    return (num * mod_exp(denom,p-2,p))%p



def lucas_binom(n,k,p):
    
    np = get_base_digits(n,p)
    kp = get_base_digits(k,p)
    
    binom = 1
    for i in range(len(np)-1,-1,-1):
        ni = np[i]
        ki = 0
        if i < len(kp):
            ki = kp[i]
        binom = (binom * fermat_binom_advanced(ni,ki,p)) % p
    return binom

mod = 10**9+7

n,m,k = map(int,input().split())
temp = 0
for d in range(1,m):
	temp += d*n*n*(m-d)
for d in range(1,n):
	temp += d*m*m*(n-d)
l = lucas_binom(n*m-2,k-2,mod)
print(temp*l%mod)