BZOJ-3782 上学路线(dp+Lucas定理+CRT)

题目描述

  矩形网格起点为 \((0,0)\)，终点为 \((n,m)\)，只能向右或向上走，有 \(t\) 个坏点不能走，求方案数。

  数据范围：\(1\leq n,m\leq10^{10},0\leq t\leq 200,p=10000003\) 或 \(p=1019663265\)。

分析

  把 $ (n,m)$ 也当成坏点，编号即 \(t+1\)，把所有坏点排序。

  设 \(dp[i]\) 为到达第 \(i\) 个坏点的合法方案数（途中不经过其他坏点），答案即 \(dp[t+1]\)。但是很难直接求 \(dp[i]\)，考虑用总方案数减去经过坏点的个数的方案数计算 $dp[i] $。

  状态转移方程为：

\[dp[i]=\dbinom{x_i+y_i}{x_i}-\sum_{j=1}^{i-1}dp[j]·\dbinom{(x_i+y_i)-(x_j+y_j)}{x_i-x_j} \]

  \(p=1000003\)，直接用 $\text{Lucas} $ 定理算组合数。

  \(p=1019663265=3\times 5\times 6793\times 10007\)，不是一个质数，且 \(p\) 太大不能用 \(\text{exLucas}\)，分别求出 \(\dbinom{n}{m}\) 在模 \(3\)，模 \(5\)，模 \(6793\)，模 \(10007\) 下的值 \(a_1,a_2,a_3,a_4\)，最后用中国剩余定理求解线性同余方程组：

\[\begin{cases}x\equiv a_1\pmod {3}\\ x\equiv a_2\pmod {5}\\ x\equiv a_3\pmod {6793}\\ x\equiv a_4\pmod {10007}\end{cases} \]

代码

#include<bits/stdc++.h>
using namespace std;
long long n,m,t,p;
struct Point
{
    long long x,y;
}P[20010];
long long dp[20010];
bool cmp(Point A,Point B)
{
    if(A.x==B.x)
        return A.y<B.y;
    return A.x<B.x;
}
long long quick_pow(long long a,long long b,long long p)
{
    long long ans=1;
    while(b)
    {
        if(b&1)
            ans=ans*a%p;
        a=a*a%p;
        b>>=1;
    }
    return ans%p;
}
const int N=1e6;
struct mod1
{
    const int mod=1000003;
    long long fac[N+10],inv[N+10];
    void init()
    {
        fac[0]=inv[0]=1;
        for(int i=1;i<=N;i++)
            fac[i]=fac[i-1]*i%mod;
        inv[N]=quick_pow(fac[N],mod-2,mod);
        for(int i=N-1;i>=1;i--)
            inv[i]=inv[i+1]*(i+1)%mod;
    }
    long long C(long long n,long long m)
    {
        if(m==0)
            return 1;
        if(m>n)
            return 0;
        return fac[n]*inv[m]%mod*inv[n-m]%mod;
    }
    long long Lucas(long long n,long long m)
    {
        if(m==0)
            return 1;
        if(m>n)
            return 0;
        return C(n%mod,m%mod)*Lucas(n/mod,m/mod)%mod;
    }
}A;
struct mod2
{
    const long long mod=1ll*1019663265;
    long long prime[5]={0,3,5,6793,10007};
    long long fac[5][10010],inv[5][10010];
    long long exgcd(long long a,long long b,long long &x,long long &y)
    {
        if(b==0)
        {
            x=1;
            y=0;
            return a;
        }
        long long gcd=exgcd(b,a%b,x,y);
        int temp=x;x=y;y=temp-a/b*y;
        return gcd;
    }
    void init()
    {
        for(int op=1;op<=4;op++)
        {
            int _N=prime[op]-1;
            fac[op][0]=inv[op][0]=1;
            for(int i=1;i<=_N;i++)
                fac[op][i]=fac[op][i-1]*i%prime[op];
            inv[op][_N]=quick_pow(fac[op][_N],prime[op]-2,prime[op]);
            for(int i=_N-1;i>=1;i--)
                inv[op][i]=inv[op][i+1]*(i+1)%mod;
        }
    }
    long long C(long long n,long long m,long long op)
    {
        if(m==0)
            return 1;
        if(m>n)
            return 0;
        return fac[op][n]*inv[op][n-m]%prime[op]*inv[op][m]%prime[op];
    }
    long long Lucas(long long n,long long m,long long op)
    {
        if(m==0)
            return 1;
        if(m>n)
            return 0;
        return C(n%prime[op],m%prime[op],op)*Lucas(n/prime[op],m/prime[op],op)%prime[op];
    }
    long long a[5],b[5];
    long long solve(long long n,long long m)
    {
        long long ans=0,lcm=1,x,y;
        for(int i=1;i<=4;i++)
            lcm=lcm*prime[i];
        for(int i=1;i<=4;i++)
            a[i]=Lucas(n,m,i);
        for(int i=1;i<=4;i++)
            b[i]=prime[i];
        for(int i=1;i<=4;i++)
            a[i]=(a[i]%b[i]+b[i])%b[i];
        for(int i=1;i<=4;i++)
        {
            long long temp=lcm/prime[i];
            exgcd(temp,b[i],x,y);
            x=(x%b[i]+b[i])%b[i];
            ans=(ans+temp*x%lcm*a[i]%lcm)%lcm;
        }
        return (ans+lcm)%lcm;
    }
}B;
int main()
{
    cin>>n>>m>>t>>p;
    if(p==1000003)
        A.init();
    if(p==1019663265)
        B.init();
    for(int i=1;i<=t;i++)
        scanf("%lld %lld",&P[i].x,&P[i].y);
    P[t+1]={n,m};
    sort(P+1,P+2+t,cmp);
    for(int i=1;i<=t+1;i++)
    {
        if(p==1000003)      dp[i]=A.Lucas(P[i].x+P[i].y,P[i].x);
        if(p==1019663265)   dp[i]=B.solve(P[i].x+P[i].y,P[i].x);
        //cout<<P[i].x<<" "<<P[i].y<<" "<<B.solve(P[i].x+P[i].y,P[i].x)<<endl;
        for(int j=1;j<=i-1;j++)
        {
            if(p==1000003)
                dp[i]=(dp[i]-dp[j]*A.Lucas((P[i].x+P[i].y)-(P[j].x+P[j].y),P[i].x-P[j].x)%p+p)%p;
            if(p==1019663265)
                dp[i]=(dp[i]-dp[j]*B.solve((P[i].x+P[i].y)-(P[j].x+P[j].y),P[i].x-P[j].x)%p+p)%p;
        }
    }
    cout<<dp[t+1]<<endl;
    return 0;
}