Bzoj3782:上学路线:组合数学+Lucas+CRT+DP

题目链接：上学路线

设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);
}