标签:namespace 数组 bzoj 一个 turn freopen scanf 容斥 中国剩余定理
权限题的话qwq简单说一下题面吧
就是要从\((0,0)\)走到\((n,m)\)然后中间有\(T\)个点是坏的(不能走),然后每次只能往上或者往右走,求方案数\(mod\ P\)的值,\(1<=n,m<=10^{10}\),\(0<=T<=200\),\(或p=1000003或1019663265\)
?
先不管范围
首先先看如果\(T=0\)的话怎么求,那么直接是\(\binom {n+m} n\)或者\(\binom {n+m} m\)就好了
具体解释的话就是,枚举往右(或者上)走的是哪几步
?
但是如果现在\(T>0\)呢?
最简单粗暴的想法就是。。大力容斥,从总的里面减掉走到某个坏点然后再走到终点的,那现在的问题就是怎么求经过坏点的方案数
我们可以枚举经过的第一个坏点,那么经过坏点的方案数就是
\[
\sum(从(0,0)走到这个坏点,且中途不经过其他坏点的方案数)*(这个坏点走到当前点的方案数)
\]
会发现前面的部分跟我们要求的问题是一样的,而后面的部分可以直接用组合数求出
?
为了方便,我们给点编一个号,第\(i\)个点的坐标为\((x_i,y_i)\)
设\(f[i]\)表示从\((0,0)\)走到\((x_i,y_i)\)的方案数,那么可以得到:
\[
f[i]=\binom {x_i+y_i}{x_i}-\sum_{j在i左下且是坏的} f[j]*\binom {(x_i-x_j)+(y_i-y_j)}{x_i-x_j}
\]
具体写的时候,我们暴力将\((n,m)\)这个点加进去作为第\(T+1\)个点计算,那么最后的答案就是\(f[T+1]\)
然而。。这并没有结束
现在我们加上范围和模数,分两类情况讨论一下:
1、如果\(p=1000003\)
这个是一个质数,那直接Lucas定理就好啦,比较开心
2、如果\(p=1019663265\)
很明显不是一个质数,\(p=1019663265=3*5*6793*10007\)
? 那么这个时候我们还需要加上一个中国剩余定理才能用Lucas。。。
大概就是我们先用Lucas定理算出在模\(,3,5,6793,10007\)下的答案,然后用中国剩余定理CRT合并
具体的话这里不展开讲了,可以传送一下到 -->这里(施工中。。。可能晚点才会传上来qwq)
然后这题就愉快地做完啦
一些实现的小细节:
1、注意当\(p=1019663265\)时是超了long long的。。
2、为了方便后面的实现建议一开始先将点排个序
3、开个数组存模数什么的写的比较爽
代码大概长这个样子
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const int N=1000010;
struct xxx{
ll x,y;
friend bool operator < (xxx x,xxx y)
{return x.x==y.x?x.y<y.y:x.x<y.x;}
}a[210];
int fac[5][N],invfac[5][N],Mod[5],rec[5];
ll f[210];
ll n,m,T,p;
int Id;
ll Lucas(ll n,ll m);
void prework();
ll ex_gcd(ll a,ll b,ll &x,ll &y);
ll CRT(int *w,int *a,int n);
void solve(int n,int id);
ll C(int n,int m,int id){return n<m?0:1LL*fac[id][n]*invfac[id][m]%Mod[id]*invfac[id][n-m]%Mod[id];}
ll get_C(int n,int m,int id);
ll ksm(ll x,int y,int id);
int main(){
#ifndef ONLINE_JUDGE
freopen("a.in","r",stdin);
#endif
scanf("%lld%lld%lld%lld",&n,&m,&T,&p);
for (int i=1;i<=T;++i)
scanf("%lld%lld",&a[i].x,&a[i].y);
sort(a+1,a+1+T);
if (a[T].x!=n||a[T].y!=m) a[++T].x=n,a[T].y=m;
prework();
solve(T,Id);
}
ll Lucas(ll n,ll m,int id){
if (n<Mod[id]&&m<Mod[id]) return C(n,m,id);
return 1LL*C(n%Mod[id],m%Mod[id],id)*Lucas(n/Mod[id],m/Mod[id],id)%Mod[id];
}
void prework(){
Mod[0]=1000003; Mod[1]=3; Mod[2]=5; Mod[3]=6793; Mod[4]=10007;
if (p==Mod[0]) Id=0;
else Id=1;
for (int i=0;i<5;++i) fac[i][0]=1;
for (int j=0;j<5;++j)
for (int i=1;i<Mod[j];++i)
fac[j][i]=1LL*fac[j][i-1]*i%Mod[j];
for (int i=0;i<5;++i) invfac[i][Mod[i]-1]=ksm(fac[i][Mod[i]-1],Mod[i]-2,i);
for (int j=0;j<5;++j)
for (int i=Mod[j]-2;i>=0;--i)
invfac[j][i]=1LL*invfac[j][i+1]*(i+1)%Mod[j];
}
ll ex_gcd(ll a,ll b,ll &x,ll &y){
if (b==0){x=1; y=0; return a;}
ll gcd=ex_gcd(b,a%b,x,y),tmp=x;
x=y; y=tmp-a/b*y;
return gcd;
}
ll CRT(int *w,int *a,int n){
ll x,y,ret=0,mod=1,tmp;
for (int i=1;i<=n;++i) mod*=w[i];
for (int i=1;i<=n;++i){
tmp=mod/w[i];
ex_gcd(w[i],tmp,x,y);
ret=(ret+y*tmp*a[i])%mod;
}
return (ret+mod)%mod;
}
ll get_C(ll n,ll m,int id){
if (!n) return 1;
if (n<m) return 0;
if (id==0) return Lucas(n,m,id);
for (int i=1;i<5;++i)
rec[i]=Lucas(n,m,i);
return CRT(Mod,rec,4);
}
void solve(int n,int id){
for (int i=1;i<=n;++i){
f[i]=get_C(a[i].x+a[i].y,a[i].x,id);
for (int j=1;j<i;++j)
if (a[j].x<=a[i].x&&a[j].y<=a[i].y)
f[i]=(f[i]+p-f[j]*get_C(a[i].x-a[j].x+a[i].y-a[j].y,a[i].x-a[j].x,id)%p)%p;
}
printf("%lld\n",f[n]);
}
ll ksm(ll x,int y,int id){
ll ret=1,base=x;
for (;y;y>>=1,base=base*base%Mod[id])
if (y&1) ret=ret*base%Mod[id];
return ret;
}
标签:namespace 数组 bzoj 一个 turn freopen scanf 容斥 中国剩余定理
原文地址:https://www.cnblogs.com/yoyoball/p/9235889.html