stablegradients/pie-gem5-pairs huggingface dataset:
// Slower C++ solution:
#include <bits/stdc++.h>

using namespace std;



typedef long long lint;

const int maxn = 2005, mod = 1e9 + 7;



int n, k, m, fac[maxn * maxn], ifac[maxn * maxn], f[maxn][maxn];



inline int gi()

{

	char c = getchar();

	while (c < '0' || c > '9') c = getchar();

	int sum = 0;

	while ('0' <= c && c <= '9') sum = sum * 10 + c - 48, c = getchar();

	return sum;

}



inline void inc(int &a, int b) {a += b; if (a >= mod) a -= mod;}

inline int C(int n, int m)

{

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

}



int main()

{

	n = gi(); k = gi();

	if (k == 1) return puts("1"), 0;



	m = n * k;

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

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

	for (int i = 2; i <= m; ++i) ifac[i] = (lint)(mod - mod / i) * ifac[mod % i] % mod;

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



	f[1][0] = 1;

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

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

			if (i > j) inc(f[i][j], f[i - 1][j]);

			if (j) inc(f[i][j], (lint)f[i][j - 1] * C(m - i - (j - 1) * (k - 1) - 1, k - 2) % mod);

		}



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

	

	return 0;

}


// Faster C++ solution for the same problem (measured speedup: 3.653043815885434x):
#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;

}
// Slower C++ solution:
#include <bits/stdc++.h>

using namespace std;



#define ull unsigned long long



const int MaxN(2003);

const int MaxL(4000003);

const int Mod(1e9 + 7);



int F[MaxN][MaxN];

int fac[MaxL], rfac[MaxL];



int fexp(int x, int k)

{

	int res = 1;

	for(; k; k >>= 1, x = (ull) x * x % Mod)

		if(k & 1) res = (ull) res * x % Mod;

	return res;

}



inline int C(int n, int m)

{

	return n < m ? 0 : (ull) fac[n] * rfac[m] % Mod * rfac[n - m] % Mod;

}



int main()

{

	int n, K;

#ifdef zxp

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

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

#endif

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

	if(K == 1) return puts("1"), 0;

	for(int i = fac[0] = 1; i <= n * K; i++)

		fac[i] = (ull) fac[i - 1] * i % Mod;

	rfac[n * K] = fexp(fac[n * K], Mod - 2);

	for(int i = n * K; i; i--)

		rfac[i - 1] = (ull) rfac[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) F[i][j] = (F[i][j] + (ull) (n - j + 1) * C((K - 1) * (n - j + 1) + n - i - 1, K - 2) % Mod * F[i][j - 1]) % Mod;

		}

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

	return 0;

}


// Faster C++ solution for the same problem (measured speedup: 2.3678667857524536x):
#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;

}
// Slower C++ solution:
#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;

}




// Faster C++ solution for the same problem (measured speedup: 2.404718181270421x):
#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;

}
// Slower C++ solution:
#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;

}


// Faster C++ solution for the same problem (measured speedup: 3.26831997397089x):
#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;

}
// Slower C++ solution:
#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;

}


// Faster C++ solution for the same problem (measured speedup: 3.2693397251009877x):
#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;

}
// Slower C++ solution:
#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]);

}


// Faster C++ solution for the same problem (measured speedup: 2.4510641994145024x):
#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;

}
// Slower C++ solution:
#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 ; 

}

// Faster C++ solution for the same problem (measured speedup: 5.281711678228875x):
#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;

}
// Slower C++ solution:
#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);

}

// Faster C++ solution for the same problem (measured speedup: 2.227654987580682x):
#include<cstdio>

#define mod 1000000007

#define maxn 2005

int n,k,fact[maxn*maxn],inv[maxn*maxn],dp[maxn][maxn];

void init()

{

	fact[0]=1;

	for(int i=1;i<=n*k;i++) fact[i]=1ll*fact[i-1]*i%mod;

	inv[1]=1;

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

	inv[0]=1;

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

}

int main()

{

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

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

	init();

	dp[0][0]=1;

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

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

		{

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

			if(j)

			{

				int x=n*k-i-(j-1)*(k-1)-1;

				dp[i][j]=(dp[i][j]+1ll*dp[i][j-1]*fact[x]%mod*inv[k-2]%mod*inv[x-k+2]%mod)%mod;

			}

		}

	int ans=1ll*dp[n][n]*fact[n]%mod;

	printf("%d\n",ans);

}
// Slower C++ solution:
#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;

}


// Faster C++ solution for the same problem (measured speedup: 2.2022739422205895x):
#include<cstdio>

#define mod 1000000007

#define maxn 2005

int n,k,fact[maxn*maxn],inv[maxn*maxn],dp[maxn][maxn];

void init()

{

	fact[0]=1;

	for(int i=1;i<=n*k;i++) fact[i]=1ll*fact[i-1]*i%mod;

	inv[1]=1;

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

	inv[0]=1;

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

}

int main()

{

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

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

	init();

	dp[0][0]=1;

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

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

		{

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

			if(j)

			{

				int x=n*k-i-(j-1)*(k-1)-1;

				dp[i][j]=(dp[i][j]+1ll*dp[i][j-1]*fact[x]%mod*inv[k-2]%mod*inv[x-k+2]%mod)%mod;

			}

		}

	int ans=1ll*dp[n][n]*fact[n]%mod;

	printf("%d\n",ans);

}
// Slower C++ solution:
#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);

}

// Faster C++ solution for the same problem (measured speedup: 3.09436171518897x):
#include<cstdio>

const int N=2007,P=1000000007;

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

void inc(int&a,int b){a+=b-P,a+=a>>31&P;}

void mul(int&a,int b){a=1ll*a*b%P;}

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

int pow(int a,int b){int r=1;for(;b;b>>=1,mul(a,a))if(b&1)mul(r,a);return r;}

int main()

{

    scanf("%d%d",&n,&k),m=n*k,f[0][0]=1;

    if(k==1) return puts("1"),0;

    fac[0]=1;for(int i=1;i<=m;++i)mul(fac[i]=fac[i-1],i);

    ifac[m]=pow(fac[m],P-2);for(int i=m;i;--i)mul(ifac[i-1]=ifac[i],i);

    for(int i=0,z=k-1;i<=n;++i,z+=k-1)

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

	    inc(f[i+1][j],1ll*f[i][j]*C(z+j-1,k-2)%P),inc(f[i][j+1],f[i][j]);

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

}
// Slower C++ solution:
#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;

}


// Faster C++ solution for the same problem (measured speedup: 3.0381773968277823x):
#include<cstdio>

const int N=2007,P=1000000007;

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

void inc(int&a,int b){a+=b-P,a+=a>>31&P;}

void mul(int&a,int b){a=1ll*a*b%P;}

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

int pow(int a,int b){int r=1;for(;b;b>>=1,mul(a,a))if(b&1)mul(r,a);return r;}

int main()

{

    scanf("%d%d",&n,&k),m=n*k,f[0][0]=1;

    if(k==1) return puts("1"),0;

    fac[0]=1;for(int i=1;i<=m;++i)mul(fac[i]=fac[i-1],i);

    ifac[m]=pow(fac[m],P-2);for(int i=m;i;--i)mul(ifac[i-1]=ifac[i],i);

    for(int i=0,z=k-1;i<=n;++i,z+=k-1)

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

	    inc(f[i+1][j],1ll*f[i][j]*C(z+j-1,k-2)%P),inc(f[i][j+1],f[i][j]);

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

}