h, w, K, MOD = map(int, input().split())

hw = h * w
powK = [pow(i, K, MOD) for i in range(hw + 1)]
pow2 = [1] * (hw + 1)
for i in range(1, hw + 1):
    pow2[i] = pow2[i - 1] * 2 % MOD

A = [[-1] * w for _ in range(h)]
didj = [(0, 1), (1, 0), (0, -1), (-1, 0)]
ans = 0


def sol(c1, c2, si, sj):
    ii = -1
    jj = -1
    for i in range(h):
        for j in range(w):
            if A[i][j] != -1:
                continue

            for di, dj in didj:
                ni = i + di
                nj = j + dj
                if ni < 0 or ni >= h or nj < 0 or nj >= w:
                    continue
                if A[ni][nj] == 1:
                    ii = i
                    jj = j
                    break

        if ii != -1:
            break

    if ii == -1:
        global ans
        ans += powK[c1] * pow2[h * w - c2] % MOD
    else:
        if si * w + sj < ii * w + jj:
            A[ii][jj] = 1
            sol(c1 + 1, c2 + 1, si, sj)
        A[ii][jj] = 0
        sol(c1, c2 + 1, si, sj)
        A[ii][jj] = -1


for si in range(h):
    for sj in range(w):
        A[si][sj] = 1
        sol(1, 1, si, sj)
        A[si][sj] = -1

ans %= MOD
ans *= pow(pow2[h * w], MOD - 2, MOD)
ans %= MOD
inv = 0
for i in range(1, h * w + 1):
    inv += pow(i, MOD - 2, MOD)

ans *= inv
print(ans % MOD)