Description：

一个网格图，有一些点不能走，只能往下或往右走，问从 $(0,0)$ 到 $(n,m)$ 的方案数。

$1\leqslant n,m\leqslant 10^9,mod=1000003/1029663265$

Solution：

首先网格图上只向两个方向走从 $(0,0)$ 到 $(n,m)$ 的方案数是 $C(n+m,n)$ ，那么我们可以先求出所有方案数，再减掉不合法的方案数，首先经过这个店的方案数可以由上述公式前后合并一下求得，但是这样肯定会算重，即有一条路径经过了多于 $1$ 个坏点，那么我们把所有这样的路径按照第一个经过的坏点分类，对于每一个点，枚举所有横纵坐标都小于等于这个点的，然后再计算一下先到那个点，然后到这个点的方案数，即 $f[j]\times C(s[i].x-s[j].x+s[i].y-s[j].y,s[i].x-s[j].x)$ ，把这些减掉即可，由于 $f[i]$ 表示的就是 $i$ 作为第一个坏点的方案数，所以这样一定是不重不漏的。

当 $mod=1000003$ 时，可以直接做， $1029663265=3\times 5\times6793\times10007$ ，可以拆开后用卢卡斯定理计算组合数最后用中国剩余定理合并。

中国剩余定理：

方程组： $$ x\mod m_1=a_1\\ x\mod m_2=a_2\\ \vdots\\ x\mod m_n=a_n\\ $$ 的通解为： $$ ans=kM+\sum_{i=1}^na_it_iM_i $$ 其中： $$ M=\prod m_i\\ M_i=M/m_i\\ t_i= M_i^{-1}(\text{mod }m_i)\\ $$

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

typedef long long ll;

ll n,m;

int t,p;

#define MAXN 210

struct block

{

	int x,y;

}s[MAXN];

bool cmp(block a,block b){return (a.x == b.x ? a.y < b.y : a.x < b.x);}

int mod;

int power(int a,int b,int mod)

{

	int res = 1;

	while(b > 0)

	{

		if(b & 1)res = 1ll * res * a % mod;

		a = 1ll * a * a % mod;

		b = b >> 1;

	}

	return res;

}

int d[5] = {0,3,5,6793,10007};

int fac[1000010],inv[1000010];

int C(int n,int m,int mod){return 1ll * fac[n] * inv[m] % mod * inv[n - m] % mod;}

int lucas(ll n,ll m,int mod)

{

	if(m > n)return 0;

	if(m == 0)return 1;

	if(n < mod)return C(n,m,mod);

	return 1ll * lucas(n / mod,m / mod,mod) * C(n % mod,m % mod,mod) % mod;

}

int cnt(ll n,ll m,int mod)

{

	return lucas(n + m,m,mod);

}

int f[MAXN];

int calc(int mod)

{

	fac[0] = 1;

	for(int i = 1;i < mod;++i)fac[i] = 1ll * fac[i - 1] * i % mod;

	inv[mod - 1] = power(fac[mod - 1],mod - 2,mod);inv[0] = 1;

	for(int i = mod - 2;i >= 1;--i)inv[i] = 1ll * (i + 1) * inv[i + 1] % mod;

	for(int i = 1;i <= t;++i)

	{

		f[i] = cnt(s[i].x,s[i].y,mod);

		for(int j = 1;j < i;++j)

		{

			if(s[j].x <= s[i].x && s[j].y <= s[i].y)

			{

				f[i] = (f[i] - 1ll * cnt(s[i].x - s[j].x,s[i].y - s[j].y,mod) * f[j] % mod + mod) % mod;

			}

		}

	}

	int ans = cnt(n,m,mod);

	for(int i = 1;i <= t;++i)

	{

		ans = (ans - 1ll * f[i] * cnt(n - s[i].x,m - s[i].y,mod) % mod + mod) % mod;

	}

	return ans;

}

int main()

{

	scanf("%lld%lld%d%d",&n,&m,&t,&p);

	for(int i = 1;i <= t;++i)

	{

		scanf("%d%d",&s[i].x,&s[i].y);

	}

	sort(s + 1,s + 1 + t,cmp);

	if(p == 1000003)

	{

		cout << calc(1000003) << endl;

	}

	else

	{

		int ans = 0;

		for(int i = 1;i <= 4;++i)

		{

			ans = (ans + 1ll * calc(d[i]) * power(p / d[i],d[i] - 2,d[i]) % p * (p / d[i]) % p) % p;

		}

		cout << ans << endl;

	}

	return 0;

}