#include <iostream>
#include <stdio.h>
#include <string.h>
#include <map>
#include <unordered_map>
#include <algorithm>
#include <vector>
#include <queue>
#include <set>
#include <math.h>

using namespace std;

#define x first
#define y second
#define vi vector<int>
#define pb push_back
#define mp make_pair
#define pii pair<int, int>
#define all(x) x.begin(),x.end()
#define SZ(x) int(x.size())
#define rep(i, a, b) for(int i = a; i < b; i++)
#define per(i, a, b) for(int i = b-1; i >= a; i--)
#define DBG(x) cerr << (#x) << "=" << x << "\n";
#define ll long long
#define inf 1000000007
#define mod 998244353//1000000007
#define N 101

template<typename U, typename V> void Min(U &a, const V & b){if(a > b)a = b;}
template<typename U, typename V> void Max(U &a, const V & b){if(a < b)a = b;}
template<typename U, typename V> void add(U &a, const V & b){a = (a + b) % mod;}
template<typename U> U gcd(U a, U b){
    if(a == 0)return b;
    if(b == 0)return a;
    if(a >= b)return gcd(a % b, b);
    else return gcd(a, b % a);
}
int pow(int a, int b){
    int ans = 1;
    while(b){
        if(b & 1)ans = 1LL * ans * a % mod;
        a = 1LL * a * a % mod;
        b >>= 1;
    }
    return ans;
}
int pow(int a, int b, int c){
    int ans = 1;
    while(b){
        if(b & 1)ans = 1LL * ans * a % c;
        a = 1LL * a * a % c;
        b >>= 1;
    }
    return ans;
}

int dp[2][N][N], trans1[N][N][N], trans2[N][N][N], inv[N];
int main(){
    int T, i, j, k, m, n, K;
    int mo;
    scanf("%d%d%d%d", &n, &m, &K, &mo);
    inv[1] = 1;
    rep(i, 2, N){
        inv[i] = mo - mo / i * (ll)inv[mo % i] % mo;
        if(inv[i] == mo)inv[i] = 0;
    }
    rep(k, 1, K+1){
        rep(j, 0, m+1){
            trans1[k][j][0] = 1;
            if(j == m)trans1[k][j][1] = 0;
            else{
                ll w = pow(k, m-j, mo) - pow(k-1, m-j, mo);
                if(w < 0)w += mo;
                w = w * pow(K-k+1, j, mo) % mo;
                rep(i, 1, n+1)trans1[k][j][i] = trans1[k][j][i-1] * w % mo * inv[i] % mo;
            }
        }
    }
    rep(k, 1, K+1){
        rep(i, 0, n+1){
            trans2[k][i][0] = 1;
            if(i == 0)trans2[k][i][1] = 0;
            else{
                ll w = pow(K-k+1, i, mo) - pow(K-k, i, mo);
                if(w < 0)w += mo;
                w = w * pow(k, n-i, mo) % mo;
                rep(j, 1, m+1)trans2[k][i][j] = trans2[k][i][j-1] * w % mo * inv[j] % mo;
            }
        }
    }
    int q = 0;
    dp[0][0][0] = 1;
    rep(k, 1, K+1){
        memset(dp[q^1], 0, sizeof(dp[0]));
        rep(i, 0, n+1){
            rep(j, 0, m+1){
                rep(x, 0, n-i+1){
                    // dp[q^1][i][j] += c[n-i][x] * (K-k+1)^(j * x) * (k^(m-j) - (k-1)^(m-j))^x * inv[x] * dp[q][i][j]
                    dp[q^1][i+x][j] = (dp[q^1][i+x][j] + (ll)dp[q][i][j] * trans1[k][j][x] % mo) % mo;
                }
            }
        }
        q ^= 1;
        memset(dp[q^1], 0, sizeof(dp[0]));
        rep(i, 0, n+1){
            rep(j, 0, m+1){
                rep(x, 0, m-j+1){
                    // dp[q^1][i][j] += c[m-j][x] * k^((n-i)*x) * ((K-k+1)^i - (K-k)^i)^x * inv[x] * dp[q][i][j]
                    dp[q^1][i][j+x] = (dp[q^1][i][j+x] + (ll)dp[q][i][j] * trans2[k][i][x] % mo) % mo;
                }
            }
        }
        q ^= 1;
    }
    int ans = dp[q][n][m];
    rep(i, 1, n+1)ans = (ll)ans * i % mo;
    rep(i, 1, m+1)ans = (ll)ans * i % mo;
    printf("%d\n", ans);
}
