题解 P4478 【[BJWC2018]上学路线】

· · 题解

考虑 dp。

我们先把 (n,m) 也当做障碍点，然后把所有的障碍点按 x 坐标为第一关键字，y 坐标为第二关键字排序。

然后设 f_i 表示到达第 i 个障碍点的合法总方案数（途中不经过障碍点）。显然，答案就是 f_{t+1}，也就是到达 (n,m) 的总方案数。

至于为什么要先排序，是因为我们要保证当处理 f_i 时，能转移到 f_i 的所有 f_j 都已经处理完了。

显然有状态转移方程：（其中 x_i 表示第 i 个障碍点的 x 坐标，y_i 同理）

f_i=\binom{x_i+y_i}{x_i}-\sum_{j=1}^{i-1}f_j\times\binom{(x_i+y_i)-(x_j+y_j)}{x_i-x_j}

首先，\dbinom{x_i+y_i}{x_i} 是从 (0,0) 到 (x_i,y_i) 的总方案（不论是否合法）。

证明：从 (0,0) 到 (x_i,y_i) 一共需要 x_i+y_i 步，其中需要横着走 x_i 步，竖着走 y_i 步，那么我们就枚举在哪几步横着走，也就是 \dbinom{x_i+y_i}{x_i}。

然后 \dbinom{(x_i+y_i)-(x_j+y_j)}{x_i-x_j} 是从 (x_j,y_j) 到 (x_i,y_i) 的总方案（不论是否合法），证明同理。

然后证明那些不合法的方案就是 \sum_{j=1}^{i-1}f_j\times\dbinom{(x_i+y_i)-(x_j+y_j)}{x_i-x_j}。

设某一条不合法路径上的第一个障碍点是第 first 个障碍点。

那么 f_j\times\dbinom{(x_i+y_i)-(x_j+y_j)}{x_i-x_j} 表示的就是那些 first=j 的不合法路径，因为这种路径从 (0,0) 到 (x_j,y_j) 是没有障碍点的，正好对应 f_j。

这么解释的话，就容易证明 \sum_{j=1}^{i-1}f_j\times\dbinom{(x_i+y_i)-(x_j+y_j)}{x_i-x_j} 和所有的不合法路径是一一对应的了。

那么这个状态转移方程就是正确的了。

现在的问题就是 \dbinom{n}{m}\bmod P 应该怎么求。

显然，当 P=1000003\in prime 时，这个可以用 Lucas 定理直接算。

但是当 P=1019663265=3\times5\times6793\times10007 时，就不能直接用 Lucas 算了。

这个需要用中国剩余定理来做，不会的请先去做一下 【模板】中国剩余定理(CRT)/曹冲养猪。

中国剩余定理的结论大概是这个东西：

对于一个方程组：

x\equiv a_1 \pmod {p_1}\\ x\equiv a_2 \pmod {p_2}\\ \cdots\\ x\equiv a_k \pmod {p_k} \end{cases}

其中 p_1,p_2,\dots,p_k 两两互质。

设 M=\prod\limits_{i=1}^{k}p_i，m_i=\dfrac{M}{p_i}，t_i 是 m_i 模 p_i 意义下的逆元。那么 x 的一个特殊解就是 x_0=\sum\limits_{i=1}^{k}a_im_it_i。

然后又因为 x 一定能表示成 a\times M+b 的形式，所以 x 的最小非负数解就是 x_0\bmod M。

这个题中，未知数 x 就是 \dbinom{n}{m}，我们要求 x\bmod P，然后我们设 p_1=3、p_2=5、p_3=6793、p_4=10007，那么 M 就是这题中的 P。

然后我们还知道 a_1=x\bmod p_1、a_2=x\bmod p_2、a_3=x\bmod p_3、a_4=x\bmod p_4 的值（这些都可以用 Lucas 定理求出来）。

那么我们通过中国剩余定理就可以把 x\bmod P 算出来了。

具体代码如下：

#include<bits/stdc++.h>

#define T 210
#define ll long long

using namespace std;

struct Point
{
    int x,y;
    Point(){};
    Point(int a,int b){x=a,y=b;}
}q[T];

int n,m,num,P;
ll f[T];

ll poww(ll a,ll b,ll p)
{
    ll ans=1;
    while(b)
    {
        if(b&1) ans=ans*a%p;
        a=a*a%p;
        b>>=1;
    }
    return ans;
}

namespace mod1//p=1019663265的情况
{
    const ll M=1019663265ll,p[5]={0,3,5,6793,10007};
    ll m[5],t[5],fac[5][10010],inv[5][10010];
    ll a[5];
    void exgcd(ll a,ll b,ll &x,ll &y)
    {
        if(!b)
        {
            x=1,y=0;
            return;
        }
        exgcd(b,a%b,x,y);
        ll tmp=x;
        x=y;
        y=tmp-a/b*y;
    }
    void init()
    {
        for(int i=1;i<=4;i++)
        {
            m[i]=M/p[i];
            ll x,y;
            exgcd(m[i],p[i],x,y);//通过exgcd求逆元
            t[i]=(x%p[i]+p[i])%p[i];
        }
        for(int i=1;i<=4;i++)
        {
            fac[i][0]=1;
            for(int j=1;j<p[i];j++)
                fac[i][j]=fac[i][j-1]*j%p[i];
        }
        for(int i=1;i<=4;i++)
            for(int j=0;j<p[i];j++)
                inv[i][j]=poww(fac[i][j],p[i]-2,p[i]);
    }
    ll C(ll n,ll m,int k)
    {
        if(n<m) return 0;
        if(!m) return 1;
        return fac[k][n]*inv[k][n-m]%p[k]*inv[k][m]%p[k];
    }
    ll Lucas(ll n,ll m,int k)
    {
        if(n<m) return 0;
        if(!m) return 1;
        return C(n%p[k],m%p[k],k)*Lucas(n/p[k],m/p[k],k)%p[k];
    }
    ll solve(ll n,ll mm)
    {
        ll x=0;
        for(int i=1;i<=4;i++)
        {
            a[i]=Lucas(n,mm,i);
            x=(x+a[i]*m[i]%M*t[i]%M)%M;//求x的解
        }
        return x;
    }
}

namespace mod2//p=1000003的情况
{
    const ll p=1000003;
    ll fac[1000005],inv[1000005];
    void init()
    {
        fac[0]=1;
        for(int i=1;i<p;i++)
            fac[i]=fac[i-1]*i%p;
        for(int i=0;i<p;i++)
            inv[i]=poww(fac[i],p-2,p);
    }
    ll C(ll n,ll m)
    {
        if(n<m) return 0;
        if(!m) return 1;
        return fac[n]*inv[n-m]%p*inv[m]%p;
    }
    ll Lucas(ll n,ll m)
    {
        if(n<m) return 0;
        if(!m) return 1;
        return C(n%p,m%p)*Lucas(n/p,m/p)%p;
    }
}

bool cmp(Point a,Point b)
{
    if(a.x==b.x) return a.y<b.y;
    return a.x<b.x;
}

int main()
{
    scanf("%d%d%d%d",&n,&m,&num,&P);
    if(P==1000003) mod2::init();
    else mod1::init();
    for(int i=1;i<=num;i++)
        scanf("%d%d",&q[i].x,&q[i].y);
    q[num+1]=Point(n,m);
    sort(q+1,q+num+2,cmp);
    for(int i=1;i<=num+1;i++)
    {
        if(P==1000003) f[i]=mod2::Lucas(q[i].x+q[i].y,q[i].x);
        else f[i]=mod1::solve(q[i].x+q[i].y,q[i].x);
        for(int j=1;j<i;j++)
        {
            if(q[j].x<=q[i].x&&q[j].y<=q[i].y)
            {
                if(P==1000003) f[i]=(f[i]-f[j]*(mod2::Lucas(q[i].x+q[i].y-q[j].x-q[j].y,q[i].x-q[j].x))%P+P)%P;
                else f[i]=(f[i]-f[j]*(mod1::solve(q[i].x+q[i].y-q[j].x-q[j].y,q[i].x-q[j].x))%P+P)%P;
            }
        }
    }
    printf("%lld\n",f[num+1]);
    return 0;
}