bzoj3782上学路线（Lucas+CRT+容斥DP+组合计数）

传送门：https://www.lydsy.com/JudgeOnline/problem.php?id=3782

有部分分的传送门：https://www.luogu.org/problemnew/show/P4478

看到标题开始还以为是AHOI的小雪和小可可……

题解：乍一看会40pts：测试点1、2：n,m<=1000的直接O(nm)DP；测试点3、4：没有障碍物直接C(n+m,n)，然后p=1e6+3是质数可以直接取模。

想了几分钟会60pts：测试点5、6：模数可以拆成几个不超过1e5的质数的乘积，直接算出C(n+m,n)对每个质数的模数，然后CRT合并一下就行了。

不会CRT的左转，我原来也是看这个博客学的：https://blog.csdn.net/niiick/article/details/80229217

其实满分也很可做，容斥一下就行了：把障碍物按x从小到大，x相同按y从小到大排序，然后f[i]表示不经过前(i-1)个障碍物但经过第i个障碍物的方案，然后增加最后一个点为(n,m)，然后可以计算f[i]=C(x[i]+y[i],y[i])+Σf[j]C(x[i]-x[j]+y[i]-y[j],x[i]-x[j])，其中j满足x[j]<=x[i]&&y[j]<=y[i]，这个计算由于T<=200算组合数+CRT合并也不会超时，复杂度O(T^2log)。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e6+7;
struct node{ll x,y;}a[N];
ll n,m,P,f[N],p[8],fac[6][N],inv[6][N],mul[8],imul[8],g[8];
int num,tp;
bool cmp(node a,node b){return a.x==b.x?a.y<b.y:a.x<b.x;}
ll qpow(ll a,ll b,ll p)
{
    ll ret=1;
    while(b)
    {
        if(b&1)ret=ret*a%p;
        a=a*a%p,b>>=1;
    }
    return ret;
}
ll c(ll a,ll b,int i)
{
    if(a<b)return 0;
    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];
}
ll C(ll a,ll b)
{
    if(!tp)return c(a,b,0);
    ll ret=0;
    for(int i=1;i<=4;i++)g[i]=c(a,b,i);
    for(int i=1;i<=4;i++)ret=(ret+g[i]*mul[i]%P*imul[i]%P)%P;
    return ret;
}
int main()
{
    scanf("%lld%lld%d%lld",&n,&m,&num,&P);
    for(int i=1;i<=num;i++)scanf("%lld%lld",&a[i].x,&a[i].y);
    a[++num]=(node){n,m};
    sort(a+1,a+num+1,cmp);
    if(P==1e6+3)p[0]=1e6+3;else p[1]=3,p[2]=5,p[3]=6793,p[4]=10007,tp=1;
    if(tp)
    {
        for(int i=1;i<=4;i++)
        {
            mul[i]=P/p[i],imul[i]=qpow(mul[i],p[i]-2,p[i]);
            fac[i][0]=1;for(int j=1;j<p[i];j++)fac[i][j]=fac[i][j-1]*j%p[i];
            inv[i][p[i]-1]=qpow(fac[i][p[i]-1],p[i]-2,p[i]);
            for(int j=p[i]-1;j;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]=qpow(fac[0][P-1],P-2,P);
        for(int i=P-1;i;i--)inv[0][i-1]=inv[0][i]*i%P;
    }
    for(int i=1;i<=num;i++)
    {
        f[i]=C(a[i].x+a[i].y,a[i].x);
        for(int j=1;j<i;j++)if(a[j].x<=a[i].x&&a[j].y<=a[i].y)
        f[i]=(f[i]-f[j]*C(a[i].x-a[j].x+a[i].y-a[j].y,a[i].x-a[j].x)%P+P)%P;
    }
    printf("%lld",f[num]);
}