P4478 [BJWC2018] 上学路线

P4478 [BJWC2018] 上学路线

思路

看到有 \(T\) 个障碍点并且障碍点非常少，所以想到暴力容斥，用总方案减去不合法的方案数。

只要经过任意一个障碍点，就是一个不合法方案。所以令 \(dp_i\) 表示只经过障碍点 \(i\) 的方案数。

设当前障碍点为 \((x_i,y_i)\)，如果不管途中是否经过其他障碍点，那么走到该点的总方案数为 \(C_{x_i+y_i}^{x_i}\)。
但是其中还包含了经过其他障碍点的方案数。

设还经过了障碍点 \((x_j,y_j)\)，满足 \(x_j \le x_i\) 且 \(y_j \le y_i\)，那么多统计的方案数为 \(dp_j \times C_{x_i-x_j+y_i-y_j}^{x_i-x_j}\) 减掉即可。
一个点能从其左下角的每个点转移过来，所以先把点按从左到右、从上到下排好序，才能保证转移是正确的。

我们发现 \(n\) 和 \(m\) 都非常大，所以要用 Lucas 定理，\(C_n^m \equiv C_{n\%P}^{m\%P} \times C_{n/P}^{m/P} \pmod P\)。
然而 Lucas 定理只适用于 \(P\) 为质数，如果 \(P\) 不是质数怎么办呢。
我们可以把 \(P\) 分解质因数，对于每个质因子都列出一个方程，得到同余方程组，然后用中国剩余定理合并求解就可以了。

我们发现 \(P = 1019663265=3 \times 5 \times 6793 \times 10007\)，所以不用繁琐的对组合数拆开取模，直接做就可以了。

代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<climits>
#include<cstdlib>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
template <typename T>
inline void read(T&x){// 快读
	int w = 0;x = 0;
	char ch = getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') w = 1;
		ch = getchar();
	}
	while(ch>='0' && ch<='9'){
		x = (x<<1)+(x<<3)+(ch^48);
		ch = getchar();
	}
	if(w) x = ~x+1;
}
template <typename T,typename...Args>
inline void read(T&t,Args&...args){
	read(t);read(args...);
}
template <typename T>
inline T Min(T x,T y){ return (x < y ? x : y); }
template <typename T>
inline T Max(T x,T y){ return (x > y ? x : y); }
template <typename T>
inline T Abs(T x){ return (x < 0 ? ~x+1 : x); }
const int N = 210,M = 1e6+10;
ll prime[4] = {3,5,6793,10007};
ll T,ans[N],P; int cnt;
struct node{
    ll x,y;
    inline int operator < (const node&G) const{
        if(x^G.x) return x < G.x;
        return y < G.y;
    }
}g[N];
ll mul[4][M],inv[4][M];
inline ll quick_pow(ll x,ll y,int t){
    ll res = 1;
    while(y){
        if(y&1) (res *= x)%=prime[t];
        (x *= x)%=prime[t];
        y >>= 1;
    }
    return res;
}
inline ll C(ll n,ll m,int t){// 组合数
    if(m>n) return 0;
    return mul[t][n]*inv[t][n-m]%prime[t]*inv[t][m]%prime[t];
}
inline ll Lucas(ll n,ll m,int t){// Lucas
    if(!m) return 1;
    return C(n%prime[t],m%prime[t],t)*Lucas(n/prime[t],m/prime[t],t)%prime[t];
}
inline ll Gcd(ll x,ll y){
    if(!x || !y) return x|y;
    ll xz = __builtin_ctzll(x);
    ll yz = __builtin_ctzll(y);
    ll z = Min(xz,yz);
    ll diff;
    y >>= yz;
    while(x){
        x >>= xz;
        diff = x-y;
        xz = __builtin_ctzll(diff);
        y = Min(x,y);
        x = Abs(diff);
    }
    return y << z;
}
inline void Exgcd(ll x,ll y,ll &Ex,ll &Ey){// 扩欧
    if(!y){
        Ex = 1;
        Ey = 0;
        return ;
    }
    Exgcd(y,x%y,Ey,Ex);
    Ey = Ey-x/y*Ex;
}
inline ll ExCRT(ll n,ll m){// ExCRT 合并方程
    ll r1 = Lucas(n,m,0);
    ll m1 = prime[0];
    ll r2,m2,Ex,Ey;
    for(int t=1;t<=cnt;++t){
        r2 = Lucas(n,m,t);
        m2 = prime[t];
        ll d = Gcd(m1,m2);
        // if((r2-r1)%d) exit(0);
        Exgcd(m1/d,m2/d,Ex,Ey);

        Ex = ((r2-r1)/d*Ex%(m2/d)+(m2/d))%(m2/d);
        r1 = ((r1+m1*Ex)%(m1/d*m2)+(m1/d*m2))%(m1/d*m2);
        m1 = m1/d*m2;
    }
    return r1;
}
int main(){
	// freopen("in.in","r",stdin);
    // freopen("out.out","w",stdout);

	ll n,m;read(n,m,T,P);
    for(int i=1;i<=T;++i) read(g[i].x,g[i].y);
    g[++T] = {n,m};
    sort(g+1,g+1+T);// 排好序
    if(P==1000003) prime[0] = 1000003,cnt = 0;
    else cnt = 3;
    for(int t=0;t<=cnt;++t){
        mul[t][0] = 1;
        for(ll i=1;i<prime[t];++i) mul[t][i] = mul[t][i-1]*i%prime[t];
        inv[t][prime[t]-1] = quick_pow(mul[t][prime[t]-1],prime[t]-2,t);
        for(ll i=prime[t]-1;i;--i) inv[t][i-1] = inv[t][i]*i%prime[t];
    }
    for(int i=1;i<=T;++i){
        ans[i] = ExCRT(g[i].x+g[i].y,g[i].x);
        for(int j=1;j<i;++j){
            if(g[j].x<=g[i].x && g[j].y<=g[i].y){// 左下角的点
                ans[i] = (ans[i]-ans[j]*ExCRT(g[i].x+g[i].y-g[j].x-g[j].y,g[i].x-g[j].x)%P+P)%P;
            }
        }
    }
    printf("%lld",ans[T]);

	fclose(stdin);
	fclose(stdout);
	return 0;
}