hdu 4089 Activation 期望dp

hdu 4089 Activation //期望dp

---

题意

给你一个队列，最开始有n个人，你前面有m-1个人。
现在每个时间点都会发生一些事情。
1. 队首的人有p1的概率继续停留在队首。
2. 队首的人有p2的概率离开队首，进入队尾。
3. 队首的人有p3的概率离开队列。
4. 队列有p4的概率会崩溃。
求队列崩溃时，你前面的人数小于等于k-1的概率。
n,m,k < 2000.

---

简要思路

令$dp[i][j]$表示队列中人数为$i$个，排在第$j$位的期望次数。
那么$ans={\sum }_{i=1}^n{\sum }_{j=1}^{k}dp[i][j]\cdot p4$
考虑转移
$dp[i][j]=dp[i+1][j+1]\cdot p3+dp[i][j+1]\cdot p2+dp[i][j]\cdot p1$

特别的对$i=j$有
$dp[i][i]=dp[i+1][i+1]\cdot p3+dp[i][i]\cdot p1+dp[i][1]\ast p2$

对$i=n,j=m$时
$dp[n][m]=dp[n][m+1]\cdot p2+dp[n][m]\cdot p1+1$
至此得到一个$n^2$个方程。
对其高斯消元即可解出。
但是$n\le 2000$ 那么这样的复杂度为$n^6$ 。肯定不行。
注意到，上述转移式子$i=n$ 时
$dp[n][j]=dp[n][j+1]\cdot p2+dp[n][j]\cdot p1$

$dp[n][j]=\frac{dp[n][j+1]\cdot p2}{1-p1}$
展开：
$$\begin{gathered} dp[n][1]=dp[n][2]\cdot p2/(1-p1) \\ dp[n][2]=dp[n][3]\cdot p2/(1-p1) \\ dp[n][3]=dp[n][4]\cdot p2/(1-p1) \\ ... \\ dp[n][n]=dp[n][1]\cdot p2/(1-p1) \end{gathered}$$
可以手动模拟在$O(n)$的时间内解出$dp[n][i]$。

那么 $dp[n-1][j]=dp[n][[j+1]\cdot p3+dp[n-1][j+1]\cdot p2+dp[n-1][j]\cdot p1$
$dp[n][j+1]$ 已经解出。那么同样的可以上述方法解出$dp[n-1][i]$。
到这里此问题已经解出来啦，总的复杂度为$O(n^2)$。

---

注意：
$n=1$的情况要特殊处理。
注意$p1=1,p2=0$时的处理。

---

代码

#include <bits/stdc++.h>

using namespace std;
#define fi first
#define se second
#define SQ(x) ((x)*(x))
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
const ll inf = 1e18;
const int maxn = 2e3 + 100;
const double eps = 1e-9;
double tp[maxn];
double G(double a[], double b[], double c[], int n) {
    for (int i = 1; i <= n; ++i) tp[i] = c[i];
    double a0 = a[1], b0 = b[1], c0 = c[1];
    for (int i = 2; i <= n - 1; ++i) {
        double temp = -b0 / a[i];
        b0 = b[i] * temp, c0 = c0 + c[i] * temp;
    }
    double temp = -a0 / a[n];
    c0 = (c[n] * temp + c0) / (b[n] * temp + b0);
    c[n] = c0;
    for (int i = n - 1; i >= 1; --i) {
        c[i] = (c[i] - b[i] * c[i + 1]) / a[i];
    }
    for (int i = 1; i < n; ++i) {
        if (fabs(b[i]) < eps) c[i] = tp[i] / a[i];
    }
    if (fabs(a[n]) < eps) c[n] = tp[n] / b[n];
}
double dp[2][maxn], a[maxn], b[maxn], c[maxn];

int main() {
//    freopen("1.in", "r", stdin);
//    freopen("1.out", "w", stdout);
    int n, m, k;
    double p1, p2, p3, p4;
    cout << fixed << setprecision(5);
    while (cin >> n >> m >> k >> p1 >> p2 >> p3 >> p4) {
        int rt = 0;
        memset(dp[rt], 0, sizeof dp[rt]);
        if (fabs(p1 - 1) < eps || fabs(p1 + p2 - 1) < eps) {
            cout << 0.0 << '\n';
            continue;
        }
        double ans = 0;
        k = min(k, n);
        for (int i = n; i >= 2; --i) {
            rt ^= 1;
            memset(dp[rt], 0, sizeof dp[rt]);
            for (int j = 1; j <= i; ++j) {
                a[j] = 1, b[j] = -p2 / (1 - p1), c[j] = dp[rt ^ 1][j + 1] * p3 / (1 - p1);
            }
            if (i == n) c[m] += 1 / (1 - p1);
            swap(a[i], b[i]);
            G(a, b, c, i);
            for (int j = 1; j <= i; ++j) dp[rt][j] = c[j];
            for (int j = 1; j <= k; ++j) ans += dp[rt][j];
        }
        if (n != 1) dp[rt][1] = dp[rt][2] * p3 / (1 - p1 - p2);
        else dp[rt][1] = 1 / (1 - p1 - p2);
        ans += dp[rt][1];
        ans *= p4;
        if (fabs(p4) < eps) ans = 0;
        ans = abs(ans);
        cout << ans << '\n';
    }
    return 0;
}