Problem:
   Problem Statement Snuke loves colorful balls. He has a total of N×K balls, K in each of his favorite N colors. The colors are numbered 1 through N .  He will arrange all of the balls in a row from left to right, in arbitrary order. Then, for each of the N colors, he will paint the leftmost ball of that color into color 0 , a color different from any of the N original colors.  After painting, how many sequences of the colors of the balls are possible? Find this number modulo 10^9+7 .      Constraints  1≤N,K≤2,000         Input The input is given from Standard Input in the following format:  N  K       Output Print the number of the possible sequences of the colors of the balls after painting, modulo 10^9+7 .        Sample Input 1 2 2      Sample Output 1 4  The following 4 sequences are possible:   (0,1,0,2)  (0,0,1,2)  (0,2,0,1)  (0,0,2,1)        Sample Input 2 3 1      Sample Output 2 1  The following 1 sequence is possible:   (0,0,0)        Sample Input 3 2 3      Sample Output 3 14       Sample Input 4 2000 2000      Sample Output 4 546381702

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])%mod * g2[n-r] % mod

mod = 10**9+7 
N = 4*10**6
g1 = [1]*(N+1) 
g2 = [1]*(N+1) 
inverse = [1]*(N+1) 

for i in range( 2, N + 1 ):
    g1[i]=( ( g1[i-1] * i ) % mod )
    inverse[i]=( ( -inverse[mod % i] * (mod//i) ) % mod )
    g2[i]=( (g2[i-1] * inverse[i]) % mod )

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

if K==1:
    print(1)
    exit()

dp=[0 for i in range(N+1)]
dp[N]=1

for i in range(N,-1,-1):
    for j in range(i-1,-1,-1):
        n=(N-(j+1))*(K-1)+N-i
        dp[j]+=(cmb(n+K-2,K-2,mod)*dp[j+1])%mod
        dp[j]%=mod

print((dp[0]*g1[N])%mod)

Solution (Python):
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

"""
・0を除く左端：N,N-1,...の順に並ぶとする。
・0は左から順に、z_N,z_{N-1},...と並ぶ（もとの正体）
・小さい数から順に処理する。n以下の数とz_n,...,z_1を並べたとき、先頭の0の個数ごとに値を持つ
"""

import numpy as np

N,K = map(int,read().split())

if K == 1:
    print(1)
    exit()

MOD = 10 ** 9 + 7

def cumprod(A, MOD = MOD):
    L = len(A); Lsq = int(L**.5+1)
    A = np.resize(A, Lsq**2).reshape(Lsq,Lsq)
    for n in range(1,Lsq):
        A[:,n] *= A[:,n-1]; A[:,n] %= MOD
    for n in range(1,Lsq):
        A[n] *= A[n-1,-1]; A[n] %= MOD
    return A.ravel()[:L]

def make_fact(U, MOD = MOD):
    x = np.arange(U, dtype = np.int64); x[0] = 1
    fact = cumprod(x, MOD)
    x = np.arange(U, 0, -1, dtype=np.int64); x[0] = pow(int(fact[-1]), MOD-2, MOD)
    fact_inv = cumprod(x, MOD)[::-1]
    return fact,fact_inv

U = (N + 10) * (K + 10)
fact, fact_inv = make_fact(U)


dp = np.array([0,1],np.int64)
for n in range(2,N+1):
    prev = dp
    dp = np.zeros(n+1,np.int64)
    """
    dp[i] には、prev[i-1,i,i+1,...]から遷移してくる。和をとって定数倍。
    累積和をいじったあと定数倍すればよい
    """
    S = prev.sum() % MOD
    np.cumsum(prev,out=prev); prev %= MOD
    dp[1] = S
    dp[2:] = S - prev[:-1]
    """
    n の挿入方法を乗ずる。
    [z^i][n][右側：小さい側が(n-1)K-i+1個、nがK-2個]
    """
    coef = fact[n*(K-1)-1:n*K][::-1].copy()
    coef *= fact_inv[K-2]; coef %= MOD
    coef *= fact_inv[(n-1)*(K-1):(n-1)*K+2][::-1]; coef %= MOD
    dp *= coef; dp %= MOD

answer = dp.sum() % MOD
answer *= fact[N]; answer %= MOD
print(answer)

Solution (Python):
def pow_mod(MOD, a, n):
    ret = 1
    while n != 0:
        if n % 2 == 1: ret = ret * a % MOD
        a = a * a % MOD
        n //= 2
    return ret

def fac(n):
    MOD = 1000000007
    ret = 1
    for i in range(1, n + 1):
        ret = ret * i % MOD
    return ret

def fac_gen(first, last):
    MOD = 1000000007
    ret = fac(first)
    yield ret
    for i in range(first + 1, last):
        ret = ret * i % MOD
        yield ret

def build_comb_list(f, m, K):
    MOD = 1000000007
    v = pow_mod(MOD, f[m - K * 2] * f[0], MOD - 2)
    for i in range(m - K, -1, -1):
        f[i] *= v
        v = v * i % MOD

def solve_dp(n, k, comb):
    MOD = 1000000007
    dp = [1] * (n + 1)
    for i in range(2, n + 1):
        ci = (i - 1) * (k - 1)
        v = 0
        for j, c in zip(range(i), comb[ci:ci + i]):
            v = (v + dp[j] * c) % MOD
            dp[j] = v
        dp[i] = v
    return dp[n] * fac(n) % MOD

def solve(n, k):
    if n == 1 or k == 1: return 1
    m = n * k
    comb = list(fac_gen(k - 2, m - 1))
    build_comb_list(comb, m - 2, k - 2)
    return solve_dp(n, k, comb)


n, k = (int(s) for s in input().split(' '))
print(solve(n, k))