bzoj 1266 [AHOI2006] 上学路线

//Achen
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<vector>
#include<queue>
#include<ctime>
#include<cmath>
const int N=207,M=1019663265;
typedef long long LL;
using namespace std;
LL n,m,p,f[5][N],pr[5][1000007],inv[5][1000007],b[10],T,mi[10]={0,3,5,6793,10007};

struct node {
    LL x,y;
    friend bool operator <(const node&A,const node&B) {
        return A.x<B.x||(A.x==B.x&&A.y<B.y);
    }
}pp[N];

template<typename T> void read(T &x) {
    char ch=getchar(); x=0; T f=1;
    while(ch!=‘-‘&&(ch<‘0‘||ch>‘9‘)) ch=getchar();
    if(ch==‘-‘) f=-1,ch=getchar();
    for(;ch>=‘0‘&&ch<=‘9‘;ch=getchar()) x=x*10+ch-‘0‘; x*=f;
}

LL lucas(LL n,LL m,int k) {
    if(m>n) return 0;
    if(n<mi[k]) return pr[k][n]*inv[k][m]%mi[k]*inv[k][n-m]%mi[k];
    return lucas(n%mi[k],m%mi[k],k)*lucas(n/mi[k],m/mi[k],k)%mi[k];
}

void exgcd(LL a,LL b,LL &x,LL &y) {
    if(!b) { x=1; y=0; return; }
    exgcd(b,a%b,y,x); y-=a/b*x;
}

int main() {
    read(n); read(m); read(T); read(p);
    for(int i=1;i<=T;i++) {
        read(pp[i].x); read(pp[i].y);
    }
    pp[T+1].x=n; pp[T+1].y=m; n=T+1;
    sort(pp+1,pp+n+1);
    if(p==1000003) {
        pr[1][0]=1; mi[1]=p;
        for(int i=1;i<=mi[1];i++) pr[1][i]=pr[1][i-1]*i%p;
        inv[1][0]=inv[1][1]=1;
        for(int i=2;i<=mi[1];i++) inv[1][i]=(p-p/i*inv[1][p%i]%p)%p;
        for(int i=2;i<=mi[1];i++) (inv[1][i]*=inv[1][i-1])%=p;
    }
    else for(int k=1;k<=4;k++) {
        pr[k][0]=1; inv[k][0]=inv[k][1]=1;
        for(int i=1;i<=mi[k];i++) pr[k][i]=pr[k][i-1]*i%mi[k];
        for(int i=2;i<=mi[k];i++) inv[k][i]=(mi[k]-mi[k]/i*inv[k][mi[k]%i]%mi[k])%mi[k];
        for(int i=2;i<=mi[k];i++) (inv[k][i]*=inv[k][i-1])%=mi[k];
    }
    for(int k=1;k<=4;k++) {
        for(int i=1;i<=n;i++) {
            f[k][i]=lucas(pp[i].x+pp[i].y,pp[i].x,k);
            for(int j=1;j<i;j++) if(pp[j].x<=pp[i].x&&pp[j].y<=pp[i].y) {
                f[k][i]-=f[k][j]*lucas(pp[i].x-pp[j].x+pp[i].y-pp[j].y,pp[i].x-pp[j].x,k)%mi[k];
                if(f[k][i]<0) f[k][i]+=mi[k];
                if(f[k][i]>=mi[k]) f[k][i]%=mi[k];
            }
        }
        if(p==1000003) {
            printf("%lld\n",f[1][n]);
            return 0;
        }
    }
    for(int i=1;i<=4;i++) 
        b[i]=f[i][n];
    LL res=0;
    for(int i=1;i<=4;i++) {
        LL tp=M/mi[i],x,y;
        exgcd(mi[i],tp,x,y);
        if(y<0) y+=M; if(y>M) y%=M;
        res=(res+y*tp%M*b[i]%M)%M;
    }
    printf("%lld\n",res);
    return 0;
}