题目描述

小B 所在的城市的道路构成了一个方形网格，它的西南角为(0,0)，东北角为(N,M)。

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

输入格式

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

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

输出格式

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

输入输出样例

输入

3 4 3 1019663265
3 0
1 1
2 2

输出

8

Solution

先来看看数据范围，就发现可以骗到分。

\(\text{40pts:}\)

测试点1、2：\(n,m≤1000\)，直接 \(O(nm)DP\)。

测试点3、4：没有施工路口，直接 \(C_{n+m}^n\) 求总方案数，然后因为 \(P\) 是质数，所以用逆元取模即可。

\(\text{60pts:}\)

测试点5、6：因为 \(P\) 不是质数，所以把 \(P\) 分解成几个质数的乘积，分别算出 \(C_{n+m}^n\) 对每个质数取模的结果，然后用 中国剩余定理\((CRT)\) 合并求解（中国剩余定理下面会讲）。

以上就是考场上的骗分思路。

---

\(\text{100pts:}\)

把每个施工路口的坐标存进 \(pos\) 结构体数组，再把 \((n,m)\) 存进 \(pos_{t+1},t++\)，然后将 \(pos\) 按 \(x\) 从小到大排序，如果 \(x\) 相等，就按 \(y\) 从小到大排序。

用 \(f_i\) 表示从 \((0,0)\) 到 \((pos_i.x,pos_i.y)\) 并且不经过其他施工路口的方案数，然后大力推式子：

\[f_i=C_{pos_i.x+pos_i.y}^{pos_i.x}-\sum_{j=1}^{t}f_jC_{pos_i.x-pos_j.x+pos_i.y-pos_j.y}^{pos_i.x-pos_j.x}\text{(符合要求的j)} \]

其中，

第一项表示从 \((0,0)\) 到 \((pos_i.x,pos_i.y)\) 的总方案数。

第二项表示所有所有符合要求的 \(j\) 从 \((0,0)\) 到 \((pos_j.x,pos_j.y)\) 且不经过其他施工路口的方案数 × 从\((pos_j.x,pos_j.y)\) 到 \((pos_i.x,pos_i.y)\) 的方案数的和。

显然，\(j\) 的要求是 \(pos_j.x<=pos_i.x\&\&pos_j.y<=pos_i.y\)。

由于已经将 \(pos\) 数组排了序，所以式子中 \(j\) 的范围可以缩小到 \(i-1\) ：

\[f_i=C_{pos_i.x+pos_i.y}^{pos_i.x}-\sum_{j=1}^{i-1}f_jC_{pos_i.x-pos_j.x+pos_i.y-pos_j.y}^{pos_i.x-pos_j.x}(pos_j.x<=pos_i.x\&\&pos_j.y<=pos_i.y) \]

举个例子（样例）：

已知 \(f_1=2\) 求 \(f_2:\)

\(f_2=C_4^2-f_1×C_2^1\)

\(\quad=6-2×2\)

\(\quad=2\)

其余同理。

最后输出 \(f_t\) 即可。

至此，本题的主体思路已经讲完了。

---

接下来考虑范围。

因为 \(P\) 不一定是质数，所以不能直接用逆元取模。

先考虑 \(P\) 为质数的时候：

直接用 \(Lucas\) 即可。

简单介绍一下：

\(C_n^m=C_{n\bmod{p}}^{m\bmod{p}}×C_{n/p}^{m/p}\bmod{p}\)

如果 \(P\) 不是质数：

如何处理上面已经提到了，就是把 \(P\) 分解成几个质数的积，存进 \(p\) 数组，再用 \(CRT\) 合并。

先介绍一下 \(CRT:\)

设 \(m_1,m_2,...,m_n\) 是两两互质的整数，\(m=∏_{i=1}^nm_i,M_i=m/m_i,t_i是线性方程M_it_i≡1(mod\ m_i)\) 的一个解。对于任意的 \(n\) 个整数 \(a_1,a_2,...,a_n,\) 方程组

\(\begin{cases} x≡a_1\ (\bmod{m_1})&\\ x≡a_2\ (\bmod{m_2})&\\ ...\ &\\ x≡a_n\ (\bmod\ m_n) \end{cases}\)

有整数解，解为 \(x=\sum_{i=1}^na_iM_it_i\)。

证明请左转自行百度。

已知 \(P=1019663265=3×5×6793×10007\)

让 \(p_1=3,p_2=5,p_3=6793,p_4=10007\)，显然它们两两互质，符合 \(CRT\) 中的 \(m\) 数组。

\(P=M\)

\(mul_{ij}=M_j\bmod{p_i}\)

\(invm_{ij}=t_j\bmod{p_i}=mul_{ij}\bmod{p_i}\)的逆元(因为 \(p\) 为质数，所以可以用逆元)

\(g_i=x\bmod{p_i}\)

原式 \(:x≡invm_{ij}\ (\bmod\ p_i)\)

\(x=\sum_{i=1}^n(g_i×mul_{ij}\bmod{p_i}×invm_{ij}\bmod{p_i})\)

还要预处理一下阶乘和阶乘的逆元，在求组合数时需要用到。

\(fac_{ij}=j!\bmod{p_i}\)

\(inv_{ij}=fac_{ij}\bmod{p_i}\) 的逆元

$C_n^m=\frac{n!}{m!(n-m)!}\bmod{p_i}=fac_n×inv_m\bmod{p_i}×inv_{n-m}\bmod{p_i}\ $(这里就需要记录下来 \(i\))

对每个 \(p_i\) 都要取模，就可以用 \(CRT\)。

\(x=C_n^m\)

原方程组 \(=>\)

\(\begin{cases} C_n^m≡g_1\ (\bmod\ p_1)&\\ C_n^m≡g_2\ (\bmod\ p_2)&\\ C_n^m≡g_3\ (\bmod\ p_3)\ &\\ C_n^m≡g_4\ (\bmod\ p_4) \end{cases}\)

这样就可以求出 \(C_n^m\bmod{P}\) 了。

然后再用 \(DP\) 转移出最后的结果。

时间复杂度 ：\(O(\text{能过})\) 开玩笑的

应该是 \(O(T^2logn)\)

Code

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#define int long long

using namespace std;

const int N=1e6+7;
struct node
{
	int x,y;
}pos[210];		//每个施工路口的坐标 
int n,m,t,P,tp;
int f[210],g[5],p[5],mul[5],fac[5][N],invm[5],inv[5][N];

inline int read()	//快读 
{
	int x=0;
	char c=getchar();
	while(c<'0'||c>'9') c=getchar();
	while(c>='0'&&c<='9') x=x*10+c-'0',c=getchar();
	return x;
}

inline bool cmp(node a,node b)	//排序 
{
	return a.x==b.x?a.y<b.y:a.x<b.x;
}

inline int power(int a,int b,int p)		//快速幂a^b%p 
{
	int ret=1;
	a%=p;
	while(b)
	{
		if(b&1) ret=ret*a%p;
		a=a*a%p;
		b>>=1;
	}
	return ret;
}

inline int C(int a,int b,int i)		//组合数
{
    if(a<b) return 0;
    if(!b||a==b) return 1;
    if(a<p[i]&&b<p[i]) return fac[i][a]*inv[i][b]%p[i]*inv[i][a-b]%p[i];
    return C(a%p[i],b%p[i],i)*C(a/p[i],b/p[i],i)%p[i];	//卢卡斯定理 
}

inline int CRT(int a,int b)		//中国剩余定理 
{
	if(!tp) return C(a,b,0);	//P为质数时直接返回C(a,b)%p[0] 
	int ret=0;
	for(int i=1; i<=4; i++) g[i]=C(a,b,i);	//g[i]=C(a,b)%p[i] 
	for(int i=1; i<=4; i++) ret=(ret+g[i]*mul[i]%P*invm[i]%P)%P;
	return ret;
}

signed main()
{
	n=read(),m=read(),t=read(),P=read();
	for(int i=1; i<=t; i++)
		pos[i].x=read(),pos[i].y=read();
	pos[++t]=(node){n,m};		//将(n,m)加到pos[t+1] 
	sort(pos+1,pos+1+t,cmp);	//排序 
	if(P==1e6+3) p[0]=P;		//如果P是质数，存到p[0] 
	else  p[1]=3,p[2]=5,p[3]=6793,p[4]=10007,tp=1;	//否则，分解成4个质数，并且tp=1,表示P不是质数 
	
	//预处理 
	if(tp)	//如果P不是质数 
	{
		for(int i=1; i<=4; i++)
		{
			mul[i]=P/p[i];		//CRT中的m数组 
			invm[i]=power(mul[i],p[i]-2,p[i]);	//mul[i]%p[i]的逆元 
			fac[i][0]=1;		//0!%p[i] 
			for(int j=1; j<p[i]; j++)
				fac[i][j]=fac[i][j-1]*j%p[i];	//j!%p[i] 
			inv[i][p[i]-1]=power(fac[i][p[i]-1],p[i]-2,p[i]);	//(p[i]-1)!%p[i]的逆元 
			for(int j=p[i]-1; j>=1; j--)
				inv[i][j-1]=inv[i][j]*j%p[i];	//逆元的递推 
		}
	}
	else
	{
		//同上，只是少一层循环，因为只有一个质数 
		fac[0][0]=1;
		for(int i=1; i<P; i++)
			fac[0][i]=fac[0][i-1]*i%P;
		inv[0][P-1]=power(fac[0][P-1],P-2,P);
		for(int i=P-1; i>=1; i--)
			inv[0][i-1]=inv[0][i]*i%P;
	}
	for(int i=1; i<=t; i++)
	{
		f[i]=CRT(pos[i].x+pos[i].y,pos[i].x);	//从(0,0)到(pos[i].x,pos[i].y)的总方案数 
		//减去经过其他施工路口的方案数 
		for(int j=1; j<i; j++)
			if(pos[j].x<=pos[i].x&&pos[j].y<=pos[i].y)
				f[i]=(f[i]-f[j]*CRT(pos[i].x-pos[j].x+pos[i].y-pos[j].y,pos[i].x-pos[j].x)%P+P)%P;
	}
	printf("%lld\n",f[t]);
	return 0;
}