[BJOI2018]治療之雨

參考博客

https://www.luogu.org/blog/ShadowassIIXVIIIIV/solution-p4457

洛谷 P4457

題目大意

\(T\) 組數據

場上有 \(m\) 個生命值無上下限的隨從與 \(1\) 個生命值上限為 \(n\) 下限為 \(0\) 當前為 \(p\) 的英雄，現在每次隨機選擇一個生命值不為下限的人物，將其生命值 \(+1\) ，然後重復 \(k\) 次，每次將隨機選擇一個生命值不為下限的人物將其生命值 \(-1\) ，問期望多少次後生命值歸 \(0\) ，如果不行則輸出 \(-1\) ，答案對 \(1000000007\) 取模，保證分母不為模數的倍數，不存在 \(n=p=k=1,m=0\)

數據範圍

\(1 \le T \le 100, 1 \le p \le n \le 1500, 0 \le m,k \le 1000000000\)

時空限制

4000ms, 512MB

分析

我研究了一個 \(n^2\) 的叠代方法，但是它讓我沒有興趣寫出來......

我們設 \(P(y)\) 為在 \(k\) 次減 \(1\) 之後，將英雄血量從 \(x\) 變成 \(x - y > 0\) 的概率
\[ P(y) = (\dfrac 1{1 + m})^y(\dfrac m{1 + m})^{k-y} \binom ky \]
設 \(f(x)\) 為 \(p = x\) 時的答案
\[ \begin{cases} f(0) = 0 \f(n) = \sum P(y) f(n - y) \f(x) = \dfrac 1{1 + m} P(0) f(x + 1) + \sum \left[ \dfrac 1{1 + m}P(y + 1) + \dfrac m{1 + m }P(y)\right] f(x - y) \end{cases} \]

那麽 \(O(n^2)\) 該怎麽辦呢（雖然我覺得 \(n^2\) 能過也很神奇）

我們叠代就可以了

觀察一下高斯消元矩陣
\[ \left[ \begin{matrix} 1 & 1 & 0 & 0 & 0 & \cdots & 0 & 0 \1 & 1 & 1 & 0 & 0 & \cdots & 0 & 0\\vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \1 & 1 & 1 & 1 & 1 & \cdots & 1 & 1 \end{matrix} \right] \]

所以我們消一消

\[ \left[ \begin{matrix} 1 & 1 & 0 & 0 & 0 & \cdots & 0 & 0 \0 & 1 & 1 & 0 & 0 & \cdots & 0 & 0\\vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \0 & 0 & 0 & 0 & 0 & \cdots & 0 & 1 \end{matrix} \right] \]

於是就是 \(n^2\) 了

Code

#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
inline char nc() {
    static char buf[100000], *l = buf, *r = buf;
    return l==r&&(r=(l=buf)+fread(buf,1,100000,stdin),l==r)?EOF:*l++;
}
template<class T> void read(T &x) {
    x = 0; int f = 1, ch = nc();
    while(ch<‘0‘||ch>‘9‘){if(ch==‘-‘)f=-1;ch=nc();}
    while(ch>=‘0‘&&ch<=‘9‘){x=x*10-‘0‘+ch;ch=nc();}
    x *= f; 
} 
typedef long long ll;
const int mod = 1000000007;
const int maxn = 1500 + 5;
int T, n, p, m, k, p1, p2;
int P[maxn], e[maxn][maxn];
inline void sub(int &x, int y) {
    x -= y; if(x < 0) x += mod;
}
inline void add(int &x, int y) {
    x += y; if(x >= mod) x -= mod;
}
ll qpow(ll x, ll y) {
    ll re = 1;
    while(y) {
        if(y & 1) re = re * x % mod;
        x = x * x % mod, y >>= 1;
    }
    return re;
}
void init() {
    memset(P, 0, sizeof(P));
    p1 = 1 * qpow(1 + m, mod - 2);
    p2 = m * qpow(1 + m, mod - 2) % mod;
    int C = 1;
    for(int i = 0; i <= min(n, k); ++i) {
        P[i] = qpow(p1, i) * qpow(p2, k - i) % mod * C % mod;
        C = (ll)C * (k - i) % mod * qpow(i + 1, mod - 2) % mod;
    }
}
void Gauss() {
    for(int i = 1; i < n; ++i) {
        int r = qpow(e[i][i], mod - 2);
        for(int j = i + 1; j <= n; ++j) {
            int t = (ll)e[j][i] * r % mod; e[j][i] = 0;
            sub(e[j][i + 1], (ll)e[i][i + 1] * t % mod);
            sub(e[j][n + 1], (ll)e[i][n + 1] * t % mod);
        }
    }
    e[n][n + 1] = (ll)e[n][n + 1] * qpow(e[n][n], mod - 2) % mod, e[n][n] = 1;
    for(int i = n - 1; i >= 1; --i) {
        sub(e[i][n + 1], (ll)e[i][i + 1] * e[i + 1][n + 1] % mod), e[i][i + 1] = 0;
        e[i][n + 1] = (ll)e[i][n + 1] * qpow(e[i][i], mod - 2) % mod, e[i][i] = 1;
    } 
}
int main() {
//  freopen("testdata.in", "r", stdin);
    read(T);
    for(int kase = 1; kase <= T; ++kase) {
        read(n), read(p), read(m), read(k);
        if(k == 0 || (m == 0 && k == 1)) {
            puts("-1");
            continue;
        }
        init();
        memset(e, 0, sizeof(e));
        for(int x = 1; x < n; ++x) {
            for(int y = 0; y < x; ++y) {
                e[x][x - y] = mod - ((ll)p1 * P[y + 1] % mod + (ll)p2 * P[y] % mod) % mod;
            }
            e[x][x + 1] = mod - (ll)p1 * P[0] % mod;
            add(e[x][x], 1), e[x][n + 1] = 1;
        }
        for(int y = 1; y < n; ++y) {
            e[n][n - y] = mod - P[y];
        }
        e[n][n] = ((1 - P[0]) + mod) % mod, e[n][n + 1] = 1;
        Gauss(); printf("%d\n", e[p][n + 1]);
    }
    return 0;
} 

總結

特殊矩陣的高斯消元可以優化

[洛谷 P4457] [BJOI2018]治療之雨