hdu 4089 Activation //期望dp
题意
给你一个队列,最开始有n个人,你前面有m-1个人。
现在每个时间点都会发生一些事情。
1. 队首的人有p1的概率继续停留在队首。
2. 队首的人有p2的概率离开队首,进入队尾。
3. 队首的人有p3的概率离开队列。
4. 队列有p4的概率会崩溃。
求队列崩溃时,你前面的人数小于等于k-1的概率。
n,m,k < 2000.
简要思路
令
d
p
[
i
]
[
j
]
dp[i][j]
dp[i][j]表示队列中人数为
i
i
i个,排在第
j
j
j位的期望次数。
那么
a
n
s
=
∑
i
=
1
n
∑
j
=
1
k
d
p
[
i
]
[
j
]
⋅
p
4
ans = \sum_{i=1}^n\sum_{j =1}^kdp[i][j] \cdot p4
ans=∑i=1n∑j=1kdp[i][j]⋅p4
考虑转移
d
p
[
i
]
[
j
]
=
d
p
[
i
+
1
]
[
j
+
1
]
⋅
p
3
+
d
p
[
i
]
[
j
+
1
]
⋅
p
2
+
d
p
[
i
]
[
j
]
⋅
p
1
dp[i][j]= dp[i+1][j+1]\cdot p3 + dp[i][j+1]\cdot p2 + dp[i][j]\cdot p1
dp[i][j]=dp[i+1][j+1]⋅p3+dp[i][j+1]⋅p2+dp[i][j]⋅p1
特别的对
i
=
j
i=j
i=j有
d
p
[
i
]
[
i
]
=
d
p
[
i
+
1
]
[
i
+
1
]
⋅
p
3
+
d
p
[
i
]
[
i
]
⋅
p
1
+
d
p
[
i
]
[
1
]
∗
p
2
dp[i][i]= dp[i+1][i+1]\cdot p3 + dp[i][i]\cdot p1 + dp[i][1]*p2
dp[i][i]=dp[i+1][i+1]⋅p3+dp[i][i]⋅p1+dp[i][1]∗p2
对
i
=
n
,
j
=
m
i=n,j=m
i=n,j=m时
d
p
[
n
]
[
m
]
=
d
p
[
n
]
[
m
+
1
]
⋅
p
2
+
d
p
[
n
]
[
m
]
⋅
p
1
+
1
dp[n][m] = dp[n][m+1]\cdot p2 + dp[n][m]\cdot p1 +1
dp[n][m]=dp[n][m+1]⋅p2+dp[n][m]⋅p1+1
至此得到一个
n
2
n^2
n2个方程。
对其高斯消元即可解出。
但是
n
≤
2000
n\leq 2000
n≤2000 那么这样的复杂度为
n
6
n^6
n6 。肯定不行。
注意到,上述转移式子
i
=
n
i=n
i=n 时
d
p
[
n
]
[
j
]
=
d
p
[
n
]
[
j
+
1
]
⋅
p
2
+
d
p
[
n
]
[
j
]
⋅
p
1
dp[n][j] = dp[n][j+1]\cdot p2 + dp[n][j]\cdot p1
dp[n][j]=dp[n][j+1]⋅p2+dp[n][j]⋅p1
d
p
[
n
]
[
j
]
=
d
p
[
n
]
[
j
+
1
]
⋅
p
2
1
−
p
1
dp[n][j] = \frac{dp[n][j+1]\cdot p2}{1-p1}
dp[n][j]=1−p1dp[n][j+1]⋅p2
展开:
d
p
[
n
]
[
1
]
=
d
p
[
n
]
[
2
]
⋅
p
2
/
(
1
−
p
1
)
d
p
[
n
]
[
2
]
=
d
p
[
n
]
[
3
]
⋅
p
2
/
(
1
−
p
1
)
d
p
[
n
]
[
3
]
=
d
p
[
n
]
[
4
]
⋅
p
2
/
(
1
−
p
1
)
.
.
.
d
p
[
n
]
[
n
]
=
d
p
[
n
]
[
1
]
⋅
p
2
/
(
1
−
p
1
)
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)\\
dp[n][1]=dp[n][2]⋅p2/(1−p1)dp[n][2]=dp[n][3]⋅p2/(1−p1)dp[n][3]=dp[n][4]⋅p2/(1−p1)...dp[n][n]=dp[n][1]⋅p2/(1−p1)
可以手动模拟在
O
(
n
)
O(n)
O(n)的时间内解出
d
p
[
n
]
[
i
]
dp[n][i]
dp[n][i]。
那么
d
p
[
n
−
1
]
[
j
]
=
d
p
[
n
]
[
[
j
+
1
]
⋅
p
3
+
d
p
[
n
−
1
]
[
j
+
1
]
⋅
p
2
+
d
p
[
n
−
1
]
[
j
]
⋅
p
1
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−1][j]=dp[n][[j+1]⋅p3+dp[n−1][j+1]⋅p2+dp[n−1][j]⋅p1
d
p
[
n
]
[
j
+
1
]
dp[n][j+1]
dp[n][j+1] 已经解出。那么同样的可以上述方法解出
d
p
[
n
−
1
]
[
i
]
dp[n-1][i]
dp[n−1][i]。
到这里此问题已经解出来啦,总的复杂度为
O
(
n
2
)
O(n^2)
O(n2)。
注意:
n
=
1
n=1
n=1的情况要特殊处理。
注意
p
1
=
1
,
p
2
=
0
p1=1,p2=0
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;
}