【BZOJ5292】[BJOI2018]治疗之雨（高斯消元）

题面

题解

设\(f[i]\)表示剩余\(i\)点生命时的期望死亡的次数。

考虑打\(k\)次下来脸上被打了\(i\)下的概率：\(\displaystyle \frac{{k\choose i}m^{k-i}}{(m+1)^k}\)。

\(m=0\)时全部打脸上了，直接判掉。

设\(P[i][j]\)表示\(i\)点血量奶完后再被打一轮下来变成\(j\)点血的概率，这个很容易算出来。

那么我们可以列出和\(f[i]\)相关的\(n+1\)个方程。

比如说：

\[\begin{cases}
f[0]=0\\
f[1]=P[1][0]*f[0]+P[1][1]*f[1]+P[1][2]*f[2]+1\\
...\\
f[n]=P[n][0]*f[0]+P[n][1]*f[1]+P[n][2]*f[2]+...+P[n][n]*f[n]+1
\end{cases}\]

考虑怎么解这个玩意，首先肯定可以高斯消元三方解决。

这个矩阵观察后发现类似于一个下三角矩阵，那么每次用第\(i\)行消掉\(i-1\)行，这样子就是一个下三角了，然后从上往下就可以直接求解。

还有一种方法就是移项后不难发现\(f[1]\)只和\(f[2]\)相关，因此可以设\(f[1]=x\)，这样子其他所有数都可以通过\(x\)给表示出来，最后带到\(f[n]\)的方程里就可以解出\(x\)。

#include<iostream>

#include<cstdio>

using namespace std;

#define ll long long

#define MAX 1600

#define MOD 1000000007

void add(int &x,int y){x+=y;if(x>=MOD)x-=MOD;}

inline int read()

{

	int x=0;bool t=false;char ch=getchar();

	while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();

	if(ch=='-')t=true,ch=getchar();

	while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();

	return t?-x:x;

}

int fpow(int a,int b){int s=1;while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}return s;}

int n,p,m,k;

int P[MAX],a[MAX][MAX];

int A[MAX],B[MAX];

int main()

{

	int T=read();

	while(T--)

	{

		n=read();p=read();m=read();k=read();

		if(!k){puts("-1");continue;}

		if(m==0)

		{

			int cnt=0;if(k==1){puts("-1");continue;}

			while(p>0)++cnt,p=min(n,p+1)-k;

			printf("%d\n",cnt);continue;

		}

		for(int i=0;i<=n;++i)P[i]=0;

		int invm1=fpow(m+1,MOD-2);P[0]=1;

		for(int i=1;i<n&&i<=k;++i)P[i]=1ll*P[i-1]*fpow(i,MOD-2)%MOD*(k-i+1)%MOD;

		for(int i=0;i<n&&i<=k;++i)P[i]=1ll*P[i]*fpow(m,k-i)%MOD*fpow(fpow(m+1,k),MOD-2)%MOD;

		if(k>=n){P[n]=1;for(int i=0;i<n;++i)P[n]=(P[n]+MOD-P[i])%MOD;}

		for(int i=0;i<=n;++i)

			for(int j=0;j<=n;++j)a[i][j]=0;

		for(int i=1;i<n;++i)

			for(int j=1;j<=i+1;++j)

				add(a[i][j],(1ll*invm1*P[i-j+1]+1ll*m*invm1%MOD*P[i-j])%MOD);

		for(int i=1;i<=n;++i)add(a[n][i],P[n-i]);

		for(int i=0;i<=n;++i)A[i]=B[i]=0;

		A[1]=1;B[1]=0;

		for(int i=1;i<n;++i)

		{

			int sa=A[i],sb=(MOD-1+B[i])%MOD;

			for(int j=1;j<=i;++j)sa=(sa+MOD-1ll*a[i][j]*A[j]%MOD)%MOD;

			for(int j=1;j<=i;++j)sb=(sb+MOD-1ll*a[i][j]*B[j]%MOD)%MOD;

			int v=fpow(a[i][i+1],MOD-2);

			A[i+1]=1ll*sa*v%MOD;B[i+1]=1ll*sb*v%MOD;

		}

		int sa=A[n],sb=(MOD-B[n]+1)%MOD;

		for(int i=1;i<=n;++i)sa=(sa+MOD-1ll*a[n][i]*A[i]%MOD)%MOD;

		for(int i=1;i<=n;++i)sb=(sb+1ll*a[n][i]*B[i])%MOD;

		int x=1ll*sb*fpow(sa,MOD-2)%MOD;

		int ans=(1ll*A[p]*x+B[p])%MOD;

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

		continue;

	}

}