stablegradients/pie-gem5-pairs huggingface dataset:
// 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.1752911211632435x):
#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 IL inline

#define RG register

#define LL long long

#define Upd(a) (a)>=mod?(a-=mod):1

using namespace std;

template <class T>

IL char gc () {

	static char buf[1<<6], *p1=buf, *p2=buf;

	return (p1==p2) && (p2=(p1=buf)+fread(buf,1,1<<6,stdin), p1==p2) ?EOF:*p1++;

}

template <typename T>

IL void read (T&data) {

	data=0;

	RG char ch=0;

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

	while (ch<='9' && ch>='0') data=(data<<1)+(data<<3)+(ch&15), ch=getchar();

}

const int mod = 1e9+7;

const int _ =2011;

int n, k;

int dp[2][_], jc[_*_], jcn[_*_];

IL int poww (RG int a, RG int b) {

	RG int ret=1, Base=a;

	for (;b;b>>=1, Base=1LL*Base*Base%mod) if (b&1) ret=1LL*ret*Base%mod;

	return ret;

}

IL int C (RG int a, RG int b) {

	if (a<b) return 0;

	return 1LL*jc[a]*jcn[b]%mod*jcn[a-b]%mod;

}

int main (/*int aa, int bb*/) {

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

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

	//n=aa, k=bb;

	read(n), read(k);

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

	jc[0]=jcn[0]=1;

	RG int i, j, pre=1, cur=0;

	for (i=1; i<=n*k; ++i) jc[i]=1LL*jc[i-1]*i%mod, jcn[i]=poww(jc[i], mod-2);

	dp[0][0]=1;

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

		swap(cur, pre);

		for (j=0; j<min(i,n+1); ++j) dp[cur][j]=0;

		for (j=1; j<=min(i,n); ++j) {

			if ((i-j)%2) continue;

			if (n*2-i+1<j) break;

			RG int a=(i+j-2)>>1, b=(i-j)>>1;

			dp[cur][j]=1LL*dp[pre][j-1]*C(a*(k-1)+b+k-1-1,k-1-1)%mod;

			dp[cur][j]+=dp[pre][j+1];

			Upd(dp[cur][j]);

		}

		if (cur%2==0) dp[cur][0]=dp[pre][1];

	}

	//return 1LL*dp[cur][0]*jc[n]%mod;

	printf("%lld",1LL*dp[cur][0]*jc[n]%mod);

}


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

}


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

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

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

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

}


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

}
// Slower C++ solution:


#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";

  }

}

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


#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]);

  }

}

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

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

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

int main(){int h,b,i,j;std::string r,a;while(std::cin>>r>>a,r!="0"&&a!="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);}}

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

int main(){int h,b,i,j;std::string r,a;while(std::cin>>r>>a,r!="0"&&a!="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);}}

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

// Faster C++ solution for the same problem (measured speedup: 5.197925748857062x):
#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);}}
// Slower C++ solution:
#include <bits/stdc++.h>

using namespace std;



int h, w, n;

map<pair<int,int> ,int> mp;

set<pair<int,int> > s[10];



void add(int x,int y) {

	int lx = max(1, x - 2), ly = max(1, y - 2);

	int rx = min(x, h - 2), ry = min(y, w - 2);

	for (x = lx; x <= rx; ++x) {

		for (y = ly; y <= ry; ++y) {

			if (mp[{x, y}]) s[mp[{x, y}]].erase({x, y});

			mp[{x, y}]++;

			s[mp[{x, y}]].insert({x, y});

		}

	}

}



int main() {

	ios_base::sync_with_stdio(false);

	cin >> h >> w >> n;

	while (n--) {

		int x, y; cin >> x >> y;

		add(x, y);

	}

	long long sum = 1LL * (h - 2) * (w - 2);

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

		sum -= (int)s[i].size();

	}

	cout << sum << '\n';

	for (int i = 1; i <= 9; ++i) cout << (int)s[i].size() << '\n';

}

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

#include<algorithm>

#include<map>

using namespace std;

int n,m,q;

long long ans[15];

map<int,map<int,int> >mp;

int main()

{

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

	while(q--)

	{

		int x,y;

		scanf("%d%d",&x,&y);

		for(int i=max(2,x-1);i<=min(n-1,x+1);i++)

			for(int j=max(2,y-1);j<=min(m-1,y+1);j++)

			{

				ans[mp[i][j]]-=1ll;

				mp[i][j]++;

				ans[mp[i][j]]+=1ll;

			}

	}

	ans[0]+=1ll*(n-2)*(m-2);

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

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

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



using namespace std;



int h, w;

int n;



int a[100005], b[100005];

long long c[10];

#define x first

#define y second

signed main() {

    ios::sync_with_stdio(false); cin.tie(0);

    cin >> h >> w >> n;

    memset(c, 0ll, sizeof c);

    set< pair<int, int> > av;

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

        cin >> a[i] >> b[i];

        av.insert({a[i], b[i]});

    }

    set< pair<int, int> > canxet;

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

        for (int j = -2; j <= 2; j++) {

            for (int k = -2; k <= 2; k++) {

                if (a[i] + j >= 1 && a[i]+j <= h-2 && b[i] + k >= 1 && b[i] + k <= w-2) {

                    canxet.insert({a[i]+j, b[i]+k});

                }

            }

        }

    }

    for (set< pair<int, int> >::iterator i = canxet.begin(); i != canxet.end(); i++) {

        int cnt = 0;

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

            for (int k = 0; k <= 2; k++) {

                if (av.find({(*i).x + j, (*i).y + k}) != av.end()) cnt++;

            }

        }

        c[cnt]++;

    }

    c[0] = 1ll*(w-2)*(h-2);

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

        c[0] -= c[i];

    }

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

        cout << c[i] << endl;

    }

}


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

#include<algorithm>

#include<map>

using namespace std;

int n,m,q;

long long ans[15];

map<int,map<int,int> >mp;

int main()

{

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

	while(q--)

	{

		int x,y;

		scanf("%d%d",&x,&y);

		for(int i=max(2,x-1);i<=min(n-1,x+1);i++)

			for(int j=max(2,y-1);j<=min(m-1,y+1);j++)

			{

				ans[mp[i][j]]-=1ll;

				mp[i][j]++;

				ans[mp[i][j]]+=1ll;

			}

	}

	ans[0]+=1ll*(n-2)*(m-2);

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

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

}