- 有一个体力上限为(n)，初始体力为(p)的英雄和(m)个体力无限、体力无上限的随从。
- 每轮操作会随机选中(1)个己方单位增加(1)点体力值，然后随机选择(k)次，每次选中(1)个己方单位扣除(1)点体力值（可重复选择）。重复执行直至英雄体力小于等于(0)。
- 求期望操作轮数。
- 数据组数(le10^3,ple nle1500,m,kle10^9)
一个前置定义
设(P_i)表示一轮操作的(k)次选择中英雄被选中(i)次的概率。
那么选中英雄的选择有(C_k^i)种组合方案，除了这(i)个以外每个选择都有(m)种选法，因此就有：
[P_i=frac{C_k^i imes m^{k-i}}{(m+1)^k}
]
然而(k)的范围是(10^9)，显然不可能预处理组合数。而有用的(P_i)实际上只有(iin[0,n])。
因此我们考虑把这个式子展开以递推，得到：
[P_0=frac{m^n}{(m+1)^n}\
P_i=P_{i-1} imes frac{k-(i-1)}{i} imes frac1m(ige 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}) imes f_j
]
特殊地，对于(j=i+1)的时候有(f_i=frac1{m+1}P_0 imes f_{i+1})，对于(i=n)的时候有(f_n=P_{n-j} imes f_j)。
显然这时候的(DP)方程是成环的，套路地考虑高斯消元。
然而，我们发现，由于(nle1500)，(O(n^3))会直接(T)飞！
特殊的高斯消元
实际上，这题的矩阵是非常特殊的。
容易发现，对于第(i)行，初始只有第(1sim i+1)个元素有值。
当我们消到第(i)行的时候，第(1sim 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
",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
",v[p]);//高斯消元后输出f[p]
}return 0;
}