#include <algorithm>  
#include <iostream>  
#include <sstream>  
#include <string>  
#include <cstring>
#include <vector>  
#include <queue>  
#include <set>  
#include <map>  
#include <cstdio>  
#include <cstdlib>  
#include <cctype>  
#include <cmath>  
#include <list>  
#include <cassert>
#include <ctime>
#include <climits>
using namespace std;  

#define PB push_back  
#define MP make_pair  
#define SZ(v) ((int)(v).size())  
#define FOR(i,a,b) for(int i=(a);i<(b);++i)  
#define REP(i,n) FOR(i,0,n)  
#define FORE(i,a,b) for(int i=(a);i<=(b);++i)  
#define REPE(i,n) FORE(i,0,n)  
#define FORSZ(i,a,v) FOR(i,a,SZ(v))  
#define REPSZ(i,v) REP(i,SZ(v))  
typedef long long ll;
typedef unsigned long long ull;
ll gcd(ll a,ll b) { return b==0?a:gcd(b,a%b); }

const int MAXH=100;
const int MAXW=100;
const int MAXDIM=100;
const int MAXVAL=100;

int h,w,mxval,mod;
void inc(int &a,int b) { if((a+=b)>=mod) a-=mod; }

// dp[nrow][ncol] = # ways of:
// * choosing nrow rows and ncol cols in the grid
// * fixing the values of the intersection of those columns and rows, such that each column has at least one cell < curvalue
// * fixing the contribution for each cell in the union of those columns, such that each of those cells has a valid contribution < curvalue
// tdp[nrow][ncol] = similar to dp[nrow][ncol], but now the contribution can be <= curvalue
int dp[MAXH+1][MAXW+1]; 
int tdp[MAXH+1][MAXW+1];

int pw[MAXVAL+1][MAXW+1];
int choose[MAXDIM+1][MAXDIM+1];

int solve() {
	REPE(val,mxval) { pw[val][0]=1; FORE(y,1,max(h,w)) pw[val][y]=(ll)pw[val][y-1]*val%mod; }
	REPE(i,MAXDIM) { choose[i][0]=choose[i][i]=1; FOR(j,1,i) choose[i][j]=(choose[i-1][j-1]+choose[i-1][j])%mod; }

	REPE(x,h) REPE(y,w) dp[x][y]=x==0&&y==0?1:0;
	FORE(val,1,mxval) {
		REPE(x,h) REPE(y,w) tdp[x][y]=0;
		REPE(x,h) REPE(y,w) {
			int eachrowcontrib=pw[val][w-y]-pw[val-1][w-y]; if(eachrowcontrib<0) eachrowcontrib+=mod;
			int eachrowvalue=pw[mxval-val+1][y];
			int cur=dp[x][y],mlt=(ll)eachrowcontrib*eachrowvalue%mod;
			REPE(dx,h-x) {
				int now=(ll)choose[h-x][dx]*cur%mod;
				inc(tdp[x+dx][y],now);
				cur=(ll)cur*mlt%mod;
			}
		}
		//printf("tdp after %d:\n",val); REPE(x,h) { REPE(y,w) { if(y!=0) printf(" "); printf("%d",tdp[x][y]); } puts(""); }
		REPE(x,h) REPE(y,w) dp[x][y]=0;
		REPE(x,h) REPE(y,w) {
			int eachcolcontrib=pw[val][h-x];
			int eachcolvalue=pw[mxval-val+1][x]-pw[mxval-val][x]; if(eachcolvalue<0) eachcolvalue+=mod;
			int cur=tdp[x][y],mlt=(ll)eachcolcontrib*eachcolvalue%mod;
			REPE(dy,w-y) {
				int now=(ll)choose[w-y][dy]*cur%mod;
				inc(dp[x][y+dy],now);
				cur=(ll)cur*mlt%mod;
			}
		}
		//printf("dp after %d:\n",val); REPE(x,h) { REPE(y,w) { if(y!=0) printf(" "); printf("%d",dp[x][y]); } puts(""); }
	}
	return dp[h][w];
}

void run() {
	scanf("%d%d%d%d",&h,&w,&mxval,&mod);
	printf("%d\n",solve());
}

int main() {
	run();
	return 0;
}
