题目链接:上学路线
设dp[i]为到了第i个坏点且是第一次到达坏点的方案数
那么dp[i]=C(x[i]+y[i],x[i])-∑C(x[i]+y[i]-x[j]-y[j],x[i]-x[j])dp[j](x[j]<=x[i]&&y[j]<=y[i])
对于mod=100003Lucas直接上
对于mod=1019663265分解质因数得到3*5*6973*10007然后每个因数分别Lucas后CRT合并
#include<cstdio>
#include<cstdlib>
#include<iostream>
#include<algorithm>
#define ll long long
using namespace std;
const int maxp=1000010;
const int maxn=1010;
ll n,m,T,mod,dp[maxn];
struct point{ll x,y;}p[maxn];
bool cmp(const point &a,const point &b){
return a.x<b.x||(a.x==b.x&&a.y<b.y);
}
namespace Type1{
ll jc[maxp],inv[maxp];
void worker(){
jc[0]=1ll; inv[0]=inv[1]=1ll;
for (int i=1;i<mod;++i) jc[i]=jc[i-1]*i%mod;
for (int i=2;i<mod;++i) inv[i]=(mod-mod/i)*inv[mod%i]%mod;
for (int i=1;i<mod;++i) inv[i]=inv[i]*inv[i-1]%mod;
}
ll C(ll n,ll m){
if (n<m) return 0;
return jc[n]*inv[m]*inv[n-m]%mod;
}
ll Lucas(ll n,ll m){
if (!m) return 1;
return C(n%mod,m%mod)*Lucas(n/mod,m/mod)%mod;
}
}
namespace Type2{
ll jc[4][10010],inv[4][10010];
int p[]={3,5,6793,10007};
void worker(){
for (int k=0;k<4;++k){
int tmp=p[k];
jc[k][0]=1ll; inv[k][0]=inv[k][1]=1ll;
for (int i=1;i<tmp;++i) jc[k][i]=jc[k][i-1]*i%mod;
for (int i=2;i<tmp;++i) inv[k][i]=(tmp-tmp/i)*inv[k][tmp%i]%tmp;
for (int i=1;i<tmp;++i) inv[k][i]=inv[k][i]*inv[k][i-1]%tmp;
}
}
ll C(ll n,ll m,int t){
if (n<m) return 0;
int tmp=p[t];
return jc[t][n]*inv[t][m]%tmp*inv[t][n-m]%tmp;
}
ll Lucas(ll n,ll m,int t){
int tmp=p[t];
if (!m) return 1;
return C(n%tmp,m%tmp,t)*Lucas(n/tmp,m/tmp,t)%tmp;
}
void exgcd(ll a,ll b,ll &x,ll &y){
ll tmp;
if (!b){x=1;y=0;return;}
exgcd(b,a%b,x,y);
tmp=x; x=y; y=tmp-a/b*y;
}
ll Crt(ll n,ll m){
ll a[4]={0},ret=0;
for (int i=0;i<4;++i) a[i]=Lucas(n,m,i);
for (int i=0;i<4;++i){
ll tmp=mod/p[i],x,y;
exgcd(tmp,p[i],x,y);
x=(x+p[i])%p[i];
ret+=a[i]*x%mod*tmp%mod;
}return ret;
}
}
ll calc(ll n,ll m){
if (mod==1000003) return Type1::Lucas(n,m);
else return Type2::Crt(n,m);
}
int main(){
scanf("%lld%lld%lld%lld",&n,&m,&T,&mod);
for (int i=1;i<=T;++i) scanf("%lld%lld",&p[i].x,&p[i].y);
p[++T].x=n; p[T].y=m;
sort(p+1,p+T+1,cmp);
if (mod==1000003) Type1::worker();
if (mod==1019663265) Type2::worker();
for (int i=1;i<=T;++i){
ll tmp=p[i].x+p[i].y;
dp[i]=calc(tmp,p[i].x);
for (int j=1;j<i;++j)
if (p[j].y<=p[i].y) dp[i]=(dp[i]-1ll*dp[j]*calc(tmp-p[j].x-p[j].y,p[i].x-p[j].x)%mod)%mod;
}printf("%lld",(dp[T]+mod)%mod);
}