bzoj 3782: 上学路线

这和51nod 1486 那题是差不多的，只不过这里过了一个合数的取模，拆成素数然后用中国剩余定理合并起来就好

（说的好轻巧啊）

/*#include<bits/stdc++.h>
#define LL long long
using namespace std;
inline LL ra()
{
    LL x=0; char ch=getchar();
    while (ch<'0' || ch>'9') ch=getchar();
    while (ch>='0' && ch<='9') {x=x*10+ch-'0'; ch=getchar();}
    return x;
}

LL p[]={3LL,5LL,6793LL,10007LL};

struct node{
    LL x,y;
}a[500];
bool cmp(node a, node b) {return a.x==b.x?a.y<b.y:a.x<b.x;}

const int maxn=1000005;
const int MAXN=10010;

LL fac[maxn],inv[maxn];
LL FAC[4][MAXN],INV[4][MAXN];
LL n,m,mod,T;

void pre()
{
    if (mod==1000003)
    {
        fac[0]=1;
        for (int i=1; i<maxn; i++) fac[i]=fac[i-1]*i%mod; 
        inv[0]=inv[1]=1;
        for (int i=2; i<maxn; i++) inv[i]=(mod-mod/i)*inv[mod%i]%mod;
        for (int i=2; i<maxn; i++) inv[i]=inv[i]*inv[i-1]%mod;
        return;
    }
    for (int j=0; j<4; j++)
    {
        FAC[j][0]=1;
        for (int i=1; i<MAXN; i++) FAC[j][i]=FAC[j][i-1]*i%p[j];
        INV[j][0]=INV[j][1]=1;
        for (int i=2; i<MAXN; i++) INV[j][i]=(p[j]-p[j]/i)*INV[j][p[j]%i]%p[j];
        for (int i=2; i<MAXN; i++) INV[j][i]=INV[j][i-1]*INV[j][i]%p[j];
    }
//    cout<<FAC[0][2]*INV[0][1]*INV[0][1];
}

LL lucas1(LL n, LL m)
{
    if (n<m) return 0;
    if (n<mod && m<mod) return fac[n]*inv[m]%mod*inv[n-m]%mod;
    return lucas1(n/mod,m/mod)*lucas1(n%mod,m%mod)%mod;
}

LL lucas2(LL n, LL m, int opt)
{
    if (n<m) return 0;// cout<<n<<' '<<m<<' '<<FAC[opt][n]*INV[opt][m]%p[opt]*INV[opt][n-m]%p[opt]; system("pause");
    if (n<p[opt] && m<p[opt]) return FAC[opt][n]*INV[opt][m]%p[opt]*INV[opt][n-m]%p[opt];
    return lucas2(n/p[opt],m/p[opt],opt)*lucas2(n%p[opt],m%p[opt],opt)%p[opt];
}
LL ksm(LL x, int p, int mod)
{
    LL sum=1;
    for (;p;p>>=1,x=x*x%mod)
        if (p&1) sum=sum*x%mod;
    return sum;
}
LL C(LL n, LL m)
{
//    printf("%d %d %d
",n,m,lucas1(n,m));
    if (mod==1000003) return lucas1(n,m);
    else{
        LL ans=0;
        for (int i=0; i<4; i++)
            ans+=lucas2(n,m,i)*ksm(mod/p[i],p[i]-2,p[i])%mod*(mod/p[i])%mod,ans=(ans+mod)%mod;
        return ans;
    }
}

LL f[MAXN];

int main()
{
    n=ra(); m=ra(); T=ra(); mod=ra(); pre();
    for (int i=1; i<=T; i++) a[i].x=ra(),a[i].y=ra();
    a[++T].x=n; a[T].y=m;
    sort(a+1,a+T+1,cmp);
    for (int i=1; i<=T; i++)
    {
        f[i]=C(a[i].x+a[i].y,a[i].x);// cout<<a[i].x<<"  "<<a[i].y<<"   "<<f[i]<<endl;
        for (int j=1; j<i; j++)
            if (a[j].y<=a[i].y)
            f[i]=(f[i]-C(a[i].x-a[j].x+a[i].y-a[j].y,a[i].x-a[j].x)*f[j]%mod+mod)%mod;
    }
    cout<<f[T]<<endl;
    return 0;
}*/

#include<bits/stdc++.h>
#define LL long long
using namespace std;
struct node{
    LL x,y;
}a[500];
bool cmp(node a, node b){
    return a.x==b.x?a.y<b.y:a.x<b.x;
}
LL n,m,P;  
LL p[10];  
LL fac[5][1000010],inv[5][1000010];  
LL f[1010];  
LL g[10],sum[10],INV[10];  
int num;  
bool w=0;  
LL ksm(LL x, LL y, LL p)
{
    LL sum=1;
    for (;y;y>>=1,x=x*x%P)
        if (y&1LL) sum=sum*x%P;
    return sum;
}
LL c(LL x, LL y, int i)
{
    if (x<p[i] && y<p[i]) return fac[i][y]*inv[i][x]%p[i]*inv[i][y-x]%p[i];
    else return c(x%p[i],y%p[i],i)*c(x/p[i],y/p[i],i)%p[i]; 
}
LL C(LL y, LL x)
{
    if (!w) return c(x,y,0);
    else
    {
        for (int i=1; i<=4; i++) g[i]=c(x,y,i);
        LL ans=0;
        for (int i=1; i<=4; i++) ans=(ans+(g[i]*sum[i])%P*INV[i]%P)%P;
        return ans;
    }
}
int main()
{
    scanf("%lld%lld%lld%lld",&n,&m,&num,&P);
    for (int i=1; i<=num; i++) scanf("%lld%lld",&a[i].x,&a[i].y);
    a[++num].x=n; a[num].y=m;
    sort(a+1,a+num+1,cmp);
    if (P!=1000003) p[1]=3,p[2]=5,p[3]=6793,p[4]=10007,w=1; else p[0]=1000003;
    if (w)
    {
        for (int j=1; j<=4; j++)
        {
            fac[j][0]=1; sum[j]=P/p[j]; INV[j]=ksm(sum[j],p[j]-2,p[j]);
            for (int i=1; i<p[j]; i++) fac[j][i]=fac[j][i-1]*i%p[j];
            inv[j][p[j]-1]=ksm(fac[j][p[j]-1],p[j]-2,p[j]);
            for (int i=p[j]-2; i>=0; i--) inv[j][i]=inv[j][i+1]*(i+1)%p[j];
        }
    }
    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]=ksm(fac[0][P-1],P-2,P);
        for (int i=P-2; i>=0; i--) inv[0][i]=inv[0][i+1]*(i+1)%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[i].y>=a[j].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;
    }
    cout<<f[num];
    return 0;
}