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>
#pragma GCC optimize ("O2,unroll-loops")
//#pragma GCC optimize("no-stack-protector,fast-math")
//#pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,tune=native")

using namespace std;
typedef long long ll;
typedef long double ld;
typedef pair<int, int> pii;
typedef pair<pii, int> piii;
typedef pair<ll, ll> pll;
#define debug(x) cerr<<#x<<'='<<(x)<<endl;
#define debugp(x) cerr<<#x<<"= {"<<(x.first)<<", "<<(x.second)<<"}"<<endl;
#define debug2(x, y) cerr<<"{"<<#x<<", "<<#y<<"} = {"<<(x)<<", "<<(y)<<"}"<<endl;
#define debugv(v) {cerr<<#v<<" : ";for (auto x:v) cerr<<x<<' ';cerr<<endl;}
#define all(x) x.begin(), x.end()
#define pb push_back
#define kill(x) return cout<<x<<'\n', 0;

const int N=105;

ll n, m, k, u, v, x, y, t, a, b, ans, mod;
ll C[N][N];
ll tav[N][N];
ll dp0[N][N][N]; // b=a
ll dp1[N][N][N]; // b=a+1

int main(){
	ios_base::sync_with_stdio(false);cin.tie(0);cout.tie(0);
	//freopen("input.txt", "r", stdin);
	//freopen("output.txt", "w", stdout);
	cin>>n>>m>>k>>mod;
	for (int i=0; i<N; i++){
		C[i][0]=C[i][i]=1;
		for (int j=1; j<i; j++) C[i][j]=(C[i-1][j-1] + C[i-1][j])%mod;
	}
	for (int x=0; x<N; x++){
		tav[x][0]=1;
		for (int y=1; y<N; y++) tav[x][y]=tav[x][y-1]*x%mod; // :|
	}
	
	dp0[0][0][0]=1;
	for (int a=0; a<k; a++){
		// b=a
		for (int x=0; x<=n; x++) for (int y=0; y<=m; y++) if (dp0[a][x][y]){
			ll val=tav[k-a][y]*(tav[a+1][m-y]-tav[a][m-y])%mod, mul=1;
			for (int xx=x; xx<=n; xx++){
				dp1[a][xx][y]=(dp1[a][xx][y] + dp0[a][x][y]*mul%mod*C[n-x][xx-x])%mod;
				mul=mul*val%mod;
			}
		}
		
		// b=a+1
		for (int x=0; x<=n; x++) for (int y=0; y<=m; y++) if (dp1[a][x][y]){
			ll val=tav[a+1][n-x]*(tav[k-a][x]-tav[k-a-1][x])%mod, mul=1;
			for (int yy=y; yy<=m; yy++){
				dp0[a+1][x][yy]=(dp0[a+1][x][yy] + dp1[a][x][y]*mul%mod*C[m-y][yy-y])%mod;
				mul=mul*val%mod;
			}
		}
	}
	ans=dp0[k][n][m];
	if (ans<0) ans+=mod;
	cout<<ans<<"\n";
	
	return 0;
}
/*
2 2 1 998244353
1 2 2 998244353
2 1 2 998244353

*/