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):
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
 
list=input().split(" ")
n=int(list[0])
m=int(list[1])
k=int(list[2])
 
niko1=int((n*(n+1)*(2*n+1)/6 - (n+1)*n*(n+1)/2 + n*n*(n+1)/2)/2)
niko2=int((m*(m+1)*(2*m+1)/6 - (m+1)*m*(m+1)/2 + m*m*(m+1)/2)/2)
niko1=niko1%(10**9+7)
niko2=niko2%(10**9+7)
niko1=niko1*m
niko2=niko2*n
niko1=niko1%(10**9+7)
niko2=niko2%(10**9+7)
niko1=niko1*m
niko2=niko2*n
niko1=niko1%(10**9+7)
niko2=niko2%(10**9+7)
niko=niko1+niko2
kosu=cmb(m*n-2, k-2)
kosu=kosu%(10**9+7)
kotae=int(niko*kosu)
kotae=kotae%(10**9+7)

print(kotae)

Solution (Python):
import sys
range = xrange
input = raw_input

MOD = 10**9+7

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

def nck(n,k):
    a = 1
    b = 1
    c = 1
    for i in range(1,n+1):
        a = a*i%MOD
    for i in range(1,k+1):
        b = b*i%MOD
    for i in range(1,n-k+1):
        c = c*i%MOD
    return a*modinv(b*c%MOD)%MOD
        
n,m,k = [int(x) for x in input().split()]

inv = modinv(n*m)
m1inv = modinv(n*m-1)

ans = 0

def summer(x):
    return x*(x-1)//2

p = (k*inv%MOD)*(m1inv*(k-1)%MOD)%MOD

xfac = n*n*p%MOD
for xi in range(m):
    ans += xfac * summer(m-xi)%MOD

yfac = m*m*p%MOD
for yi in range(n):
    ans += yfac * summer(n-yi)%MOD




        
        
        
        
        
        
        
        
        
        
        
        
print nck(n*m,k)*ans%MOD

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

tn = 0
tm = 0
for n in range(N):
  tn = (tn + n*(n+1)//2)%mod
tn = tn*M*M %mod
for m in range(M):
  tm = (tm + m*(m+1)//2)%mod
tm = tm*N*N %mod


class combinationmod():
  def __init__(self,MAX, mod):
    self.mod =mod
    self.fac=[0]*(MAX+1)
    self.finv=[0]*(MAX+1)
    self.inv=[0]*(MAX+1)
    self.fac[0]=self.fac[1]=1
    self.finv[0]=self.finv[1]=1
    self.inv[1]=1
    for i in range(2,MAX+1):
      self.fac[i] = self.fac[i-1]*i%mod
      self.inv[i] = (-self.inv[mod%i]*(mod//i))%mod
      self.finv[i] = self.finv[i-1]* self.inv[i] %mod

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

c = combinationmod(200000,int(1e9+7))


ans = (tn + tm)*c.nCk(N*M-2, K-2) %mod
print(ans)

Solution (Python):





import sys

n,m,k = [int(i) for i in input().split()]
MOD=10**9+7
SIZE = 2*(10**5)
inv = [0]*(SIZE+1)
fac = [0]*(SIZE+1)
finv = [0]*(SIZE+1)
inv[1] = 1
fac[0] = fac[1] = 1
finv[0] = finv[1] = 1
for i in range(2,SIZE+1):
    inv[i] = MOD -(MOD//i)*inv[MOD%i]%MOD
    fac[i] = fac[i-1]*i%MOD
    finv[i]= finv[i-1]*inv[i]%MOD

def choose(n,r): 
    if 0 <= r <= n:
        return (fac[n]*finv[r]%MOD)*finv[n-r]%MOD
    else:
        return 0


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

for i in range(m):
    ans1 += i*(m-i)
    ans1 %= MOD
ans1 *= n
ans1 %= MOD
ans1 *= n
ans1 %= MOD

ans = (ans0+ans1)%MOD


print(ans*choose(n*m-2,k-2)%MOD)

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

MOD = int(1e9 + 7)
NM = N * M

def su(n):
  return n * (n + 1) // 2
class comb():
  F = [1, 1]
  Fi = [1, 1]
  I = [0, 1]
  def __init__(self, num, mod):
    self.MOD = mod
    for i in range(2, num + 1):
      self.F.append((self.F[-1] * i) % mod)
      self.I.append(mod - self.I[mod % i] * (mod // i) % mod)
      self.Fi.append(self.Fi[-1] * self.I[i] % mod)
  def com(self, n, k):
    if n < k: return 0
    if n < 0 or k < 0: return 0
    return self.F[n] * (self.Fi[k] * self.Fi[n - k] % self.MOD) % self.MOD

co = comb(NM, MOD)
times = co.com(NM - 2, K - 2)
Ans = 0
for i in range(N):
  for j in range(M):
    k = (su(i) + su(N - i - 1)) * M
    k = k + (su(j) + su(M - j - 1)) * N
    k = k * times
    Ans += k

print(Ans // 2 % MOD)

Solution (Python):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,copy,functools
import random

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return list(map(int, sys.stdin.readline().split()))
def LLI(): return [list(map(int, l.split())) for l in sys.stdin.readlines()]
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 LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)
def pe(s): return print(str(s), file=sys.stderr)

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


def comb(n, b):
    if b > n - b:
        b = n - b
    r = 1
    for k in range(n, n-b, -1):
        r = r * k % mod
    d = 1
    for k in range(1, b+1):
        d = d * k % mod
    return r * inv(d) % mod

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

    r = comb(n*m,k)
    cc = comb(n*m,2)
    tt = 0
    m2 = m * m
    n2 = n * n
    for i in range(n):
        v = (n-i) * (n-i-1) // 2
        tt += v * m2 % mod
        tt %= mod
    for i in range(m):
        v = (m-i) * (m-i-1) // 2
        tt += v * n2 % mod
        tt %= mod

    rr = r * k * (k-1) // 2 * tt
    rr = rr * inv(cc) % mod

    return rr%mod


print(main())

Solution (Python):
import sys
from collections import deque
from collections import defaultdict
import math
sys.setrecursionlimit(20000000)
input = sys.stdin.readline
n,m,k = map(int,input().split())
mod = 10**9 + 7
fact = [1]*(n*m+10000+1)
rfact = [1]*(n*m+10000+1)
for i in range(n*m+10000):
    fact[i+1] = rr = ((i+1) * fact[i]) % mod
def comb(N, K,mod):
    return fact[N] * pow(fact[K], mod-2, mod) * pow(fact[N-K], mod-2, mod) % mod
ans = 0
ni = comb(n*m-2,k-2,mod)
for i in range(n):
    for j in range(m):
        ko = (1+j)*j*n//2
        ko += (m-j)*(m-1-j)*n//2
        ko += (1+i)*i*m//2
        ko += (n-i)*(n-1-i)*m//2
        ans += ni*ko
        ans %= mod
print(ans*pow(2,mod-2,mod)%mod)

Solution (Python):
from collections import defaultdict,deque
import sys,heapq,bisect,math,itertools,string,queue,datetime
sys.setrecursionlimit(10**8)
INF = float('inf')
mod = 10**9+7
eps = 10**-7
def inpl(): return list(map(int, input().split()))
def inpl_str(): return list(input().split())


H,W,K = inpl()

MAX = H*W+5
fac = [1]*(MAX+1)
for i in range(1,MAX+1):
	fac[i] = (fac[i-1]*i)%mod

gyakugen = [1]*(MAX+1)
gyakugen[MAX] = pow(fac[MAX],mod-2,mod)
for i in range(MAX,0,-1):
	gyakugen[i-1] = (gyakugen[i]*i)%mod

def Comb(n,k):
	return (fac[n]*gyakugen[k]*gyakugen[n-k])%mod

ans = 0
ruiseki = [0]
tmp = 0
for i in range(1,max(H,W)+5):
    tmp += i
    ruiseki.append(tmp)

for y in range(H):
    ans += (ruiseki[y] + ruiseki[H-y-1])*W*W

for x in range(W):
    ans += (ruiseki[x] + ruiseki[W-x-1])*H*H

ans //= 2
ans *= Comb(H*W-2,K-2)

print(ans%mod)

Solution (Python):
N, M, K = map(int, input().split())
P = 10**9+7
fact = [1]
for i in range(N*M):
    fact.append(fact[-1]*(i+1)%P)
def inv(n):
    return pow(n, P-2, P)
def C(a, b):
    return fact[a]*inv(fact[a-b]*fact[b])%P
print((M**2*(N**3-N)+N**2*(M**3-M))//6*C(M*N-2, K-2)%P)

Solution (Python):
n, m, k = list(map(int, input().split()))
MOD = 10 ** 9 + 7

t = 0
pre = [0]
for i in range(1, max(n, m) + 1):
    t += i
    pre.append(t)


def calc_base(n, m, pre):
    base = 0
    for i in range(n // 2):
        base = base + (pre[i] + pre[n - i - 1]) * 2 % MOD
    if n % 2 == 1:
        i = n // 2
        base = (base + pre[i] * 2) % MOD
    return base * m * m % MOD


bi = calc_base(n, m, pre)
bj = calc_base(m, n, pre)


def prepare(n, MOD):
    f = 1
    for m in range(1, n + 1):
        f *= m
        f %= MOD
    fn = f

    inv = pow(f, MOD - 2, MOD)
    invs = [1] * (n + 1)
    invs[n] = inv
    for m in range(n, 1, -1):
        inv *= m
        inv %= MOD
        invs[m - 1] = inv

    return fn, invs


nm = n * m - 2
f, inv = prepare(nm, MOD)
print((bi + bj) * f * inv[k - 2] % MOD * inv[nm - k + 2] * inv[2] % MOD)

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

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

inv = [0]*(k+3)
inv[1] = 1
for i in range(2,k+2):
  
  inv[i] = ((Mod-Mod//i) * inv[Mod%i]) % Mod


s = (n*n*(m+1)*m*(m-1) + m*m*(n+1)*n*(n-1))//6
for i in range(k-2):
   s = (s*(n*m-2-i)*inv[i+1])%Mod
print(s)

Solution (Python):
import sys

stdin = sys.stdin

mod = 1000000007

ni = lambda: int(ns())
na = lambda: list(map(int, stdin.readline().split()))
ns = lambda: stdin.readline().rstrip()  

n, m, k = na()

def invl(a, mod):
    b = mod
    p = 1; q = 0
    while b > 0:
        c = a // b
        a, b = b, a%b
        p, q = q, p-c*q
    return p + mod if p < 0 else p

num = 1
den = 1
for u in range(k-2):
    num = num * (n*m-2-u) % mod
for u in range(k-2):
    den = den * (1+u) % mod
C = num * invl(den, mod) % mod






v = n*(n+1)*(n-1)//3*m*m + m*(m+1)*(m-1)//3*n*n
v *= C
print(v*invl(2,mod)%mod)

Solution (Python):
from collections import defaultdict as ddic, deque, Counter
from itertools import permutations, combinations, product
import bisect, heapq

rr = raw_input
rri = lambda: int(raw_input())
rrm = lambda: map(int, raw_input().split())
MOD = 10**9 + 7

def binom(n, k):
    numer = 1
    for i in xrange(1, n+1):
        numer *= i
        numer %= MOD
    denom = 1
    for i in xrange(1, k+1):
        denom *= i
        denom %= MOD
    denom2 = 1
    for i in xrange(1, n-k+1):
        denom2 *= i
        denom2 %= MOD
    ans = numer * pow(denom, MOD-2, MOD) % MOD
    ans *= pow(denom2, MOD-2, MOD)
    ans %= MOD
    return ans

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

def flat(R, C):
    if C > R: R,C = C,R
    
    
    
    ans = 0
    for c in xrange(C):
        left = c
        right = C - 1 - c
        base = sumToN(left) + sumToN(right)
        for r in xrange(R):
            up, down = r, R-1-r
            base2 = sumToN(up) + sumToN(down)
            ans += base * R + base2 * C
    ans /= 2
    return ans % MOD


def solve(R, C, K):
    N = R * C
    return flat(R, C) * binom(N - 2, K - 2) % MOD

print solve(*rrm())

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

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

U = 2*10**5
fact = [1] * (U+1)
fact_inv = [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)
  x %= MOD
  return x

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

MM = M*M%MOD
NN = N*N%MOD
answer = MM*F(N)%MOD + NN*F(M)%MOD
answer %= MOD
answer *= comb(N*M-2,K-2)
answer %= MOD
print(answer)