在洛谷上发了份题解，链接：https://www.luogu.org/blog/user22151/solution-p4457
我不是很会矩阵做这道题(不会卡常)。
但是我们可以用本质差不多的方法，但是不用矩阵，可能常数会小一点？反正目前洛谷开O2不开O2、bz都是最快的……
用$dp[i]$表示还剩$i$点血的时候，期望还要有多少轮操作。
用$f[x]$表示$K$次中有$i$次打到英雄的概率。
$f[x]=C(K,x)(frac{1}{m+1})^x(frac{m}{m+1})^{(k-x)}$
首先对于$m=0$的情况特判。
对于$m!=0$情况，我们发现有
$dp[0]=0$
$dp[1]=a_{1,2}dp[2]+a_{1,1}dp[1]+1$
...
$dp[i]=a_{i,i+1}dp[i+1]+a_{i,i}dp[i]+...+a_{i,1}dp[1]+1$
...
$dp[n]=a_{n,n}dp[n]+a_{n,n-1}dp[n-1]+...+a_{n,1}dp[1]+1$
其中$a_{i,j}$表示英雄$i$点血在一轮操作后变成$j$点血的概率，这个可以根据$f$算出来，不再赘述。
所以说我们可以得到$n$个等式，第$i$个式子是$dp[i]$等于多少多少
$n$个未知数，$n$个等式(方程)，很优秀的样子。
但是我们发现，第$i$个等式右边有$dp[i+1]$，所以不能直接递推，感觉很恼火。
我们把第$i$个式子右边的$dp[i+1]$放到等式左边来，然后把$dp[i]$放到等式右边，那么这个式子就变成了$dp[i+1]$等于多少多少。
注意第$n$个式子例外，因为和$dp[n+1]$无关，仍是$dp[n]$等于多少多少。
所以这样我们可以用第$i$个式子算$dp[i+1]$？
但是我们算不出$dp[1]$呀，也没法直接往后递推。
那我们先设$dp[1]=X$，那么我们可以根据第$1$~$n-1$个式子，把$dp[i]$表示成$AX+B$的形式，递推算出$dp[n]=A_1X+B_1$，然后根据第$n$个式子算出$dp[n]=A_2X+B_2$。
那么$A_1X+B_1=A_2X+B_2$，解一元一次方程直接算出$X$。
注意特判解不出$X$的情况输出-1。
//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(register int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(register int i=(a);i>=(b);--i)
const int maxn=3000+7;
const ll mod=1e9+7;
ll Td,n,P,m,K,X;
char cc; ll ff;
template<typename T>void read(T& aa) {
aa=0;cc=getchar();ff=1;
while((cc<'0'||cc>'9')&&cc!='-') cc=getchar();
if(cc=='-') ff=-1,cc=getchar();
while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
aa*=ff;
}
ll qp(ll x,ll k) {
if(x<=1) return x;
ll rs=1;
while(k) {
if(k&1) rs=rs*x%mod;
k>>=1; x=x*x%mod;
}
return rs;
}
ll mo(ll x) {if(x>=mod) x-=mod; return x;}
struct T{
ll x,y;
T(ll x=0.,ll y=0.):x(x),y(y){}
inline T operator + (const T& b) const{return T(mo(x+b.x),mo(y+b.y));}
inline T operator - (const T& b) const{return T(mo(x-b.x+mod),mo(y-b.y+mod));}
inline T operator * (const T&b) const{return T((x*b.y+y*b.x)%mod,y*b.y%mod);}
}dp[maxn],f[maxn],o,r,t,inv,finv;
ll work() {
if(K==0) return -1;
if(K==1&&n!=1) return -1;
dp[0]=T(0,0);
ll x; T p;
For(i,1,n-1) {
dp[i]=T(0,1)+dp[max(i+1-K,(ll)0)];
}
dp[n]=T(0,1)+dp[max(n-K,(ll)0)];
return dp[P].y;
}
ll solve() {
if(P==0) return 0;
if(n==1&&K==0) return -1;
ll x; T p;
o=T(0,m+1); r=T(0,m); inv=T(0,qp(qp(m+1,K),mod-2));
f[0]=T(0,1); x=min(K,n);
For(i,1,x) f[i]=f[i-1]*T(0,(K-i+1)*qp(i,mod-2)%mod);
For(i,0,x) f[i]=f[i]*inv*T(0,qp(m,K-i));
if(m==0) return work();
finv=T(0,qp(f[0].y,mod-2));
if(f[0].y==1||f[0].y==0) return -1;
// cerr<<clock()<<"
";
//f[i]: probability of get i of K , finv : inv of f[0]
dp[0]=T(0,0);
dp[1]=T(1,0);
//dp[i+1]=((m+1)(dp[i]-1)-m*sum(f[j]*dp[i-j])-sum(f[j]*dp[i-j+1]))/f[0]
For(i,1,n-1) {
dp[i+1]=o*dp[i]-o;
x=min((ll)i,K);
p=T(0,0);
For(j,0,x) p=p+f[j]*dp[i-j];
dp[i+1]=dp[i+1]-p*r;
x=min((ll)i+1,K);
For(j,1,x) dp[i+1]=dp[i+1]-f[j]*dp[i-j+1];
dp[i+1]=dp[i+1]*finv;
}
// cerr<<clock()<<"
";
//(1-f[0])*dp[n]=sum(f[j]*dp[n-j])+1;
t=T(0,1); x=min(n,K);
For(j,1,x) t=t+f[j]*dp[n-j];
t=t*T(0,qp(mo(mod+1-f[0].y),mod-2));
t=t-dp[n];
if(t.x==0&&dp[P].x!=0) return -1;
if(t.x) X=(mod-1)*t.y%mod*qp(t.x,mod-2)%mod;
return (X*dp[P].x%mod+dp[P].y)%mod;
}
int main() {
read(Td);
while(Td--) {
read(n); read(P); read(m); read(K);
printf("%lld
",solve());
}
// cerr<<clock()<<"
";
return 0;
}