洛谷P4457/loj#2513 [BJOI2018]治疗之雨（高斯消元+概率期望）

题面

传送门(loj)

传送门(洛谷)

题解

模拟赛的时候只想出了高斯消元然后死活不知道怎么继续……结果正解居然就是高斯消元卡常？

首先有个比较难受的地方是它一个回合可能不止扣一滴血……我们得算出(P_i)表示一回合扣(i)滴血的概率，为

[P_i={{kchoose i}m^{k-i}over (m+1)^k} ]

所以这个柿子啥意思？

我们可以把(k)次扣血看成一个长度为(k)的序列，每个序列有(m+1)种选择方法，于是总的选法就是((m+1)^k)。我们要钦定它扣(i)滴血，就是令其中(i)个数强制选择(1)，方案数为({kchoose i})。剩下的数都不能选(1)，所以方案数为(m^{k-i})

有了这个(P_i)我们就可以考虑找关系了

然而这里还有一个很讨厌的地方就是它时不时会被奶一口……

我们记(f_i)表示当血量为(i)时被干掉的期望回合数，(g_i)表示一回合内打出伤害大于等于(i)的概率，然后考虑这东西该怎么转移

[f_i={mover m+1}left(g_i+sum_{k=0}^{i-1}P_k(f_{i-k}+1) ight)+{1over m+1}left(g_{i+1}+sum_{k=0}^{i}P_k(f_{i+1-k}+1) ight) ]

边界条件为

[f_n=g_n+sum_{k=0}^{n-1}P_k(f_{n-k}+1) ]

所以这柿子是啥？

({mover m+1})表示没有被奶到的概率，那么我们枚举它被(A)了几下。如果它一回合被干掉了，那么期望局数为(1)，否则(A)了(k)下之后血量会到(i-k)，这一部分的期望步数就是(f_{i-k}+1)。后面那个就是如果被奶了之后的情况。顺便因为满血的时候是不可能被奶的，所以(f_n)要特殊考虑

然而这柿子一点都不好看，特别是(g_i)很不爽，那就继续推倒

[f_i={mover m+1}left(g_i+sum_{k=0}^{i-1}P_k+sum_{k=0}^{i-1}P_kf_{i-k} ight)+{1over m+1}left(g_{i+1}+sum_{k=0}^{i}P_k+sum_{k=0}^{i}P_kf_{i+1-k} ight) ]

[f_i={mover m+1}left(1+sum_{k=0}^{i-1}P_kf_{i-k} ight)+{1over m+1}left(1+sum_{k=0}^{i}P_kf_{i+1-k} ight) ]

同理有(f_n=1+sum_{k=0}^{n-1}P_kf_{n-k})

那么我们就可以愉快地递推……等会儿这咋推啊……转移好像都成环了啊……

那么我们就把它看成一个方程组来高斯消元求解吧

哈？(n=1500)你让我高斯消元？

我们来好好观察一下这个方程组，(f_i)所对应的方程只与(f_1,...f_{i+1})的值有关，也就是说每一行的对角线右边只有在它右边一位的那个系数不为(0)

因为我们高斯消元的时候是拿自己这行去减下面的，所以每一行中只有(2)个系数要去和下面的相减，复杂度就能化到(O(n^2))了

虽然我还是不知道这个复杂度是怎么过去的

//minamoto
#include<bits/stdc++.h>
#define R register
#define fp(i,a,b) for(R int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(R int i=a,I=b-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
template<class T>inline bool cmin(T&a,const T&b){return a>b?a=b,1:0;}
using namespace std;
char buf[1<<21],*p1=buf,*p2=buf;
inline char getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}
int read(){
    R int res,f=1;R char ch;
    while((ch=getc())>'9'||ch<'0')(ch=='-')&&(f=-1);
    for(res=ch-'0';(ch=getc())>='0'&&ch<='9';res=res*10+ch-'0');
    return res*f;
}
char sr[1<<21],z[20];int K=-1,Z=0;
inline void Ot(){fwrite(sr,1,K+1,stdout),K=-1;}
void print(R int x){
    if(K>1<<20)Ot();if(x<0)sr[++K]='-',x=-x;
    while(z[++Z]=x%10+48,x/=10);
    while(sr[++K]=z[Z],--Z);sr[++K]='
';
}
const int N=1505,P=1e9+7;
inline int add(R int x,R int y){return x+y>=P?x+y-P:x+y;}
inline int dec(R int x,R int y){return x-y<0?x-y+P:x-y;}
inline int mul(R int x,R int y){return 1ll*x*y-1ll*x*y/P*P;}
int ksm(R int x,R int y){
	R int res=1;
	for(;y;y>>=1,x=mul(x,x))if(y&1)res=mul(res,x);
	return res;
}
int a[N][N],ans[N],s[N],inv[N],g[N];
int n,m,p,k,T,qwq,aa,bb,tmp,res;
void init(int n){
	inv[0]=inv[1]=1;
	fp(i,2,n)inv[i]=mul(P-P/i,inv[P%i]);
}
void Gauss(int n){
	fp(i,1,n-1){
		int t=ksm(a[i][i],P-2);
		a[i][i]=1,a[i][n]=mul(a[i][n],t);
		if(i!=n-1)a[i][i+1]=mul(a[i][i+1],t);
		fp(j,i+1,n-1){
			t=a[j][i],a[j][i]=0;
			a[j][i+1]=dec(a[j][i+1],mul(t,a[i][i+1])),
			a[j][n]=dec(a[j][n],mul(t,a[i][n]));
		}
	}
    ans[n-1]=a[n-1][n];
    fd(i,n-2,1)ans[i]=dec(a[i][n],mul(a[i][i+1],ans[i+1]));
}
int main(){
//	freopen("testdata.in","r",stdin);
	T=read();
	init(N-5);
	while(T--){
		n=read(),p=read(),m=read(),k=read();
		if(!k||(!m&&k==1)){puts("-1");continue;}
		if(!m){
			while(p>0){if(p<n)++p;p-=k,++res;}
			printf("%d
",res),res=0;continue;
		}
		qwq=ksm(m+1,k),qwq=ksm(qwq,P-2)%P,tmp=1;
		fp(i,0,n)s[i]=0;
		fp(i,0,min(n,k)){
			s[i]=1ll*tmp*ksm(m,k-i)%P*qwq%P;
			tmp=1ll*tmp*inv[i+1]%P*(k-i)%P;
		}
		memset(a,0,sizeof(a));
		bb=ksm(m+1,P-2),aa=mul(m,bb);
		fp(i,1,n-1){
			++a[i][i],++a[i][n+1];
			fp(j,0,i-1)a[i][i-j]=dec(a[i][i-j],mul(s[j],aa));
			fp(j,0,i)a[i][i+1-j]=dec(a[i][i+1-j],mul(s[j],bb));
		}
		++a[n][n],++a[n][n+1];
		fp(j,0,n-1)a[n][n-j]=dec(a[n][n-j],s[j]);
		Gauss(n+1);
		printf("%d
",ans[p]);
	}
	return 0;
}