stablegradients/pie-gem5-bysrc huggingface dataset:
// Original C++ solution (slower):
#include <iostream>

#include <cstdio>

using namespace std;



const int MOD(1000000007);

int n, k;

int fac[4000000], inv_fac[4000000], dp[2001][2001];





int Power(long long a, int b)

{

	long long ret(1ll);

	

	while (b) {

		b & 1 && (ret *= a);

		ret >= MOD && (ret %= MOD);

		a *= a;

		a >= MOD && (a %= MOD);

		b >>= 1;

	}

	return ret;

}





int main()

{

	long long lt;

	int t;

	

	scanf("%d%d", &n, &k);

	if (k == 1) {

		printf("1\n");

		return 0;

	}

	fac[0] = 1;

	t = n * k;

	for (int i = 1; i < t; ++i) {

		lt = (long long)i * fac[i - 1];

		fac[i] = lt < MOD ? lt : lt % MOD; 

	}

	inv_fac[t - 1] = Power((long long)fac[t - 1], MOD - 2);

	for (int i = t - 2; i; --i) {

		lt = (long long)(i + 1) * inv_fac[i + 1];

		inv_fac[i] = lt < MOD ? lt : lt % MOD;

	}

	inv_fac[0] = 1;

	dp[0][0] = 1;

	t = -2;

	for (int i = 1; i <= n; ++i) {

		t += k - 1;

		for (int j = 0; j <= i; ++j) {

			++t;

			lt = (long long)fac[t] * inv_fac[t - k + 2];

			lt >= MOD && (lt %= MOD);

			lt *= (long long)inv_fac[k - 2];

			lt >= MOD && (lt %= MOD);

			lt *= (long long)dp[i - 1][j];

			j && (lt += dp[i][j - 1]);

			dp[i][j] = lt < MOD ? lt : lt % MOD;

		}

		t -= i + 1;

	}

	lt = (long long)dp[n][n] * fac[n] % MOD;

	printf("%lld\n", lt);

	return 0;

}




// Optimized solution to the same problem (faster):
#include <cstdio>

#include <algorithm>

using namespace std;

const int MAXN=2010;

const int MOD=1e9+7;

int n,k;

long long fac[MAXN*MAXN],inv[MAXN*MAXN],f[MAXN][MAXN];

long long C(long long x,long long y)

{

	return fac[x]*inv[y]%MOD*inv[x-y]%MOD;

}

long long qpow(long long x,long long y)

{

	long long tmp=1;

	for(;y;y>>=1,x=x*x%MOD)if(y&1)tmp=tmp*x%MOD;

	return tmp;

}

int main()

{

//	freopen("test.in","r",stdin);

//	freopen("test.out","w",stdout);

	scanf("%d%d",&n,&k);

	if(k==1)

	{

		printf("1");

		return 0;

	}

	fac[0]=1;

	for(int i=1;i<=n*k;i++)fac[i]=fac[i-1]*i%MOD;

	inv[n*k]=qpow(fac[n*k],MOD-2);

	for(int i=n*k;i>=1;i--)inv[i-1]=inv[i]*i%MOD;

	f[0][0]=1;

	for(int i=1;i<=n;i++)

		for(int j=0;j<=i;j++)

		{

			f[i][j]=f[i-1][j];

			if(!j)continue;

			f[i][j]=(f[i][j]+f[i][j-1]*(n-j+1)%MOD*C(n-i+(n-j+1)*(k-1)-1,k-2))%MOD;

		}

	printf("%lld\n",f[n][n]);

	return 0;

}

// Fastest known solution for this problem:
// Original C++ solution (slower):
#include<bits/stdc++.h>

using namespace std;

#define LL long long

#define int long long

#define Mod 1000000007ll

#define pow ksm

inline int ksm(int a,int x){

	LL ans=1,w=a;

	while(x){

		if(x&1)ans=ans*w%Mod;

		x>>=1;

		w=w*w%Mod;

	}

	return ans;

}

#define MN 4000005

#define mn 2005

int fac[MN],dfac[MN],n,k,f[mn][mn];

inline int C(int n,int m){

	return 1ll*fac[m]*dfac[n]%Mod*dfac[m-n]%Mod;

}

void init(int N){

	fac[0]=1;

	for(int i=1;i<=N;++i)

		fac[i]=fac[i-1]*i%Mod;

	dfac[N]=ksm(fac[N],Mod-2);

	for(int i=N-1;i>=0;--i){

		dfac[i]=(LL)dfac[i+1]*(LL)(i+1)%Mod;

	}

}

signed main(){

	scanf("%lld%lld",&n,&k);

	if(k==1){

		puts("1");

		return 0;

	}

	init(n*k);

	f[0][0]=1;

	for(int i=1;i<=n;++i)

		for(int j=0;j<=i;++j){

			f[i][j]=(f[i][j]+f[i-1][j])%Mod;

			if(j)f[i][j]=(f[i][j]+(f[i][j-1]*(n-j+1)%Mod*C(k-2,n-i+(k-1ll)*(n-j+1ll)-1ll)+Mod))%Mod;	

		}

	printf("%lld",f[n][n]);

	return 0;

}


// Optimized solution to the same problem (faster):
#include <cstdio>



typedef long long ll;

const int maxn=2000+10,mod=1000000007;

ll fac[maxn*maxn],ifac[maxn*maxn],f[maxn][maxn];



ll power(ll a,ll b){

	ll ret=1;

	while(b){

		if(b&1)ret=ret*a%mod;

		a=a*a%mod;

		b>>=1;

	}

	return ret;

}



ll C(ll n,ll m){

	if(n<m)return 0;

	return fac[n]*ifac[m]%mod*ifac[n-m]%mod;

}



int main(){

	// freopen("in.txt","r",stdin);

	int n,m;scanf("%d%d",&n,&m);

	if(m==1){printf("%d\n",1);return 0;}

	fac[0]=ifac[0]=1;

	for(int i=1;i<=n*m;++i)fac[i]=fac[i-1]*i%mod;

	ifac[n*m]=power(fac[n*m],mod-2);

	for(int i=n*m-1;i>=1;--i)ifac[i]=ifac[i+1]*(i+1)%mod;

	f[0][0]=1;

	for(int i=0;i<=n;++i)

		for(int j=0;j<=i;++j){

			if(i>j)f[i][j]+=f[i-1][j];

			if(j)f[i][j]+=f[i][j-1]*C(n*m-i-(j-1)*(m-1)-1,m-2)%mod;

			f[i][j]%=mod;

			// printf("%d %d %lld\n",i,j,f[i][j]);

		}

	printf("%lld\n",f[n][n]*fac[n]%mod);

	return 0;

}

// Fastest known solution for this problem:
// Original C++ solution (slower):
#include <bits/stdc++.h>

#define rep(i, x, y) for (int i = x; i <= y; i++)

using namespace std;



const int N = 2e3 + 10, mod = 1e9 + 7;

typedef long long ll;

ll n, K;

ll fac[N * N], f[N][N], inv[N * N];



ll quick_pow(ll a, ll b) {

    ll ret = 1;

    for (; b; b >>= 1) {

        if (b & 1) ret = ret * a % mod;

        a = a * a % mod;

    } return ret;

}



void pre() {

    fac[0] = inv[0] = 1;

    rep(i, 1, n * K) fac[i] = fac[i - 1] * i % mod;

    inv[n * K] = quick_pow(fac[n * K], mod - 2);

    for (int i = n * K - 1; i >= 1; i--) inv[i] = inv[i + 1] * (i + 1) % mod;

}



ll C(ll n, ll m) { return fac[n] * inv[m] % mod * inv[n - m] % mod; }



ll upd(ll x, ll y) { return (x + y) % mod; }



int main() {

    cin >> n >> K;

    if (K == 1) { puts("1"); return 0; }

    pre();

    rep(i, 1, n) f[i][0] = 1;

    rep(i, 1, n)

        rep(j, 1, i) {

            f[i][j] = upd(f[i - 1][j], (n - j + 1) * f[i][j - 1] % mod * C((n - j + 1) * (K - 1) + n - i - 1, K - 2) % mod);

        }

    printf("%lld\n", f[n][n]);

    return 0;

}


// Optimized solution to the same problem (faster):
#include <cstdio>



#define N 2001

#define P 4000001

#define M 1000000007LL



long long p[2][P], f[N][N];

int n, k, t;



long long sel(int n, int m)

{

    return n >= m ? p[0][n] * p[1][m] % M * p[1][n - m] % M : 0;

}



int main()

{

    scanf("%d%d", &n, &k);

    t = n * k;



    if (k == 1)

    {

        puts("1");

        return 0;

    }



    p[0][0] = 1;

    for (int i = 1; i <= t; i++)

        p[0][i] = p[0][i - 1] * i % M;



    p[1][t] = 1;

    for (long long x = p[0][t], e = M - 2; e; e >>= 1, x = x * x % M)

        if (e & 1)

            p[1][t] = p[1][t] * x % M;



    for (int i = t; i; i--)

        p[1][i - 1] = p[1][i] * i % M;



    f[0][0] = 1;

    for (int i = 0; i <= n; i++)

    {

        for (int j = i; j <= n; j++)

        {

            f[i][j + 1] = (f[i][j + 1] + f[i][j]) % M;

            f[i + 1][j] = (f[i + 1][j] + f[i][j] * sel(t - i * (k - 1) - j - 1, k - 2) % M) % M;

        }

    }



    printf("%lld\n", f[n][n] * p[0][n] % M);

    return 0;

}


// Fastest known solution for this problem:
// Original C++ solution (slower):
#include <cstdio>

#include <cstring>

#include <iostream>

#include <algorithm>

using namespace std;

#define ll long long

#define mod 1000000007

#define N 2005

#define M 4000005

void Add(int &x,int y)

{

	x += y;

	if (x >= mod)

		x -= mod;

}

int n, k, m, f[N][N];

int fac[M], inv[M];

int C(int n, int m){return (ll)fac[n] * inv[m] % mod * inv[n - m] % mod;}

int Pow(int a, int k)

{

	int ret = 1;

	while (k)

	{

		if (k & 1)

			ret = (ll)ret * a % mod;

		a = (ll)a * a % mod;

		k >>= 1;

	}

	return ret;

}

int main()

{

	scanf("%d%d", &n, &k);

	m = n * k;

	if (k == 1)

	{

		puts("1");

		return 0;

	}

	fac[0] = 1;

	for (int i = 1; i <= m; i++)

		fac[i] = (ll)fac[i - 1] * i % mod;

	inv[m] = Pow(fac[m], mod - 2);

	for (int i = m; i >= 1; i--)

		inv[i - 1] = (ll)inv[i] * i % mod;

	f[0][0] = 1;

	for (int i = 1; i <= n; i++)

		for (int j = 0; j <= i; j++)

		{

			Add(f[i][j], f[i - 1][j]);

			if(j)

				Add(f[i][j], (ll)f[i][j - 1] * (n - j + 1) % mod * C((n - i) + (n - j + 1) * (k - 1) - 1, k - 2) % mod);

		}

	printf("%d\n", f[n][n]);

}


// Optimized solution to the same problem (faster):
#include <cstdio>



#define N 2001

#define P 4000001

#define M 1000000007LL



long long p[2][P], f[N][N];

int n, k, t;



long long sel(int n, int m)

{

    return n >= m ? p[0][n] * p[1][m] % M * p[1][n - m] % M : 0;

}



int main()

{

    scanf("%d%d", &n, &k);

    t = n * k;



    if (k == 1)

    {

        puts("1");

        return 0;

    }



    p[0][0] = 1;

    for (int i = 1; i <= t; i++)

        p[0][i] = p[0][i - 1] * i % M;



    p[1][t] = 1;

    for (long long x = p[0][t], e = M - 2; e; e >>= 1, x = x * x % M)

        if (e & 1)

            p[1][t] = p[1][t] * x % M;



    for (int i = t; i; i--)

        p[1][i - 1] = p[1][i] * i % M;



    f[0][0] = 1;

    for (int i = 0; i <= n; i++)

    {

        for (int j = i; j <= n; j++)

        {

            f[i][j + 1] = (f[i][j + 1] + f[i][j]) % M;

            f[i + 1][j] = (f[i + 1][j] + f[i][j] * sel(t - i * (k - 1) - j - 1, k - 2) % M) % M;

        }

    }



    printf("%lld\n", f[n][n] * p[0][n] % M);

    return 0;

}


// Fastest known solution for this problem:
// Original C++ solution (slower):
#include <cstdio>

#include <iostream>



#define MAXN 2050

#define ll long long

#define Mod 1000000007



using namespace std ; int i, j ; int N, M ; 

ll Fac[MAXN * MAXN], Inv[MAXN * MAXN]/*Co[MAXN][MAXN], */, dp[MAXN][MAXN] ;



ll expow(ll a, int b){

	ll res = 1 ; 

	while (b){

		if (b & 1) res = res * a % Mod ;

		a = a * a % Mod, b = b >> 1 ;

	}

	return res ;

} 

ll Co(int n, int m){

	if (n < m) return 0 ; 

	return Fac[n] * Inv[m] % Mod * Inv[n - m] % Mod ; 

}

int main(){

	cin >> N >> M ; /*Co[0][0] = 1 ; */	 Inv[0] = 1 ;

	for (i = 0 ; i <= N ; ++ i) dp[i][0] = 1 ; Fac[0] = 1 ;

	for (i = 1 ; i <= N * M ; ++ i) Fac[i] = Fac[i - 1] * i % Mod ;

	for (i = 1 ; i <= N * M ; ++ i) Inv[i] = expow(Fac[i], Mod - 2) ;  

	/*for (i = 1 ; i <= N ; ++ i){

        Co[i][0] = 1 ;

        for (j = 1 ; j <= N ; ++ j) Co[i][j] = (Co[i - 1][j - 1] + Co[i - 1][j]) % Mod ;

    }*/

	for (i = 1 ; i <= N ; ++ i)

		for (j = 1 ; j <= i ; ++ j)

			dp[i][j] = (dp[i - 1][j]  + dp[i][j - 1] * (N - j + 1) % Mod * Co(N * M - i - (j - 1) * (M - 1) - 1, M - 2) % Mod) % Mod ;

	cout << (M <= 1 ? 1 : dp[N][N]) << endl ; return 0 ; 

}

// Optimized solution to the same problem (faster):
#include <cstdio>

#include <algorithm>

#include <cctype>

using namespace std;

inline int read(int f = 1, int x = 0, char ch = ' ')

{

    while(!isdigit(ch = getchar())) if(ch == '-') f = -1;

    while(isdigit(ch)) x = x*10+ch-'0', ch = getchar();

    return f*x;

}

const int N = 2e3+5, M = 4e6+5, P = 1e9+7;

int n, m, f[N][N], fac[M], ifac[M], inv[M], ans;

int C(int n, int m) { return 1ll*fac[n]*ifac[n-m]%P*ifac[m]%P; }

int main()

{

    n = read(), m = read(); 

    inv[1] = ifac[0] = fac[0] = 1; 

    for(int i = 2; i <= n*m; ++i) inv[i] = 1ll*(P-P/i)*inv[P%i]%P;

    for(int i = 1; i <= n*m; ++i) fac[i] = 1ll*i*fac[i-1]%P, ifac[i] = 1ll*inv[i]*ifac[i-1]%P;

    f[0][0] = 1;

    for(int i = 1; i <= n; ++i)

        for(int j = 0; j <= i; ++j)

        {

            f[i][j] = f[i-1][j];

            if(j) f[i][j] = (f[i][j]+1ll*f[i][j-1]*(n-j+1)%P*C(n*m-i-(j-1)*(m-1)-1, m-2)%P)%P;

        }

    printf("%d\n", m == 1?1:f[n][n]);

    return 0;

}

// Fastest known solution for this problem:
// Original C++ solution (slower):
#include<iostream>

#include<cstring>

#include<cstdio>

#include<algorithm>

#define rep(i,a,b) for(int i=a;i<=b;i++)

#define mod 1000000007

#define N 2010

#define lint long long

#define debug(x) cerr<<#x<<"="<<x

#define sp <<" "

#define ln <<endl

using namespace std;

int dp[N][N],f[N*N],fi[N*N];

inline int fast_pow(int x,int k,int ans=1)

{	for(;k;k>>=1,x=(lint)x*x%mod) (k&1)?ans=(lint)ans*x%mod:0;return ans; }

inline int prelude(int n)

{

	for(int i=f[0]=1;i<=n;i++) f[i]=(lint)f[i-1]*i%mod;

	fi[n]=fast_pow(f[n],mod-2);

	for(int i=n-1;i>=0;i--) fi[i]=fi[i+1]*(i+1ll)%mod;

	return 0;

}

inline int C(int n,int m)

{

	if(m>n||m<0||n<0) return 0;

	return (lint)f[n]*fi[m]%mod*fi[n-m]%mod;

}

int main()

{

	int n,k;scanf("%d%d",&n,&k);

	dp[0][0]=1,prelude(n*k);

	if(k==1) return !printf("1\n");

	rep(i,0,n-1) rep(j,i,n) if(dp[i][j])

		(i<j?(dp[i+1][j]+=dp[i][j])%=mod:0),

		(dp[i][j+1]+=(lint)dp[i][j]*C((j+1)*(k-1)+i-1,k-2)%mod)%=mod;

//	rep(i,0,n) rep(j,i,n) debug(i)sp,debug(j)sp,debug(dp[i][j])ln;

	return !printf("%lld\n",(lint)dp[n][n]*f[n]%mod);

}

// Optimized solution to the same problem (faster):
#include <cstdio>



const int mod = 1e9 + 7;

const int MAXN = 2005, MAXS = MAXN * MAXN;



template<typename _T>

void read( _T &x )

{

	x = 0;char s = getchar();int f = 1;

	while( s > '9' || s < '0' ){if( s == '-' ) f = -1; s = getchar();}

	while( s >= '0' && s <= '9' ){x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar();}

	x *= f;

}



template<typename _T>

void write( _T x )

{

	if( x < 0 ){ putchar( '-' ); x = ( ~ x ) + 1; }

	if( 9 < x ){ write( x / 10 ); }

	putchar( x % 10 + '0' );

}



int f[MAXN][MAXN];

int fac[MAXS], ifac[MAXS];

int N, K;



int qkpow( int, int );

int inv( const int a ) { return qkpow( a, mod - 2 ); }

void sub( int &x, const int v ) { x = ( x < v ? x - v + mod : x - v ); }

void add( int &x, const int v ) { x = ( x + v >= mod ? x + v - mod : x + v ); }



int qkpow( int base, int indx )

{

	int ret = 1;

	while( indx )

	{

		if( indx & 1 ) ret = 1ll * ret * base % mod;

		base = 1ll * base * base % mod, indx >>= 1;

	}

	return ret;

}



void init( const int siz )

{

	fac[0] = ifac[0] = 1;

	for( int i = 1 ; i <= siz ; i ++ ) fac[i] = 1ll * fac[i - 1] * i % mod;

	ifac[siz] = inv( fac[siz] );

	for( int i = siz - 1 ; ~ i ; i -- ) ifac[i] = 1ll * ifac[i + 1] * ( i + 1 ) % mod;

}



int C( const int n, const int m ) { return 1ll * fac[n] * ifac[m] % mod * ifac[n - m] % mod; }



int main()

{

	read( N ), read( K );

	if( K == 1 ) { puts( "1" ); return 0; }

	init( N * K );

	for( int i = 0 ; i <= N ; i ++ ) f[i][0] = 1;

	for( int i = 1 ; i <= N ; i ++ )

		for( int j = 1 ; j <= i ; j ++ )

		{

			f[i][j] = f[i - 1][j];

			int t = N * K - i - ( K - 1 ) * ( j - 1 ) - 1;

			add( f[i][j], 1ll * f[i][j - 1] * ( N - j + 1 ) % mod * C( t, K - 2 ) % mod );

		}

	write( f[N][N] ), putchar( '\n' );

	return 0;

}

// Fastest known solution for this problem:
// Original C++ solution (slower):
#include<bits/stdc++.h>

#define For(i,j,k) for(int i=j;i<=k;++i)

#define Forr(i,j,k) for(int i=j;i>=k;--i)

#define ll long long 

using namespace std;

const int N = 2010, INF = 0x3f3f3f3f, Mod = 1e9 + 7;

template<class T>void read(T &x){

	x=0;char c=getchar();

	while(!isdigit(c))c=getchar();

	while( isdigit(c))x=x*10+c-48,c=getchar();

}

inline void file(){

	freopen("test.in","r",stdin);

	freopen("test.out","w",stdout);

}

int n, k, fac[N * N], ifac[N * N];

ll qpow(ll a, ll b){

	ll ret=1;

	for(;b;b>>=1,a=a*a%Mod)if(b&1)ret=ret*a%Mod;

	return ret;

}

void init(){

	read(n), read(k);

	fac[0] = 1;

	For(i, 1, n * k)fac[i] = 1ll * fac[i - 1] * i % Mod;

	ifac[n * k] = qpow(fac[n * k], Mod - 2);

	Forr(i, n * k, 1)ifac[i - 1] = 1ll * ifac[i] * i % Mod;

}

inline ll C(int a, int b){

	if(a < 0 || b < 0 || a < b)return 0;

	return fac[a] * 1ll * ifac[b] % Mod * ifac[a - b] % Mod;

}

int dp[N][N];

void solve(){

	dp[0][0] = 1;

	For(i, 0, n)

		For(j, i, n){

			if(i) dp[i][j] = dp[i - 1][j];

			if(j) dp[i][j] = (dp[i][j] + dp[i][j - 1] * C(i + j * (k - 1) - 1, k - 2)) % Mod;

		}

	if(k != 1)

		printf("%lld\n", 1ll * dp[n][n] * fac[n] % Mod);

	else puts("1");

}

int main(){

	//file();

	init();

	solve();

	return 0;

}


// Optimized solution to the same problem (faster):
#include <cstdio>



const int mod = 1e9 + 7;

const int MAXN = 2005, MAXS = MAXN * MAXN;



template<typename _T>

void read( _T &x )

{

	x = 0;char s = getchar();int f = 1;

	while( s > '9' || s < '0' ){if( s == '-' ) f = -1; s = getchar();}

	while( s >= '0' && s <= '9' ){x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar();}

	x *= f;

}



template<typename _T>

void write( _T x )

{

	if( x < 0 ){ putchar( '-' ); x = ( ~ x ) + 1; }

	if( 9 < x ){ write( x / 10 ); }

	putchar( x % 10 + '0' );

}



int f[MAXN][MAXN];

int fac[MAXS], ifac[MAXS];

int N, K;



int qkpow( int, int );

int inv( const int a ) { return qkpow( a, mod - 2 ); }

void sub( int &x, const int v ) { x = ( x < v ? x - v + mod : x - v ); }

void add( int &x, const int v ) { x = ( x + v >= mod ? x + v - mod : x + v ); }



int qkpow( int base, int indx )

{

	int ret = 1;

	while( indx )

	{

		if( indx & 1 ) ret = 1ll * ret * base % mod;

		base = 1ll * base * base % mod, indx >>= 1;

	}

	return ret;

}



void init( const int siz )

{

	fac[0] = ifac[0] = 1;

	for( int i = 1 ; i <= siz ; i ++ ) fac[i] = 1ll * fac[i - 1] * i % mod;

	ifac[siz] = inv( fac[siz] );

	for( int i = siz - 1 ; ~ i ; i -- ) ifac[i] = 1ll * ifac[i + 1] * ( i + 1 ) % mod;

}



int C( const int n, const int m ) { return 1ll * fac[n] * ifac[m] % mod * ifac[n - m] % mod; }



int main()

{

	read( N ), read( K );

	if( K == 1 ) { puts( "1" ); return 0; }

	init( N * K );

	for( int i = 0 ; i <= N ; i ++ ) f[i][0] = 1;

	for( int i = 1 ; i <= N ; i ++ )

		for( int j = 1 ; j <= i ; j ++ )

		{

			f[i][j] = f[i - 1][j];

			int t = N * K - i - ( K - 1 ) * ( j - 1 ) - 1;

			add( f[i][j], 1ll * f[i][j - 1] * ( N - j + 1 ) % mod * C( t, K - 2 ) % mod );

		}

	write( f[N][N] ), putchar( '\n' );

	return 0;

}

// Fastest known solution for this problem:
// Original C++ solution (slower):
#include <bits/stdc++.h>

#define ll long long

using namespace std;

const int maxn=2000+10;

const int p=1e9+7;

int n,k,f[maxn][maxn],farc[maxn*maxn],inv[maxn*maxn];



int C(int n,int m){

	return (ll)farc[n]*inv[m]%p*(ll)inv[n-m]%p;

}



int main()

{

	scanf("%d%d",&n,&k);k--;

	if(!k){printf("1\n");return 0;}

	farc[0]=farc[1]=inv[0]=inv[1]=1;

	for(int i=2;i<=n*k+n+k;i++){

		farc[i]=(ll)farc[i-1]*i%p;

		inv[i]=(ll)(p-p/i)*inv[p%i]%p;

	}

	for(int i=2;i<=n*k+n+k;i++) inv[i]=(ll)inv[i]*inv[i-1]%p;

	f[0][0]=1;

	for(int i=0;i<=n;i++)

		for(int j=i;j<=n;j++){

			if(i>0) f[i][j]=(f[i][j]+f[i-1][j])%p;

			if(j>0) f[i][j]=(f[i][j]+(ll)f[i][j-1]*C(i+(j-1)*k+k-1,k-1)%p)%p;

		}

	printf("%d\n",(ll)f[n][n]*farc[n]%p);

	return 0;

}

// Optimized solution to the same problem (faster):
#include <cstdio>



const int mod = 1e9 + 7;

const int MAXN = 2005, MAXS = MAXN * MAXN;



template<typename _T>

void read( _T &x )

{

	x = 0;char s = getchar();int f = 1;

	while( s > '9' || s < '0' ){if( s == '-' ) f = -1; s = getchar();}

	while( s >= '0' && s <= '9' ){x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar();}

	x *= f;

}



template<typename _T>

void write( _T x )

{

	if( x < 0 ){ putchar( '-' ); x = ( ~ x ) + 1; }

	if( 9 < x ){ write( x / 10 ); }

	putchar( x % 10 + '0' );

}



int f[MAXN][MAXN];

int fac[MAXS], ifac[MAXS];

int N, K;



int qkpow( int, int );

int inv( const int a ) { return qkpow( a, mod - 2 ); }

void sub( int &x, const int v ) { x = ( x < v ? x - v + mod : x - v ); }

void add( int &x, const int v ) { x = ( x + v >= mod ? x + v - mod : x + v ); }



int qkpow( int base, int indx )

{

	int ret = 1;

	while( indx )

	{

		if( indx & 1 ) ret = 1ll * ret * base % mod;

		base = 1ll * base * base % mod, indx >>= 1;

	}

	return ret;

}



void init( const int siz )

{

	fac[0] = ifac[0] = 1;

	for( int i = 1 ; i <= siz ; i ++ ) fac[i] = 1ll * fac[i - 1] * i % mod;

	ifac[siz] = inv( fac[siz] );

	for( int i = siz - 1 ; ~ i ; i -- ) ifac[i] = 1ll * ifac[i + 1] * ( i + 1 ) % mod;

}



int C( const int n, const int m ) { return 1ll * fac[n] * ifac[m] % mod * ifac[n - m] % mod; }



int main()

{

	read( N ), read( K );

	if( K == 1 ) { puts( "1" ); return 0; }

	init( N * K );

	for( int i = 0 ; i <= N ; i ++ ) f[i][0] = 1;

	for( int i = 1 ; i <= N ; i ++ )

		for( int j = 1 ; j <= i ; j ++ )

		{

			f[i][j] = f[i - 1][j];

			int t = N * K - i - ( K - 1 ) * ( j - 1 ) - 1;

			add( f[i][j], 1ll * f[i][j - 1] * ( N - j + 1 ) % mod * C( t, K - 2 ) % mod );

		}

	write( f[N][N] ), putchar( '\n' );

	return 0;

}

// Fastest known solution for this problem:
// Original C++ solution (slower):
#include<bits/stdc++.h>

#define int long long

using namespace std;

const int maxN = 2e3 + 13, mod = 1e9 + 7;

int n, k, dp[maxN][maxN], fac[maxN * maxN], rfac[maxN * maxN];



int pw (int a, int b) {

    int ret = 1;

    for (; b; b >>= 1, a = 1ll * a * a % mod)

	if (b & 1)

	    ret = 1ll * ret * a % mod;

    return ret;

}



int c (int a, int b) {

    return 1ll * fac[b] * rfac[a] % mod * rfac[b - a] % mod;

}



int32_t main () {

    cin >> n >> k;

    dp[0][0] = 1;

    fac[0] = rfac[0] = 1;

    for (int i = 1; i <= n * k; i++) {

	fac[i] = 1ll * fac[i - 1] * i % mod;

	rfac[i] = 1ll * rfac[i - 1] * pw(i, mod - 2) % mod;

    }



    for (int i = 1; i <= n; i++)

	for (int j = 0; j <= i; j++) {

	    if (j)

		dp[i][j] = dp[i][j - 1];

	    dp[i][j] = (dp[i][j] + (1ll * dp[i - 1][j] * c(k - 2, i * (k - 1) + j - 1) % mod)) % mod;

	    // cout << i << " " << j << " " << dp[i][j] << endl;

	}

    //cout << dp[n][n] << endl;

    if (k == 1)

	return cout << 1 << endl, 0;

    cout << 1ll * dp[n][n] * fac[n] % mod << endl;

}


// Optimized solution to the same problem (faster):
#include <cstdio>



const int mod = 1e9 + 7;

const int MAXN = 2005, MAXS = MAXN * MAXN;



template<typename _T>

void read( _T &x )

{

	x = 0;char s = getchar();int f = 1;

	while( s > '9' || s < '0' ){if( s == '-' ) f = -1; s = getchar();}

	while( s >= '0' && s <= '9' ){x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar();}

	x *= f;

}



template<typename _T>

void write( _T x )

{

	if( x < 0 ){ putchar( '-' ); x = ( ~ x ) + 1; }

	if( 9 < x ){ write( x / 10 ); }

	putchar( x % 10 + '0' );

}



int f[MAXN][MAXN];

int fac[MAXS], ifac[MAXS];

int N, K;



int qkpow( int, int );

int inv( const int a ) { return qkpow( a, mod - 2 ); }

void sub( int &x, const int v ) { x = ( x < v ? x - v + mod : x - v ); }

void add( int &x, const int v ) { x = ( x + v >= mod ? x + v - mod : x + v ); }



int qkpow( int base, int indx )

{

	int ret = 1;

	while( indx )

	{

		if( indx & 1 ) ret = 1ll * ret * base % mod;

		base = 1ll * base * base % mod, indx >>= 1;

	}

	return ret;

}



void init( const int siz )

{

	fac[0] = ifac[0] = 1;

	for( int i = 1 ; i <= siz ; i ++ ) fac[i] = 1ll * fac[i - 1] * i % mod;

	ifac[siz] = inv( fac[siz] );

	for( int i = siz - 1 ; ~ i ; i -- ) ifac[i] = 1ll * ifac[i + 1] * ( i + 1 ) % mod;

}



int C( const int n, const int m ) { return 1ll * fac[n] * ifac[m] % mod * ifac[n - m] % mod; }



int main()

{

	read( N ), read( K );

	if( K == 1 ) { puts( "1" ); return 0; }

	init( N * K );

	for( int i = 0 ; i <= N ; i ++ ) f[i][0] = 1;

	for( int i = 1 ; i <= N ; i ++ )

		for( int j = 1 ; j <= i ; j ++ )

		{

			f[i][j] = f[i - 1][j];

			int t = N * K - i - ( K - 1 ) * ( j - 1 ) - 1;

			add( f[i][j], 1ll * f[i][j - 1] * ( N - j + 1 ) % mod * C( t, K - 2 ) % mod );

		}

	write( f[N][N] ), putchar( '\n' );

	return 0;

}

// Fastest known solution for this problem:
// Original C++ solution (slower):


#include<iostream>

int main(){

  int f[2],i,j;

  std::string r,a;

  while(std::cin>>r>>a,r!="0"&&a!="0"){

    for(f[0]=f[1]=i=0;i<4;i++)for(j=4;j-->0;)if(a[i]==r[j])f[i!=j]++;

    std::cout<<f[0]<<" "<<f[1]<<"\n";

  }

}

// Optimized solution to the same problem (faster):
#include<cstdio>

int main(){int h,b,i,j;char r[5],a[5];while(scanf("%s%s",r,a),r[1]!='\0'&&a[1]!='\0'){for(h=b=i=0;i<4;i++)for(j=4;j-->0;)if(a[i]==r[j])i==j?h++:b++;printf("%d %d\n",h,b);}}

// Fastest known solution for this problem:
// Original C++ solution (slower):


#include<iostream>

int main(){

  int f[2],i,j;

  std::string r,a;

  while(std::cin>>r>>a,r!="0"&&a!="0"){

    for(f[0]=f[1]=i=0;i<4;i++)for(j=4;j-->0;)if(a[i]==r[j])f[i!=j]++;

    printf("%d %d\n",f[0],f[1]);

  }

}

// Optimized solution to the same problem (faster):
#include<cstdio>

int main(){int h,b,i,j;char r[5],a[5];while(scanf("%s%s",r,a),r[1]!='\0'&&a[1]!='\0'){for(h=b=i=0;i<4;i++)for(j=4;j-->0;)if(a[i]==r[j])i==j?h++:b++;printf("%d %d\n",h,b);}}

// Fastest known solution for this problem: