[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} \]
那么我们高斯消元即可。。。然而 \(1500\) 怎么 \(n^3\)

那么 \(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;
} 

总结

特殊矩阵的高斯消元可以优化