1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
#include <bits/stdc++.h>

#define ll long long
#define ld long double
#define endl '\n'
#define st first
#define nd second
#define pb push_back
#define sz(x) (int)(x).size()
#define all(x) (x).begin(), (x).end()
#define FOR(i,l,r) for(int i=(l);i<=(r);i++)
#define ROF(i,r,l) for(int i=(r);i>=(l);i--)
using namespace std;
ll infl=1000000000000000007;
int inf=1000000007;
ll mod;

#ifdef DEBUG
auto&operator<<(auto&o,pair<auto,auto>p){return o<<"("<<p.first<<", "<<p.second<<")";}
auto operator<<(auto&o,auto x)->decltype(x.end(),o){o<<"{";int i=0;for(auto e:x)o<<","+!i++<<e;return o<<"}";}
#define debug(X...)cerr<<"["#X"]: ",[](auto...$){((cerr<<$<<"; "),...)<<endl;}(X)
#else
#define debug(...){}
#endif

const int N=1000007;

ll pot(ll x,ll p) {ll res=1;for(;p;p>>=1) {if(p&1) res=res*x%mod; x=x*x%mod;} return res;}

ll fac[N];
ll inv[N];
int n;

ll calc(vector<int>t)
{
	debug(t);
	int a=sz(t),S=0;
	FOR(i,0,a-1) S+=t[i];
	ll res=fac[S];
	FOR(i,0,a-1)
	{
		int it=a-1;
		FOR(j,1,t[i])
		{
			while(t[it]<j) it--;
			res=res*inv[t[i]-j+it-i+1]%mod;
		}
	}	
	return res;
}

map<vector<int>,ll>mem;

int cnt;

ll get(vector<int>t)
{
	ll X=mem[t];
	if(X) return X;
	cnt++;
	ll res=0;
	int a=sz(t);
	bool nie=0;
	FOR(i,1,a-1) if(t[i]<t[i-1]) nie=1;
	if(nie)
	{
		bool is=0;
		FOR(i,0,a-1)
		{
			if(t[i]>0&&(i==0||t[i-1]<t[i]))
			{
				t[i]--;
				res+=get(t);
				t[i]++;
				is=1;
			}
		}
		if(!is) res=1;
		res%=mod;
	}
	else 
	{
		vector<int>w=t;
		reverse(all(w));
		res=calc(w);
	}
	return mem[t]=res;
}

signed main()
{
    ios_base::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    int a,b,c;
    cin>>a>>b>>c>>mod;
    int t=a+b-1-c;
    n=a*b-t*(t+1)/2;
    fac[0]=1;
    for(ll i=1;i<=n;i++) fac[i]=fac[i-1]*i%mod;
    FOR(i,1,n) inv[i]=pot(i,mod-2);
    vector<int>r(a);
    FOR(i,0,a-1) r[i]=b;
    FOR(i,1,t) r[a-1-(t-i)]-=i;
    debug(r);
    ll ans=calc(r);
    if(t==0) ans=ans*ans%mod;
    else ans=ans*get(r)%mod;
    cout<<n<<" "<<ans<<endl;
    
    return 0;
}