Problem: For each of the K^{NM} ways to write an integer between 1 and K (inclusive) in every square in a square grid with N rows and M columns, find the value defined below, then compute the sum of all those K^{NM} values, modulo D.

* For each of the NM squares, find the minimum among the N+M-1 integers written in the square's row or the square's column. The value defined for the grid is the product of all these NM values.

Constraints

* 1 \leq N,M,K \leq 100
* 10^8 \leq D \leq 10^9
* N,M,K, and D are integers.
* D is prime.

Input

Input is given from Standard Input in the following format:


N M K D


Output

Print the sum of the K^{NM} values, modulo D.

Examples

Input

2 2 2 998244353


Output

35


Input

2 3 4 998244353


Output

127090


Input

31 41 59 998244353


Output

827794103

Solution (in CPP):

Solution:
// ※※※ 解答不能 ※※※
// ecnerwal氏
// https://atcoder.jp/contests/agc039/submissions/7876684
#include<bits/stdc++.h>
using namespace std;
namespace fft{

#define rep(i, a, b) for(int i = a; i < (b); ++i)
#define trav(a, x) for(auto& a : x)
#define all(x) x.begin(), x.end()
#define sz(x) (int)(x).size()
typedef long long ll;
typedef pair<int, int> pii;
typedef vector<int> vi;
// FFT
using dbl = double;
struct num{ /// start-hash
    dbl x, y;
    num(dbl x_ = 0, dbl y_ = 0) : x(x_), y(y_) { }
};
inline num operator+(num a, num b) { return num(a.x + b.x, a.y + b.y); }
inline num operator-(num a, num b) { return num(a.x - b.x, a.y - b.y); }
inline num operator*(num a, num b) { return num(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x); }
inline num conj(num a) { return num(a.x, -a.y); }
inline num inv(num a) { dbl n = (a.x*a.x+a.y*a.y); return num(a.x/n,-a.y/n); }
/// end-hash

using vn = vector<num>;
vi rev({0, 1});
vn rt(2, num(1)), fa, fb;

inline void init(int n) { /// start-hash
    if (n <= sz(rt)) return;
    rev.resize(n);
    rep(i,0,n) rev[i] = (rev[i>>1] | ((i&1)*n)) >> 1;
    rt.reserve(n);
    for (int k = sz(rt); k < n; k *= 2) {
        rt.resize(2*k);
        double a=M_PI/k; num z(cos(a),sin(a)); // FFT
        rep(i,k/2,k) rt[2*i] = rt[i], rt[2*i+1] = rt[i]*z;
    }
} /// end-hash

inline void fft(vector<num> &a, int n) { /// start-hash
    init(n);
    int s = __builtin_ctz(sz(rev)/n);
    rep(i,0,n) if (i < rev[i]>>s) swap(a[i], a[rev[i]>>s]);
    for (int k = 1; k < n; k *= 2)
        for (int i = 0; i < n; i += 2 * k) rep(j,0,k) {
            num t = rt[j+k] * a[i+j+k];
            a[i+j+k] = a[i+j] - t;
            a[i+j] = a[i+j] + t;
        }
} /// end-hash

// Complex/NTT
vn multiply(vn a, vn b) { /// start-hash
    int s = sz(a) + sz(b) - 1;
    if (s <= 0) return {};
    int L = s > 1 ? 32 - __builtin_clz(s-1) : 0, n = 1 << L;
    a.resize(n), b.resize(n);
    fft(a, n);
    fft(b, n);
    num d = inv(num(n));
    rep(i,0,n) a[i] = a[i] * b[i] * d;
    reverse(a.begin()+1, a.end());
    fft(a, n);
    a.resize(s);
    return a;
} /// end-hash

// Complex/NTT power-series inverse
// Doubles b as b[:n] = (2 - a[:n] * b[:n/2]) * b[:n/2]
vn inverse(const vn& a) { /// start-hash
    if (a.empty()) return {};
    vn b({inv(a[0])});
    b.reserve(2*a.size());
    while (sz(b) < sz(a)) {
        int n = 2*sz(b);
        b.resize(2*n, 0);
        if (sz(fa) < 2*n) fa.resize(2*n);
        fill(fa.begin(), fa.begin()+2*n, 0);
        copy(a.begin(), a.begin()+min(n,sz(a)), fa.begin());
        fft(b, 2*n);
        fft(fa, 2*n);
        num d = inv(num(2*n));
        rep(i, 0, 2*n) b[i] = b[i] * (2 - fa[i] * b[i]) * d;
        reverse(b.begin()+1, b.end());
        fft(b, 2*n);
        b.resize(n);
    }
    b.resize(a.size());
    return b;
} /// end-hash

// Integer multiply mod m (num = complex) /// start-hash
vi multiply_mod(const vi& a, const vi& b, int m) {
    int s = sz(a) + sz(b) - 1;
    if (s <= 0) return {};
    int L = s > 1 ? 32 - __builtin_clz(s-1) : 0, n = 1 << L;
    if (sz(fa) < n) fa.resize(n);
    if (sz(fb) < n) fb.resize(n);

    rep(i,0,sz(a)) fa[i] = num(a[i] & ((1<<15)-1), a[i] >> 15);
    fill(fa.begin()+sz(a), fa.begin() + n, 0);
    rep(i,0,sz(b)) fb[i] = num(b[i] & ((1<<15)-1), b[i] >> 15);
    fill(fb.begin()+sz(b), fb.begin() + n, 0);

    fft(fa, n);
    fft(fb, n);
    double r0 = 0.5 / n; // 1/2n
    rep(i,0,n/2+1) {
        int j = (n-i)&(n-1);
        num g0 = (fb[i] + conj(fb[j])) * r0;
        num g1 = (fb[i] - conj(fb[j])) * r0;
        swap(g1.x, g1.y); g1.y *= -1;
        if (j != i) {
            swap(fa[j], fa[i]);
            fb[j] = fa[j] * g1;
            fa[j] = fa[j] * g0;
        }
        fb[i] = fa[i] * conj(g1);
        fa[i] = fa[i] * conj(g0);
    }
    fft(fa, n);
    fft(fb, n);
    vi r(s);
    rep(i,0,s) r[i] = int((ll(fa[i].x+0.5)
                + (ll(fa[i].y+0.5) % m << 15)
                + (ll(fb[i].x+0.5) % m << 15)
                + (ll(fb[i].y+0.5) % m << 30)) % m);
    return r;
} /// end-hash

} // namespace fft

template <int& MOD_> struct modnum {
    static constexpr int& MOD = MOD_;
    //static_assert(MOD_ > 0, "MOD must be positive");

private:
    using ll = long long;
    int v;
    static int minv(int a, int m) {
        a %= m;
        assert(a);
        return a == 1 ? 1 : int(m - ll(minv(m, a)) * ll(m) / a);
    }

public:
    modnum() : v(0) {}
    modnum(ll v_) : v(int(v_ % MOD)) { if (v < 0) v += MOD; }
    explicit operator int() const { return v; }
    friend std::ostream& operator << (std::ostream& out, const modnum& n) { return out << int(n); }
    friend std::istream& operator >> (std::istream& in, modnum& n) { ll v_; in >> v_; n = modnum(v_); return in; }

    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 inv() const {
        modnum res;
        res.v = minv(v, MOD);
        return res;
    }
    modnum neg() const {
        modnum res;
        res.v = v ? MOD-v : 0;
        return res;
    }

    modnum operator- () const {
        return neg();
    }
    modnum operator+ () const {
        return modnum(*this);
    }

    modnum& operator ++ () {
        v ++;
        if (v == MOD) v = 0;
        return *this;
    }
    modnum& operator -- () {
        if (v == 0) v = MOD;
        v --;
        return *this;
    }
    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.inv();
    }

    friend modnum operator ++ (modnum& a, int) { modnum r = a; ++a; return r; }
    friend modnum operator -- (modnum& a, int) { modnum r = a; --a; return r; }
    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; }
};

int MOD;
using num = modnum<MOD>;

vector<num> multiply(const vector<num>& a, const vector<num>& b) {
    vector<int> ai(a.begin(), a.end());
    vector<int> bi(b.begin(), b.end());
    vector<int> r = fft::multiply_mod(ai, bi, MOD);
    return vector<num>(r.begin(), r.end());
}

const int MAXN = 120;
const int MAXM = 120;
int N, M, K;
num dp[MAXN][MAXM];

const int MAXV = 1.1e3;
const int V = 1010;
num fact[MAXV];
num ifact[MAXV];
num prepow[MAXV][MAXV];

int main(){
    ios_base::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    cin >> N >> M >> K >> MOD;
    fact[0] = 1;
    for(int i = 1; i <= V; i++) fact[i] = fact[i - 1] * num(i);
    ifact[V] = fact[V].inv();
    for(int i = V; i > 0; i--) ifact[i - 1] = ifact[i] * num(i);
    for(int b = 0; b <= V; b++){
        prepow[b][0] = 1;
        for(int e = 1; e <= V; e++) prepow[b][e] = prepow[b][e - 1] * num(b);
    }
    for(int i = 0; i <= N; i++) for(int j = 0; j <= M; j++) dp[i][j] = 0;
    dp[0][0] = fact[N] * fact[M];
    for(int k = 0; k < 2 * K; k++){
        if(k % 2 == 0){
            // increment row
            for(int j = 0; j <= M; j++){
                num rowFac = (prepow[k / 2 + 1][M - j] - prepow[k / 2][M - j]) * prepow[K - k / 2][j];
                num powRowFac = 1;
                vector<num> polyFac(N + 1);
                for(int di = 0; di <= N; di++, powRowFac *= rowFac){
                    polyFac[di] = powRowFac * ifact[di];
                }
                vector<num> polyDP(N + 1);
                for(int i = 0; i <= N; i++) polyDP[i] = dp[i][j];
                vector<num> polyNDP = multiply(polyDP, polyFac);
                for(int i = 0; i <= N; i++) dp[i][j] = polyNDP[i];
            }
        }else{
            // increment col
            for(int i = 0; i <= N; i++){
                num colFac = (prepow[K - k / 2][i] - prepow[K - 1 - k / 2][i]) * prepow[k / 2 + 1][N - i];
                num powColFac = 1;
                vector<num> polyFac(M + 1);
                for(int dj = 0; dj <= M; dj++, powColFac *= colFac) polyFac[dj] = powColFac * ifact[dj];
                vector<num> polyDP(M + 1);
                for(int j = 0; j <= M; j++) polyDP[j] = dp[i][j];
                vector<num> polyNDP = multiply(polyDP, polyFac);
                for(int j = 0; j <= M; j++) dp[i][j] = polyNDP[j];
            }
        }
    }
    cout << dp[N][M] << '\n';
    return 0;
}