In an NxM grid with each cell's dimension being 1x1, there will be (N+1) x (M+1) cross points. 
You task is to count the number of ways (S) of choosing K different points from these cross points such that all of them lie on a straight line and at least one of the cross points lies on the border.

Input Format  

A single line containing 3 integers N, M & K separated by a single space. 

Output Format  

A single integer denoting the number of ways (S) modulo 1000000007

Constraints

0 < N, M <= 3000 

2 <= K <= max(N, M) + 1

Sample Input  

2 2 3

Sample Output  

8

Explanation  

If you imagine a grid of the first quadrant of the co-ordinate system. Then, we have, 8 such 3 points of which at least 1 point on the borders. 

(0,0), (0,1), (0,2) 

(1,0), (1,1), (1,2) 

(2,0), (2,1), (2,2) 

(0,0), (1,0), (2,0) 

(0,1), (1,1), (2,1) 

(0,2), (1,2), (2,2) 

(0,0), (1,1), (2,2) and 

(0,2), (1,1), (2,0)

Difficulty: MEDIUM_HARD
Tags: Combinatorics

Starter code:

Solutions:
from time import time
from itertools import compress

def rwh_primes_iter(n):
	half = n >> 1
	sieve = bytearray([True]) * half
	for i in range(3, int(n ** 0.5) + 1, 2):
		if sieve[i >> 1]:
			sieve[i * i >> 1::i] = bytearray((n - i * i - 1) // (2 * i) + 1)
	return [2, *compress(range(3, n, 2), sieve[1:])]
mod = 7 + 10 ** 9
N_max = 10 ** 4
div = [[] for _ in range(N_max + 1)]
for d in range(1, int(N_max ** 0.5) + 1):
	div[d * d].append(d)
	for m in range(d * (d + 1), N_max + 1, d):
		div[m].extend([d, m // d])
for i in range(1, N_max + 1):
	div[i].sort()
primes = rwh_primes_iter(N_max)
mu = [1] * (N_max + 1)
for p in primes:
	for n in range(p, N_max + 1, p):
		mu[n] = -mu[n]
	step = p ** 2
	if step > N_max:
		continue
	for n in range(step, N_max + 1, step):
		mu[n] = 0
fac = [1] * (N_max + 1)
invfac = [1] * (N_max + 1)
for i in range(2, N_max + 1):
	fac[i] = fac[i - 1] * i % mod
invfac[N_max] = pow(fac[N_max], mod - 2, mod)
for i in reversed(range(2, N_max)):
	invfac[i] = invfac[i + 1] * (i + 1) % mod

def gcd(a, b):
	while b:
		(a, b) = (b, a % b)
	return a

def binomial(n, k):
	if n < k or k < 0:
		return 0
	return fac[n] * invfac[k] * invfac[n - k] % mod

def sum_floor(X, Y, d):
	L = X * d // Y
	U = [(x * Y - 1) // d for x in range(L + 2)]
	ans = 0
	for x in range(1, L + 1):
		if U[x + 1] >= X:
			ans += x * (X - U[x])
			break
		ans += x * (U[x + 1] - U[x])
	return ans

def sum_bin_floor(X, Y, d, p):
	U = [0] + [Y * d // y for y in range(1, Y * d + 2)]

def F(N, M, K):
	ans = 0
	for x in range(1, M):
		tmp = 0
		if (M - 1) // x < K:
			continue
		for d in div[x]:
			tmp += mu[x // d] * sum_floor(N - 1, M - 1, d)
		ans = (ans + binomial((M - 1) // x, K) * tmp) % mod
	for n in range(1, N):
		for x in range(1, n + 1):
			if n // x < K:
				continue
			tmp = 0
			for d in div[x]:
				tmp += mu[x // d] * ((M - 1) * d // n)
			ans = (ans + binomial(n // x, K) * tmp) % mod
	for n in range(1, N):
		ans = (ans - binomial(gcd(n, M - 1), K)) % mod
	return ans

def g(n, m, K):
	ans = 0
	for d in range(1, min(n, m) + 1):
		tmp = 0
		for b in range(1, m // d + 1):
			tmp = (tmp + binomial(m // (b * d), K) * min(n // d, n * b // m)) % mod
		ans = (ans + mu[d] * tmp) % mod
	return ans

def f(n, m, K):
	if n == 0 or m == 0:
		return 0
	d = gcd(n, m)
	return (g(n, m, K) + g(m, n, K) - binomial(d, K)) % mod

def f_1_2_bf(N, M, K):
	ans = binomial(N, K)
	for n in range(1, N + 1):
		for m in range(1, M):
			d = gcd(n, m)
			ans = (ans + binomial(d - 1, K - 1)) % mod
	return ans

def f_1_2(N, M, K):
	(X, Y) = (N, M - 1)
	L = min(X, Y)
	ans = binomial(N, K)
	for g in range(1, L + 1):
		tmp = 0
		for (d, g_d) in enumerate(range(g, L + 1, g)):
			if mu[d + 1] == 0:
				continue
			tmp += mu[d + 1] * (X // g_d) * (Y // g_d)
		ans = (ans + tmp * binomial(g - 1, K - 1)) % mod
	return ans

def f_1_3(N, M, K):
	ans = 0
	for n in range(N - 1, N + 1):
		ans = (ans + binomial(gcd(n, M) - 1, K - 1) * (N - n + 1)) % mod
	for n in range(1, N - 1):
		ans = (ans + binomial(gcd(n, M) - 1, K - 1) * 2 * (N - n)) % mod
	return ans

def solve(N, M, K):
	ans = 0
	ans = (F(N, M, K) + F(M, N, K)) % mod
	ans = (4 * ans - 2 * f(N - 1, M - 1, K) - 2 * f(M - 1, N - 1, K)) % mod
	ans = (ans + 2 * (N - 1) * binomial(M - 1, K) + 2 * (M - 1) * binomial(N - 1, K)) % mod
	ans = (ans + 2 * (binomial(N, K + 1) + binomial(M, K + 1))) % mod
	ans = (ans + 2 * f_1_2(N, M, K)) % mod
	ans = (ans + 2 * f_1_2(M, N, K)) % mod
	ans = (ans + f_1_3(N, M, K) + f_1_3(M, N, K)) % mod
	ans = (ans + (M - 1) * binomial(N - 1, K - 1) + (N - 1) * binomial(M - 1, K - 1)) % mod
	return ans
(N, M, K) = map(int, input().strip().split())
start = time()
print(solve(N, M, K - 1))


from time import time
from itertools import compress

def rwh_primes_iter(n):
	half = n >> 1
	sieve = bytearray([True]) * half
	for i in range(3, int(n ** 0.5) + 1, 2):
		if sieve[i >> 1]:
			sieve[i * i >> 1::i] = bytearray((n - i * i - 1) // (2 * i) + 1)
	return [2, *compress(range(3, n, 2), sieve[1:])]
mod = 7 + 10 ** 9
N_max = 10 ** 4
div = [[] for _ in range(N_max + 1)]
for d in range(1, int(N_max ** 0.5) + 1):
	div[d * d].append(d)
	for m in range(d * (d + 1), N_max + 1, d):
		div[m].extend([d, m // d])
for i in range(1, N_max + 1):
	div[i].sort()
primes = rwh_primes_iter(N_max)
mu = [1] * (N_max + 1)
for p in primes:
	for n in range(p, N_max + 1, p):
		mu[n] = -mu[n]
	step = p ** 2
	if step > N_max:
		continue
	for n in range(step, N_max + 1, step):
		mu[n] = 0
fac = [1] * (N_max + 1)
invfac = [1] * (N_max + 1)
for i in range(2, N_max + 1):
	fac[i] = fac[i - 1] * i % mod
invfac[N_max] = pow(fac[N_max], mod - 2, mod)
for i in reversed(range(2, N_max)):
	invfac[i] = invfac[i + 1] * (i + 1) % mod

def gcd(a, b):
	while b:
		(a, b) = (b, a % b)
	return a

def binomial(n, k):
	if n < k or k < 0:
		return 0
	return fac[n] * invfac[k] * invfac[n - k] % mod

def sum_floor(X, Y, d):
	L = X * d // Y
	U = [(x * Y - 1) // d for x in range(L + 2)]
	ans = 0
	for x in range(1, L + 1):
		if U[x + 1] >= X:
			ans += x * (X - U[x])
			break
		ans += x * (U[x + 1] - U[x])
	return ans

def sum_bin_floor(X, Y, d, p):
	U = [0] + [Y * d // y for y in range(1, Y * d + 2)]

def F(N, M, K):
	ans = 0
	for x in range(1, M):
		tmp = 0
		if (M - 1) // x < K:
			continue
		for d in div[x]:
			tmp += mu[x // d] * sum_floor(N - 1, M - 1, d)
		ans = (ans + binomial((M - 1) // x, K) * tmp) % mod
	for n in range(1, N):
		for x in range(1, n + 1):
			if n // x < K:
				continue
			tmp = 0
			for d in div[x]:
				tmp += mu[x // d] * ((M - 1) * d // n)
			ans = (ans + binomial(n // x, K) * tmp) % mod
	for n in range(1, N):
		ans = (ans - binomial(gcd(n, M - 1), K)) % mod
	return ans

def g(n, m, K):
	ans = 0
	for d in range(1, min(n, m) + 1):
		tmp = 0
		for b in range(1, m // d + 1):
			tmp = (tmp + binomial(m // (b * d), K) * min(n // d, n * b // m)) % mod
		ans = (ans + mu[d] * tmp) % mod
	return ans

def f(n, m, K):
	if n == 0 or m == 0:
		return 0
	d = gcd(n, m)
	return (g(n, m, K) + g(m, n, K) - binomial(d, K)) % mod

def f_1_2_bf(N, M, K):
	ans = binomial(N, K)
	for n in range(1, N + 1):
		for m in range(1, M):
			d = gcd(n, m)
			ans = (ans + binomial(d - 1, K - 1)) % mod
	return ans

def f_1_2(N, M, K):
	(X, Y) = (N, M - 1)
	L = min(X, Y)
	ans = binomial(N, K)
	for g in range(1, L + 1):
		tmp = 0
		for (d, g_d) in enumerate(range(g, L + 1, g)):
			if mu[d + 1] == 0:
				continue
			tmp += mu[d + 1] * (X // g_d) * (Y // g_d)
		ans = (ans + tmp * binomial(g - 1, K - 1)) % mod
	return ans

def f_1_3(N, M, K):
	ans = 0
	for n in range(N - 1, N + 1):
		ans = (ans + binomial(gcd(n, M) - 1, K - 1) * (N - n + 1)) % mod
	for n in range(1, N - 1):
		ans = (ans + binomial(gcd(n, M) - 1, K - 1) * 2 * (N - n)) % mod
	return ans

def solve(N, M, K):
	ans = 0
	ans = (F(N, M, K) + F(M, N, K)) % mod
	ans = (4 * ans - 2 * f(N - 1, M - 1, K) - 2 * f(M - 1, N - 1, K)) % mod
	ans = (ans + 2 * (N - 1) * binomial(M - 1, K) + 2 * (M - 1) * binomial(N - 1, K)) % mod
	ans = (ans + 2 * (binomial(N, K + 1) + binomial(M, K + 1))) % mod
	ans = (ans + 2 * f_1_2(N, M, K)) % mod
	ans = (ans + 2 * f_1_2(M, N, K)) % mod
	ans = (ans + f_1_3(N, M, K) + f_1_3(M, N, K)) % mod
	ans = (ans + (M - 1) * binomial(N - 1, K - 1) + (N - 1) * binomial(M - 1, K - 1)) % mod
	return ans
(N, M, K) = map(int, input().strip().split())
start = time()
print(solve(N, M, K - 1))


from time import time
from itertools import compress

def rwh_primes_iter(n):
	half = n >> 1
	sieve = bytearray([True]) * half
	for i in range(3, int(n ** 0.5) + 1, 2):
		if sieve[i >> 1]:
			sieve[i * i >> 1::i] = bytearray((n - i * i - 1) // (2 * i) + 1)
	return [2, *compress(range(3, n, 2), sieve[1:])]
mod = 7 + 10 ** 9
N_max = 10 ** 4
div = [[] for _ in range(N_max + 1)]
for d in range(1, int(N_max ** 0.5) + 1):
	div[d * d].append(d)
	for m in range(d * (d + 1), N_max + 1, d):
		div[m].extend([d, m // d])
for i in range(1, N_max + 1):
	div[i].sort()
primes = rwh_primes_iter(N_max)
mu = [1] * (N_max + 1)
for p in primes:
	for n in range(p, N_max + 1, p):
		mu[n] = -mu[n]
	step = p ** 2
	if step > N_max:
		continue
	for n in range(step, N_max + 1, step):
		mu[n] = 0
fac = [1] * (N_max + 1)
invfac = [1] * (N_max + 1)
for i in range(2, N_max + 1):
	fac[i] = fac[i - 1] * i % mod
invfac[N_max] = pow(fac[N_max], mod - 2, mod)
for i in reversed(range(2, N_max)):
	invfac[i] = invfac[i + 1] * (i + 1) % mod

def gcd(a, b):
	while b:
		(a, b) = (b, a % b)
	return a

def binomial(n, k):
	if n < k or k < 0:
		return 0
	return fac[n] * invfac[k] * invfac[n - k] % mod

def sum_floor(X, Y, d):
	L = X * d // Y
	U = [(x * Y - 1) // d for x in range(L + 2)]
	ans = 0
	for x in range(1, L + 1):
		if U[x + 1] >= X:
			ans += x * (X - U[x])
			break
		ans += x * (U[x + 1] - U[x])
	return ans

def sum_bin_floor(X, Y, d, p):
	U = [0] + [Y * d // y for y in range(1, Y * d + 2)]

def F(N, M, K):
	ans = 0
	for x in range(1, M):
		tmp = 0
		if (M - 1) // x < K:
			continue
		for d in div[x]:
			tmp += mu[x // d] * sum_floor(N - 1, M - 1, d)
		ans = (ans + binomial((M - 1) // x, K) * tmp) % mod
	for n in range(1, N):
		for x in range(1, n + 1):
			if n // x < K:
				continue
			tmp = 0
			for d in div[x]:
				tmp += mu[x // d] * ((M - 1) * d // n)
			ans = (ans + binomial(n // x, K) * tmp) % mod
	for n in range(1, N):
		ans = (ans - binomial(gcd(n, M - 1), K)) % mod
	return ans

def g(n, m, K):
	ans = 0
	for d in range(1, min(n, m) + 1):
		tmp = 0
		for b in range(1, m // d + 1):
			tmp = (tmp + binomial(m // (b * d), K) * min(n // d, n * b // m)) % mod
		ans = (ans + mu[d] * tmp) % mod
	return ans

def f(n, m, K):
	if n == 0 or m == 0:
		return 0
	d = gcd(n, m)
	return (g(n, m, K) + g(m, n, K) - binomial(d, K)) % mod

def f_1_2_bf(N, M, K):
	ans = binomial(N, K)
	for n in range(1, N + 1):
		for m in range(1, M):
			d = gcd(n, m)
			ans = (ans + binomial(d - 1, K - 1)) % mod
	return ans

def f_1_2(N, M, K):
	(X, Y) = (N, M - 1)
	L = min(X, Y)
	ans = binomial(N, K)
	for g in range(1, L + 1):
		tmp = 0
		for (d, g_d) in enumerate(range(g, L + 1, g)):
			if mu[d + 1] == 0:
				continue
			tmp += mu[d + 1] * (X // g_d) * (Y // g_d)
		ans = (ans + tmp * binomial(g - 1, K - 1)) % mod
	return ans

def f_1_3(N, M, K):
	ans = 0
	for n in range(N - 1, N + 1):
		ans = (ans + binomial(gcd(n, M) - 1, K - 1) * (N - n + 1)) % mod
	for n in range(1, N - 1):
		ans = (ans + binomial(gcd(n, M) - 1, K - 1) * 2 * (N - n)) % mod
	return ans

def solve(N, M, K):
	ans = 0
	ans = (F(N, M, K) + F(M, N, K)) % mod
	ans = (4 * ans - 2 * f(N - 1, M - 1, K) - 2 * f(M - 1, N - 1, K)) % mod
	ans = (ans + 2 * (N - 1) * binomial(M - 1, K) + 2 * (M - 1) * binomial(N - 1, K)) % mod
	ans = (ans + 2 * (binomial(N, K + 1) + binomial(M, K + 1))) % mod
	ans = (ans + 2 * f_1_2(N, M, K)) % mod
	ans = (ans + 2 * f_1_2(M, N, K)) % mod
	ans = (ans + f_1_3(N, M, K) + f_1_3(M, N, K)) % mod
	ans = (ans + (M - 1) * binomial(N - 1, K - 1) + (N - 1) * binomial(M - 1, K - 1)) % mod
	return ans
(N, M, K) = map(int, input().strip().split())
start = time()
print(solve(N, M, K - 1))


from time import time
from itertools import compress

def rwh_primes_iter(n):
	half = n >> 1
	sieve = bytearray([True]) * half
	for i in range(3, int(n ** 0.5) + 1, 2):
		if sieve[i >> 1]:
			sieve[i * i >> 1::i] = bytearray((n - i * i - 1) // (2 * i) + 1)
	return [2, *compress(range(3, n, 2), sieve[1:])]
mod = 7 + 10 ** 9
N_max = 10 ** 4
div = [[] for _ in range(N_max + 1)]
for d in range(1, int(N_max ** 0.5) + 1):
	div[d * d].append(d)
	for m in range(d * (d + 1), N_max + 1, d):
		div[m].extend([d, m // d])
for i in range(1, N_max + 1):
	div[i].sort()
primes = rwh_primes_iter(N_max)
mu = [1] * (N_max + 1)
for p in primes:
	for n in range(p, N_max + 1, p):
		mu[n] = -mu[n]
	step = p ** 2
	if step > N_max:
		continue
	for n in range(step, N_max + 1, step):
		mu[n] = 0
fac = [1] * (N_max + 1)
invfac = [1] * (N_max + 1)
for i in range(2, N_max + 1):
	fac[i] = fac[i - 1] * i % mod
invfac[N_max] = pow(fac[N_max], mod - 2, mod)
for i in reversed(range(2, N_max)):
	invfac[i] = invfac[i + 1] * (i + 1) % mod

def gcd(a, b):
	while b:
		(a, b) = (b, a % b)
	return a

def binomial(n, k):
	if n < k or k < 0:
		return 0
	return fac[n] * invfac[k] * invfac[n - k] % mod

def sum_floor(X, Y, d):
	L = X * d // Y
	U = [(x * Y - 1) // d for x in range(L + 2)]
	ans = 0
	for x in range(1, L + 1):
		if U[x + 1] >= X:
			ans += x * (X - U[x])
			break
		ans += x * (U[x + 1] - U[x])
	return ans

def sum_bin_floor(X, Y, d, p):
	U = [0] + [Y * d // y for y in range(1, Y * d + 2)]

def F(N, M, K):
	ans = 0
	for x in range(1, M):
		tmp = 0
		if (M - 1) // x < K:
			continue
		for d in div[x]:
			tmp += mu[x // d] * sum_floor(N - 1, M - 1, d)
		ans = (ans + binomial((M - 1) // x, K) * tmp) % mod
	for n in range(1, N):
		for x in range(1, n + 1):
			if n // x < K:
				continue
			tmp = 0
			for d in div[x]:
				tmp += mu[x // d] * ((M - 1) * d // n)
			ans = (ans + binomial(n // x, K) * tmp) % mod
	for n in range(1, N):
		ans = (ans - binomial(gcd(n, M - 1), K)) % mod
	return ans

def g(n, m, K):
	ans = 0
	for d in range(1, min(n, m) + 1):
		tmp = 0
		for b in range(1, m // d + 1):
			tmp = (tmp + binomial(m // (b * d), K) * min(n // d, n * b // m)) % mod
		ans = (ans + mu[d] * tmp) % mod
	return ans

def f(n, m, K):
	if n == 0 or m == 0:
		return 0
	d = gcd(n, m)
	return (g(n, m, K) + g(m, n, K) - binomial(d, K)) % mod

def f_1_2_bf(N, M, K):
	ans = binomial(N, K)
	for n in range(1, N + 1):
		for m in range(1, M):
			d = gcd(n, m)
			ans = (ans + binomial(d - 1, K - 1)) % mod
	return ans

def f_1_2(N, M, K):
	(X, Y) = (N, M - 1)
	L = min(X, Y)
	ans = binomial(N, K)
	for g in range(1, L + 1):
		tmp = 0
		for (d, g_d) in enumerate(range(g, L + 1, g)):
			if mu[d + 1] == 0:
				continue
			tmp += mu[d + 1] * (X // g_d) * (Y // g_d)
		ans = (ans + tmp * binomial(g - 1, K - 1)) % mod
	return ans

def f_1_3(N, M, K):
	ans = 0
	for n in range(N - 1, N + 1):
		ans = (ans + binomial(gcd(n, M) - 1, K - 1) * (N - n + 1)) % mod
	for n in range(1, N - 1):
		ans = (ans + binomial(gcd(n, M) - 1, K - 1) * 2 * (N - n)) % mod
	return ans

def solve(N, M, K):
	ans = 0
	ans = (F(N, M, K) + F(M, N, K)) % mod
	ans = (4 * ans - 2 * f(N - 1, M - 1, K) - 2 * f(M - 1, N - 1, K)) % mod
	ans = (ans + 2 * (N - 1) * binomial(M - 1, K) + 2 * (M - 1) * binomial(N - 1, K)) % mod
	ans = (ans + 2 * (binomial(N, K + 1) + binomial(M, K + 1))) % mod
	ans = (ans + 2 * f_1_2(N, M, K)) % mod
	ans = (ans + 2 * f_1_2(M, N, K)) % mod
	ans = (ans + f_1_3(N, M, K) + f_1_3(M, N, K)) % mod
	ans = (ans + (M - 1) * binomial(N - 1, K - 1) + (N - 1) * binomial(M - 1, K - 1)) % mod
	return ans
(N, M, K) = map(int, input().strip().split())
start = time()
print(solve(N, M, K - 1))


from time import time
from itertools import compress

def rwh_primes_iter(n):
	half = n >> 1
	sieve = bytearray([True]) * half
	for i in range(3, int(n ** 0.5) + 1, 2):
		if sieve[i >> 1]:
			sieve[i * i >> 1::i] = bytearray((n - i * i - 1) // (2 * i) + 1)
	return [2, *compress(range(3, n, 2), sieve[1:])]
mod = 7 + 10 ** 9
N_max = 10 ** 4
div = [[] for _ in range(N_max + 1)]
for d in range(1, int(N_max ** 0.5) + 1):
	div[d * d].append(d)
	for m in range(d * (d + 1), N_max + 1, d):
		div[m].extend([d, m // d])
for i in range(1, N_max + 1):
	div[i].sort()
primes = rwh_primes_iter(N_max)
mu = [1] * (N_max + 1)
for p in primes:
	for n in range(p, N_max + 1, p):
		mu[n] = -mu[n]
	step = p ** 2
	if step > N_max:
		continue
	for n in range(step, N_max + 1, step):
		mu[n] = 0
fac = [1] * (N_max + 1)
invfac = [1] * (N_max + 1)
for i in range(2, N_max + 1):
	fac[i] = fac[i - 1] * i % mod
invfac[N_max] = pow(fac[N_max], mod - 2, mod)
for i in reversed(range(2, N_max)):
	invfac[i] = invfac[i + 1] * (i + 1) % mod

def gcd(a, b):
	while b:
		(a, b) = (b, a % b)
	return a

def binomial(n, k):
	if n < k or k < 0:
		return 0
	return fac[n] * invfac[k] * invfac[n - k] % mod

def sum_floor(X, Y, d):
	L = X * d // Y
	U = [(x * Y - 1) // d for x in range(L + 2)]
	ans = 0
	for x in range(1, L + 1):
		if U[x + 1] >= X:
			ans += x * (X - U[x])
			break
		ans += x * (U[x + 1] - U[x])
	return ans

def sum_bin_floor(X, Y, d, p):
	U = [0] + [Y * d // y for y in range(1, Y * d + 2)]

def F(N, M, K):
	ans = 0
	for x in range(1, M):
		tmp = 0
		if (M - 1) // x < K:
			continue
		for d in div[x]:
			tmp += mu[x // d] * sum_floor(N - 1, M - 1, d)
		ans = (ans + binomial((M - 1) // x, K) * tmp) % mod
	for n in range(1, N):
		for x in range(1, n + 1):
			if n // x < K:
				continue
			tmp = 0
			for d in div[x]:
				tmp += mu[x // d] * ((M - 1) * d // n)
			ans = (ans + binomial(n // x, K) * tmp) % mod
	for n in range(1, N):
		ans = (ans - binomial(gcd(n, M - 1), K)) % mod
	return ans

def g(n, m, K):
	ans = 0
	for d in range(1, min(n, m) + 1):
		tmp = 0
		for b in range(1, m // d + 1):
			tmp = (tmp + binomial(m // (b * d), K) * min(n // d, n * b // m)) % mod
		ans = (ans + mu[d] * tmp) % mod
	return ans

def f(n, m, K):
	if n == 0 or m == 0:
		return 0
	d = gcd(n, m)
	return (g(n, m, K) + g(m, n, K) - binomial(d, K)) % mod

def f_1_2_bf(N, M, K):
	ans = binomial(N, K)
	for n in range(1, N + 1):
		for m in range(1, M):
			d = gcd(n, m)
			ans = (ans + binomial(d - 1, K - 1)) % mod
	return ans

def f_1_2(N, M, K):
	(X, Y) = (N, M - 1)
	L = min(X, Y)
	ans = binomial(N, K)
	for g in range(1, L + 1):
		tmp = 0
		for (d, g_d) in enumerate(range(g, L + 1, g)):
			if mu[d + 1] == 0:
				continue
			tmp += mu[d + 1] * (X // g_d) * (Y // g_d)
		ans = (ans + tmp * binomial(g - 1, K - 1)) % mod
	return ans

def f_1_3(N, M, K):
	ans = 0
	for n in range(N - 1, N + 1):
		ans = (ans + binomial(gcd(n, M) - 1, K - 1) * (N - n + 1)) % mod
	for n in range(1, N - 1):
		ans = (ans + binomial(gcd(n, M) - 1, K - 1) * 2 * (N - n)) % mod
	return ans

def solve(N, M, K):
	ans = 0
	ans = (F(N, M, K) + F(M, N, K)) % mod
	ans = (4 * ans - 2 * f(N - 1, M - 1, K) - 2 * f(M - 1, N - 1, K)) % mod
	ans = (ans + 2 * (N - 1) * binomial(M - 1, K) + 2 * (M - 1) * binomial(N - 1, K)) % mod
	ans = (ans + 2 * (binomial(N, K + 1) + binomial(M, K + 1))) % mod
	ans = (ans + 2 * f_1_2(N, M, K)) % mod
	ans = (ans + 2 * f_1_2(M, N, K)) % mod
	ans = (ans + f_1_3(N, M, K) + f_1_3(M, N, K)) % mod
	ans = (ans + (M - 1) * binomial(N - 1, K - 1) + (N - 1) * binomial(M - 1, K - 1)) % mod
	return ans
(N, M, K) = map(int, input().strip().split())
start = time()
print(solve(N, M, K - 1))