BZOJ3782上学路线 组合数+容斥+CRT

Posted 2020-10-02 CQzhangyu

【BZOJ3782】上学路线

Description

小C所在的城市的道路构成了一个方形网格，它的西南角为(0,0)，东北角为(N,M)。小C家住在西南角，学校在东北角。现在有T个路口进行施工，小C不能通过这些路口。小C喜欢走最短的路径到达目的地，因此他每天上学时都只会向东或北行走；而小C又喜欢走不同的路径，因此他问你按照他走最短路径的规则，他可以选择的不同的上学路线有多少条。由于答案可能很大，所以小C只需要让你求出路径数mod P的值。

Input

第一行，四个整数N、M、T、P。

接下来的T行，每行两个整数，表示施工的路口的坐标。

Output

一行，一个整数，路径数mod P的值。

Sample Input

3 4 3 1019663265
3 0
1 1
2 2

Sample Output

8

HINT

1<=N,M<=10^10

0<=T<=200

p=1000003或p=1019663265

题解：从(0,0)走到(n,m)的总方案数=C(n+m,n)。

依旧考虑容斥，先将点排序，用f[i]表示从(0,0)走到(x[i],y[i])，途中不经过其它障碍的方案数，那么如果j在i的左下方，则f[i]-=f[j]*(从j走到i的方案数)。

然而1019663265不是质数？分解质因数的1019663265=3*5*6793*10007，分别求解，再用中国剩余定理合并即可。

EXGCD还能写错~

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
ll mod,ans;
ll n,m,P;
int T;
ll f[210],jc[1000010],jcc[1000010],ine[1000010];
struct node
{
	ll x,y;
}p[210];
bool cmp(node a,node b)
{
	return (a.x==b.x)?(a.y<b.y):(a.x<b.x);
}
ll C(ll a,ll b)
{
	if(a<b)	return 0;
	if(!b)	return 1;
	if(a<mod&&b<mod)	return jc[a]*jcc[a-b]%mod*jcc[b]%mod;
	return C(a%mod,b%mod)*C(a/mod,b/mod)%mod;
}
ll calc()
{
	memset(f,0,sizeof(f));
	int i,j;
	jc[0]=jcc[0]=1;
	ine[1]=1;
	for(i=2;i<mod;i++)	ine[i]=(mod-(mod/i)*ine[mod%i]%mod)%mod;
	for(i=1;i<mod;i++)	jc[i]=jc[i-1]*i%mod,jcc[i]=jcc[i-1]*ine[i]%mod;
	for(i=1;i<=T;i++)
	{
		f[i]=C(p[i].x+p[i].y,p[i].x);
		for(j=1;j<i;j++)
			if(p[i].x>=p[j].x&&p[i].y>=p[j].y)	f[i]=(f[i]-f[j]*C(p[i].x-p[j].x+p[i].y-p[j].y,p[i].x-p[j].x)%mod+mod)%mod;
	}
	return f[T];
}
void exgcd(ll a,ll b,ll &x,ll &y)
{
	if(!b)
	{
		x=1,y=0;
		return ;
	}
	exgcd(b,a%b,x,y);
	ll tmp=x;
	x=y,y=tmp-a/b*x;
}
ll work(ll a,ll c,ll b)
{
	ll x,y;
	exgcd(a,b,x,y),x=(x+b)%b;
	return x*a%P*c%P;
}
int main()
{
	scanf("%lld%lld%d%lld",&n,&m,&T,&P);
	int i;
	for(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(P==1000003)
	{
		mod=P,printf("%lld",calc());
		return 0;
	}
	ll a1,a2,a3,a4;
	mod=3,a1=calc();
	mod=5,a2=calc();
	mod=6793,a3=calc();
	mod=10007,a4=calc();
	ans=(ans+work(P/3,a1,3))%P;
	ans=(ans+work(P/5,a2,5))%P;
	ans=(ans+work(P/6793,a3,6793))%P;
	ans=(ans+work(P/10007,a4,10007))%P;
	printf("%lld",(ans+P)%P);
	return 0;
}//3 4 3 1000003 3 0 1 1 2 2