#include <bits/stdc++.h>
using namespace std;
using ll = long long;
ll modinv(ll a, ll m) {
  assert(m > 0);
  if (m == 1) return 0;
  a %= m;
  if (a < 0) a += m;
  assert(a != 0);
  if (a == 1) return 1;
  return m - modinv(m, a) * m / a;
}
int MOD;
struct modnum {
 private:
  int v;

 public:
  modnum() : v(0) {}
  modnum(ll v_) : v(int(v_ % MOD)) {
    if (v < 0) v += MOD;
  }
  explicit operator int() const { return v; }
  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 operator~() const {
    modnum res;
    res.v = modinv(v, MOD);
    return res;
  }
  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); }
  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;
  }
};
using num = modnum;
num powmod(num x, int a) {
  if (a == 0) return 1;
  if (a & 1) return x * powmod(x, a - 1);
  return powmod(x * x, a / 2);
}
vector<num> fact;
vector<num> ifact;
void init() {
  fact = {1};
  for (int i = 1; i < 100000; i++) fact.push_back(i * fact[i - 1]);
  for (num x : fact) ifact.push_back(1 / x);
}
num ncr(int n, int k) {
  if (k < 0 || k > n) return 0;
  return fact[n] * ifact[k] * ifact[n - k];
}
int n, m, k;
void trans(vector<vector<num> >& A, num c) {
  vector<num> pows(n * m + 1);
  pows[0] = 1;
  for (int i = 1; i < pows.size(); i++) {
    pows[i] = pows[i - 1] * c;
  }
  for (int i = 0; i <= n; i++) {
    for (int j = 0; j <= m; j++) {
      A[i][j] *= pows[i * j];
    }
  }
}
int main() {
  ios_base::sync_with_stdio(0), cin.tie(0), cout.tie(0);
  cin >> n >> m >> k;
  cin >> MOD;
  init();
  vector<vector<num> > DP(n + 1, vector<num>(m + 1, 0));
  vector<vector<num> > DP2(n + 1, vector<num>(m + 1, 0));
  vector<num> par(n * m + 1, 0);
  for (int i = 0; i <= n * m; i++) {
    if (i & 1) {
      par[i] = -1;
    } else {
      par[i] = 1;
    }
  }
  DP[0][0] = 1;
  for (int c = 1; c < k; c++) {
    for (int x = 0; x <= n; x++) {
      for (int y = 0; y <= m; y++) {
        DP[x][y] /= powmod(c, n * m - (n - x) * (m - y));
      }
    }
    trans(DP, num(1) / num(k - c));
    for (int x = 0; x <= n; x++) {
      for (int y = 0; y <= m; y++) {
        DP[x][y] *= ifact[x] * ifact[y];
      }
    }
    DP2.assign(n + 1, vector<num>(m + 1, 0));
    for (int x = 0; x <= n; x++) {
      for (int y = 0; y <= m; y++) {
        for (int a = 0; x + a <= n; a++) {
          DP2[x + a][y] += DP[x][y] * ifact[a] * par[a];
        }
      }
    }
    DP.assign(n + 1, vector<num>(m + 1, 0));
    for (int x = 0; x <= n; x++) {
      for (int y = 0; y <= m; y++) {
        for (int b = 0; y + b <= m; b++) {
          DP[x][y + b] += DP2[x][y] * ifact[b] * par[b];
        }
      }
    }
    trans(DP, num(k - c) / num(k - c + 1));
    DP2.assign(n + 1, vector<num>(m + 1, 0));
    for (int x = 0; x <= n; x++) {
      for (int y = 0; y <= m; y++) {
        for (int a = 0; x + a <= n; a++) {
          DP2[x + a][y] += DP[x][y] * ifact[a];
        }
      }
    }
    DP.assign(n + 1, vector<num>(m + 1, 0));
    for (int x = 0; x <= n; x++) {
      for (int y = 0; y <= m; y++) {
        for (int b = 0; y + b <= m; b++) {
          DP[x][y + b] += DP2[x][y] * ifact[b];
        }
      }
    }
    trans(DP, num(k - c + 1));
    for (int x = 0; x <= n; x++) {
      for (int y = 0; y <= m; y++) {
        DP[x][y] *= fact[x] * fact[y];
      }
    }
    for (int x = 0; x <= n; x++) {
      for (int y = 0; y <= m; y++) {
        DP[x][y] *= powmod(c, n * m - (n - x) * (m - y));
      }
    }
  }
  num ans = 0;
  for (int i = 0; i <= n; i++) {
    for (int j = 0; j <= m; j++) {
      ans += DP[i][j] * ncr(n, i) * ncr(m, j) * powmod(k, (n - i) * (m - j));
    }
  }
  cout << (int)ans << '\n';
}
