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 <bits/stdc++.h>
#define clr(x) memset(x,0,sizeof x)
#define For(i,a,b) for (int i=(a);i<=(b);i++)
#define Fod(i,b,a) for (int i=(b);i>=(a);i--)
#define pb(x) push_back(x)
#define mp(x,y) make_pair(x,y)
#define fi first
#define se second
#define outval(x) cerr<<#x" = "<<x<<endl
#define outtag(x) cerr<<"-----------------"#x"-----------------\n"
#define outarr(a,L,R) cerr<<#a"["<<L<<".."<<R<<"] = ";\
                    For(_x,L,R) cerr<<a[_x]<<" ";cerr<<endl;
using namespace std;
typedef long long LL;
typedef unsigned long long ULL;
typedef pair <int,int> pii;
LL read(){
    LL x=0,f=0;
    char ch=getchar();
    while (!isdigit(ch))
        f=ch=='-',ch=getchar();
    while (isdigit(ch))
        x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
    return f?-x:x;
}
const int N=105;
int n,m,k,mod;
int Pow(int x,int y){
	int ans=1;
	for (;y;y>>=1,x=(LL)x*x%mod)
		if (y&1)
			ans=(LL)ans*x%mod;
	return ans;
}
void Add(int &x,int y){
	if ((x+=y)>=mod)
		x-=mod;
}
void Del(int &x,int y){
	if ((x-=y)<0)
		x+=mod;
}
int Add(int x){
	return x>=mod?x-mod:x;
}
int Del(int x){
	return x<0?x+mod:x;
}
int Fac[N],Inv[N];
int pw[N][N*N];
void prework(){
	int n=N-1;
	Fac[0]=1;
	For(i,1,n)
		Fac[i]=(LL)Fac[i-1]*i%mod;
	Inv[n]=Pow(Fac[n],mod-2);
	Fod(i,n,1)
		Inv[i-1]=(LL)Inv[i]*i%mod;
	For(i,0,n){
		pw[i][0]=1;
		For(j,1,n*n)
			pw[i][j]=(LL)pw[i][j-1]*i%mod;
	}
}
int C(int n,int m){
	if (m<0||m>n)
		return 0;
	return (LL)Fac[n]*Inv[m]%mod*Inv[n-m]%mod;
}
int dp[4][N][N][N];
int f[4][N][N];
int main(){
	n=read(),m=read(),k=read(),mod=read();
	prework();
	if (n<m)
		swap(n,m);
	clr(dp);
	dp[0][1][0][0]=1;
	For(i,1,k){
		For(b,0,m)
			For(j,0,n){
				f[0][b][j]=(LL)Inv[j]*pw[k-i+1][b*j]%mod*pw[i][(m-b)*j]%mod;
				f[2][b][j]=(LL)Inv[j]*pw[k-i][b*j]%mod*pw[i][(m-b)*j]%mod;
			}
		For(a,0,n)
			For(j,0,m){
				f[1][a][j]=(LL)Inv[j]*pw[k-i+1][a*j]%mod*pw[i][(n-a)*j]%mod;
				f[3][a][j]=(LL)Inv[j]*pw[k-i][a*j]%mod*pw[i][(n-a)*j]%mod;
			}
		For(a,0,n)
			For(b,0,m)
				For(t,0,3){
					int v=dp[t][i][a][b];
					if (!v)
						continue;
					int T=(t+1)%4;
					if (t==0){
						For(j,0,n-a)
							Add(dp[T][i][a+j][b],(LL)v*f[0][b][j]%mod);
					}
					else if (t==1){
						For(j,0,m-b)
							Add(dp[T][i][a][b+j],(LL)v*f[1][a][j]%mod);
					}
					else if (t==2){
						For(j,0,n-a){
							int tmp=(LL)v*f[2][b][j]%mod;
							if (j&1)
								Del(dp[T][i][a+j][b],tmp);
							else
								Add(dp[T][i][a+j][b],tmp);
						}
					}
					else if (t==3){
						For(j,0,m-b){
							int tmp=(LL)v*f[3][a][j]%mod;
							if (j&1)
								Del(dp[T][i+1][a][b+j],tmp);
							else
								Add(dp[T][i+1][a][b+j],tmp);
						}
					}
				}
	}
	int ans=(LL)dp[0][k+1][n][m]*Fac[n]%mod*Fac[m]%mod;
	cout<<ans<<endl;
	return 0;
}