【洛谷4457】[BJOI2018] 治疗之雨（高斯消元好题）

点此看题面

• 有一个体力上限为\(n\)，初始体力为\(p\)的英雄和\(m\)个体力无限、体力无上限的随从。
• 每轮操作会随机选中\(1\)个己方单位增加\(1\)点体力值，然后随机选择\(k\)次，每次选中\(1\)个己方单位扣除\(1\)点体力值（可重复选择）。重复执行直至英雄体力小于等于\(0\)。
• 求期望操作轮数。
• 数据组数\(\le10^3,p\le n\le1500,m,k\le10^9\)

一个前置定义

设\(P_i\)表示一轮操作的\(k\)次选择中英雄被选中\(i\)次的概率。

那么选中英雄的选择有\(C_k^i\)种组合方案，除了这\(i\)个以外每个选择都有\(m\)种选法，因此就有：

\[P_i=\frac{C_k^i\times m^{k-i}}{(m+1)^k} \]

然而\(k\)的范围是\(10^9\)，显然不可能预处理组合数。而有用的\(P_i\)实际上只有\(i\in[0,n]\)。

因此我们考虑把这个式子展开以递推，得到：

\[P_0=\frac{m^n}{(m+1)^n}\\ P_i=P_{i-1}\times \frac{k-(i-1)}{i}\times \frac1m(i\ge 1) \]

这样一来就能预处理出\(P\)数组了，它在之后的解题过程中会起到非常关键的作用。

概率\(DP\)

我们设\(f_i\)表示英雄体力值达到\(i\)的期望轮数。

考虑\(j\)对于\(i\)的贡献，有\(\frac m{m+1}\)的概率没有选中\(i\)恢复体力，有\(\frac1{m+1}\)的概率选中\(i\)恢复体力，因此有转移方程：

\[f_i=(\frac{m}{m+1}P_{i-j}+\frac1{m+1}P_{i-j+1})\times f_j \]

特殊地，对于\(j=i+1\)的时候有\(f_i=\frac1{m+1}P_0\times f_{i+1}\)，对于\(i=n\)的时候有\(f_n=P_{n-j}\times f_j\)。

显然这时候的\(DP\)方程是成环的，套路地考虑高斯消元。

然而，我们发现，由于\(n\le1500\)，\(O(n^3)\)会直接\(T\)飞！

特殊的高斯消元

实际上，这题的矩阵是非常特殊的。

容易发现，对于第\(i\)行，初始只有第\(1\sim i+1\)个元素有值。

当我们消到第\(i\)行的时候，第\(1\sim i-1\)个元素都已被消去，就只剩\(i\)和\(i+1\)两个元素，那么这时候高斯消元枚举列的一重循环的复杂度就变成了常数级别。

而最后倒枚的过程中，实际上只有第\(i-1\)行能有第\(i\)个元素，因此我们只用去消上一行就可以了。

因此总复杂度\(O(n^2)\)。

代码：\(O(Tn^2)\)

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 1500
#define X 1000000007
using namespace std;
int n,m,p,k,P[N+5],C[N+5][N+5];
I int QP(RI x,RI y) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
I int BF(RI x) {RI t=0;W(x>0) x^n&&++x,x-=k,++t;return t;}//暴力
int a[N+5][N+5],v[N+5];I void Gauss()//高斯消元
{
	RI i,j,t;for(i=1;i<=n;++i)//枚举行
	{
		t=QP(a[i][i],X-2),a[i][i]=1,a[i][i+1]=1LL*a[i][i+1]*t%X,v[i]=1LL*v[i]*t%X;//把a[i][i]变成1，方便后续操作
		for(j=i+1;j<=n;++j) v[j]=(v[j]-1LL*a[j][i]*v[i]%X+X)%X,//枚举后面的行
			a[j][i+1]=(a[j][i+1]-1LL*a[j][i]*a[i][i+1]%X+X)%X,a[j][i]=0;//常数级别枚举列消元
	}
	for(i=n;i^1;--i) v[i-1]=(v[i-1]-1LL*a[i-1][i]*v[i]%X+X)%X,a[i-1][i]=0;//最后倒枚一遍
}
int main()
{
	RI Tt,i,j,Inv;scanf("%d",&Tt);W(Tt--)
	{
		if(scanf("%d%d%d%d",&n,&p,&m,&k),!k) {puts("-1");continue;}//特判k=0
		if(!m) {k==1?puts("-1"):printf("%d\n",BF(p));continue;}//特判m=0
		for(Inv=QP(m+1,X-2),P[0]=1LL*QP(m,k)*QP(Inv,k)%X,
			i=1;i<=n;++i) P[i]=i<=k?1LL*P[i-1]*QP(i,X-2)%X*(k-i+1)%X*QP(m,X-2)%X:0;//递推求出P[i]
		for(i=1;i^n;++i) for(a[i][i+1]=1LL*Inv*P[0]%X,//预处理前i行第i+1项的系数
			j=1;j<=i;++j) a[i][j]=(1LL*P[i-j]*m+P[i-j+1])%X*Inv%X;//预处理前i行第1~i项的系数
		for(j=1;j<=n;++j) a[n][j]=P[n-j];//预处理第n行的系数
		for(i=1;i<=n;++i) !a[i][i]--&&(a[i][i]+=X),v[i]=X-1;//每行两个特殊项的系数
		Gauss(),printf("%d\n",v[p]);//高斯消元后输出f[p]
	}return 0;
}