「BJOI2018」治疗之雨

传送门

Description

有\(m+1\)个数，第一个数为\(p\)，每轮：选一个数\(+1\)，再依次选\(k\)个数\(-1\)

要求如果第一个数\(=N\)，不能选它\(+1\)，如果第一个数\(=0\)，不能选它\(-1\)

如果没有可选的数，跳过该次选择

问使得第一个数\(=0\)的期望步数

\(N\le1500\)，\(Case\le10\)

Solution

设\(f_i\)表示当第一个数为\(i\)时期望多少轮变为\(0\)

\[f_i=1+\sum_{j=1}^{i+1}p_{i,j}f_j,1\le i<n
\\
f_n=1+\sum_{j=1}^np_{i,j}f_j
\]

其中\(p_{i,j}\)表示一轮将第一个数从\(i\)变为\(j\)的概率

对第一个式子进行移项

\[p_{i,i+1}f_{i+1}=-1-\sum_{j=1}^{i-1}f_j+(1-p_{i,i})f_i,2\le i \le n
\]

问题在于不知道\(f_1\)，所以可以先把\(f_1\)设为\(x\)，推出\(f_n=a_1x+b_1\)

然后根据第二个式子，得到\(f_n=a_1x+b_2\)

最后解出\(f_1\)，求得\(f_p\)

Code

#include<bits/stdc++.h>

#define ll long long

#define dbg1(x) cerr<<#x<<"="<<(x)<<" "

#define dbg2(x) cerr<<#x<<"="<<(x)<<"\n"

#define dbg3(x) cerr<<#x<<"\n"

using namespace std;

#define reg register

inline int read()

{

    int x=0,f=1;char ch=getchar();

    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}

    while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}

    return x*f;

}

const int MN=1505,P=1e9+7;

int Mul(int x,int y){return (1ll*x*y)%P;}

int Add(int x,int y){return (x+y)%P;}

int fp(int x,int y=P-2){int r=1;for(;y;y>>=1,x=Mul(x,x))if(y&1)r=Mul(r,x);return r;}

struct d{

	int x,y;

	d(int x=0,int y=0):x(x),y(y){}

	d operator+(const d&o)const{return d(Add(x,o.x),Add(y,o.y));}

	d operator-(const d&o)const{return d(Add(x,P-o.x),Add(y,P-o.y));}

	d operator*(const int o)const{return d(Mul(x,o),Mul(y,o));}

}f[MN],or_;

int N,p,M,K,invM_,invM,a[MN][MN],g[MN],b[MN],inv[MN],X;

int solve()

{

	reg int i,j;

	if(K==0)return -1;

	if(M==0)

	{

		if(K==1&&N>1) return -1;

		g[0]=0;for(i=1;i<=N;++i)g[i]=g[max(min(N,i+1)-K,0)]+1;

		return g[p];

	}

	memset(b,0,sizeof b);memset(a,0,sizeof a);

	invM_=fp(M+1);invM=fp(M);b[0]=Mul(fp(M,K),fp(invM_,K));

	for(i=1;i<=N&&i<=K;++i)b[i]=Mul(Mul(b[i-1],invM),Mul(inv[i],K-i+1));

	for(i=1;i<N;++i)for(j=1;j<=i+1;++j)if(i-K<=j)

		a[i][j]=Add(Mul(invM_,b[i+1-j]),Mul(invM_,Mul(M,b[i-j])));

	for(i=1;i<=N;++i)a[N][i]=b[N-i];

	for(f[1]=d(1,0),i=2;i<=N;++i)

	{

		f[i]=d(0,P-1);

		for(j=1;j<i-1;++j)f[i]=f[i]-(f[j]*a[i-1][j]);

		f[i]=f[i]+f[i-1]*((1-a[i-1][i-1]+P)%P);

		f[i]=f[i]*fp(a[i-1][i]);

	}

	for(or_=d(0,1),i=1;i<N;++i)or_=or_+f[i]*a[N][i];or_=or_*fp(Add(1,P-a[N][N]));

	or_=or_-f[N];if(!or_.x&&f[p].x) return -1;

	X=Mul(or_.y,fp(P-or_.x));

	return Add(Mul(f[p].x,X),f[p].y);

}

int main()

{

	int cas=read();inv[0]=inv[1]=1;

	for(int i=2;i<=1500;++i)inv[i]=Mul((P-P/i),inv[P%i]);

	while(cas--)

	{

		N=read(),p=read(),M=read(),K=read();

		printf("%d\n",solve());

	}

	return 0;

}

Blog来自PaperCloud，未经允许，请勿转载，TKS！

巴特西

「BJOI2018」治疗之雨

传送门

Description

Solution

Code

最新文章

热门文章