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:
// I need some holiday
// But I don't wish that everyday is holiday
// Because I didn't bear the pain
// Because it wasn't a true freedom
//         ——『イヤホンと蝉時雨』·Orangestar

#include <bits/stdc++.h>

using namespace std;

typedef long long LL;
typedef unsigned int UI;
typedef pair<int, int> pii;

#define DEBUG(...) fprintf(stderr, __VA_ARGS__)
#define MP make_pair
#define fst first
#define snd second
#define SZ(u) ((int) (u).size())
#define ALL(u) (u).begin(), (u).end()

inline void proc_status()
{
	ifstream t("/proc/self/status");
	cerr << string(istreambuf_iterator<char>(t), istreambuf_iterator<char>()) << endl;
}

template<typename T> inline int chkmin(T &a, const T &b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline int chkmax(T &a, const T &b) { return b > a ? a = b, 1 : 0; }
template<typename T> inline T sqr(const T &val) { return val * val; }

namespace fastIO
{
	const int MAX_BUFFER_SIZE = 1 << 16;

	char buffer[MAX_BUFFER_SIZE], *cur = buffer, *ed = buffer;

	inline char getc() { return cur == ed && (ed = (cur = buffer) + fread(buffer, 1, MAX_BUFFER_SIZE, stdin), cur == ed) ? EOF : *cur++; }
}
using fastIO::getc;

template<typename T> inline T read()
{
	register T sum(0), fg(1);
	register char ch(getc());
	for(; !isdigit(ch); ch = getc()) if(ch == '-') fg = -1;
	for(;  isdigit(ch); ch = getc()) sum = sum * 10 - '0' + ch;
	return sum * fg;
}

inline char *read_str(char *s)
{
	register char ch(getc());
	while(!isgraph(ch)) ch = getc();
	for(; isgraph(ch); ch = getc()) *s++ = ch;
	*s = '\0';
	return s;
}

// END tpl.

int MOD;

inline int fpm(int x, int y)
{
	int res = 1;
	for(; y; y >>= 1, x = (LL) x * x % MOD) if(y & 1) res = (LL) res * x % MOD;
	return res;
}

inline int inv(int x) { assert(x); return fpm(x, MOD - 2); }

namespace MATH
{
	const int n = 100;

	int fac[n + 5], ifac[n + 5];

	inline void init()
	{
		fac[0] = 1;
		for(int i = 1; i <= n; ++i) fac[i] = (LL) fac[i - 1] * i % MOD;
		ifac[n] = inv(fac[n]);
		for(int i = n - 1; i >= 0; --i) ifac[i] = (LL) ifac[i + 1] * (i + 1) % MOD;
	}
}
using MATH::fac;
using MATH::ifac;

const int MAXN = 100;

int n, m, K;

inline void input()
{
	n = read<int>(), m = read<int>(), K = read<int>(), MOD = read<int>();
}

inline void solve()
{
	static int pw[MAXN + 5][MAXN * MAXN + 5];

	for(int i = 0; i <= K; ++i)
	{
		pw[i][0] = 1;
		for(int j = 1; j <= n * m; ++j) pw[i][j] = (LL) pw[i][j - 1] * i % MOD;
	}

	static int f[MAXN + 5][MAXN + 5], g[MAXN + 5][MAXN + 5];

	f[0][0] = 1;
	for(int D = 1; D <= K; ++D)
	{
		static int coef[MAXN + 5][MAXN + 5];

		swap(f, g), memset(f, 0, sizeof f);
		for(int j = 0; j <= n; ++j) for(int k = 0; k <= m; ++k)
			coef[j][k] = (LL) (j % 2 == 0 ? +1 : -1) * ifac[j] * pw[D - 1][j * (m - k)] % MOD * pw[K - D + 1][j * k] % MOD;

		for(int i = 0; i <= n; ++i) for(int j = 0; j <= m; ++j)
			for(int k = 0; k <= i; ++k) (f[i][j] += (LL) g[i - k][j] * coef[k][j] % MOD) %= MOD;

		swap(f, g), memset(f, 0, sizeof f);
		for(int j = 0; j <= m; ++j) for(int k = 0; k <= n; ++k)
			coef[j][k] = (LL) (j % 2 == 0 ? +1 : -1) * ifac[j] * pw[D - 1][j * (n - k)] % MOD * pw[K - D + 1][j * k] % MOD;

		for(int i = 0; i <= n; ++i) for(int j = 0; j <= m; ++j)
			for(int k = 0; k <= j; ++k) (f[i][j] += (LL) g[i][j - k] * coef[k][i] % MOD) %= MOD;

		swap(f, g), memset(f, 0, sizeof f);
		for(int j = 0; j <= n; ++j) for(int k = 0; k <= m; ++k)
			coef[j][k] = (LL) ifac[j] * pw[D][j * (m - k)] % MOD * pw[K - D + 1][j * k] % MOD;

		for(int i = 0; i <= n; ++i) for(int j = 0; j <= m; ++j)
			for(int k = 0; k <= i; ++k) (f[i][j] += (LL) g[i - k][j] * coef[k][j] % MOD) %= MOD;

		swap(f, g), memset(f, 0, sizeof f);
		for(int j = 0; j <= m; ++j) for(int k = 0; k <= n; ++k)
			coef[j][k] = (LL) ifac[j] * pw[D][j * (n - k)] % MOD * pw[K - D + 1][j * k] % MOD;

		for(int i = 0; i <= n; ++i) for(int j = 0; j <= m; ++j)
			for(int k = 0; k <= j; ++k) (f[i][j] += (LL) g[i][j - k] * coef[k][i] % MOD) %= MOD;
	}
	int ans = (LL) f[n][m] * fac[n] % MOD * fac[m] % MOD;
	printf("%d\n", (ans + MOD) % MOD);
}

int main()
{
#ifdef K_ON
	freopen("F.in", "r", stdin);
	freopen("F.out", "w", stdout);
#endif

	input();
	MATH::init();
	solve();

	return 0;
}