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;
typedef bool boolean;

#define ll long long

void exgcd(int a, int b, int& x, int& y) {
	if (!b) {
		x = 1, y = 0;
	} else {
		exgcd(b, a % b, y, x);
		y -= (a / b) * x;
	}
}

int inv(int a, int n) {
	int x, y;
	exgcd(a, n, x, y);
	return (x < 0) ? (x + n) : (x);
}

int Mod;

typedef class Z {
	public:
		int v;

		Z() : v(0) {	}
		Z(int x) : v(x){	}
		Z(ll x) : v(x % Mod) {	}

		friend Z operator + (const Z& a, const Z& b) {
			int x;
			return Z(((x = a.v + b.v) >= Mod) ? (x - Mod) : (x));
		}
		friend Z operator - (const Z& a, const Z& b) {
			int x;
			return Z(((x = a.v - b.v) < 0) ? (x + Mod) : (x));
		}
		friend Z operator * (const Z& a, const Z& b) {
			return Z(a.v * 1ll * b.v);
		}
		friend Z operator ~(const Z& a) {
			return inv(a.v, Mod);
		}
		friend Z operator - (const Z& a) {
			return Z(0) - a;
		}
		Z& operator += (Z b) {
			return *this = *this + b;
		}
		Z& operator -= (Z b) {
			return *this = *this - b;
		}
		Z& operator *= (Z b) {
			return *this = *this * b;
		}
		friend boolean operator == (const Z& a, const Z& b) {
			return a.v == b.v;
		} 
} Z;

Z qpow(Z a, int p) {
	Z rt = Z(1), pa = a;
	for ( ; p; p >>= 1, pa = pa * pa) {
		if (p & 1) {
			rt = rt * pa;
		}
	}
	return rt;
}

typedef Z Zi;

const int N = 102;

int n, m, K;
Zi comb[N][N];
Zi f[2][N][N];
Zi coef1[N * N], ceof2[N];

int main() {
	scanf("%d%d%d%d", &n, &m, &K, &Mod);
	comb[0][0] = 1;
	int qaq = max(n, m);
	for (int i = 1; i <= qaq; i++) {
		comb[i][0] = comb[i][i] = 1;
		for (int j = 1; j < i; j++) {
			comb[i][j] = comb[i - 1][j - 1] + comb[i - 1][j];
		}
	}
	int cur = 0;
	f[cur][0][0] = 1;
	for (int k = 1; k <= K; k++) {
		memset(f[cur ^= 1], 0, sizeof(f[0]));
		coef1[0] = 1, coef1[1] = K - k + 1;
		for (int i = 2, _ = n * m; i <= _; i++)
			coef1[i] = coef1[i - 1] * coef1[1];
		for (int i = 0; i <= n; i++) {
			for (int j = 0; j <= m; j++) {
				Zi v = f[cur ^ 1][i][j];
				if (!v.v)
					continue;
				Zi y = qpow(k, m - j) - qpow(k - 1, m - j), pw = 1;
				for (int x = 0; i + x <= n; x++, pw *= y)
					f[cur][i + x][j] += v * comb[n - i][x] * coef1[x * j] * pw;
			}
		}
		memset(f[cur ^= 1], 0, sizeof(f[0]));
		coef1[0] = 1, coef1[1] = k;
		for (int i = 2, _ = n * m; i <= _; i++)
			coef1[i] = coef1[i - 1] * coef1[1];
		for (int i = 0; i <= n; i++) {
			Zi y = qpow(K - k + 1, i) - qpow(K - k, i);
			for (int j = 0; j <= m; j++) {
				Zi v = f[cur ^ 1][i][j];
				if (!v.v)
					continue;
				Zi pw = 1;
				for (int x = 0; j + x <= m; x++, pw *= y)
					f[cur][i][j + x] += v * comb[m - j][x] * pw * coef1[(n - i) * x];
			}
		}
	}
	Zi ans = f[cur][n][m];
	printf("%d\n", ans.v);
	return 0;
}