题意
初始时玩家有 \(p\) 滴血，满血 \(n\) 滴，每个回合会进行如下操作：
-
若当前还没有满血，则以 \(\frac{1}{m+1}\) 的概率增加一滴血；
-
\(k\) 次判定，每次以 \(\frac{1}{m+1}\) 的概率扣除一滴血，当血量减少为 \(0\) 时游戏结束。
求玩家能存活的期望回合数。
\(1 \leq p \leq n \leq 1500\) ，\(0 \leq m, k \leq 10^9\)。
思路
设一个回合扣除 \(i\) 滴血的概率为 \(P_i(i \leq k)\)，由数学选修三的知识可知，这其实是一个伯努利试验的模型，可以得到表达式：
\[P_i=\binom{n}{i} (\frac{1}{{m+1}})^i (\frac{m}{{m+1}})^{n-i} \]由于：
\[\dfrac{P_i}{P_{i-1}}=\dfrac{\binom{n}{i} (\frac{1}{{m+1}})^i (\frac{m}{{m+1}})^{n-i}}{\binom{n}{i-1} (\frac{1}{{m+1}})^{i-1} (\frac{m}{{m+1}})^{n-i+1}} \]可以得到递推式：
\[P_i=P_{i-1} \times \frac{1}{m} \times \frac{k-i+1}{i} \]由于一个回合最多一次将满血全部扣完，所以这部分可以 \(O(n)\) 预处理得到。
设 \(f_i\) 为有 \(i\) 滴血的期望存活回合数，当 \(i\in [1,n-1]\) 时，简单分类讨论不难得到表达式：
\[f_i=\sum_{j=1}^{i}f_j (\frac{1}{m+1} P_{i+1-j}+\frac{m}{m+1}P_{i-j})+\frac{1}{m+1}P_0 f_{i+1}+1 \]其中当 \(i=n\) 时，由于此时不能再加血，所以此时的表达式为：
\[f_n=\sum_{j=1}^{n} P_{n-j} f_j+1 \]这里一共有 \(n\) 个方程以及 \(n\) 个位置数，高斯消元可以得出答案。然而直接消元的复杂度为 \(O(n^3)\)，需要考虑优化。
列出方程组，当 \(a_{i,j}=1\) 时表示这一项的系数不为 \(0\)，反之则表示这一项不存在，可以得到：
\[\begin{bmatrix} 1& 1& 0& 0& 0& \ldots 1&\\ 1& 1& 1& 0& 0& \ldots 1&\\ 1& 1& 1& 1& 0& \ldots 1&\\ \ldots \\ 1& 1& 1& 1& 1& \ldots 1&\\ \end{bmatrix}\]发现这个方程的形式很特殊，每次往下消元的时候只需要带入两个值，消去这一行第一个非空的系数，最终得到的形式是这样的：
\[\begin{bmatrix} 1& 1& 0& 0& 0& \ldots 1&\\ 0& 1& 1& 0& 0& \ldots 1&\\ 0& 0& 1& 1& 0& \ldots 1&\\ \ldots \\ 0& 0& 0& 0& 1& \ldots 1&\\ \end{bmatrix}\]最后一行就是 \(x=A\) 的形式了，那么只需要向上回带就可以解出所有的未知数，最终的复杂度为 \(O(n^2)\)。
最后，还需要特判一下 \(k=0\) 和 \(m=0\) 的情况。
code：
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define LL long long
const int N=2023,mod=1e9+7;
int n,p,m,k;LL P[N],a[N][N];//初始p滴血，满血n滴，1/(m+1)概率加血，1/(m+1)概率扣血->k次
int mul(int a,int b){int res=1;while(b) ((b&1)&&(res=1ll*res*a%mod,0)),a=1ll*a*a%mod,b>>=1;return res;}
void Gauss()
{
for(int i=1;i<=n;i++)
{
LL inv=mul(a[i][i],mod-2);a[i][i]=1;a[i][n+1]=a[i][n+1]*inv%mod;
if(i!=n) a[i][i+1]=a[i][i+1]*inv%mod;
for(int j=i+1;j<=n;j++)
{
LL mult=a[j][i];a[j][i]=0;a[j][i+1]=(a[j][i+1]-mult*a[i][i+1]%mod+mod)%mod;
a[j][n+1]=(a[j][n+1]-mult*a[i][n+1]%mod+mod)%mod;
}
}
for(int i=n;i>1;i--) a[i-1][n+1]=(a[i-1][n+1]-a[i-1][i]*a[i][n+1]%mod+mod)%mod,a[i-1][i]=0;
}
void solve()
{
scanf("%d%d%d%d",&n,&p,&m,&k);
if(k==0) return puts("-1"),void();
if(m==0)
{
if(k==1) return puts("-1"),void();int res=0;
while(p>0){if(p<n) p++;p-=k;res++;} return printf("%d\n",res),void();
}
for(int i=1;i<N;i++) P[i]=0;
LL invm=mul(m,mod-2),invm1=mul(m+1,mod-2);P[0]=mul(1ll*m*mul(m+1,mod-2)%mod,k);
for(int i=1;i<=min(n+1,k);i++) P[i]=P[i-1]*invm%mod*mul(i,mod-2)%mod*(k-i+1)%mod;
for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++) a[i][j]=0;
for(int i=1;i<n;i++) for(int j=1;j<=i;j++) a[i][j]=(P[i-j]*m+P[i-j+1])%mod*invm1%mod;
for(int i=1;i<n;i++) a[i][i+1]=P[0]*invm1%mod;
for(int i=1;i<n;i++) a[i][i]=(a[i][i]-1+mod)%mod;
for(int i=1;i<=n;i++) a[i][n+1]=mod-1;
for(int i=1;i<=n;i++) a[n][i]=P[n-i];a[n][n]=(a[n][n]-1+mod)%mod;
Gauss();printf("%lld\n",(a[p][n+1]%mod+mod)%mod);
}
int main()
{
// freopen("in.txt","r",stdin);
int T;scanf("%d",&T);while(T--) solve();
return 0;
}