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:

//#define USEPB_DS
#define USETR1
#define CPPELEVEN
#define GPP

/*
 * temp.cpp
 *
 *  Created on: 2012-7-18
 *      Author: BSBandme
 */
//#pragma comment(linker, "/STACK:1024000000,1024000000")

#include <iostream>
#include <fstream>
#include <string.h>
#include <cstdio>
#include <algorithm>
#include <string>
#include <vector>
#include <queue>
#include <cassert>
#include <list>
#include <iomanip>
#include <math.h>
#include <deque>
#include <utility>
#include <map>
#include <set>
#include <bitset>
#include <numeric>
#include <climits>
#include <cctype>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <functional>
#include <sstream>
#include <tuple>

using namespace std;

#ifndef CPPELEVEN
#ifdef USETR1
#include <tr1/unordered_map>
#include <tr1/unordered_set>
using namespace tr1;
#endif
#else
#include <unordered_map>
#include <unordered_set>
#endif

#ifdef USEPB_DS
#include <ext/pb_ds/priority_queue.hpp>
#include <ext/pb_ds/assoc_container.hpp>
using namespace __gnu_pbds;
// binomial_heap_tag, rc_binomial_heap_tag, thin_heap_tag, binary_heap_tag
typedef __gnu_pbds::priority_queue<int, greater<int>, pairing_heap_tag> pq_type;
// splay_tree_tag, ov_tree_tag
typedef tree <int, null_type, less <int>, rb_tree_tag, tree_order_statistics_node_update> tree_type;
#endif

#define mpr make_pair
typedef unsigned int ui;
typedef unsigned long long ull;
typedef long long ll;
typedef pair <int, int> pii;
typedef pair <ll, ll> pll;
typedef pair <double, double> pdd;
typedef vector <int> vi;
typedef vector <ll> vll;
typedef vector <double> vd;
typedef vector <string> vs;
typedef map <string, int> mpsi;
typedef map <double, int> mpdi;
typedef map <int, int> mpii;

const double pi = acos(0.0) * 2.0;
const long double eps = 1e-10;
const int step[8][2] = {{-1, 0}, {0, 1}, {1, 0}, {0, -1}, {-1, 1}, {1, 1}, {1, -1}, {-1, -1}};

template <class T> inline T abs1(T a) {return a < 0 ? -a : a;}

#ifndef CPPELEVEN
template <class T> inline T max1(T a, T b) { return b < a ? a : b; }
template <class T> inline T max1(T a, T b, T c) { return max1(max1(a, b), c); }
template <class T> inline T max1(T a, T b, T c, T d) { return max1(max1(a, b, c), d); }
template <class T> inline T max1(T a, T b, T c, T d, T e) { return max1(max1(a, b, c, d), e); }
template <class T> inline T min1(T a, T b) { return a < b ? a : b; }
template <class T> inline T min1(T a, T b, T c) { return min1(min1(a, b), c); }
template <class T> inline T min1(T a, T b, T c, T d) { return min1(min1(a, b, c), d); }
template <class T> inline T min1(T a, T b, T c, T d, T e) { return min1(min1(a, b, c, d), e); }
#else
template <typename t, typename t1>
t min1(t a, t1 b) { return a < b ? a : b; }
template <typename t, typename... arg>
t min1(t a, arg... arr) { return min1(a, min1(arr...)); }
template <typename t, typename t1>
t max1(t a, t1 b) { return a > b ? a : b; }
template <typename t, typename... arg>
t max1(t a, arg... arr) { return max1(a, max1(arr...)); }
#endif

inline int jud(double a, double b){
	if(abs(a) < eps && abs(b) < eps) return 0;
	else if(abs1(a - b) / max(abs1(a), abs1(b)) < eps) return 0;
	if(a < b) return -1;
	return 1;
}
template <typename t> inline int jud(t a, t b){
	if(a < b) return -1;
	if(a == b) return 0;
	return 1;
}

// f_lb == 1代表返回相同的一串的左边界，f_small == 1代表返回如果没有寻找的值返回小的数
template <typename it, typename t1>
inline int find(t1 val, it a, int na, bool f_small = 1, bool f_lb = 1){
	if(na == 0) return 0;
	int be = 0, en = na - 1;
	if(*a <= *(a + na - 1)){
		if(f_lb == 0) while(be < en){
			int mid = (be + en + 1) / 2;
			if(jud(*(a + mid), val) != 1) be = mid;
			else en = mid - 1;
		}else while(be < en){
			int mid = (be + en) / 2;
			if(jud(*(a + mid), val) != -1) en = mid;
			else be = mid + 1;
		}
		if(f_small && jud(*(a + be), val) == 1) be--;
		if(!f_small && jud(*(a + be), val) == -1) be++;
	} else {
		if(f_lb) while(be < en){
			int mid = (be + en + 1) / 2;
			if(jud(*(a + mid), val) != -1) be = mid;
			else en = mid - 1;
		}else while(be < en){
			int mid = (be + en) / 2;
			if(jud(*(a + mid), val) != 1) en = mid;
			else be = mid + 1;
		}
		if(!f_small && jud(*(a + be), val) == -1) be--;
		if(f_small && jud(*(a + be), val) == 1) be++;
	}
	return be;
}

template <class T> inline T lowb(T num) {return num & (-num); }
#ifdef GPP
inline int bitnum(ui nValue) { return __builtin_popcount(nValue); }
inline int bitnum(int nValue) { return __builtin_popcount(nValue); }
inline int bitnum(ull nValue) { return __builtin_popcount(nValue) + __builtin_popcount(nValue >> 32); }
inline int bitnum(ll nValue) { return __builtin_popcount(nValue) + __builtin_popcount(nValue >> 32); }
inline int bitmaxl(ui a) { if(a == 0) return 0; return 32 - __builtin_clz(a); }
inline int bitmaxl(int a) { if(a == 0) return 0; return 32 - __builtin_clz(a); }
inline int bitmaxl(ull a) { int temp = a >> 32; if(temp) return 32 - __builtin_clz(temp) + 32; return bitmaxl(int(a)); }
inline int bitmaxl(ll a) { int temp = a >> 32; if(temp) return 32 - __builtin_clz(temp) + 32; return bitmaxl(int(a)); }
#else
#endif

long long pow(long long n, long long m, long long mod = 0){
	if(m < 0) return 0;
	long long ans = 1;
	long long k = n;
	while(m){
		if(m & 1) {
			ans *= k;
			if(mod) ans %= mod;
		}
		k *= k;
		if(mod) k %= mod;
		m >>= 1;
	}
	return ans;
}

#define  MOD 1000000007
template <class t1, class t2>
inline void add(t1 &a, t2 b, int mod = -1) {
	if(mod == -1) mod = MOD;
	a += b;
	while(a >= mod) a -= mod;
	while(a < 0) a += mod;
}
template <class t>
void output1(t arr) {
	for(int i = 0; i < (int)arr.size(); i++)
		cerr << arr[i] << ' ';
	cerr << endl;
}
template <class t>
void output2(t arr) {
	for(int i = 0; i < (int)arr.size(); i++)
		output1(arr[i]);
}

//....................密..........封..........线..........下..........禁..........止..........hack...............................................

const int maxn = 110;
ll dp[maxn][maxn], dp1[maxn][maxn];
ll c[maxn][maxn], fac[maxn], invfac[maxn], p[maxn][maxn];
int n, m, pk, mod;

int main() {

//............................不要再忘了检查maxn大小了！！！！BSBandme你个SB！！！！...................................................

	ios_base::sync_with_stdio(0);
	#ifdef DEBUG //......................................................................................................
	freopen("input.txt", "r", stdin);
	#endif //...........................................................................................................

	scanf("%d%d%d%d%d", &n, &m, &pk, &mod);

	for (int i = 0; i < maxn; i++) {
		c[i][0] = 1;
		for (int j = 1; j <= i; j++) {
			c[i][j] = c[i - 1][j] + c[i - 1][j - 1];
			c[i][j] %= mod;
		}
	}
	fac[0] = 1; invfac[0] = 1;
	for (int i = 1; i < maxn; i++) {
		fac[i] = fac[i - 1] * i % mod;
		invfac[i] = pow(fac[i], mod - 2, mod);
	}
	for (int i = 0; i < maxn; i++) {
		p[i][0] = 1;
		for (int j = 1; j <= i; j++) {
			p[i][j] = p[i][j - 1] * (i - j + 1) % mod;
		}
	}

	dp[0][0] = 1;
	for (int i = 1; i <= pk; i++) {
//		for (int j = n; j >= 0; j--) for (int k = m; k >= 0; k--) if(dp[j][k]) {
//			for (int rj = j; rj <= n; rj++) for (int rk = k; rk <= m; rk++) {
//				if (rj == 0 && rk != 0) continue;
//				if (rj != 0 && rk == 0) continue;
//				if (rj == j && rk == k) continue;
//				int all = (rj - j) * k + (rk - k) * j + (rj - j) * (rk - k);
//				int multi = (rj - j) * (m - k) + (rk - k) * (n - j) - (rj - j) * (rk - k);
//				multi = pow(i, multi, mod) * c[n - j][rj - j] % mod * c[m - k][rk - k] % mod;
//				for (int dj = 0; dj <= rj - j; dj++) for (int dk = 0; dk <= rk - k; dk++) {
//					int flag = 1;
//					if ((dj + dk) & 1) flag = -1;
//					int rd = rj * dk + rk * dj - dj * dk;
//					dp[rj][rk] += dp[j][k] * pow(pk - i + 1, all - rd, mod) % mod * pow(pk - i, rd, mod) % mod * flag * multi % mod * c[rj - j][dj] % mod * c[rk - k][dk] % mod;
//					dp[rj][rk] %= mod;
//				}
//			}
//		}
		memset(dp1, 0, sizeof(dp1));
		for (int j = n; j >= 0; j--) for (int k = m; k >= 0; k--) if (dp[j][k]) {
			ll multi = pow(pk - i + 1, k, mod) * (pow(i, m - k, mod) - pow(i - 1, m - k, mod)) % mod;
			for (ll t = j, rm = 1; t <= n; t++) {
				(dp1[t][k] += dp[j][k] * rm % mod * c[n - j][t - j]) %= mod;
				rm = rm * multi % mod;
			}
		}
		memset(dp, 0, sizeof(dp));
		for (int j = n; j >= 0; j--) for (int k = m; k >= 0; k--) if (dp1[j][k]) {
			ll multi = (pow(pk - i + 1, j, mod) - pow(pk - i, j, mod))* pow(i, n - j, mod) % mod;
			for (ll t = k, rm = 1; t <= m; t++) {
				(dp[j][t] += dp1[j][k] * rm % mod * c[m - k][t - k]) %= mod;
				rm = rm * multi % mod;
			}
		}
	}

	cout << (dp[n][m] + mod) % mod << endl;


	return 0;
}