考虑 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)/曹冲养猪。
中国剩余定理的结论大概是这个东西:
对于一个方程组:
\(\begin{cases} 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;
}
标签:dbinom,crt,BJWC2018,int,bmod,times,ll,Lucas
From: https://www.cnblogs.com/ez-lcw/p/16837097.html