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;

template <int& MOD_> struct modnum {
	static constexpr int& MOD = MOD_;
	//static_assert(MOD_ > 0, "MOD must be positive");

private:
	using ll = long long;

	int v;

	static int minv(int a, int m) {
		a %= m;
		assert(a);
		return a == 1 ? 1 : int(m - ll(minv(m, a)) * ll(m) / a);
	}

public:

	modnum() : v(0) {}
	modnum(ll v_) : v(int(v_ % MOD)) { if (v < 0) v += MOD; }
	explicit operator int() const { return v; }
	friend std::ostream& operator << (std::ostream& out, const modnum& n) { return out << int(n); }
	friend std::istream& operator >> (std::istream& in, modnum& n) { ll v_; in >> v_; n = modnum(v_); return in; }

	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 inv() const {
		modnum res;
		res.v = minv(v, MOD);
		return res;
	}
	modnum neg() const {
		modnum res;
		res.v = v ? MOD-v : 0;
		return res;
	}

	modnum operator- () const {
		return neg();
	}
	modnum operator+ () const {
		return modnum(*this);
	}

	modnum& operator ++ () {
		v ++;
		if (v == MOD) v = 0;
		return *this;
	}
	modnum& operator -- () {
		if (v == 0) v = MOD;
		v --;
		return *this;
	}
	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.inv();
	}

	friend modnum operator ++ (modnum& a, int) { modnum r = a; ++a; return r; }
	friend modnum operator -- (modnum& a, int) { modnum r = a; --a; return r; }
	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; }
};

int MOD;
using num = modnum<MOD>;

const int MAXN = 120;
const int MAXM = 120;
int N, M, K;
num dp[MAXN][MAXM];
num ndp[MAXN][MAXM];

const int MAXV = 1.1e3;
const int V = 1010;
num fact[MAXV];
num ifact[MAXV];
num prepow[MAXV][MAXV];

num facs[MAXV][MAXV];

int main() {
	ios_base::sync_with_stdio(0), cin.tie(0), cout.tie(0);
	cin >> N >> M >> K >> MOD;

	fact[0] = 1;
	for (int i = 1; i <= V; i++) {
		fact[i] = fact[i-1] * num(i);
	}
	ifact[V] = fact[V].inv();
	for (int i = V; i > 0; i--) {
		ifact[i-1] = ifact[i] * num(i);
	}

	for (int b = 0; b <= V; b++) {
		prepow[b][0] = 1;
		for (int e = 1; e <= V; e++) {
			prepow[b][e] = prepow[b][e-1] * num(b);
		}
	}

	for (int i = 0; i <= N; i++) {
		for (int j = 0; j <= M; j++) {
			dp[i][j] = 0;
		}
	}
	dp[0][0] = fact[N] * fact[M];

	for (int k = 0; k < 2 * K; k++) {
		for (int i = 0; i <= N; i++) {
			for (int j = 0; j <= M; j++) {
				ndp[i][j] = 0;
			}
		}

		if (k % 2 == 0) {
			// increment row
			for (int j = 0; j <= M; j++) {
				num rowFac = (prepow[k/2+1][M-j] - prepow[k/2][M-j]) * prepow[K-k/2][j];
				num powRowFac = 1;
				for (int di = 0; di <= N; di++, powRowFac *= rowFac) {
					facs[di][j] = powRowFac * ifact[di];
				}
			}
			for (int i = 0; i <= N; i++) {
				for (int ni = i; ni <= N; ni++) {
					for (int j = 0; j <= M; j++) {
						ndp[ni][j] += dp[i][j] * facs[ni-i][j];
					}
				}
			}
		} else {
			// increment col
			for (int i = 0; i <= N; i++) {
				num colFac = (prepow[K-k/2][i] - prepow[K-1-k/2][i]) * prepow[k/2+1][N-i];
				for (int j = 0; j <= M; j++) {
					num powColFac = 1;
					for (int nj = j; nj <= M; nj++, powColFac *= colFac) {
						num val = dp[i][j];
						val *= powColFac;
						val *= ifact[nj-j];
						ndp[i][nj] += val;
					}
				}
			}
		}
		for (int i = 0; i <= N; i++) {
			for (int j = 0; j <= M; j++) {
				dp[i][j] = ndp[i][j];
			}
		}
	}
	cout << dp[N][M] << '\n';

	return 0;
}