Problem: For each of the K^{NM} ways to write an integer between 1 and K (inclusive) in every square in a square grid with N rows and M columns, find the value defined below, then compute the sum of all those K^{NM} values, modulo D.

* For each of the NM squares, find the minimum among the N+M-1 integers written in the square's row or the square's column. The value defined for the grid is the product of all these NM values.

Constraints

* 1 \leq N,M,K \leq 100
* 10^8 \leq D \leq 10^9
* N,M,K, and D are integers.
* D is prime.

Input

Input is given from Standard Input in the following format:


N M K D


Output

Print the sum of the K^{NM} values, modulo D.

Examples

Input

2 2 2 998244353


Output

35


Input

2 3 4 998244353


Output

127090


Input

31 41 59 998244353


Output

827794103

Solution (in CPP):

Solution:
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<bitset>
#include<set>
#define N (110)
#define inf (0x7f7f7f7f)
#define rg register int
#define Label puts("NAIVE")
#define spa print(' ')
#define ent print('\n')
#define rand() (((rand())<<(15))^(rand()))
#define file(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
typedef long double ld;
typedef long long LL;
typedef unsigned long long ull;
using namespace std;
namespace fastIO1{
	inline char read(){
		static const int IN_LEN=1000000;
		static char buf[IN_LEN],*s,*t;
		return (s==t?t=(s=buf)+fread(buf,1,IN_LEN,stdin),(s==t?-1:*s++):*s++);
	}
	template<class T>
	inline void read(T &x){
		static bool iosig;
		static char c;
		for(iosig=false,c=read();!isdigit(c);c=read()){
			if(c=='-')iosig=true;
			if(c==-1)return;
		}
		for(x=0;isdigit(c);c=read())x=((x+(x<<2))<<1)+(c^'0');
		if(iosig)x=-x;
	}
	inline char readc(char &c){
		for(c=read();!isalpha(c)&&!isdigit(c);c=read())
		if(c==-1)return 0;
	}
	const int OUT_LEN = 10000000;
	char obuf[OUT_LEN],*ooh=obuf;
	inline void print(char c){
		if(ooh==obuf+OUT_LEN)fwrite(obuf,1,OUT_LEN,stdout),ooh=obuf;
		*ooh++=c;
	}
	template<class T>
	inline void print(T x){
		static int buf[30],cnt;
		if(x==0)print('0');
		else{
			if(x<0)print('-'),x=-x;
			for(cnt=0;x;x/=10)buf[++cnt]=x%10+48;
			while(cnt)print((char)buf[cnt--]);
		}
	}
	inline void flush(){fwrite(obuf,1,ooh-obuf,stdout);}
}
namespace fastIO2{
	template<class T>
	inline void read(T &x){
		static bool iosig;
		static char c;
		for(iosig=false,c=getchar();!isdigit(c);c=getchar()){
			if(c=='-')iosig=true;
			if(c==-1)return;
		}
		for(x=0;isdigit(c);c=getchar())x=((x+(x<<2))<<1)+(c^'0');
		if(iosig)x=-x;
	}
}
using namespace fastIO1;
int n,m,K,P;
int pw[N][N*N],f[2][N][N],C[N][N],g[N][N]; 
void Init(int n){
	for(int i=0;i<=n;i++){
		pw[i][0]=1;
		for(int j=1;j<=n*n;j++)pw[i][j]=1ll*pw[i][j-1]*i%P;
//		cout<<pw[i][5]<<endl; 
	}
	C[0][0]=1;
	for(int i=1;i<=n;i++){
		C[i][0]=1;
		for(int j=1;j<=i;j++)
		C[i][j]=(C[i-1][j]+C[i-1][j-1])%P;
		
//		cout<<C[i][10]<<endl; 
	}
}
int main(){
	read(n),read(m),read(K),read(P);
	Init(100); int nw=0; f[0][0][0]=1;
	for(int w=1;w<=K;w++){
		nw^=1,memset(f[nw],0,sizeof(f[nw]));
		for(int j=0;j<=m;j++)
		for(int i=0;i<=n;i++)
		g[j][i]=1ll*pw[w][i*(m-j)]*pw[K-w+1][i*j]%P;
		for(int i=0;i<=n;i++)
		for(int j=0;j<=m;j++)
		if(f[nw^1][i][j]){
			for(int t=0;t+i<=n;t++)
			(f[nw][i+t][j]+=1ll*f[nw^1][i][j]*C[i+t][i]%P*g[j][t]%P)%=P
			//,cout<<f[nw][i+t][j]<<endl
			;
		}
		
		nw^=1,memset(f[nw],0,sizeof(f[nw]));
		for(int i=0;i<=n;i++)
		for(int j=0;j<=m;j++)
		g[i][j]=1ll*pw[w][j*(n-i)]*pw[K-w+1][i*j]%P;
		for(int i=0;i<=n;i++)
		for(int j=0;j<=m;j++)
		if(f[nw^1][i][j]){
			for(int t=0;t+j<=m;t++)
			(f[nw][i][j+t]+=1ll*f[nw^1][i][j]*C[j+t][j]%P*g[i][t]%P)%=P;
//			cout<<i<<" "<<j+t<<" "<<f[nw][i][j+t]<<endl;
		}
		
		nw^=1,memset(f[nw],0,sizeof(f[nw]));
		for(int j=0;j<=m;j++)
		for(int i=0;i<=n;i++)
		g[j][i]=1ll*pw[w][i*(m-j)]*pw[K-w][i*j]%P;
		for(int i=0;i<=n;i++)
		for(int j=0;j<=m;j++)
		if(f[nw^1][i][j]){
			for(int t=0;t+i<=n;t++)
			(f[nw][i+t][j]+=(t&1?-1ll:1ll)*f[nw^1][i][j]*C[i+t][i]%P*g[j][t]%P)%=P;;
//			cout<<i+t<<" "<<j<<" "<<f[nw][i+t][j]<<" "<<f[nw^1][i+t][j]<<endl;
		}
		
		nw^=1,memset(f[nw],0,sizeof(f[nw]));
		for(int i=0;i<=n;i++)
		for(int j=0;j<=m;j++)
		g[i][j]=1ll*pw[w][j*(n-i)]*pw[K-w][i*j]%P;
		for(int i=0;i<=n;i++)
		for(int j=0;j<=m;j++)
		if(f[nw^1][i][j]){
			for(int t=0;t+j<=m;t++)
			(f[nw][i][j+t]+=(t&1?-1ll:1ll)*f[nw^1][i][j]*C[j+t][j]%P*g[i][t]%P)%=P;
//			cout<<i<<" "<<j+t<<" "<<f[nw][i][j+t]<<endl;
		}
	}
	printf("%d\n",(f[nw][n][m]+P)%P);
}