bzoj 3782 上学路线

3782: 上学路线

Time Limit: 10 Sec   Memory Limit: 128 MB
Submit: 261   Solved: 106
[ Submit][ Status][ Discuss]

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

Source

【分析】

这道题可以DP

dp[i]表示从起点走到编号为i的坏点且不经过其它坏点的路径方案数。

那么dp[i]=C(x[i]+y[i],x[i]) - Σ(j,x[j]<=x[i],y[j]<=y[i]) C(x[i]-x[j]+y[i]-y[j],x[i]-x[j])*dp[j]

对于p为质数，可以用Lucas定理解决

对于p不是质数，可以用Lucas+CRT解决

【代码】

//bzoj 3782 上学路线
#include<bits/stdc++.h>
#define ll long long
#define M(a) memset(a,0,sizeof a)
#define fo(i,j,k) for(int i=j;i<=k;i++)
using namespace std;
const int mxn=1000005;
int T,P;//mod;
ll n,m,dp[205];
struct point {ll x,y;} p[205];
ll prime[5]={0,3,5,6793,10007};
inline bool comp(point u,point v)
{
	return (u.x==v.x)?u.y<v.y:u.x<v.x;
}
struct work1
{
	int fac[mxn],inv[mxn];
	inline void init()
	{
		fac[0]=inv[0]=inv[1]=1;
		fo(i,1,P) fac[i]=(ll)fac[i-1]*i%P;
		fo(i,2,P) inv[i]=(ll)(P-P/i)*inv[P%i]%P;
		fo(i,1,P) inv[i]=(ll)inv[i]*inv[i-1]%P;
	}
	inline ll C(ll n,ll m)
	{
		if(n<m) return 0;
		if(n<P && m<P)
		  return (ll)fac[n]*inv[m]%P*inv[n-m]%P;
		return C(n/P,m/P)*C(n%P,m%P)%P;
	}
}S;
struct work2
{
	int fac[5][mxn],inv[5][mxn];
	inline void init()
	{
		fo(k,1,4)
		{
			ll mod=prime[k];
			fac[k][0]=inv[k][0]=inv[k][1]=1;
			fo(i,1,mod) fac[k][i]=(ll)fac[k][i-1]*i%mod;
			fo(i,2,mod) inv[k][i]=(ll)(mod-mod/i)*inv[k][mod%i]%mod;
			fo(i,1,mod) inv[k][i]=(ll)inv[k][i]*inv[k][i-1]%mod;
		}
	}
	inline ll C(int k,ll n,ll m,ll mod)
	{
		if(n<m) return 0;
		if(n<mod && m<mod)
		  return (ll)fac[k][n]*inv[k][m]%mod*inv[k][n-m]%mod;
		return C(k,n/mod,m/mod,mod)*C(k,n%mod,m%mod,mod)%mod;
	}
	inline void exgcd(ll a,ll b,ll &x,ll &y)
	{
		if(b==0) {x=1,y=0;return;}
		exgcd(b,a%b,y,x);
		y-=(a/b)*x;
	}
	inline ll CRT(ll n,ll m)
	{
		ll t,x,y,ans=0,M=P;
		fo(k,1,4)
		{
			t=M/prime[k];
			exgcd(t,prime[k],x,y);
			ans=(ans+t*x*C(k,n,m,prime[k])%M)%M;
		}
		return (ans+M)%M;
	}
}W;
inline ll C(ll n,ll m)
{
	if(P==1000003) return S.C(n,m);
	else return W.CRT(n,m);
}
int main()
{
	scanf("%lld%lld%d%d",&n,&m,&T,&P);
	fo(i,1,T) scanf("%lld%lld",&p[i].x,&p[i].y);
	p[++T]=(point){n,m};
	sort(p+1,p+T+1,comp);
	if(P==1000003) S.init();
	else W.init();
	fo(i,1,T)
	{
		dp[i]=C(p[i].x+p[i].y,p[i].x);
		fo(j,1,i-1) if(p[j].y<=p[i].y)
		  dp[i]=(dp[i]-C(p[i].x-p[j].x+p[i].y-p[j].y,p[i].x-p[j].x)*dp[j]%P)%P;
		dp[i]=(dp[i]+P)%P;
	}
	printf("%lld\n",dp[T]);
	return 0;
}