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.730824754425739x):
#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;

}
// 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.432215594845853x):
#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;

}
// 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.4700684586175576x):
#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;

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

}

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



int f[2001][2001],n,k,invmul[4000001],_mul[4000001];

const int P=1000000007;

inline int mul(const int &a,const int &b){return 1ll*a*b%P;}

inline int add(int a,const int &b){a+=b;return(a>=P)?a-P:a;}

inline int sub(int a,const int &b){a-=b;return (a<0)?a+P:a;}

inline int C(const int &n,const int &m){return mul(mul(_mul[n],invmul[m]),invmul[n-m]);}

int main(){

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

	if(k==1){

		puts("1");

		return 0;

	}

	_mul[0]=_mul[1]=invmul[0]=invmul[1]=1;

	for(int i=2;i<=n*k;i++)_mul[i]=mul(_mul[i-1],i);

	for(int i=2;i<=n*k;i++)invmul[i]=mul(P-P/i,invmul[P%i]);

	for(int i=2;i<=n*k;i++)invmul[i]=mul(invmul[i-1],invmul[i]);

	f[0][0]=1;

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

		f[i][0]=1;

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

			f[i][j]=add((j<i)?f[i-1][j]:0,mul(f[i][j-1],C(n*k-i-(j-1)*(k-1)-1,k-2)));

	}

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

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

#define LL long long

#define RG register

#define R RG int

using namespace std;

const LL N=2009,M=N*N,YL=1e9+7;

int f[N][N],fac[M],inv[M];

inline int fm(RG LL b){

	RG LL a=1;

	for(R k=YL-2;k;k>>=1,b=b*b%YL)

		if(k&1)a=a*b%YL;

	return a;

}

inline int C(R n,R m){

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

	return(LL)fac[n]*inv[m]%YL*inv[n-m]%YL;

}

int main(){

	R n,k,nk;

	scanf("%d%d",&n,&k);nk=n*k;

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

	for(R i=fac[0]=1;i<=nk;++i)

		fac[i]=fac[i-1]*(LL)i%YL;

	inv[nk]=fm(fac[nk]);

	for(R i=nk;i;--i)

		inv[i-1]=inv[i]*(LL)i%YL;

	f[0][0]=1;

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

		for(R j=f[i][0]=1;j<=i;++j)

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

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

	return 0;

}

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



int f[2001][2001],n,k,invmul[4000001],_mul[4000001];

const int P=1000000007;

inline int mul(const int &a,const int &b){return 1ll*a*b%P;}

inline int add(int a,const int &b){a+=b;return(a>=P)?a-P:a;}

inline int sub(int a,const int &b){a-=b;return (a<0)?a+P:a;}

inline int C(const int &n,const int &m){return mul(mul(_mul[n],invmul[m]),invmul[n-m]);}

int main(){

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

	if(k==1){

		puts("1");

		return 0;

	}

	_mul[0]=_mul[1]=invmul[0]=invmul[1]=1;

	for(int i=2;i<=n*k;i++)_mul[i]=mul(_mul[i-1],i);

	for(int i=2;i<=n*k;i++)invmul[i]=mul(P-P/i,invmul[P%i]);

	for(int i=2;i<=n*k;i++)invmul[i]=mul(invmul[i-1],invmul[i]);

	f[0][0]=1;

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

		f[i][0]=1;

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

			f[i][j]=add((j<i)?f[i-1][j]:0,mul(f[i][j-1],C(n*k-i-(j-1)*(k-1)-1,k-2)));

	}

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

}
// Slower C++ solution:
#include<iostream>

#include<algorithm>

#include<cstdio>

#include<cstring>

#define Mod 1000000007

using namespace std;

long long f[2010][2010];

int n,k;

long long fac[4000010],inv[4000010];

long long fpow(long long a,long long b)

{

	long long ans=1,t=a;

	while(b)

	{

		if(b&1)ans=ans*t%Mod;

		t=t*t%Mod;

		b>>=1;

	}

	return ans;

}

void getfac()

{

	fac[0]=1;

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

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

	inv[n*k]=fpow(fac[n*k],Mod-2);

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

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

	return ;

}

long long C(int a,int b)

{

	return fac[a]*inv[b]%Mod*inv[a-b]%Mod;

}

int main()

{

	cin>>n>>k;

	if(k==1){

		cout<<1<<"\n";

		return 0;

	}

	getfac();

	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]+f[i][j-1]*(n-j+1)%Mod*C(n-i+(k-1)*(n-j+1)-1,k-2)%Mod)%Mod;

		}

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

	return 0;

}

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



int f[2001][2001],n,k,invmul[4000001],_mul[4000001];

const int P=1000000007;

inline int mul(const int &a,const int &b){return 1ll*a*b%P;}

inline int add(int a,const int &b){a+=b;return(a>=P)?a-P:a;}

inline int sub(int a,const int &b){a-=b;return (a<0)?a+P:a;}

inline int C(const int &n,const int &m){return mul(mul(_mul[n],invmul[m]),invmul[n-m]);}

int main(){

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

	if(k==1){

		puts("1");

		return 0;

	}

	_mul[0]=_mul[1]=invmul[0]=invmul[1]=1;

	for(int i=2;i<=n*k;i++)_mul[i]=mul(_mul[i-1],i);

	for(int i=2;i<=n*k;i++)invmul[i]=mul(P-P/i,invmul[P%i]);

	for(int i=2;i<=n*k;i++)invmul[i]=mul(invmul[i-1],invmul[i]);

	f[0][0]=1;

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

		f[i][0]=1;

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

			f[i][j]=add((j<i)?f[i-1][j]:0,mul(f[i][j-1],C(n*k-i-(j-1)*(k-1)-1,k-2)));

	}

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

}
// 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: 2.60668920245542x):
#include <cstdio>



int f[2001][2001],n,k,invmul[4000001],_mul[4000001];

const int P=1000000007;

inline int mul(const int &a,const int &b){return 1ll*a*b%P;}

inline int add(int a,const int &b){a+=b;return(a>=P)?a-P:a;}

inline int sub(int a,const int &b){a-=b;return (a<0)?a+P:a;}

inline int C(const int &n,const int &m){return mul(mul(_mul[n],invmul[m]),invmul[n-m]);}

int main(){

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

	if(k==1){

		puts("1");

		return 0;

	}

	_mul[0]=_mul[1]=invmul[0]=invmul[1]=1;

	for(int i=2;i<=n*k;i++)_mul[i]=mul(_mul[i-1],i);

	for(int i=2;i<=n*k;i++)invmul[i]=mul(P-P/i,invmul[P%i]);

	for(int i=2;i<=n*k;i++)invmul[i]=mul(invmul[i-1],invmul[i]);

	f[0][0]=1;

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

		f[i][0]=1;

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

			f[i][j]=add((j<i)?f[i-1][j]:0,mul(f[i][j-1],C(n*k-i-(j-1)*(k-1)-1,k-2)));

	}

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

}
// 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.1901753559974595x):
#include <cstdio>



int f[2001][2001],n,k,invmul[4000001],_mul[4000001];

const int P=1000000007;

inline int mul(const int &a,const int &b){return 1ll*a*b%P;}

inline int add(int a,const int &b){a+=b;return(a>=P)?a-P:a;}

inline int sub(int a,const int &b){a-=b;return (a<0)?a+P:a;}

inline int C(const int &n,const int &m){return mul(mul(_mul[n],invmul[m]),invmul[n-m]);}

int main(){

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

	if(k==1){

		puts("1");

		return 0;

	}

	_mul[0]=_mul[1]=invmul[0]=invmul[1]=1;

	for(int i=2;i<=n*k;i++)_mul[i]=mul(_mul[i-1],i);

	for(int i=2;i<=n*k;i++)invmul[i]=mul(P-P/i,invmul[P%i]);

	for(int i=2;i<=n*k;i++)invmul[i]=mul(invmul[i-1],invmul[i]);

	f[0][0]=1;

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

		f[i][0]=1;

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

			f[i][j]=add((j<i)?f[i-1][j]:0,mul(f[i][j-1],C(n*k-i-(j-1)*(k-1)-1,k-2)));

	}

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

}
// 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.213044520463744x):
#include <cstdio>



int f[2001][2001],n,k,invmul[4000001],_mul[4000001];

const int P=1000000007;

inline int mul(const int &a,const int &b){return 1ll*a*b%P;}

inline int add(int a,const int &b){a+=b;return(a>=P)?a-P:a;}

inline int sub(int a,const int &b){a-=b;return (a<0)?a+P:a;}

inline int C(const int &n,const int &m){return mul(mul(_mul[n],invmul[m]),invmul[n-m]);}

int main(){

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

	if(k==1){

		puts("1");

		return 0;

	}

	_mul[0]=_mul[1]=invmul[0]=invmul[1]=1;

	for(int i=2;i<=n*k;i++)_mul[i]=mul(_mul[i-1],i);

	for(int i=2;i<=n*k;i++)invmul[i]=mul(P-P/i,invmul[P%i]);

	for(int i=2;i<=n*k;i++)invmul[i]=mul(invmul[i-1],invmul[i]);

	f[0][0]=1;

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

		f[i][0]=1;

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

			f[i][j]=add((j<i)?f[i-1][j]:0,mul(f[i][j-1],C(n*k-i-(j-1)*(k-1)-1,k-2)));

	}

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

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

}

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

#define RG register

#define R RG int

using namespace std;

const LL N=2009,M=N*N,YL=1e9+7;

int f[N][N],fac[M],inv[M];

inline int fm(RG LL b){

	RG LL a=1;

	for(R k=YL-2;k;k>>=1,b=b*b%YL)

		if(k&1)a=a*b%YL;

	return a;

}

inline int C(R n,R m){

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

	return(LL)fac[n]*inv[m]%YL*inv[n-m]%YL;

}

int main(){

	R n,k,nk;

	scanf("%d%d",&n,&k);nk=n*k;

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

	for(R i=fac[0]=1;i<=nk;++i)

		fac[i]=fac[i-1]*(LL)i%YL;

	inv[nk]=fm(fac[nk]);

	for(R i=nk;i;--i)

		inv[i-1]=inv[i]*(LL)i%YL;

	f[0][0]=1;

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

		for(R j=f[i][0]=1;j<=i;++j)

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

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

	return 0;

}

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

#include<cstdio>

#include<cstring>

#define Mod 1000000007

using namespace std;

long long f[2010][2010];

int n,k;

long long fac[4000010],inv[4000010];

long long fpow(long long a,long long b)

{

	long long ans=1,t=a;

	while(b)

	{

		if(b&1)ans=ans*t%Mod;

		t=t*t%Mod;

		b>>=1;

	}

	return ans;

}

void getfac()

{

	fac[0]=1;

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

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

	inv[n*k]=fpow(fac[n*k],Mod-2);

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

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

	return ;

}

long long C(int a,int b)

{

	return fac[a]*inv[b]%Mod*inv[a-b]%Mod;

}

int main()

{

	cin>>n>>k;

	if(k==1){

		cout<<1<<"\n";

		return 0;

	}

	getfac();

	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]+f[i][j-1]*(n-j+1)%Mod*C(n-i+(k-1)*(n-j+1)-1,k-2)%Mod)%Mod;

		}

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

	return 0;

}

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

}