首先这个题需要dp.dp[i]=C(x[i]+y[i],x[i])-Σdp[j]*C(x[i]-x[j]+y[i]-y[j],x[i]-x[j])(x[i]>=x[j],y[i]>=y[j])。
然后就是喜闻乐见的lucas+CRT.
#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#define maxn 1000050
using namespace std;
long long n,m,t,mod,p[],ans[],inv1[][maxn],inv2[][maxn],up,f[],a[];
struct pnt
{
long long x,y;
}pt[];
bool cmp(pnt x,pnt y)
{
if (x.y!=y.y) return x.y<y.y;
return x.x<y.x;
}
long long f_pow(long long x,long long y,long long mod)
{
long long ans=,base=x;
while (y)
{
if (y&) ans=(ans*base)%mod;
base=(base*base)%mod;
y>>=;
}
return ans;
}
void get_mod(long long type)
{
up=type;
if (type==)
{
p[]=;inv1[][]=;inv2[][]=;
for (long long i=;i<=p[]-;i++)
{
inv1[][i]=inv1[][i-]*i%p[];
inv2[][i]=f_pow(inv1[][i],p[]-,p[]);
}
}
else
{
p[]=;p[]=;p[]=;p[]=;
for (long long i=;i<=;i++)
{
inv1[i][]=;inv2[i][]=;
for (long long j=;j<=p[i]-;j++)
{
inv1[i][j]=inv1[i][j-]*j%p[i];
inv2[i][j]=f_pow(inv1[i][j],p[i]-,p[i]);
}
}
}
}
long long comb(long long n,long long m,long long p,long long type)
{
if (n<m) return ;
return inv1[type][n]*inv2[type][m]%p*inv2[type][n-m]%p;
}
long long lucas(long long n,long long m,long long p,long long type)
{
if (!m) return ;
return comb(n%p,m%p,p,type)*lucas(n/p,m/p,p,type)%p;
}
long long CRT()
{
long long ret=;
for (long long i=;i<=up*up;i++)
ret=(ret+a[i]*(mod/p[i])%mod*f_pow(mod/p[i],p[i]-,p[i])%mod)%mod;
return ret%mod;
}
long long get_C(long long n,long long m)
{
for (long long i=;i<=up*up;i++) a[i]=lucas(n,m,p[i],i);
return CRT();
}
int main()
{
scanf("%lld%lld%lld%lld",&n,&m,&t,&mod);
for (long long i=;i<=t;i++) scanf("%d%d",&pt[i].x,&pt[i].y);
if (mod==) get_mod();else get_mod();
sort(pt+,pt+t+,cmp);t++;pt[t].x=n;pt[t].y=m;
for (long long i=;i<=t;i++)
{
f[i]=get_C(pt[i].x+pt[i].y,pt[i].x);
for (long long j=;j<=i-;j++)
{
if (pt[j].x<=pt[i].x && pt[j].y<=pt[i].y)
f[i]=(f[i]-f[j]*get_C(pt[i].x-pt[j].x+pt[i].y-pt[j].y,pt[i].x-pt[j].x)%mod+mod)%mod;
}
}
printf("%lld\n",(f[t]+mod)%mod);
return ;
}