### [bzoj5292][Bjoi2018]治療之雨【概率與期望】【高斯消元】

【題目鏈接】
https://www.lydsy.com/JudgeOnline/problem.php?id=5292
【題解】
首先，如果 $k=0$$k=0$ $k=1,m=0,p>1$$k=1,m=0,p>1$時答案為-1。
記 ${f}_{i}$$f_{i}$表示當前剩余生命值為 $i$$i$的期望步數。 $nu{m}_{i}$$num_{i}$為一次行動對英雄造成 $i$$i$點傷害的概率。
那么有：
${f}_{i}\left(i$f_{i}(i
${f}_{n}=1+\underset{j=1}{\overset{i}{\sum }}{f}_{j}\ast nu{m}_{i-j}$$f_{n}=1+\sum_{j=1}^{i}f_{j}*num_{i-j}$
顯然這個方程組可以用高斯消元 $O\left({N}^{3}\right)$$O(N^3)$求解。
觀察高斯消元的矩陣，第i行只有前面i+1個數有值。所以可以把矩陣消成下三角的形式。從下往上做，每次只要與上一行做一次消元。
時間復雜度 $O\left({N}^{2}\right)$$O(N^2)$
【代碼】

# include <bits/stdc++.h>
# define ll long long
# define inf 0x3f3f3f3f
# define N 1501
using namespace std;
int tmp = 0, fh = 1; char ch = getchar();
while (ch < '0' || ch > '9') {if (ch == '-') fh = -1; ch = getchar(); }
while (ch >= '0' && ch <= '9'){tmp = tmp * 10 + ch - '0'; ch = getchar(); }
return tmp * fh;
}

const int P = 1e9 + 7;
int A[N][N], B[N], num[N], inv[N], lim, p, m, k, mul[N];
int power(int x, int y){
int i = x; x = 1;
while (y > 0){
if (y % 2 == 1) x = 1ll * x * i % P;
i = 1ll * i * i % P;
y /= 2;
}
return x;
}
int C(int n, int m){
return 1ll * mul[m] * inv[m] % P;
}
void guass(int n){
for (int i = n; i >= 2; i--){
int inv = power(A[i][i], P - 2);
if (A[i][i] == 0) continue;
int tmp = 1ll * A[i - 1][i] * inv % P;
for (int k = 1; k <= n; k++)
A[i - 1][k] = (A[i - 1][k] - 1ll * A[i][k] * tmp) % P;
B[i - 1] = (B[i - 1] - 1ll * B[i] * tmp) % P;
}
for (int i = 1; i <= n; i++){
for (int j = 1; j < i; j++) B[i] = (B[i] - 1ll * B[j] * A[i][j]) % P;
B[i] = 1ll * B[i] * power(A[i][i], P - 2) % P;
}
}
int main(){
inv[0] = 1;
for (int i = 1; i < N; i++)
inv[i] = 1ll * power(i, P - 2) * inv[i - 1] % P;
for (int opt = read(); opt--;){
if (k == 0 || (m == 0 && k == 1)){
puts("-1");
continue;
}
for (int i = 1; i <= lim; i++)
for (int j = 1; j <= lim; j++)
A[i][j] = 0;
mul[0] = 1;
for (int i = 1; i <= lim; i++)
mul[i] = 1ll * mul[i - 1] * (k - i + 1) % P;
int tot = power(m + 1, k), invtot = power(tot, P - 2), invnum = power(m + 1, P - 2);
int numm = 1ll * m * invnum % P, num1 = invnum;
if (m == 0) numm = 0, num1 = 1;
memset(num, 0, sizeof(num));
for (int i = 0; i <= min(k , lim); i++)
num[i] = 1ll * C(k, i) * power(m, k - i) % P * invtot % P;
for (int i = 1; i <= lim; i++){
B[i] = 1;
if (i != lim){
for (int j = 1; j <= i; j++)
A[i][j] = (A[i][j] - 1ll * num[i - j] * numm) % P;
for (int j = 1; j <= i + 1; j++)
A[i][j] = (A[i][j] - 1ll * num[i - j + 1] * num1) % P;
}
else {
for (int j = 1; j <= i; j++)
A[i][j] = (A[i][j] - num[i - j]) % P;
}
A[i][i] = A[i][i] + 1;
}
guass(lim);
printf("%d\n", (B[p] + P) % P);
}
return 0;
}