bzoj3782上学路线

题意:从n*m网格图的左下角走到右上角(n,m<=10^10),有t个坐标不能经过(t<=200),只能向上向右走,问有多少种不同的走法,对p取模,

p只有两种取值,1000003(质数)和1019663265(四个质数的乘积, 3*5*6793*10007)

考试的时候有部分分…

1. n,m<=1000时,O(nm)DP即可

2.t=0的时候求一个组合数就可以了,1000003:预处理阶乘及逆元,l1019663265:lucas定理求组合数之后用CRT合并.

那么考虑t=1的情况,我们只需要把总的路径条数减去经过那个障碍点的路径条数就可以了.走法=”左下角到障碍点的走法”*”障碍点到右上角的做法”

t=2时,设两个障碍点为A,B,”总的路径条数”-“经过A的路径条数”-“经过B的路径条数”算出来的答案可能偏小,如果A,B可以同时经过,那么最终答案要加上”同时经过A,B的路径条数”

那么这道题就可以用容斥来做.方便起见给障碍点从左到右从下到上排个序,记f[i][j]表示走到了第i个障碍点且包括第i个点在内经过了j个障碍点的路径条数,枚举从左下角的哪一个点k走过来即可.转移的时候乘上一个组合数表示从k到i的走法数目.

注意循环的时候在最内层枚举j只需要O(n^2)次计算组合数,最内层枚举k需要计算O(n^3)次组合数,效率相差很大.我一开始在别的地方卡了半天常数终于在bz上卡进总时限了,但加了这个优化瞬间从10s+变成2s以内…代码里有很多卡常的痕迹,其实我在考试的时候连wys那个鬼畜的循环展开都用了然而组合数计算次数太多了还是没卡过去...

#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long ll;
ll n,m;int t,mod;
int QuickPow(int a,int x){
    int ans=1;
    for(;x;x>>=1,a=a*1LL*a%mod){
        if(x&1)ans=ans*1LL*a%mod;
    }
    return ans;
}
namespace Prime{
    const int maxn=500003;
    int fac[maxn],inv[maxn];
    void init(){
        fac[0]=1;
        for(int i=1;i<maxn;++i)fac[i]=fac[i-1]*1LL*i%mod;
        inv[maxn-1]=QuickPow(fac[maxn-1],mod-2);
        for(int i=maxn-1;i>=1;--i){
            inv[i-1]=inv[i]*1LL*i%mod;
        }
    }
    int C(ll n,ll m){
        if(n<m)return 0;
        return fac[n]*1LL*inv[m]%mod*1LL*inv[n-m]%mod;
    }
    int Lucas(ll n,ll m){//printf("%lld %lld
",n,m);
        if(n<maxn)return C(n,m);
        int t1=n/maxn,t2=m/maxn;
        return Lucas(t1,t2)*1LL*C(n-maxn*1LL*t1,m-maxn*1LL*t2)%maxn;
    }
};
namespace NotPrime{
    const int maxn=12000,Mod=1019663265;
    int p[5]={3,5,6793,10007,1019663265};
    int fac[4][maxn],inv[4][maxn];
    int BUF3[4][4]={{1,0,0,0},{1,1,0,0},{1,2,1,0},{1,0,0,1}};
    int BUF5[6][6]={{1,0,0,0,0,0},{1,1,0,0,0,0},{1,2,1,0,0,0},{1,3,3,1,0,0},{1,4,1,4,1,0}};
    void init(){
        for(int i=0;i<4;++i){
            fac[i][0]=1;
            for(int j=1;j<p[i];++j)fac[i][j]=fac[i][j-1]*1LL*j%p[i];
            inv[i][p[i]-1]=QuickPow(fac[i][p[i]-1],p[i]-2);
            for(int j=p[i]-1;j>=1;--j){
                inv[i][j-1]=inv[i][j]*1LL*j%p[i];
            }
        }
    }
    int M;
    int lucas(ll n,ll m,int t){
        if(n<M){
            if(m>n)return 0;
            if(!t)return BUF3[n][m];
            if(t==1)return BUF5[n][m];
            return fac[t][n]*inv[t][m]*1LL*inv[t][n-m]%M;
        }
        int t1=n/M,t2=m/M;
        return lucas(t1,t2,t)*lucas(n-t1*1LL*M,m-t2*1LL*M,t)%M;
    }
    int res[4];
    int CRT(){
        int ans=0,lcm,Inv;
        for(int i=0;i<4;++i){
            lcm=p[4]/p[i];
            Inv=QuickPow(lcm%p[i],p[i]-2);
            ans=(ans+res[i]*1LL*lcm%mod*1LL*Inv%mod)%p[4];
        }
        return ans;
    }
};
namespace Ending{
    int f[205][205];
    struct point{
        ll x,y,sum;
        void read(){scanf("%lld%lld",&x,&y);sum=x+y;}
        bool operator <(const point &B)const{
            return (x==B.x)?(y<B.y):x<B.x;
        }
    }P[205];
    void DP(int q){
        using namespace NotPrime;
        M=p[q];
        f[0][0]=1;long long tmp;
        for(int i=1;i<=t;++i){
            for(int j=1;j<=i;++j)f[i][j]=0;
            for(int k=0;k<i;++k){
                if(P[k].y>P[i].y)continue;
                int c=lucas(P[i].sum-P[k].sum,P[i].y-P[k].y,q);
                for(int j=1;j<=i;++j){
                    f[i][j]=(f[i][j]+f[k][j-1]*1LL*c%M)%M;
                }
            }
        }
        int lim=t+1;
        for(int k=0;k<=t;++k){
            tmp=0;
            for(int i=0;i<=t;++i){
                tmp+=f[i][k]*lucas(P[lim].sum-P[i].sum,P[lim].y-P[i].y,q);   
            }
            f[lim][k]=tmp%M;
        }
        int ans=0;
        for(int k=0;k<=t;++k){
            if(k&1)ans=(ans-f[t+1][k]+M)%M;
            else ans=(ans+f[t+1][k])%M;
        }
        res[q]=ans;
    }
};
int main(){
   //freopen("path.in","r",stdin); 
    //freopen("path.out","w",stdout);
    scanf("%lld%lld%d%d",&n,&m,&t,&mod);

    using namespace Ending;
    for(int i=1;i<=t;++i){
        P[i].read();
    }
    P[0].x=0;P[0].y=0;P[t+1].x=n;P[t+1].y=m;P[t+1].sum=n+m;
    sort(P+1,P+t+1);
    if(mod==1000003){
        using namespace Prime;
        init();
        f[0][0]=1;
        for(int i=1;i<=t;++i){
            for(int k=0;k<i;++k){
                for(int j=1;j<=i;++j){
                    if(P[i].y>=P[k].y)f[i][j]=(f[i][j]+f[k][j-1]*1LL*Lucas(P[i].x-P[k].x+P[i].y-P[k].y,P[i].y-P[k].y)%mod)%mod;
                }
            }
        }
        for(int k=0;k<=t;++k){
            for(int i=0;i<=t;++i){
                f[t+1][k]=(f[t+1][k]+f[i][k]*1LL*Lucas(P[t+1].x-P[i].x+P[t+1].y-P[i].y,P[t+1].y-P[i].y)%mod)%mod;              
            }
        }
        int ans=0;
        for(int k=0;k<=t;++k){
            if(k&1)ans=(ans-f[t+1][k]+mod)%mod;
            else ans=(ans+f[t+1][k])%mod;
        }
        printf("%d
",ans);
    }else{
        using namespace NotPrime;
        init();
        for(int i=0;i<4;++i)DP(i);
        printf("%d
",CRT());
    }
    //fclose(stdin);fclose(stdout);
    return 0;
}