stablegradients/pie-gem5-pairs huggingface dataset:
// Slower C++ solution:
#include<iostream>

#include<iomanip>

#include<stack>

#include<queue>

#include<list>

#include<vector>

#include<set>

#include<map>

#include<string>

#include<algorithm>

#include<cmath>

#include<cstdio>

#include<cstring>

#include<cstdlib>

#include<ctime>

#define ll long long

#define ld double

#define inf 2001

#define infm inf*inf

#define mod ((ll)1e9+7)

#define INF (int)(1e9) 

#define pi acos(-1)

#define rd(n) {n=0;char ch;int f=0;do{ch=getchar();if(ch=='-'){f=1;}}while(ch<'0'||ch>'9');while('0'<=ch&&ch<='9'){n=(n<<1)+(n<<3)+ch-48;ch=getchar();}if(f)n=-n;}

using namespace std;



ll getinv(ll x){

	ll k=mod-2,ans=1;

	while (k){

		if (k&1){

			ans=ans*x%mod;

		}

		x=x*x%mod;

		k/=2;

	}

	return ans;

}



ll fact[infm],invf[infm];



void prework(int N){

	fact[0]=1;

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

		fact[i]=fact[i-1]*i%mod;

	}

	invf[N]=getinv(fact[N]);

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

		invf[i]=invf[i+1]*(i+1)%mod;

	}

	return;

}



ll C(int x,int y){

	return fact[x]*invf[y]%mod*invf[x-y]%mod;

}



ll cal(int x,int y){

	return C(x+y,x);

}



int n,K;

int f[inf][inf];



int main(){

	rd(n) rd(K)

	if (K==1){

		puts("1");

		return 0;

	}

	prework(n*K);

	f[0][0]=1;

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

		ll sum=0;

		for (int j=i;j>=1;j--){

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

			f[i][j]=sum*cal((i-1)*K-j+1,K-2)%mod;

		}

	}

	ll ans=0;

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

		ans=(ans+f[n][i])%mod;

	}

	printf("%lld\n",ans*fact[n]%mod);

	return 0;

}




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

#include <cmath>

#include <cstdio>

#include <cstring>



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



int pw(int x, int y) {

	int r = 1;

	for (; y; y >>= 1, x = 1ll * x * x % P)

		if (y & 1) r = 1ll * r * x % P;

	return r;

}



int fac[M], ifac[M];

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



int n, k, dp[N][N];

int main() {

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

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

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

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

	(n * k)[ifac] = pw((n * k)[fac], P - 2);

	for (int i = n * k - 1; i >= 1; i--) i[ifac] = 1ll * (i + 1)[ifac] * (i + 1) % P;

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

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

			if (i || j)

				dp[i][j] = (dp[i - 1][j] +

							(n - j + 1ll) * c(n * k - i - (j - 1) * (k - 1) - 1, k - 2) %

								P * dp[i][j - 1]) %

						   P;

			else

				dp[0][0] = 1;

	printf("%d", dp[n][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: 4.175151610495772x):
#include<cstdio>

#define mo 1000000007

int fac[4000005],inv[4000005],n,k,f[2005][2005];

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

int ksm(int a,int p)

{int ans=1;for(;p;p>>=1,a=1ll*a*a%mo)if(p&1)ans=1ll*ans*a%mo;return ans;}

int main()

{

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

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

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

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

	inv[n*k]=ksm(fac[n*k],mo-2);

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

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

	{

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

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

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

	}printf("%d\n",f[n][n]*1ll*fac[n]%mo);

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

#define mo 1000000007

int fac[4000005],inv[4000005],n,k,f[2005][2005];

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

int ksm(int a,int p)

{int ans=1;for(;p;p>>=1,a=1ll*a*a%mo)if(p&1)ans=1ll*ans*a%mo;return ans;}

int main()

{

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

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

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

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

	inv[n*k]=ksm(fac[n*k],mo-2);

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

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

	{

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

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

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

	}printf("%d\n",f[n][n]*1ll*fac[n]%mo);

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

#define mo 1000000007

int fac[4000005],inv[4000005],n,k,f[2005][2005];

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

int ksm(int a,int p)

{int ans=1;for(;p;p>>=1,a=1ll*a*a%mo)if(p&1)ans=1ll*ans*a%mo;return ans;}

int main()

{

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

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

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

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

	inv[n*k]=ksm(fac[n*k],mo-2);

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

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

	{

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

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

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

	}printf("%d\n",f[n][n]*1ll*fac[n]%mo);

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

#define mo 1000000007

int fac[4000005],inv[4000005],n,k,f[2005][2005];

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

int ksm(int a,int p)

{int ans=1;for(;p;p>>=1,a=1ll*a*a%mo)if(p&1)ans=1ll*ans*a%mo;return ans;}

int main()

{

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

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

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

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

	inv[n*k]=ksm(fac[n*k],mo-2);

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

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

	{

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

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

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

	}printf("%d\n",f[n][n]*1ll*fac[n]%mo);

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

#include<cstring>

#include<algorithm>

#define ll long long

using namespace std;

const int N=2010,M=4000010;

const ll mod=1e9+7;

ll dp[N][N],d1[M],d2[M];

int n,k,max1;

inline ll f(ll y,ll x){

	ll t=y,s=1;if(x%2)s=y;x>>=1;

	while(x){t=(t*t)%mod;if(x%2)s=(s*t)%mod;x>>=1;}

	return s;

}

inline ll c(int m,int n){

	//printf("c%d %d %lld\n",m,n,((d1[n]*d2[m])%mod*d2[n-m])%mod);

	return ((d1[n]*d2[m])%mod*d2[n-m])%mod;

}

int main(){

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

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

	d1[0]=1;d2[0]=1;

	for(int i=1;i<=max1;i++)d1[i]=(d1[i-1]*i)%mod;

	d2[max1]=f(d1[max1],mod-2);//printf("%lld ",d2[max1]);

	for(int i=max1-1;i>=1;i--)d2[i]=(d2[i+1]*(i+1))%mod;

	dp[1][0]=1;dp[1][1]=(c(k-2,max1-2)*n)%mod;//printf("1 1 %lld\n",dp[1][1]);

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

		dp[i][0]=1;

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

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

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

		}

	}

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

	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.4800398501181684x):
#include<cstdio>

#include<cstring>

#include<algorithm>

#define ll long long

using namespace std;

const int N=2010,M=4000010;

const ll mod=1e9+7;

ll dp[N][N],d1[M],d2[M];

int n,k,max1;

inline ll f(ll y,ll x){

	ll t=y,s=1;if(x%2)s=y;x>>=1;

	while(x){t=(t*t)%mod;if(x%2)s=(s*t)%mod;x>>=1;}

	return s;

}

inline ll c(int m,int n){

	//printf("c%d %d %lld\n",m,n,((d1[n]*d2[m])%mod*d2[n-m])%mod);

	return ((d1[n]*d2[m])%mod*d2[n-m])%mod;

}

int main(){

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

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

	d1[0]=1;d2[0]=1;

	for(int i=1;i<=max1;i++)d1[i]=(d1[i-1]*i)%mod;

	d2[max1]=f(d1[max1],mod-2);//printf("%lld ",d2[max1]);

	for(int i=max1-1;i>=1;i--)d2[i]=(d2[i+1]*(i+1))%mod;

	dp[1][0]=1;dp[1][1]=(c(k-2,max1-2)*n)%mod;//printf("1 1 %lld\n",dp[1][1]);

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

		dp[i][0]=1;

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

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

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

		}

	}

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

	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.3356746862740745x):
#include<cstdio>

#include<cstring>

#include<algorithm>

#define ll long long

using namespace std;

const int N=2010,M=4000010;

const ll mod=1e9+7;

ll dp[N][N],d1[M],d2[M];

int n,k,max1;

inline ll f(ll y,ll x){

	ll t=y,s=1;if(x%2)s=y;x>>=1;

	while(x){t=(t*t)%mod;if(x%2)s=(s*t)%mod;x>>=1;}

	return s;

}

inline ll c(int m,int n){

	//printf("c%d %d %lld\n",m,n,((d1[n]*d2[m])%mod*d2[n-m])%mod);

	return ((d1[n]*d2[m])%mod*d2[n-m])%mod;

}

int main(){

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

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

	d1[0]=1;d2[0]=1;

	for(int i=1;i<=max1;i++)d1[i]=(d1[i-1]*i)%mod;

	d2[max1]=f(d1[max1],mod-2);//printf("%lld ",d2[max1]);

	for(int i=max1-1;i>=1;i--)d2[i]=(d2[i+1]*(i+1))%mod;

	dp[1][0]=1;dp[1][1]=(c(k-2,max1-2)*n)%mod;//printf("1 1 %lld\n",dp[1][1]);

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

		dp[i][0]=1;

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

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

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

		}

	}

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

	return 0;

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

#include<cstring>

#include<algorithm>

#define ll long long

using namespace std;

const int N=2010,M=4000010;

const ll mod=1e9+7;

ll dp[N][N],d1[M],d2[M];

int n,k,max1;

inline ll f(ll y,ll x){

	ll t=y,s=1;if(x%2)s=y;x>>=1;

	while(x){t=(t*t)%mod;if(x%2)s=(s*t)%mod;x>>=1;}

	return s;

}

inline ll c(int m,int n){

	//printf("c%d %d %lld\n",m,n,((d1[n]*d2[m])%mod*d2[n-m])%mod);

	return ((d1[n]*d2[m])%mod*d2[n-m])%mod;

}

int main(){

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

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

	d1[0]=1;d2[0]=1;

	for(int i=1;i<=max1;i++)d1[i]=(d1[i-1]*i)%mod;

	d2[max1]=f(d1[max1],mod-2);//printf("%lld ",d2[max1]);

	for(int i=max1-1;i>=1;i--)d2[i]=(d2[i+1]*(i+1))%mod;

	dp[1][0]=1;dp[1][1]=(c(k-2,max1-2)*n)%mod;//printf("1 1 %lld\n",dp[1][1]);

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

		dp[i][0]=1;

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

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

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

		}

	}

	printf("%lld",dp[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: 3.639299024486953x):
#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

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

}