1 #include<cstdio>
 2 #include<algorithm>
 3 #define rep(i,l,r) for (ll i=(l); i<=(r); i++)
 4 typedef long long ll;
 5 using namespace std;
 6
 7 const int N=100010;
 8 ll n,m,K,P,tot,num[N],pi[N],pk[N],sm[N];
 9
10 ll ksm(ll a,ll b,ll p){
11     ll res=1;
12     for (; b; a=a*a%p,b>>=1)
13         if (b & 1) res=res*a%p;
14     return res;
15 }
16
17 void exgcd(ll a,ll b,ll &x,ll &y){
18     if (!b){ x=1; y=0; return; }
19         else exgcd(b,a%b,y,x),y-=a/b*x;
20 }
21
22 ll inv(ll a,ll p){ ll x,y; exgcd(a,p,x,y); return (x%p+p)%p; }
23
24 ll F(ll n,ll pi,ll pk){ return n ? ksm(sm[pk],n/pk,pk)*sm[n%pk]%pk*F(n/pi,pi,pk)%pk : 1; }
25
26 ll exLucas(ll n,ll m,ll pi,ll pk){
27     if (n<m) return 0;
28     sm[0]=sm[1]=1; rep(i,2,pk) sm[i]=sm[i-1]*((i%pi)?i:1)%pk;
29     ll a=F(n,pi,pk),b=F(m,pi,pk),c=F(n-m,pi,pk),k=0;
30     for (ll i=n; i; i/=pi) k+=i/pi;
31     for (ll i=m; i; i/=pi) k-=i/pi;
32     for (ll i=n-m; i; i/=pi) k-=i/pi;
33     return a*inv(b*c%pk,pk)%pk*ksm(pi,k,pk)%pk;
34 }
35
36 ll CRT(ll r[],ll m[]){
37     ll res=0;
38     rep(i,1,tot) res=(res+P/m[i]*inv(P/m[i],m[i])%P*r[i])%P;
39     return res;
40 }
41
42 void Frac(ll n){
43     tot=0;
44     for (ll i=2; i*i<=n; i++) if (n%i==0){
45         pi[++tot]=i; pk[tot]=i; n/=i;
46         while (n%i==0) n/=i,pk[tot]=pk[tot]*i;
47     }
48     if (n>1) pi[++tot]=n,pk[tot]=n;
49 }
50
51 ll Pe(ll n,ll m,ll P){ ll res=1; rep(i,1,m) res=res*(n-i+P+1)%P; return res; }
52
53 int main(){
54     scanf("%lld%lld%lld%lld",&n,&m,&K,&P);
55     Frac(P); ll res=0;
56     rep(i,1,tot){
57         ll s=1;
58         if (pi[i]==pk[i] && pi[i]>m)
59             rep(j,1,m) s=s*(K+m-j)%pi[i]*ksm(j,pi[i]-2,pi[i])%pi[i];
60         else s=exLucas(K+m-1,m,pi[i],pk[i]);
61         s=Pe(s,n,pi[i]); res=(res+P/pi[i]*inv(P/pi[i],pi[i])%P*s)%P;
62     }
63     printf("%lld\n",res);
64     return 0;
65 }