#include <bits/stdc++.h>
using namespace std;
#define rep(i,a,n) for (int i=a;i<n;i++)
#define per(i,a,n) for (int i=n-1;i>=a;i--)
#define pb push_back
#define mp make_pair
#define all(x) (x).begin(),(x).end()
#define fi first
#define se second
#define SZ(x) ((int)(x).size())
typedef vector<int> VI;
typedef long long ll;
typedef pair<int,int> PII;
mt19937 mrand(random_device{}()); 
int rnd(int x) { return mrand() % x;}
ll gcd(ll a,ll b) { return b?gcd(b,a%b):a;}
// head

const int N=110;

int n,m,k;

typedef int i32;
typedef unsigned int u32;
typedef long long i64;
typedef unsigned long long u64;
typedef __int128_t i128;
typedef __uint128_t u128;

struct Mod32 {
	Mod32():n_(0) {}
	Mod32(u32 n):n_(init(n)) {}
	static u32 init(u32 w) { return reduce(u64(w) * r2); }
	static void set_mod(u32 m) {
		mod=m; assert(mod&1);
		inv=m; rep(i,0,4) inv*=2-inv*m;
		r2=-u64(m)%m;
	}
	static u32 reduce(u64 x) {
		u32 y=u32(x>>32)-u32((u64(u32(x)*inv)*mod)>>32);
		return i32(y)<0?y+mod:y;
	}
	Mod32& operator += (Mod32 rhs) { n_+=rhs.n_-mod; if (i32(n_)<0) n_+=mod; return *this; }
	Mod32 operator + (Mod32 rhs) const { return Mod32(*this)+=rhs; }
	Mod32& operator -= (Mod32 rhs) { n_-=rhs.n_; if (i32(n_)<0) n_+=mod; return *this; }
	Mod32 operator - (Mod32 rhs) const { return Mod32(*this)-=rhs; }
	Mod32& operator *= (Mod32 rhs) { n_=reduce(u64(n_)*rhs.n_); return *this; }
	Mod32 operator * (Mod32 rhs) const { return Mod32(*this)*=rhs; }
	u32 get() const { return reduce(n_); }
	static u32 mod,inv,r2;
	u32 n_;
};
u32 Mod32::mod,Mod32::inv,Mod32::r2;
Mod32 one,zero;
ll mod;

Mod32 powmod(Mod32 a,ll b) {
	Mod32 res(1);
	for(;b;b>>=1) {
		if(b&1) res*=a;
		a*=a;
	}
	return res;
}

Mod32 fac[N],fnv[N],sfnv[N],comb[N][N],dp[N][N][N],pd[N][N],tmp[N][N];
Mod32 pw[N][N*N],ipw[N][N*N];


int main() {
//	scanf("%d%d%d%lld",&n,&m,&k,&mod);
	scanf("%d%d%d%lld",&n,&m,&k,&mod);
	Mod32::set_mod(mod); one=Mod32(1); zero=Mod32(0);
	fac[0]=one;
	rep(i,1,101) fac[i]=fac[i-1]*Mod32(i);
	rep(i,0,101) {
		fnv[i]=powmod(fac[i],mod-2);
		sfnv[i]=(i%2==0)?fnv[i]:zero-fnv[i];
	}
	comb[0][0]=Mod32(1);
	rep(i,1,101) {
		comb[i][0]=comb[i][i]=one;
		rep(j,1,i) {
			comb[i][j]=comb[i-1][j]+comb[i-1][j-1];
		}
	}
	dp[0][0][0]=fac[n]*fac[m];
	rep(i,0,101) {
		Mod32 pi(i);
		pw[i][0]=one;
		rep(j,1,10001) pw[i][j]=pw[i][j-1]*pi;
		Mod32 pv=powmod(pi,mod-2);
		ipw[i][0]=one;
		rep(j,1,10001) ipw[i][j]=ipw[i][j-1]*pv;
	}
	rep(i,0,k) {
		rep(p1,0,n+1) rep(p2,0,m+1) {
			pd[p1][p2]=zero;
			dp[i][p1][p2]=(dp[i][p1][p2]*
				pw[i+1][(n-p1)*(m-p2)]*ipw[k-i][p1*p2]);
		}
		rep(p1,0,n+1) rep(p2,0,m+1) tmp[p1][p2]=zero;
		rep(p1,0,n+1) rep(p2,0,m+1) for (int u=0;p1+u<=n;u++) {
			tmp[p1+u][p2]+=dp[i][p1][p2]*fnv[u];
		}
		rep(p1,0,n+1) rep(p2,0,m+1) for (int v=0;p2+v<=m;v++) {
			pd[p1][p2+v]+=tmp[p1][p2]*fnv[v];
		}
		rep(p1,0,n+1) rep(p2,0,m+1) {
			pd[p1][p2]=(pd[p1][p2]*pw[k-i][p1*p2]);
		}
		if (i==k-1) {
			rep(p1,0,n+1) rep(p2,0,m+1) dp[i+1][p1][p2]=pd[p1][p2];
			continue;
		}
		rep(p1,0,n+1) rep(p2,0,m+1) {
			pd[p1][p2]=pd[p1][p2]*ipw[k-i-1][p1*p2];
			tmp[p1][p2]=zero;
		}
		rep(p1,0,n+1) rep(p2,0,m+1) for (int u=0;p1+u<=n;u++) {
			tmp[p1+u][p2]+=pd[p1][p2]*sfnv[u];
		}
		rep(p1,0,n+1) rep(p2,0,m+1) for (int v=0;p2+v<=m;v++) {
			dp[i+1][p1][p2+v]+=tmp[p1][p2]*sfnv[v];
		}
		rep(p1,0,n+1) rep(p2,0,m+1) {
			dp[i+1][p1][p2]=(dp[i+1][p1][p2]*ipw[i+1][(n-p1)*(m-p2)]*
				pw[k-i-1][p1*p2]);
		}
	}
	printf("%d\n",dp[k][n][m].get());
}

Measured memory_bytes: 16128