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


【題目鏈接】
  https://www.lydsy.com/JudgeOnline/problem.php?id=5292
【題解】
  首先,如果 k = 0 k = 1 , m = 0 , p > 1 時答案為-1。
  記 f i 表示當前剩余生命值為 i 的期望步數。 n u m i 為一次行動對英雄造成 i 點傷害的概率。
  那么有:
   f i ( i < n ) = 1 + m / ( m + 1 ) ( j = 1 i f j n u m i j ) + 1 / ( m + 1 ) ( j = 1 i + 1 f j n u m i j + 1 )
   f n = 1 + j = 1 i f j n u m 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;
}

注意!

本站转载的文章为个人学习借鉴使用,本站对版权不负任何法律责任。如果侵犯了您的隐私权益,请联系我们删除。



 
粤ICP备14056181号  © 2014-2020 ITdaan.com