#include <bits/stdc++.h>
using namespace std;
template <class S, class T>
ostream& operator<<(ostream& o, const pair<S, T>& p) {
  return o << "(" << p.first << "," << p.second << ")";
}
template <class T>
ostream& operator<<(ostream& o, const vector<T>& vc) {
  o << "{";
  for (const T& v : vc) o << v << ",";
  o << "}";
  return o;
}
using ll = long long;
template <class T>
using V = vector<T>;
template <class T>
using VV = vector<vector<T>>;
constexpr ll TEN(int n) { return (n == 0) ? 1 : 10 * TEN(n - 1); }
unsigned int mod = 1;
struct ModInt {
  using uint = unsigned int;
  using ll = long long;
  using ull = unsigned long long;
  uint v;
  ModInt() : v(0) {}
  ModInt(ll _v) : v(normS(_v % mod + mod)) {}
  explicit operator bool() const { return v != 0; }
  static uint normS(const uint& x) { return (x < mod) ? x : x - mod; }
  static ModInt make(const uint& x) {
    ModInt m;
    m.v = x;
    return m;
  }
  ModInt operator+(const ModInt& b) const { return make(normS(v + b.v)); }
  ModInt operator-(const ModInt& b) const { return make(normS(v + mod - b.v)); }
  ModInt operator-() const { return make(normS(mod - v)); }
  ModInt operator*(const ModInt& b) const { return make((ull)v * b.v % mod); }
  ModInt operator/(const ModInt& b) const { return *this * b.inv(); }
  ModInt& operator+=(const ModInt& b) { return *this = *this + b; }
  ModInt& operator-=(const ModInt& b) { return *this = *this - b; }
  ModInt& operator*=(const ModInt& b) { return *this = *this * b; }
  ModInt& operator/=(const ModInt& b) { return *this = *this / b; }
  ModInt& operator++(int) { return *this = *this + 1; }
  ModInt& operator--(int) { return *this = *this - 1; }
  ll extgcd(ll a, ll b, ll& x, ll& y) const {
    ll p[] = {a, 1, 0}, q[] = {b, 0, 1};
    while (*q) {
      ll t = *p / *q;
      for (int i = 0; i < (int)(3); i++) swap(p[i] -= t * q[i], q[i]);
    }
    if (p[0] < 0)
      for (int i = 0; i < (int)(3); i++) p[i] = -p[i];
    x = p[1], y = p[2];
    return p[0];
  }
  ModInt inv() const {
    ll x, y;
    extgcd(v, mod, x, y);
    return make(normS(x + mod));
  }
  ModInt pow(ll p) const {
    if (p < 0) return inv().pow(-p);
    ModInt a = 1;
    ModInt x = *this;
    while (p) {
      if (p & 1) a *= x;
      x *= x;
      p >>= 1;
    }
    return a;
  }
  bool operator==(const ModInt& b) const { return v == b.v; }
  bool operator!=(const ModInt& b) const { return v != b.v; }
  friend istream& operator>>(istream& o, ModInt& x) {
    ll tmp;
    o >> tmp;
    x = ModInt(tmp);
    return o;
  }
  friend ostream& operator<<(ostream& o, const ModInt& x) { return o << x.v; }
};
using mint = ModInt;
V<mint> fact, ifact;
void InitFact(int N) {
  fact.resize(N);
  ifact.resize(N);
  fact[0] = 1;
  for (int i = 1; i <= (int)(N - 1); i++) fact[i] = fact[i - 1] * i;
  ifact[N - 1] = fact[N - 1].inv();
  for (int i = N - 2; i >= 0; i--) ifact[i] = ifact[i + 1] * (i + 1);
}
mint C[1000][1000];
mint pw[101][10001];
int bsr(int x) { return 31 - __builtin_clz(x); }
using D = double;
const D pi = acos(-1);
using Pc = complex<D>;
void fft(bool type, vector<Pc>& c) {
  static vector<Pc> buf[30];
  int N = c.size();
  int s = bsr(N);
  assert(1 << s == N);
  if (buf[s].empty()) {
    buf[s] = vector<Pc>(N);
    for (int i = 0; i < (int)(N); i++) buf[s][i] = polar<D>(1, i * 2 * pi / N);
  }
  vector<Pc> a(N), b(N);
  copy(begin(c), end(c), begin(a));
  for (int i = 1; i <= (int)(s); i++) {
    int W = 1 << (s - i);
    for (int y = 0; y < N / 2; y += W) {
      Pc now = buf[s][y];
      if (type) now = conj(now);
      for (int x = 0; x < (int)(W); x++) {
        auto l = a[y << 1 | x];
        auto r = now * a[y << 1 | x | W];
        b[y | x] = l + r;
        b[y | x | N >> 1] = l - r;
      }
    }
    swap(a, b);
  }
  copy(begin(a), end(a), begin(c));
}
template <class Mint>
vector<Mint> multiply_fft(const vector<Mint>& x, const vector<Mint>& y) {
  if (x.empty() || y.empty()) return {};
  const int B = 15;
  const int K = 2;
  int S = x.size() + y.size() - 1;
  int N = 1;
  while (N < S) N *= 2;
  vector<Pc> a[K], b[K];
  for (int t = 0; t < (int)(K); t++) {
    a[t] = vector<Pc>(N);
    b[t] = vector<Pc>(N);
    for (int i = 0; i < (int)(x.size()); i++)
      a[t][i] = Pc((x[i].v >> (t * B)) & ((1 << B) - 1), 0);
    for (int i = 0; i < (int)(y.size()); i++)
      b[t][i] = Pc((y[i].v >> (t * B)) & ((1 << B) - 1), 0);
    fft(false, a[t]);
    fft(false, b[t]);
  }
  vector<Mint> z(S);
  vector<Pc> c(N);
  Mint base = 1;
  for (int t = 0; t < (int)(K + K - 1); t++) {
    fill_n(begin(c), N, Pc(0, 0));
    for (int xt = 0; xt < (int)(K); xt++) {
      int yt = t - xt;
      if (0 <= yt && yt < K) {
        for (int i = 0; i < (int)(N); i++) c[i] += a[xt][i] * b[yt][i];
      }
    }
    fft(true, c);
    for (int i = 0; i < (int)(S); i++) {
      c[i] *= 1.0 / N;
      z[i] += Mint(ll(round(c[i].real()))) * base;
    }
    base *= 1 << B;
  }
  return z;
}
int main() {
  cin.tie(0);
  ios::sync_with_stdio(false);
  cout << fixed << setprecision(20);
  int H, W, K;
  cin >> H >> W >> K >> mod;
  for (int i = 0; i < (int)(101); i++) {
    pw[i][0] = 1;
    for (int j = 1; j <= (int)(10000); j++) pw[i][j] = pw[i][j - 1] * i;
  }
  InitFact(1000);
  mint cur[101][101] = {};
  mint nxt[101][101] = {};
  cur[0][0] = 1;
  for (int t = 0; t < (int)(K); t++) {
    for (int a = 0; a < (int)(H + 1); a++)
      for (int b = 0; b < (int)(W + 1); b++)
        cur[a][b] = cur[a][b] * fact[H - a] * fact[W - b] / pw[K - t][a * b] *
                    pw[t + 1][(H - a) * (W - b)];
    {
      const int B = 201;
      V<mint> f(B * B), g(B * B);
      for (int i = 0; i < (int)(H + 1); i++)
        for (int j = 0; j < (int)(W + 1); j++) {
          f[i * B + j] = ifact[i] * ifact[j];
          g[i * B + j] = cur[i][j];
        }
      auto h = multiply_fft(f, g);
      for (int a = 0; a < (int)(H + 1); a++)
        for (int b = 0; b < (int)(W + 1); b++) nxt[a][b] = h[a * B + b];
    }
    for (int c = 0; c < (int)(H + 1); c++)
      for (int d = 0; d < (int)(W + 1); d++)
        nxt[c][d] *= ifact[H - c] * ifact[W - d] * pw[K - t][c * d] /
                     pw[t + 1][(H - c) * (W - d)];
    for (int a = 0; a < (int)(H + 1); a++)
      for (int b = 0; b < (int)(W + 1); b++)
        cur[a][b] = nxt[a][b], nxt[a][b] = 0;
    if (t == K - 1) {
      for (int a = 0; a < (int)(H + 1); a++)
        for (int b = 0; b < (int)(W + 1); b++) nxt[a][b] = cur[a][b];
    } else {
      for (int a = 0; a < (int)(H + 1); a++)
        for (int b = 0; b < (int)(W + 1); b++)
          cur[a][b] = cur[a][b] * fact[H - a] * fact[W - b] /
                      pw[K - t - 1][a * b] * pw[t + 1][(H - a) * (W - b)];
      {
        const int B = 201;
        V<mint> f(B * B), g(B * B);
        for (int i = 0; i < (int)(H + 1); i++)
          for (int j = 0; j < (int)(W + 1); j++) {
            f[i * B + j] = ifact[i] * ifact[j] * ((i + j) & 1 ? -1 : 1);
            g[i * B + j] = cur[i][j];
          }
        auto h = multiply_fft(f, g);
        for (int a = 0; a < (int)(H + 1); a++)
          for (int b = 0; b < (int)(W + 1); b++) nxt[a][b] = h[a * B + b];
      }
      for (int c = 0; c < (int)(H + 1); c++)
        for (int d = 0; d < (int)(W + 1); d++)
          nxt[c][d] *= ifact[H - c] * ifact[W - d] * pw[K - t - 1][c * d] /
                       pw[t + 1][(H - c) * (W - d)];
    }
    for (int a = 0; a < (int)(H + 1); a++)
      for (int b = 0; b < (int)(W + 1); b++)
        cur[a][b] = nxt[a][b], nxt[a][b] = 0;
  }
  cout << cur[H][W] << endl;
}
