BZOJ-3782 上學路線(dp+Lucas定理+CRT)
阿新 • • 發佈:2020-12-08
題目描述
矩形網格起點為 \((0,0)\),終點為 \((n,m)\),只能向右或向上走,有 \(t\) 個壞點不能走,求方案數。
資料範圍:\(1\leq n,m\leq10^{10},0\leq t\leq 200,p=10000003\) 或 \(p=1019663265\)。
分析
把 $ (n,m)$ 也當成壞點,編號即 \(t+1\),把所有壞點排序。
設 \(dp[i]\) 為到達第 \(i\) 個壞點的合法方案數(途中不經過其他壞點),答案即 \(dp[t+1]\)。但是很難直接求 \(dp[i]\),考慮用總方案數減去經過壞點的個數的方案數計算 $dp[i] $。
狀態轉移方程為:
\(p=1000003\),直接用 $\text{Lucas} $ 定理算組合數。
\(p=1019663265=3\times 5\times 6793\times 10007\),不是一個質數,且 \(p\) 太大不能用 \(\text{exLucas}\),分別求出 \(\dbinom{n}{m}\) 在模 \(3\),模 \(5\),模 \(6793\),模 \(10007\) 下的值 \(a_1,a_2,a_3,a_4\)
程式碼
#include<bits/stdc++.h> using namespace std; long long n,m,t,p; struct Point { long long x,y; }P[20010]; long long dp[20010]; bool cmp(Point A,Point B) { if(A.x==B.x) return A.y<B.y; return A.x<B.x; } long long quick_pow(long long a,long long b,long long p) { long long ans=1; while(b) { if(b&1) ans=ans*a%p; a=a*a%p; b>>=1; } return ans%p; } const int N=1e6; struct mod1 { const int mod=1000003; long long fac[N+10],inv[N+10]; void init() { fac[0]=inv[0]=1; for(int i=1;i<=N;i++) fac[i]=fac[i-1]*i%mod; inv[N]=quick_pow(fac[N],mod-2,mod); for(int i=N-1;i>=1;i--) inv[i]=inv[i+1]*(i+1)%mod; } long long C(long long n,long long m) { if(m==0) return 1; if(m>n) return 0; return fac[n]*inv[m]%mod*inv[n-m]%mod; } long long Lucas(long long n,long long m) { if(m==0) return 1; if(m>n) return 0; return C(n%mod,m%mod)*Lucas(n/mod,m/mod)%mod; } }A; struct mod2 { const long long mod=1ll*1019663265; long long prime[5]={0,3,5,6793,10007}; long long fac[5][10010],inv[5][10010]; long long exgcd(long long a,long long b,long long &x,long long &y) { if(b==0) { x=1; y=0; return a; } long long gcd=exgcd(b,a%b,x,y); int temp=x;x=y;y=temp-a/b*y; return gcd; } void init() { for(int op=1;op<=4;op++) { int _N=prime[op]-1; fac[op][0]=inv[op][0]=1; for(int i=1;i<=_N;i++) fac[op][i]=fac[op][i-1]*i%prime[op]; inv[op][_N]=quick_pow(fac[op][_N],prime[op]-2,prime[op]); for(int i=_N-1;i>=1;i--) inv[op][i]=inv[op][i+1]*(i+1)%mod; } } long long C(long long n,long long m,long long op) { if(m==0) return 1; if(m>n) return 0; return fac[op][n]*inv[op][n-m]%prime[op]*inv[op][m]%prime[op]; } long long Lucas(long long n,long long m,long long op) { if(m==0) return 1; if(m>n) return 0; return C(n%prime[op],m%prime[op],op)*Lucas(n/prime[op],m/prime[op],op)%prime[op]; } long long a[5],b[5]; long long solve(long long n,long long m) { long long ans=0,lcm=1,x,y; for(int i=1;i<=4;i++) lcm=lcm*prime[i]; for(int i=1;i<=4;i++) a[i]=Lucas(n,m,i); for(int i=1;i<=4;i++) b[i]=prime[i]; for(int i=1;i<=4;i++) a[i]=(a[i]%b[i]+b[i])%b[i]; for(int i=1;i<=4;i++) { long long temp=lcm/prime[i]; exgcd(temp,b[i],x,y); x=(x%b[i]+b[i])%b[i]; ans=(ans+temp*x%lcm*a[i]%lcm)%lcm; } return (ans+lcm)%lcm; } }B; int main() { cin>>n>>m>>t>>p; if(p==1000003) A.init(); if(p==1019663265) B.init(); for(int i=1;i<=t;i++) scanf("%lld %lld",&P[i].x,&P[i].y); P[t+1]={n,m}; sort(P+1,P+2+t,cmp); for(int i=1;i<=t+1;i++) { if(p==1000003) dp[i]=A.Lucas(P[i].x+P[i].y,P[i].x); if(p==1019663265) dp[i]=B.solve(P[i].x+P[i].y,P[i].x); //cout<<P[i].x<<" "<<P[i].y<<" "<<B.solve(P[i].x+P[i].y,P[i].x)<<endl; for(int j=1;j<=i-1;j++) { if(p==1000003) dp[i]=(dp[i]-dp[j]*A.Lucas((P[i].x+P[i].y)-(P[j].x+P[j].y),P[i].x-P[j].x)%p+p)%p; if(p==1019663265) dp[i]=(dp[i]-dp[j]*B.solve((P[i].x+P[i].y)-(P[j].x+P[j].y),P[i].x-P[j].x)%p+p)%p; } } cout<<dp[t+1]<<endl; return 0; }