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>
using namespace std;

using ll = long long;
ll modinv(ll a, ll m) {
	assert(m > 0);
	if (m == 1) return 0;
	a %= m;
	if (a < 0) a += m;
	assert(a != 0);
	if (a == 1) return 1;
	return m - modinv(m, a) * m / a;
}

int MOD;

struct modnum {
private:
	int v;
public:

	modnum() : v(0) {}
	modnum(ll v_) : v(int(v_ % MOD)) { if (v < 0) v += MOD; }
	explicit operator int () const { return v; }
	friend bool operator == (const modnum& a, const modnum& b) { return a.v == b.v; }
	friend bool operator != (const modnum& a, const modnum& b) { return a.v != b.v; }

	modnum operator ~ () const {
		modnum res;
		res.v = modinv(v, MOD);
		return res;
	}

	modnum& operator += (const modnum& o) {
		v += o.v;
		if (v >= MOD) v -= MOD;
		return *this;
	}
	modnum& operator -= (const modnum& o) {
		v -= o.v;
		if (v < 0) v += MOD;
		return *this;
	}
	modnum& operator *= (const modnum& o) {
		v = int(ll(v) * ll(o.v) % MOD);
		return *this;
	}
	modnum& operator /= (const modnum& o) {
		return *this *= (~o);
	}

	friend modnum operator + (const modnum& a, const modnum& b) { return modnum(a) += b; }
	friend modnum operator - (const modnum& a, const modnum& b) { return modnum(a) -= b; }
	friend modnum operator * (const modnum& a, const modnum& b) { return modnum(a) *= b; }
	friend modnum operator / (const modnum& a, const modnum& b) { return modnum(a) /= b; }
};

using num = modnum;

num powmod(num x, int a){
	if(a == 0) return 1;
	if(a & 1) return x * powmod(x, a-1);
	return powmod(x * x, a / 2);
}

vector<num> fact;
vector<num> ifact;
vector<num> ipfact;

void init(){
	fact = {1};
	for(int i = 1; i < 11000; i++) fact.push_back(i * fact[i-1]);
	for(num x : fact) ifact.push_back(1 / x);	
	num f = 1;
	for(num x : ifact){
		ipfact.push_back(x * f);
		f *= -1;
	}
}

num ncr(int n, int k){
	if(k < 0 || k > n) return 0;
	return fact[n] * ifact[k] * ifact[n-k];
}

int n, m, k;
void trans(vector<vector<num> >& A, num c){
	vector<num> pows(n*m+1);
	pows[0] = 1;
	for(int i = 1; i < pows.size(); i++){
		pows[i] = pows[i-1] * c;
	}
	for(int i = 0; i <= n; i++){
		for(int j = 0; j <= m; j++){
			A[i][j] *= pows[i*j];
		}
	}
}

int main(){
	ios_base::sync_with_stdio(0), cin.tie(0), cout.tie(0);
	cin >> n >> m >> k;
	cin >> MOD;
	init();
	vector<vector<num> > DP(n+1, vector<num>(m+1, 0));
	vector<vector<num> > DP2(n+1, vector<num>(m+1, 0));

	vector<num> par(n*m + 1, 0);
	for(int i = 0; i <= n*m; i++){
		if(i & 1){
			par[i] = -1;
		} else {
			par[i] = 1;
		}
	}

	DP[0][0] = 1;

	for(int c = 1; c < k; c++){
		for(int x = 0; x <= n; x++){
			for(int y = 0; y <= m; y++){
				DP[x][y] /= powmod(c, n * m - (n - x) * (m - y));
			}
		}
		trans(DP, num(1) / num(k - c + 1));
		for(int x = 0; x <= n; x++){
			for(int y = 0; y <= m; y++){
				DP[x][y] *= ifact[x] * ifact[y];
			}
		}

		DP2.assign(n+1, vector<num>(m+1, 0));
		for(int x = 0; x <= n; x++){
			for(int y = 0; y <= m; y++){
				for(int a = 0; x + a <= n; a++){
					DP2[x + a][y] += DP[x][y] * ifact[a];
				}
			}
		}
		DP.assign(n+1, vector<num>(m+1, 0));
		for(int x = 0; x <= n; x++){
			for(int y = 0; y <= m; y++){
				for(int b = 0; y + b <= m; b++){
					DP[x][y + b] += DP2[x][y] * ifact[b];
				}
			}
		}

		trans(DP, num(k - c + 1) / num(k - c));

		DP2.assign(n+1, vector<num>(m+1, 0));
		for(int x = 0; x <= n; x++){
			for(int y = 0; y <= m; y++){
				for(int a = 0; x + a <= n; a++){
					DP2[x + a][y] += DP[x][y] * ipfact[a];
				}
			}
		}
		DP.assign(n+1, vector<num>(m+1, 0));
		for(int x = 0; x <= n; x++){
			for(int y = 0; y <= m; y++){
				for(int b = 0; y + b <= m; b++){
					DP[x][y + b] += DP2[x][y] * ipfact[b];
				}
			}
		}

		trans(DP, num(k - c));
		for(int x = 0; x <= n; x++){
			for(int y = 0; y <= m; y++){
				DP[x][y] *= fact[x] * fact[y];
			}
		}
		for(int x = 0; x <= n; x++){
			for(int y = 0; y <= m; y++){
				DP[x][y] *= powmod(c, n * m - (n - x) * (m - y));
			}
		}
	}
	num ans = 0;
	for(int i = 0; i <= n; i++){
		for(int j = 0; j <= m; j++){
			ans += DP[i][j] * ncr(n, i) * ncr(m, j) * powmod(k, (n - i) * (m - j));
		}
	}
	cout << (int)ans << '\n';
}