传送门
首先设(m)为随从个数，(k)为暗影打击装甲的个数，(p)为剩余生命值，(n)为生命上限
然后考虑每个回合受到伤害，设(A_i)为每个回合被攻击(i)次的概率
[A_i=C^{i}_{k}*（frac{1}{m+1}）^i*(frac{m}{m+1})^{k-i}\
]
然后考虑答案，我们设(F_i)为当前有(i)滴血被干掉的期望回合数，再设(g_i)为当前回合打出超过(i)滴血的概率，容易得出
[F_i=frac{1}{m+1}(g_{i+1}+sum_{j=0}^{i}A_{j}(F_{i+1-j}-1))+frac{m}{m+1}(g_i+sum_{j=0}^{i-1}A_j(F_{i-j}+1))\
F_i=frac{1}{m+1}(1+sum_{j=0}^{i}A_{j}F_{i+1-j})+frac{m}{m+1}(1+sum_{j=0}^{i-1}A_jF_{i-j})\
]
好像没有什么问题，但是还得考虑满血的情况，这个时候是回不上血的
所以最后得出来
[i!=n:\
F_i=frac{1}{m+1}(1+sum_{j=0}^{i}A_{j}F_{i+1-j})+frac{m}{m+1}(1+sum_{j=0}^{i-1}A_jF_{i-j})\
i==n:\
F_n=1+sum_{j=0}^{n-1}A_jF_{n-j}
]
然后考虑概率期望的通常做法：高斯消元
由于(n=1500)，高斯消元显然会TLE
然后我们观察一下这个方程组，发现第(i)个方程只有前(i+1)个系数不为(0)，所以考虑每次消元都只留下(i,i+1)，然后最后一个方程只有第(n)个位置上系数不为(0)，然后解出(F_n)，递推出答案就好了，那么复杂度就可以优化到(O(n^2))了
提示：bzoj太慢，mod太多会TLE，请优化
代码：
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<vector>
#include<map>
using namespace std;
void read(int &x) {
char ch; bool ok;
for(ok=0,ch=getchar(); !isdigit(ch); ch=getchar()) if(ch=='-') ok=1;
for(x=0; isdigit(ch); x=x*10+ch-'0',ch=getchar()); if(ok) x=-x;
}
#define rg register
const int maxn=1.5e3+10,mod=1e9+7;
int n,T,m,p,k,a[maxn][maxn],fac[maxn],inv[maxn],A[maxn],facc[maxn];
int mul(int x,int y){return 1ll*x*y-1ll*x*y/mod*mod;}
int mi(int a,int b)
{
int ans=1;
while(b)
{
if(b&1)ans=1ll*ans*a%mod;
b>>=1,a=1ll*a*a%mod;
}
return ans;
}
void prepare()
{
fac[0]=inv[0]=1;
for(rg int i=1;i<=1500;i++)fac[i]=1ll*fac[i-1]*i%mod;
inv[1500]=mi(fac[1500],mod-2);;
for(rg int i=1499;i;i--)inv[i]=1ll*inv[i+1]*(i+1)%mod;
}
int C(int x,int y)
{
int ans=1;
if(y<=1500)ans=1ll*ans*fac[y]*inv[y-x]%mod;
else ans=1ll*ans*facc[x]%mod;
return 1ll*ans*inv[x]%mod;
}
void gauss()
{
int y[maxn];
for(rg int i=1;i<=n;i++)
{
for(rg int j=1;j<i;j++)
{
int t=1ll*a[i][j]*y[j]%mod;
a[i][j+1]=(a[i][j+1]-mul(t,a[j][j+1])+mod)%mod;
a[i][j]=0,a[i][n+1]=(a[i][n+1]-mul(t,a[j][n+1])+mod)%mod;
}
y[i]=mi(a[i][i],mod-2);
}
int ans=mul(a[n][n+1],y[n]);
for(rg int i=n-1;i>=p;i--)ans=1ll*(a[i][n+1]-mul(a[i][i+1],ans)+mod)%mod*y[i]%mod;
printf("%d
",ans?(ans+mod)%mod:-1);
}
int main()
{
read(T),prepare();
while(T--)
{
memset(a,0,sizeof a);
read(n),read(p),read(m),read(k);int res=0;
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;
}
for(rg int i=1;i<=n;i++)a[i][i]=a[i][n+1]=mod-1;
int ny=mi(m+1,mod-2),u=mul(ny,m);memset(A,0,sizeof A);facc[0]=1;
if(k>1500)for(rg int i=1;i<=n;i++)facc[i]=mul(facc[i-1],(k-i+1));
for(rg int i=0;i<=min(k,n);i++)A[i]=mul(C(i,k),mul(mi(ny,i),mi(u,k-i)));
for(rg int i=1;i<n;i++)
{
for(rg int j=0;j<=i;j++)a[i][i+1-j]=(a[i][i+1-j]+mul(ny,A[j]))%mod;
for(rg int j=0;j<i;j++)a[i][i-j]=(a[i][i-j]+mul(u,A[j]))%mod;
}
for(rg int i=0;i<n;i++)a[n][n-i]=(a[n][n-i]+A[i])%mod;
gauss();
}
}