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

来源：互联网发布：258优化网站编辑：程序博客网时间：2024/04/30 21:34

题目链接：上学路线

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

1 0