[bzoj5292][Bjoi2018]治疗之雨【概率与期望】【高斯消元】

【题目链接】
　　https://www.lydsy.com/JudgeOnline/problem.php?id=5292
【题解】
　　首先，如果$k=0$或$k=1,m=0,p>1$时答案为-1。
　　记$f_{i}$表示当前剩余生命值为$i$的期望步数。$num_{i}$为一次行动对英雄造成$i$点伤害的概率。
　　那么有：
　　$f_{i}(i<n)=1+m/(m+1)*(\sum_{j=1}^{i}f_{j}*num_{i-j})+1/(m+1)*(\sum_{j=1}^{i+1}fj*num_{i-j+1})$
　　$f_{n}=1+\sum_{j=1}^{i}f_{j}*num_{i-j}$
　　显然这个方程组可以用高斯消元$O(N^3)$求解。
　　观察高斯消元的矩阵，第i行只有前面i+1个数有值。所以可以把矩阵消成下三角的形式。从下往上做，每次只要与上一行做一次消元。
　　时间复杂度$O(N^2)$
【代码】

# include <bits/stdc++.h>
# define    ll      long long
# define    inf     0x3f3f3f3f
# define    N       1501
using namespace std;
int read(){
    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--;){
        lim = read(), p = read(), m = read(), k = read();
        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;
}