[BZOI2018]治疗之雨

题面：Luogu
题解：期望dp+gauss消元
先吐槽一下这个题面，根本看不懂，看讨论区半天才发现，另外m个就是用来降低概率的
设\(prob_i\)表示一轮降了\(i\)点的概率

\[Prob_i=\frac{\binom{k}{i}m^{k-i}}{(m+1)^k} \]

总共\((m+1)^{k}\)种方案，选\(i\)次降一点，剩下的不降
由于\(k\)很大，所以\(\binom{k}{i}\)不能通过预处理阶乘算出，而要递推（否则会#6~#10RE）
递推方法很多，我的是

\[prob_i=\frac{1}{i!}*m^{k-i}*\frac{k!}{(k-i)!}*\frac{1}{(m+1)^k} \]

设\(Q_{x,y}\)表示一轮操作从\(x\)点变成\(y\)点

\[Q_{x,y}=\begin{split} \left\{ \begin{array}{lr} 0 & x=n \&\& y=n+1 &\\ Prob_{x-y} & x=n\\ \frac{1}{m+1}Prob_0 & y=x+1\\ \frac{m}{m+1}Prob_{x-y}+\frac{1}{m+1}Prob_{x-y+1} & \text{otherwise} \end{array} \right. \end{split} \]

设\(dp_i\)表示从\(i\)点变成\(0\)的期望步数

\[\begin{split} & dp_i=1+\sum_{j=1}^{i+1} Q_{i,j}dp_j & (1\le i<n) \\ & dp_n=1+\sum_{j=1}^nQ_{n,j}dp_j \end{split} \]

于是可以推一下式子

\[dp_{i+1}=\frac{dp_i-1-\sum_{j=1}^{i}{Q_{i,j}dp_{j}}}{Q_{i,i+1}} (1\le i<n) \]

可以发现，我们只需求出\(dp_1\)即可求出其它
然后是两种办法
1.高斯消元
矩阵大概是这样

\[\left[ \begin{matrix} 1-Q_{1,1} & -Q_{1,2} & 0 &\cdots & 0\\ -Q_{2,1} & 1-Q_{2,2} & -Q_{2,3} &\cdots & 0 \\ \vdots & \vdots & \vdots &\ddots & \vdots \\ -Q_{n-1,1} & -Q_{n-1,2} & -Q_{n-1,3} &\cdots & -Q_{n-1,n} \\ -Q_{n,1} & -Q_{n,2} & -Q_{n,3} &\cdots & 1-Q_{n,n} \end{matrix} \right]* \left[ \begin{matrix} dp_{1}&\\ dp_{2}&\\ \vdots&\\ dp_{n-1}&\\ dp_{n}& \end{matrix} \right]= \left[ \begin{matrix} 1&\\ 1&\\ \vdots&\\ 1&\\ 1& \end{matrix} \right] \]

然后可以通过一些手法把它搞成只要\(n^2\)的

2.设\(dp_1=X\)
则每一项都可以表示为\(A[i]*X+B[i]\)
然后通过\(dp_n\)的递推式也能得到一个\(dp_n=An*X+Bn\)
所以

\[A[n]*X+B[n]=An*X+Bn \\ X=\frac{Bn-B[n]}{A[n]-An} \]

然后带回去就可以了
我也是写的这种方法
还有一件令我迷惑的事：prob每次不清零会WA#1~#5，但是我写的每次应该都会重新算一遍，在没被改变的地方应该也没有被调用

#include<bits/stdc++.h>
using namespace std;
template<typename T>
inline void read(T& x)
{
	x = 0; char c = getchar();
	while (!isdigit(c)) c = getchar();
	while (isdigit(c)) x = x * 10 + c - '0', c = getchar();
}
#define maxn 1505
#define P 1000000007
int fac[maxn], inv[maxn], prob[maxn];
int A[maxn], B[maxn], An, Bn, X;
int n, p, m, k, invM;
inline int mul(int x, int y) { return 1ll * x * y % P; }
inline int qpow(int x, int y)
{
	int ans = 1;
	for (; y; x = mul(x, x), y >>= 1) if (y & 1) ans = mul(ans, x);
	return ans;
}
inline int getp(int x, int y)//Q_{x,y}
{
	if (x == n && y == n + 1) return 0;
	if (x == n) return prob[x - y];
	if (y == x + 1) return mul(invM, prob[0]);
	return (0ll + mul(mul(m, prob[x - y]), invM) + mul(invM, prob[x - y + 1])) % P;
}
inline void exgcd(int a, int b, int& x, int& y)
{
	if (!b) { x = 1, y = 0; return; }
	exgcd(b, a % b, y, x);
	y -= a / b * x;
}
inline int getinv(int a)
{
	int x, y;
	exgcd(a, P, x, y);
	x %= P;
	return (x < 0) ? x + P : x;
}
inline void init(int maxnum = 1500)
{
	fac[0] = inv[0] = 1;
	for (int i = 1; i <= maxnum; ++i) fac[i] = 1ll * fac[i - 1] * i % P;
	inv[maxnum] = getinv(fac[maxnum]);
	for (int i = maxnum - 1; i; --i) inv[i] = 1ll * inv[i + 1] * (i + 1) % P;
}
inline void getprob()
{
	invM = getinv(m + 1);
	int inv1 = qpow(invM, k);
	for (int i = 0, tp = qpow(m, k), tp2 = getinv(m), tpm = 1; i <= min(n, k); ++i)
		prob[i] = mul(mul(mul(inv[i], tpm), tp), inv1), tp = mul(tp, tp2), tpm = mul(tpm, k - i);
	//tp=m^{k-i},inv[i]=\frac{1}{i!},tpm=\frac{k!}{(k-i)!},inv1=\frac{1}{(m+1)^k}
}
inline void solve()
{
	A[1] = 1, B[1] = 0; An = 0, Bn = 1;
	for (int i = 2; i <= n; ++i)
	{
		A[i] = A[i - 1], B[i] = B[i - 1] - 1;
		for (int j = 1, tp; j < i; ++j)
		{
			tp = getp(i - 1, j);
			A[i] = (0ll + A[i] - mul(tp, A[j])) % P;
			B[i] = (0ll + B[i] - mul(tp, B[j])) % P;
		}
		if (A[i] < 0) A[i] += P;
		if (B[i] < 0) B[i] += P;
		int tp = getinv(getp(i - 1, i));
		A[i] = mul(A[i], tp), B[i] = mul(B[i], tp);
	}
	for (int i = 1, tp; i <= n; ++i)
	{
		tp = getp(n, i);
		An = (0ll + An + mul(tp, A[i])) % P;
		Bn = (0ll + Bn + mul(tp, B[i])) % P;
	}
	X = mul((Bn - B[n] + P) % P, getinv((A[n] - An + P) % P));
	printf("%d\n", (mul(A[p], X) + B[p]) % P);
	memset(prob, 0, sizeof(prob));
}
int main()
{
	int T; read(T);
	init();
	while (T--)
	{
		read(n), read(p), read(m), read(k);
		if (!k || (!m && k == 1)) { puts("-1"); continue; }//特判
		if (!m)
		{
			int ans = 0;
			while (p > 0)
			{
				if (p < n) ++p;
				p -= k, ++ans;
			}
			printf("%d\n", ans);
			continue;
		}
		getprob();
		solve();
	}
	return 0;
}