#include <bits/stdc++.h>
using namespace std;

int pw[110][11000];
int qw[110][11000];
int f[110], g[110];
int dp[2][110][110];
int ex[11000];

int a[110][110];
int b[110][110];

int fast_exponentiation(int a, int n, int p) {
	int ret = 1;
	for (; n; n >>= 1, a = 1LL * a * a % p) if (n & 1) ret = 1LL * ret * a % p;
	return ret;
}

int main() {
	ios::sync_with_stdio(0);
	cin.tie(0); cout.tie(0);
	int n, m, k, mod; cin >> n >> m >> k >> mod;
	for (int x = 0; x <= k; x++) {
		for (int j = 0; j <= n * m; j++) {
			pw[x][j] = j == 0 ? 1 : pw[x][j - 1] * 1LL * x % mod;
		}
		int y = fast_exponentiation(x, mod - 2, mod);
		for (int j = 0; j <= n * m; j++) {
			qw[x][j] = j == 0 ? 1 : qw[x][j - 1] * 1LL * y % mod;
		}
	}
	for (int i = 0; i <= max(n, m); i++) {
		f[i] = i == 0 ? 1 : f[i - 1] * 1LL * i % mod;
		g[i] = fast_exponentiation(f[i], mod - 2, mod);
	}
	int M = 2 * m + 1, N = 2 * n + 1;
	dp[0][0][0] = 1;
	int now = 0;
	for (int i = 1; i < k; i++) {
		now ^= 1;
		int alpha = (1 + qw[k - i][1]) % mod;
		for (int j = 0; j <= n * m; j++) {
			ex[j] = j == 0 ? 1 : ex[j - 1] * 1LL * alpha % mod;
		}
		for (int r = 0; r <= n; r++) {
			for (int c = 0; c <= m; c++) {
				a[r][c] = dp[!now][r][c] * 1LL * pw[i][(n - r) * (m - c)] % mod;
			}
		}
		for (int r = 0; r <= n; r++) {
			for (int c = 0; c <= m; c++) {
				b[r][c] = 0;
				for (int x = 0; x <= r; x++) {
					b[r][c] = (b[r][c] + a[x][c] * 1LL * g[r - x]) % mod;
				}
			}
		}
		for (int r = 0; r <= n; r++) {
			for (int c = 0; c <= m; c++) {
				a[r][c] = 0;
				for (int x = 0; x <= c; x++) {
					a[r][c] = (a[r][c] + b[r][x] * 1LL * g[c - x]) % mod;
				}
				a[r][c] = a[r][c] * 1LL * ex[r * c] % mod;
			}
		}
		for (int r = 0; r <= n; r++) {
			for (int c = 0; c <= m; c++) {
				b[r][c] = 0;
				for (int x = 0; x <= r; x++) {
					if (r - x & 1) b[r][c] = (b[r][c] + a[x][c] * 1LL * (mod - g[r - x])) % mod;
					else b[r][c] = (b[r][c] + a[x][c] * 1LL * (g[r - x])) % mod;
				}
			}
		}
		for (int r = 0; r <= n; r++) {
			for (int c = 0; c <= m; c++) {
				a[r][c] = 0;
				for (int x = 0; x <= c; x++) {
					if (c - x & 1) a[r][c] = (a[r][c] + b[r][x] * 1LL * (mod - g[c - x])) % mod;
					else a[r][c] = (a[r][c] + b[r][x] * 1LL * (g[c - x])) % mod;
				}
				dp[now][r][c] = a[r][c] * 1LL * qw[i][(n - r) * (m - c)] % mod;
			}
		}
	}
	int ans = 0;
	for (int r = 0; r <= n; r++) {
		for (int c = 0; c <= m; c++) {
			ans = (ans + g[n - r] * 1LL * g[m - c] % mod * dp[now][r][c] % mod * pw[k][(n - r) * (m - c)]) % mod;
		}
	}
	ans = 1LL * ans * f[n] % mod * f[m] % mod;
	cout << ans << endl;
	return 0;
}
