#include <bits/stdc++.h>
using namespace std;
using ll = long long;
template <class t, class u>
void chmax(t& first, u second) {
  if (first < second) first = second;
}
template <class t, class u>
void chmin(t& first, u second) {
  if (second < first) first = second;
}
template <class t>
using vc = vector<t>;
template <class t>
using vvc = vc<vc<t>>;
using pi = pair<int, int>;
using vi = vc<int>;
template <class t, class u>
ostream& operator<<(ostream& os, const pair<t, u>& p) {
  return os << "{" << p.first << "," << p.second << "}";
}
template <class t>
ostream& operator<<(ostream& os, const vc<t>& v) {
  os << "{";
  for (auto e : v) os << e << ",";
  return os << "}";
}
using uint = unsigned;
using ull = unsigned long long;
template <int i, class T>
void print_tuple(ostream&, const T&) {}
template <int i, class T, class H, class... Args>
void print_tuple(ostream& os, const T& t) {
  if (i) os << ",";
  os << get<i>(t);
  print_tuple<i + 1, T, Args...>(os, t);
}
template <class... Args>
ostream& operator<<(ostream& os, const tuple<Args...>& t) {
  os << "{";
  print_tuple<0, tuple<Args...>, Args...>(os, t);
  return os << "}";
}
void print(ll x, int suc = 1) {
  cout << x;
  if (suc == 1) cout << "\n";
  if (suc == 2) cout << " ";
}
ll read() {
  ll i;
  cin >> i;
  return i;
}
vi readvi(int n, int off = 0) {
  vi v(n);
  for (int i = int(0); i < int(n); i++) v[i] = read() + off;
  return v;
}
template <class T>
void print(const vector<T>& v, int suc = 1) {
  for (int i = int(0); i < int(v.size()); i++)
    print(v[i], i == int(v.size()) - 1 ? suc : 2);
}
string readString() {
  string s;
  cin >> s;
  return s;
}
template <class T>
T sq(const T& t) {
  return t * t;
}
void yes(bool ex = true) {
  cout << "Yes" << endl;
  if (ex) exit(0);
}
void no(bool ex = true) {
  cout << "No" << endl;
  if (ex) exit(0);
}
constexpr ll ten(int n) { return n == 0 ? 1 : ten(n - 1) * 10; }
const ll infLL = LLONG_MAX / 3;
const int inf = INT_MAX / 2 - 100;
int topbit(signed t) { return t == 0 ? -1 : 31 - __builtin_clz(t); }
int topbit(ll t) { return t == 0 ? -1 : 63 - __builtin_clzll(t); }
int botbit(signed first) { return first == 0 ? 32 : __builtin_ctz(first); }
int botbit(ll first) { return first == 0 ? 64 : __builtin_ctzll(first); }
int popcount(signed t) { return __builtin_popcount(t); }
int popcount(ll t) { return __builtin_popcountll(t); }
bool ispow2(int i) { return i && (i & -i) == i; }
int mask(int i) { return (int(1) << i) - 1; }
bool inc(int first, int second, int c) {
  return first <= second && second <= c;
}
template <class t>
void mkuni(vc<t>& v) {
  sort(v.begin(), v.end());
  v.erase(unique(v.begin(), v.end()), v.end());
}
ll rand_int(ll l, ll r) {
  static random_device rd;
  static mt19937_64 gen(rd());
  return uniform_int_distribution<ll>(l, r)(gen);
}
using uint = unsigned;
using ull = unsigned long long;
uint mod = 1;
struct mint {
  uint v;
  mint(ll vv = 0) { s(vv % mod + mod); }
  mint& s(uint vv) {
    v = vv < mod ? vv : vv - mod;
    return *this;
  }
  mint operator-() const { return mint() - *this; }
  mint& operator+=(const mint& rhs) { return s(v + rhs.v); }
  mint& operator-=(const mint& rhs) { return s(v + mod - rhs.v); }
  mint& operator*=(const mint& rhs) {
    v = ull(v) * rhs.v % mod;
    return *this;
  }
  mint& operator/=(const mint& rhs) { return *this *= rhs.inv(); }
  mint operator+(const mint& rhs) const { return mint(*this) += rhs; }
  mint operator-(const mint& rhs) const { return mint(*this) -= rhs; }
  mint operator*(const mint& rhs) const { return mint(*this) *= rhs; }
  mint operator/(const mint& rhs) const { return mint(*this) /= rhs; }
  mint pow(int n) const {
    mint res(1), x(*this);
    while (n) {
      if (n & 1) res *= x;
      x *= x;
      n >>= 1;
    }
    return res;
  }
  mint inv() const { return pow(mod - 2); }
  friend ostream& operator<<(ostream& os, const mint& m) { return os << m.v; }
  bool operator<(const mint& r) const { return v < r.v; }
  bool operator==(const mint& r) const { return v == r.v; }
};
const int vmax = (1 << 21) + 10;
mint fact[vmax], finv[vmax], invs[vmax];
void initfact() {
  fact[0] = 1;
  for (int i = int(1); i < int(vmax); i++) {
    fact[i] = fact[i - 1] * i;
  }
  finv[vmax - 1] = fact[vmax - 1].inv();
  for (int i = vmax - 2; i >= 0; i--) {
    finv[i] = finv[i + 1] * (i + 1);
  }
  for (int i = vmax - 1; i >= 1; i--) {
    invs[i] = finv[i] * fact[i - 1];
  }
}
mint choose(int n, int k) { return fact[n] * finv[n - k] * finv[k]; }
mint binom(int first, int second) {
  return fact[first + second] * finv[first] * finv[second];
}
mint catalan(int n) {
  return binom(n, n) - (n - 1 >= 0 ? binom(n - 1, n + 1) : 0);
}
const int nmax = 110;
mint dp[2][nmax][nmax], buf[nmax][nmax];
mint ps[nmax][nmax * nmax];
mint ips[nmax][nmax * nmax];
signed main() {
  cin.tie(0);
  ios::sync_with_stdio(0);
  cout << fixed << setprecision(20);
  int n, m, k;
  cin >> n >> m >> k;
  cin >> mod;
  initfact();
  for (int i = int(0); i < int(k + 1); i++) {
    ps[i][0] = 1;
    for (int j = int(0); j < int(n * m); j++) ps[i][j + 1] = ps[i][j] * i;
  }
  for (int i = int(0); i < int(k + 1); i++) {
    ips[i][0] = 1;
    mint w = mint(i).inv();
    for (int j = int(0); j < int(n * m); j++) ips[i][j + 1] = ips[i][j] * w;
  }
  int cur = 0;
  dp[cur][0][0] = 1;
  for (int w = int(k) - 1; w >= int(0); w--) {
    int nx = cur ^ 1;
    memset(dp[nx], 0, sizeof(dp[nx]));
    for (int ij = int(0); ij < int(n + m + 1); ij++) {
      memset(buf, 0, sizeof(buf));
      for (int i = int(0); i < int(n + 1); i++)
        for (int j = int(max(ij - i, int(0))); j < int(m + 1); j++) {
          if (i + j == ij)
            buf[i][j] += dp[cur][i][j] * finv[i] * finv[j] *
                         ps[k - w][(n - i) * (m - j)];
          if (i < n) buf[i + 1][j] += buf[i][j] * ps[w + 1][j];
          if (j < m) buf[i][j + 1] += buf[i][j] * ps[w + 1][i];
          dp[nx][i][j] += buf[i][j] * finv[i + j - ij];
        }
    }
    memset(dp[cur], 0, sizeof(dp[cur]));
    for (int ij = int(0); ij < int(n + m + 1); ij++) {
      memset(buf, 0, sizeof(buf));
      for (int i = int(0); i < int(n + 1); i++)
        for (int j = int(max(ij - i, int(0))); j < int(m + 1); j++) {
          if (i + j == ij) buf[i][j] += dp[nx][i][j];
          dp[cur][i][j] += buf[i][j] * finv[i + j - ij];
          if (i < n) buf[i + 1][j] -= buf[i][j] * ps[w][j];
          if (j < m) buf[i][j + 1] -= buf[i][j] * ps[w][i];
        }
    }
    for (int i = int(0); i < int(n + 1); i++)
      for (int j = int(0); j < int(m + 1); j++) {
        dp[cur][i][j] *= ips[k - w][(n - i) * (m - j)] * fact[i] * fact[j];
      }
  }
  cout << dp[cur][n][m] << endl;
}
