考虑 dp。
我们先把 ( n , m ) (n,m) (n,m) 也当做障碍点，然后把所有的障碍点按 x x x 坐标为第一关键字， y y y 坐标为第二关键字排序。
然后设 f i f_i fi 表示到达第 i i i 个障碍点的合法总方案数（途中不经过障碍点）。显然，答案就是 f t + 1 f_{t+1} ft+1，也就是到达 ( n , m ) (n,m) (n,m) 的总方案数。
至于为什么要先排序，是因为我们要保证当处理 f i f_i fi 时，能转移到 f i f_i fi 的所有 f j f_j fj 都已经处理完了。
显然有状态转移方程：（其中 x i x_i xi 表示第 i i i 个障碍点的 x x x 坐标， y i y_i yi 同理）
f i = ( x i + y i x i ) − ∑ j = 1 i − 1 f j × ( ( x i + y i ) − ( x j + y j ) x i − x j ) f_i=inom{x_i+y_i}{x_i}-sum_{j=1}^{i-1}f_j imesinom{(x_i+y_i)-(x_j+y_j)}{x_i-x_j} fi=(xixi+yi)−j=1∑i−1fj×(xi−xj(xi+yi)−(xj+yj))
首先， ( x i + y i x i ) dbinom{x_i+y_i}{x_i} (xixi+yi) 是从 ( 0 , 0 ) (0,0) (0,0) 到 ( x i , y i ) (x_i,y_i) (xi,yi) 的总方案（不论是否合法）。
证明：从 ( 0 , 0 ) (0,0) (0,0) 到 ( x i , y i ) (x_i,y_i) (xi,yi) 一共需要 x i + y i x_i+y_i xi+yi 步，其中需要横着走 x i x_i xi 步，竖着走 y i y_i yi 步，那么我们就枚举在哪几步横着走，也就是 ( x i + y i x i ) dbinom{x_i+y_i}{x_i} (xixi+yi)。
然后 ( ( x i + y i ) − ( x j + y j ) x i − x j ) dbinom{(x_i+y_i)-(x_j+y_j)}{x_i-x_j} (xi−xj(xi+yi)−(xj+yj)) 是从 ( x j , y j ) (x_j,y_j) (xj,yj) 到 ( x i , y i ) (x_i,y_i) (xi,yi) 的总方案（不论是否合法），证明同理。
然后证明那些不合法的方案就是 ∑ j = 1 i − 1 f j × ( ( x i + y i ) − ( x j + y j ) x i − x j ) sum_{j=1}^{i-1}f_j imesdbinom{(x_i+y_i)-(x_j+y_j)}{x_i-x_j} ∑j=1i−1fj×(xi−xj(xi+yi)−(xj+yj))。
设某一条不合法路径上的第一个障碍点是第 f i r s t first first 个障碍点。
那么 f j × ( ( x i + y i ) − ( x j + y j ) x i − x j ) f_j imesdbinom{(x_i+y_i)-(x_j+y_j)}{x_i-x_j} fj×(xi−xj(xi+yi)−(xj+yj)) 表示的就是那些 f i r s t = j first=j first=j 的不合法路径，因为这种路径从 ( 0 , 0 ) (0,0) (0,0) 到 ( x j , y j ) (x_j,y_j) (xj,yj) 是没有障碍点的，正好对应 f j f_j fj。
这么解释的话，就容易证明 ∑ j = 1 i − 1 f j × ( ( x i + y i ) − ( x j + y j ) x i − x j ) sum_{j=1}^{i-1}f_j imesdbinom{(x_i+y_i)-(x_j+y_j)}{x_i-x_j} ∑j=1i−1fj×(xi−xj(xi+yi)−(xj+yj)) 和所有的不合法路径是一一对应的了。
那么这个状态转移方程就是正确的了。
现在的问题就是 ( n m ) m o d P dbinom{n}{m}mod P (mn)modP 应该怎么求。
显然，当 P = 1000003 ∈ p r i m e P=1000003in prime P=1000003∈prime 时，这个可以用 Lucas 定理直接算。
但是当 P = 1019663265 = 3 × 5 × 6793 × 10007 P=1019663265=3 imes5 imes6793 imes10007 P=1019663265=3×5×6793×10007 时，就不能直接用 Lucas 算了。
这个需要用中国剩余定理来做，不会的请先去做一下 【模板】中国剩余定理(CRT)/曹冲养猪。
中国剩余定理的结论大概是这个东西：
对于一个方程组：
{ x ≡ a 1 ( m o d p 1 ) x ≡ a 2 ( m o d p 2 ) ⋯ x ≡ a k ( m o d p k ) egin{cases} xequiv a_1 pmod {p_1}\ xequiv a_2 pmod {p_2}\ cdots\ xequiv a_k pmod {p_k} end{cases} ⎩⎪⎪⎪⎨⎪⎪⎪⎧x≡a1(modp1)x≡a2(modp2)⋯x≡ak(modpk)
其中 p 1 , p 2 , … , p k p_1,p_2,dots,p_k p1,p2,…,pk 两两互质。
设 M = ∏ i = 1 k p i M=prodlimits_{i=1}^{k}p_i M=i=1∏kpi， m i = M p i m_i=dfrac{M}{p_i} mi=piM， t i t_i ti 是 m i m_i mi 模 p i p_i pi 意义下的逆元。那么 x x x 的一个特殊解就是 x 0 = ∑ i = 1 k a i m i t i x_0=sumlimits_{i=1}^{k}a_im_it_i x0=i=1∑kaimiti。
然后又因为 x x x 一定能表示成 a × M + b a imes M+b a×M+b 的形式，所以 x x x 的最小非负数解就是 x 0 m o d M x_0mod M x0modM。
这个题中，未知数 x x x 就是 ( n m ) dbinom{n}{m} (mn)，我们要求 x m o d P xmod P xmodP，然后我们设 p 1 = 3 p_1=3 p1=3、 p 2 = 5 p_2=5 p2=5、 p 3 = 6793 p_3=6793 p3=6793、 p 4 = 10007 p_4=10007 p4=10007，那么 M M M 就是这题中的 P P P。
然后我们还知道 a 1 = x m o d p 1 a_1=xmod p_1 a1=xmodp1、 a 2 = x m o d p 2 a_2=xmod p_2 a2=xmodp2、 a 3 = x m o d p 3 a_3=xmod p_3 a3=xmodp3、 a 4 = x m o d p 4 a_4=xmod p_4 a4=xmodp4 的值（这些都可以用 Lucas 定理求出来）。
那么我们通过中国剩余定理就可以把 x m o d P xmod P xmodP 算出来了。
具体代码如下：
#include<bits/stdc++.h>
#define T 210
#define ll long long
using namespace std;
struct Point
{
int x,y;
Point(){};
Point(int a,int b){x=a,y=b;}
}q[T];
int n,m,num,P;
ll f[T];
ll poww(ll a,ll b,ll p)
{
ll ans=1;
while(b)
{
if(b&1) ans=ans*a%p;
a=a*a%p;
b>>=1;
}
return ans;
}
namespace mod1//p=1019663265的情况
{
const ll M=1019663265ll,p[5]={0,3,5,6793,10007};
ll m[5],t[5],fac[5][10010],inv[5][10010];
ll a[5];
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*y;
}
void init()
{
for(int i=1;i<=4;i++)
{
m[i]=M/p[i];
ll x,y;
exgcd(m[i],p[i],x,y);//通过exgcd求逆元
t[i]=(x%p[i]+p[i])%p[i];
}
for(int i=1;i<=4;i++)
{
fac[i][0]=1;
for(int j=1;j<p[i];j++)
fac[i][j]=fac[i][j-1]*j%p[i];
}
for(int i=1;i<=4;i++)
for(int j=0;j<p[i];j++)
inv[i][j]=poww(fac[i][j],p[i]-2,p[i]);
}
ll C(ll n,ll m,int k)
{
if(n<m) return 0;
if(!m) return 1;
return fac[k][n]*inv[k][n-m]%p[k]*inv[k][m]%p[k];
}
ll Lucas(ll n,ll m,int k)
{
if(n<m) return 0;
if(!m) return 1;
return C(n%p[k],m%p[k],k)*Lucas(n/p[k],m/p[k],k)%p[k];
}
ll solve(ll n,ll mm)
{
ll x=0;
for(int i=1;i<=4;i++)
{
a[i]=Lucas(n,mm,i);
x=(x+a[i]*m[i]%M*t[i]%M)%M;//求x的解
}
return x;
}
}
namespace mod2//p=1000003的情况
{
const ll p=1000003;
ll fac[1000005],inv[1000005];
void init()
{
fac[0]=1;
for(int i=1;i<p;i++)
fac[i]=fac[i-1]*i%p;
for(int i=0;i<p;i++)
inv[i]=poww(fac[i],p-2,p);
}
ll C(ll n,ll m)
{
if(n<m) return 0;
if(!m) return 1;
return fac[n]*inv[n-m]%p*inv[m]%p;
}
ll Lucas(ll n,ll m)
{
if(n<m) return 0;
if(!m) return 1;
return C(n%p,m%p)*Lucas(n/p,m/p)%p;
}
}
bool cmp(Point a,Point b)
{
if(a.x==b.x) return a.y<b.y;
return a.x<b.x;
}
int main()
{
scanf("%d%d%d%d",&n,&m,&num,&P);
if(P==1000003) mod2::init();
else mod1::init();
for(int i=1;i<=num;i++)
scanf("%d%d",&q[i].x,&q[i].y);
q[num+1]=Point(n,m);
sort(q+1,q+num+2,cmp);
for(int i=1;i<=num+1;i++)
{
if(P==1000003) f[i]=mod2::Lucas(q[i].x+q[i].y,q[i].x);
else f[i]=mod1::solve(q[i].x+q[i].y,q[i].x);
for(int j=1;j<i;j++)
{
if(q[j].x<=q[i].x&&q[j].y<=q[i].y)
{
if(P==1000003) f[i]=(f[i]-f[j]*(mod2::Lucas(q[i].x+q[i].y-q[j].x-q[j].y,q[i].x-q[j].x))%P+P)%P;
else f[i]=(f[i]-f[j]*(mod1::solve(q[i].x+q[i].y-q[j].x-q[j].y,q[i].x-q[j].x))%P+P)%P;
}
}
}
printf("%lld
",f[num+1]);
return 0;
}