首页 > 其他分享 >【BJWC2018】上学路线(dp,Lucas,crt)

【BJWC2018】上学路线(dp,Lucas,crt)

时间:2022-10-28 18:46:28浏览次数:79  
标签:dbinom crt BJWC2018 int bmod times ll Lucas

考虑 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

相关文章