CodeEff/ECCO huggingface dataset:
# Original solution
import sys

def isort(a, n, g):
    global cnt
    for i in range(g, n):
        v = a[i]
        j = i - g
        while j >= 0 and a[j] > v:
            a[j + g] = a[j]
            j -= g
            cnt += 1
        a[j + g] = v

# input data
n = int(sys.stdin.readline())
a = []
for i in range(n):
    a.append(int(sys.stdin.readline()))

# determine m, g
g = [1]
m = 1
while True:
    x = 3 * g[m - 1] + 1
    if x >= n: break
    g.append(x)
    m += 1
g = g[::-1]

# 1 output m
print(m)
# 2 output g
print((' '.join(map(str, g))))

# sort it
cnt = 0
for i in range(m):
    isort(a, n, g[i])

# 3 output cnt
print(cnt)
# 4 output sorted a
print(('\n'.join(map(str, a))))

# Optimized solution
def isort(a, n, g):
    global cnt
    for i in range(g, n):
        v = a[i]
        j = i - g
        while j >= 0 and a[j] > v:
            a[j + g] = a[j]
            j -= g
            cnt += 1
        a[j + g] = v

# input data
n = int(input())
a = []
for i in range(n):
    a.append(int(input()))

# determine m, g
g = [1]
m = 1
while True:
    x = 3 * g[m - 1] + 1
    if x >= n: break
    g.append(x)
    m += 1
g = g[::-1]

# 1 output m
print(m)
# 2 output g
print((' '.join(map(str, g))))

# sort it
cnt = 0
for i in range(m):
    isort(a, n, g[i])

# 3 output cnt
print(cnt)
# 4 output sorted a
print(('\n'.join(map(str, a))))
# Original solution
swap_count = 0


def main():
  num = int(eval(input()))
  elements = []
  for i in range(0, num):
    elements.append(int(eval(input())))

  # Get intervals
  interval = 1
  intervals = [interval]
  while True:
    interval = int_f(interval)
    if interval > num:
      break
    intervals.append(interval)
  intervals.reverse()

  for interval in intervals:
    insertion_sort(elements, num, interval)

  print((len(intervals)))
  print((' '.join([str(interval) for interval in intervals])))
  print(swap_count)
  for val in elements:
    print(val)


def int_f(x):
  return 3 * x + 1


def insertion_sort(elements, num, interval=1):
  global swap_count

  i = interval
  while i < num:
    j = i
    while interval <= j:
      if elements[j - interval] <= elements[j]:
        break

      elements[j - interval], elements[j] = elements[j], elements[j - interval]
      swap_count += 1
      j -= interval
    i += 1


if __name__ == '__main__':
  main()

# Optimized solution
swap_count = 0


def main():
  num = int(eval(input()))
  elements = []
  for i in range(0, num):
    elements.append(int(eval(input())))

  # Get intervals
  interval = 1
  intervals = [interval]
  while True:
    interval = int_f(interval)
    if interval > num:
      break
    intervals.insert(0, interval)
  # intervals.reverse()

  for interval in intervals:
    insertion_sort(elements, num, interval)

  print((len(intervals)))
  print((' '.join([str(interval) for interval in intervals])))
  print(swap_count)
  for val in elements:
    print(val)


def int_f(x):
  return 3 * x + 1


def insertion_sort(elements, num, interval=1):
  global swap_count

  i = interval
  while i < num:
    j = i
    while interval <= j:
      if elements[j - interval] <= elements[j]:
        break

      elements[j - interval], elements[j] = elements[j], elements[j - interval]
      swap_count += 1
      j -= interval
    i += 1


if __name__ == '__main__':
  main()
# Original solution
import sys
import math
from collections import deque

def insertion_sort(a, n, g):
    ct = 0
    for i in range(g,n):
        v = a[i]
        j = i-g
        while j >= 0 and a[j] > v:
            a[j+g] = a[j]
            j = j-g
            ct += 1
        a[j+g] = v
    return ct

n = int(input())
a = list(map(int, sys.stdin.readlines()))
b = 701
ct= 0
g = deque([x for x in [701,301,132,57,23,10,4,1] if x <= n])
while True:
    b = math.floor(2.25*b)
    if b > n:
        break
    g.appendleft(b)

for i in g:
    ct += insertion_sort(a, n, i)

print(len(g))
print(*g, sep=" ")
print(ct)
print(*a, sep="\n")


# Optimized solution
import sys
import math

def insertion_sort(a, n, g):
    ct = 0
    for i in range(g,n):
        v = a[i]
        j = i-g
        while j >= 0 and a[j] > v:
            a[j+g] = a[j]
            j = j-g
            ct += 1
        a[j+g] = v
    return ct

n = int(input())
a = list(map(int, sys.stdin.readlines()))
b = 701
ct= 0
g = [x for x in [1,4,10,23,57,132,301,701] if x <= n]
while True:
    b = math.floor(2.25*b)
    if b > n:
        break
    g.append(b)
g = g[::-1]

for i in g:
    ct += insertion_sort(a, n, i)

print(len(g))
print(*g, sep=" ")
print(ct)
print(*a, sep="\n")
# Original solution
import sys
import math

def insertion_sort(a, n, g):
    ct = 0
    for j in range(0,n-g):
        v = a[j+g]
        while j >= 0 and a[j] > v:
            a[j+g] = a[j]
            j = j-g
            ct += 1
        a[j+g] = v
    return ct

n = int(input())
a = list(map(int, sys.stdin.readlines()))
b = 701
ct= 0
g = [x for x in [1,4,10,23,57,132,301,701] if x <= n]
while True:
    b = math.floor(2.25*b)
    if b > n:
        break
    g.append(b)
g = g[::-1]

for i in g:
    ct += insertion_sort(a, n, i)

print(len(g))
print(*g, sep=" ")
print(ct)
print(*a, sep="\n")


# Optimized solution
import sys
import math

ct= 0
def insertion_sort(a, n, g):
    global ct
    for j in range(0,n-g):
        v = a[j+g]
        while j >= 0 and a[j] > v:
            a[j+g] = a[j]
            j = j-g
            ct += 1
        a[j+g] = v

n = int(input())
a = list(map(int, sys.stdin.readlines()))
b = 701
g = [x for x in [1,4,10,23,57,132,301,701] if x <= n]
while True:
    b = math.floor(2.25*b)
    if b > n:
        break
    g.append(b)
g = g[::-1]

for i in g:
    insertion_sort(a, n, i)

print(len(g))
print(*g, sep=" ")
print(ct)
print(*a, sep="\n")
# Original solution
def insertion_sort(num_list,gap):
    global cnt 

    for i in range(gap,len(num_list)):
        evc = num_list[i]
        j   = i-gap
        while j >= 0 and num_list[j] > evc:
            num_list[j+gap] = num_list[j]
            j = j-gap
            cnt += 1
            num_list[j+gap] = evc

def shell_sort(num_list):
    global cnt
    cnt = 0
    m_list = []
    h = 1
    while h <= len(num_list):
        m_list.append(h)
        h = 3*h+1
    m_list.reverse()
    m = len(m_list)

    print(m)
    print((' '.join([str(i) for i in m_list])))
    for i in range(m):
        insertion_sort(num_list,m_list[i])

n = int(eval(input()))
num_list = [int(eval(input())) for _ in range(n)]

shell_sort(num_list)
print(cnt)
for num in num_list:
    print(num)


# Optimized solution
def insertion_sort(num_list,gap):
    global cnt
    for i in range(gap,len(num_list)):
        swap_num = num_list[i]
        j = i-gap
        while j >= 0 and num_list[j] > swap_num:
            num_list[j+gap] = num_list[j]
            j -= gap
            cnt += 1
        num_list[j+gap] = swap_num

def shell_sort(num_list):
    global cnt 
    cnt = 0
    m = 1
    m_list = []
    while len(num_list) >= m:
        m_list.append(m)
        m = 3*m+1
    m_list.reverse()
    print((len(m_list)))
    print((' '.join([str(i) for i in m_list])))

    for i in range(len(m_list)):
        insertion_sort(num_list,m_list[i])
    print(cnt)

n = int(eval(input()))
num_list = [int(eval(input())) for _ in range(n)]
shell_sort(num_list)
for num in num_list:
    print(num)
# Original solution
n = int(eval(input()))
A = [int(eval(input())) for i in range(n)]

def insertionSort(A, n, g):
 global cnt
 for i in range(g, n):
  v = A[i]
  j = i - g
  while j >= 0 and A[j] > v:
   A[j + g] = A[j]
   j -= g
   cnt +=1
  A[j + g] = v

def shellSort(A, n, m, G):
 A1 = A[:]
 global cnt
 cnt = 0
 for i in range(m):
  insertionSort(A1, n, G[i])
 return A1

def Gmake(n):
 G=[]
 k = 1
 while k <= n and len(G) < 100:
  G.append(k)
  k += 3
 G.reverse()
 return G

G = Gmake(n)
A2 = shellSort(A, n, len(G), G)
print((str(len(G))))
print((" ".join(map(str, G))))
print((str(cnt)))
for i in A2:
 print((str(i)))

# Optimized solution
n = int(eval(input()))
A = [int(eval(input())) for i in range(n)]

def insertionSort(A, n, g):
 global cnt
 for i in range(g, n):
  v = A[i]
  j = i - g
  while j >= 0 and A[j] > v:
   A[j + g] = A[j]
   j -= g
   cnt +=1
  A[j + g] = v

def shellSort(A, n, m, G):
 A1 = A[:]
 global cnt
 cnt = 0
 for i in range(m):
  insertionSort(A1, n, G[i])
 return A1

def Gmake(n):
 G=[]
 k = 1
 while k <= n and len(G) < 100:
  G.append(k)
  k = 3 * k + 1
 G.reverse()
 return G

G = Gmake(n)
A2 = shellSort(A, n, len(G), G)
print((str(len(G))))
print((" ".join(map(str, G))))
print((str(cnt)))
for i in A2:
 print((str(i)))
# Original solution
def insertionSort( nums, n, g ):
	cnt = 0
	i = g
	while i < len( nums ):
		v = nums[i]
		j = i - g
		while 0 <= j and v < nums[j]:
			nums[ j+g ] = nums[j]
			j -= g
			cnt += 1
		nums[ j+g ] = v
		i += 1
	
	return cnt


def shellSort( nums, n ):
	cnt = 0
	g = []
	val =0
	for i in range( 0, n ):
		val = 3*val+1
		if  n < val:
			break
		g.append( val )

	g.reverse( )
	m = len( g )
	print( m )
	print(( " ".join( map( str, g ) ) ))
	for i in range( 0, m ):
		cnt += insertionSort( nums, n, g[i] )
	print( cnt )


n = int( input( ) )
nums = []
for i in range( 0, n ):
	nums.append( int( input( ) ) )

shellSort( nums, n )
for val in nums:
	print( val )

# Optimized solution
def insertionSort( nums, g ):
	cnt = 0
	i = g
	n = len( nums )
	while i < n:
		v = nums[i]
		j = i - g
		while 0 <= j and v < nums[j]:
			nums[ j+g ] = nums[j]
			j -= g
			cnt += 1
		nums[ j+g ] = v
		i += 1
	return cnt


def shellSort( nums, n ):
	cnt = 0
	g = []
	val =0
	for i in range( 0, n ):
		val = 3*val+1
		if  n < val:
			break
		g.append( val )

	g.reverse( )
	m = len( g )
	print( m )
	print(( " ".join( map( str, g ) ) ))
	for i in range( 0, m ):
		cnt += insertionSort( nums, g[i] )
	print( cnt )


n = int( input( ) )
nums = []
for i in range( 0, n ):
	nums.append( int( input( ) ) )

shellSort( nums, n )
for val in nums:
	print( val )
# Original solution
def insertionSort( nums, g ):
	cnt = 0
	i = g
	n = len( nums )
	while i < n:
		v = nums[i]
		j = i - g
		while 0 <= j and v < nums[j]:
			nums[ j+g ] = nums[j]
			j -= g
			cnt += 1
		nums[ j+g ] = v
		i += 1
	return cnt


def shellSort( nums, n ):
	cnt = 0
	g = []
	val =0
	for i in range( 0, n ):
		val = 3*val+1
		if  n < val:
			break
		g.append( val )

	g.reverse( )
	m = len( g )
	print( m )
	print(( " ".join( map( str, g ) ) ))
	for i in range( 0, m ):
		cnt += insertionSort( nums, g[i] )
	print( cnt )


n = int( input( ) )
nums = []
for i in range( 0, n ):
	nums.append( int( input( ) ) )

shellSort( nums, n )
for val in nums:
	print( val )

# Optimized solution
import sys

def insertionSort( nums, n, g ):
	cnt = 0
	i = g
	while i < n:
		v = nums[i]
		j = i - g
		while 0 <= j and v < nums[j]:
			nums[ j+g ] = nums[j]
			j -= g
			cnt += 1
		nums[ j+g ] = v
		i += 1
	return cnt
		
def shellSort( nums, n ):
	g = []
	val =0
	for i in range( n ):
		val = 3*val+1
		if  n < val:
			break
		g.append( val )
	
	g.reverse( )
	m = len( g )
	print( m )
	print(( " ".join( map( str, g ) ) ))
	cnt = 0
	for i in range( m ):
		cnt += insertionSort( nums, n, g[i] )
	print( cnt )
    
n = int( sys.stdin.readline( ) )
nums = []
for i in range( n ):
	nums.append( int( sys.stdin.readline( ) ) )

cnt = 0
shellSort( nums, n )
print(( "\n".join( map( str, nums ) ) ))
# Original solution
def insertionSort( nums, n, g ):
    cnt = 0
    for i in range( g, n ):
        v = nums[i]
        j = i - g
        while j >= 0 and nums[j] > v:
            nums[ j+g ] = nums[j]
            j -= g
            cnt += 1
         
        nums[ j+g ] = v
    return cnt
         
def shellSort( nums, n ):
    g = []
    v = 1
    while v <= n:
        g.append( v )
        v = 3*v+1     
    g.reverse( )
    m = len( g )
    print( m )
    print(( " ".join( map( str, g ) ) ))  
    cnt = 0
    for i in range( m ):
        cnt += insertionSort( nums, n, g[i] )
    print( cnt )
 
 
n = int( eval(input( )) )
nums = [ int( eval(input( )) ) for i in range( n ) ]
nums.append( 0 )
shellSort( nums, n )
nums.pop( )
for v in nums:
    print( v )

# Optimized solution
import sys
 
def insertionSort( nums, n, g ):
    cnt = 0
    i = g
    while i < n:
        v = nums[i]
        j = i - g
        while 0 <= j and v < nums[j]:
            nums[ j+g ] = nums[j]
            j -= g
            cnt += 1
        nums[ j+g ] = v
        i += 1
    return cnt
         
def shellSort( nums, n ):
    g = []
    val = 1
    while val <= n:
        g.append( val )
        val = 3*val+1
     
    g.reverse( )
    m = len( g )
    print( m )
    print(( " ".join( map( str, g ) ) ))
    cnt = 0
    for i in range( m ):
        cnt += insertionSort( nums, n, g[i] )
    print( cnt )
     
n = int( sys.stdin.readline( ) )
nums = []
for i in range( n ):
    nums.append( int( sys.stdin.readline( ) ) )
 
shellSort( nums, n )
print(( "\n".join( map( str, nums ) ) ))
# Original solution
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import math


def generate_intervals(n):
    if n > 1:
        l = math.ceil(math.log(2 * n + 1, 3))
        return [(3 ** i - 1) // 2 for i in range(l - 1, 0, -1)]
    else:
        return [1]


class Sort(object):

    def __init__(self, array, num):
        self.a = array
        self.n = num
        self.cnt = 0
        self.gs = generate_intervals(self.n)
        self.m = len(self.gs)

    def insertionsort(self, g):
        for i in range(g, self.n):
            v = self.a[i]
            j = i - g
            while j >= 0 and self.a[j] > v:
                self.a[j + g] = self.a[j]
                j = j - g
                self.cnt += 1
            self.a[j + g] = v

    def shellsort(self):
        for i in range(0, self.m):
            self.insertionsort(self.gs[i])


def main():
    n = int(eval(input()))
    a = [int(eval(input())) for i in range(n)]
    s = Sort(a, n)
    s.shellsort()
    print((s.m))
    print((" ".join(map(str, s.gs))))
    print((s.cnt))
    for i in range(n):
        print((s.a[i]))


if __name__ == "__main__":
    main()

# Optimized solution
#!/usr/bin/env python3

def construct_gs(n):
    gs = []
    g = 1
    while g <= n:
        gs.append(g)
        g = 3 * g + 1
    return gs[::-1]


def shell_sort(xs, n):
    global cnt
    cnt = 0
    gs = construct_gs(n)
    m = len(gs)
    def insertion_sort(g):
        global cnt
        for i in range(g, n):
            v = xs[i]
            j = i - g
            while j >= 0 and xs[j] > v:
                xs[j + g] = xs[j]
                j = j - g
                cnt += 1
            xs[j + g] = v
    for g in gs:
        insertion_sort(g)
    return m, gs, cnt, xs


def main():
    n = int(eval(input()))
    xs = [int(eval(input())) for _ in range(n)]
    m, gs, cnt, xs = shell_sort(xs, n)
    print(m)
    print((" ".join(map(str, gs))))
    print(cnt)
    for x in xs:
        print(x)


if __name__ == '__main__':
    main()
# Original solution
n = int(eval(input()))
A = list()
for i in range(0, n):
    A.append(int(eval(input())))

G = [1]
cnt = 0

for l in range(1, 100):
    h = 3 * G[l - 1] + 1
    if(h > n):
        m = l
        break
    else:
        G.append(3*G[l-1] + 1)
for g in reversed(G):
    for t in range(g, n):
        v = A[t]
        j = t - g
        while j >= 0 and A[j] > v:
            A[j+g] = A[j]
            j = j - g
            cnt += 1
        A[j+g] = v


print(m)
print((" ".join(str(x) for x in reversed(G))))
print(cnt)

print(("\n".join(str(x) for x in A)))

# Optimized solution
n = int(eval(input()))
A = list()
for i in range(0, n):
    A.append(int(eval(input())))

G = [1]  # G[0]????????¨?????????
cnt = 0

# G[i] = 3*Gi-1 + 1
for i in range(1, 100):
    h = 3 * G[i - 1] + 1
    if(h > n):  # h???n?¶???????????????§??????
        m = i
        break
    else:
        G.append(h)


for g in reversed(G):
    for i in range(g, n):
        v = A[i]
        j = i - g
        while j >= 0 and A[j] > v:
            A[j+g] = A[j]
            j = j - g
            cnt += 1
        A[j+g] = v


print(m)
print((" ".join(str(x) for x in reversed(G))))
print(cnt)

print(("\n".join(str(x) for x in A)))
# Original solution
cnt = 0
n = int(input())
A = [int(input()) for _ in range(n)]
G = [1] if n < 5 else [1, 4]
while G[-1]*3+1 <= n/5: G.append(G[-1]*3+1)
G.reverse()

def insersionSort(g):
	global cnt
	for i in range(g, n):
		v = A[i]
		j = i - g
		while j >= 0 and A[j] > v:
			A[j+g] = A[j]
			j = j - g
			cnt += 1
		A[j+g] = v

def shellSort():
	for g in G:
		insersionSort(g)

shellSort()
print(len(G))
print(" ".join(map(str, G)))
print(cnt)
print("\n".join(map(str, A)))

# Optimized solution
cnt = 0
n = int(input())
A = [int(input()) for _ in range(n)]
G = [1] if n < 5 else [1, 4]
while G[-1]*3+1 <= n/3: G.append(G[-1]*3+1)
G.reverse()

def insersionSort(g):
	global cnt
	for i in range(g, n):
		v = A[i]
		j = i - g
		while j >= 0 and A[j] > v:
			A[j+g] = A[j]
			j = j - g
			cnt += 1
		A[j+g] = v

def shellSort():
	for g in G:
		insersionSort(g)

shellSort()
print(len(G))
print(" ".join(map(str, G)))
print(cnt)
print("\n".join(map(str, A)))
# Original solution
N, M, K = list(map(int,input().split()))
MOD = 10**9+7

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

# 0~N*M-1までの1次元上と考える
ans = 0

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

c = cmb(N*M-2,K-2)

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)


    
    
    

# Optimized solution
N, M, K = list(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

# 0~N*M-1までの1次元上と考える
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)
# Original solution
from functools import reduce
from operator import mul

MOD = 1000000007


def combination(n, r):
  r = min(n - r, r)
  if r == 0:
    return 1
  over = reduce(mul, list(range(n, n - r, -1)))
  under = reduce(mul, list(range(1, r + 1)))

  return over // under


def main():
  N, M, K = [int(i) for i in input().strip().split(' ')]
  c = combination(N * M - 2, K - 2) % MOD

  X = sum([(M - i) * i for i in range(M)])
  Y = sum([(N - i) * i for i in range(N)])
  total = X * (N ** 2) + Y * (M ** 2) % MOD
  total = total * c % MOD
  print(total)


main()

# Optimized solution
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()
# Original solution
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()

# Optimized solution
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()
# Original solution
N, M, K = [int(item) for item in input().split()]
MOD = 10**9 + 7

n = 5 * 10**5 
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)

# Optimized solution
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)