#include <bits/stdc++.h>

#define LOG(FMT...) fprintf(stderr, FMT)

using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef vector<int> vi;

// mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());

template<class T>
istream &operator>>(istream &is, vector<T> &v) {
  for (T &x : v)
    is >> x;
  return is;
}

ostream &operator<<(ostream &os, const pair<char, int> &unit) {
  return os << unit.first << "^" << unit.second;
}

template<class T>
ostream &operator<<(ostream &os, const vector<T> &v) {
  if (!v.empty()) {
    os << v.front();
    for (int i = 1; i < v.size(); ++i)
      os << ' ' << v[i];
  }
  return os;
}

const int N = 110;

int P;
int fac[N], ifac[N], pw[N][N * N];

int mpow(int x, int k) {
  int ret = 1;
  while (k) {
    if (k & 1)
      ret = ret * (ll) x % P;
    k >>= 1;
    x = x * (ll) x % P;
  }
  return ret;
}

int norm(int x) { return x >= P ? x - P : x; }

void add(int &x, int y) {
  if ((x += y) >= P)
    x -= P;
}

void sub(int &x, int y) {
  if ((x -= y) < 0)
    x += P;
}

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

int inv(int a) {
  int x, y;
  exGcd(a, P, x, y);
  return norm(x + P);
}

void prepare(int n, int k) {
  fac[0] = 1;
  for (int i = 1; i <= n; ++i) fac[i] = fac[i - 1] * (ll) i % P;
  ifac[1] = 1;
  for (int i = 2; i <= n; ++i)
    ifac[i] = -(P / i) * (ll) ifac[P % i] % P + P;
  ifac[0] = 1;
  for (int i = 1; i <= n; ++i)
    ifac[i] = ifac[i - 1] * (ll) ifac[i] % P;
  for (int i = 0; i <= k; ++i) {
    pw[i][0] = 1;
    for (int j = 1; j <= n * n; ++j)
      pw[i][j] = pw[i][j - 1] * (ll) i % P;
  }
}

int f[N][N], g[N][N];

int main() {
#ifdef LBT
  freopen("test.in", "r", stdin);
  int nol_cl = clock();
#endif
  ios::sync_with_stdio(false);
  cin.tie(nullptr);

  int n, m, k;
  cin >> n >> m >> k >> P;
  prepare(max(n, m), k);
  f[0][0] = 1;
  for (int i = 1; i <= k; ++i) {
    for (int a = 0; a <= n; ++a)
      for (int b = 0; b <= m; ++b) {
        int v = f[a][b] * (ll) mpow(i, P - 1 - a * m - b * n + a * b) % P * fac[n - a] % P * fac[m - b] % P;
        for (int c = a; c <= n; ++c)
          for (int d = b; d <= m; ++d)
            for (int j = 0; j <= c - a; ++j)
              for (int l = 0; l <= d - b; ++l)
                add(g[c][d],
                    mpow(P - 1, j + l) * (ll) ifac[j] % P * ifac[l] % P * ifac[c - a - j] % P * ifac[d - b - l] % P
                        * mpow(k - i + 1, (c - j) * (d - l) - a * b) % P * mpow(k - i, c * l + d * j - j * l) % P * v
                        % P);
      }
    for (int c = 0; c <= n; ++c)
      for (int d = 0; d <= m; ++d)
        g[c][d] = g[c][d] * (ll) mpow(i, c * m + d * n - c * d) % P * ifac[n - c] % P * ifac[m - d] % P;
    memcpy(f, g, sizeof(f));
    for (int c = 0; c <= n; ++c)
      cerr << vector<int>(g[c], g[c] + m + 1) << '\n';
    cerr << "-----------------\n";
    memset(g, 0, sizeof(g));
  }
  cout << f[n][m];

#ifdef LBT
  LOG("Time: %dms\n", int ((clock()
      -nol_cl) / (double)CLOCKS_PER_SEC * 1000));
#endif
  return 0;
}