HDU4089概率dp解题报告

HDU4089概率dp解题报告

K - Activation
Time Limit:10000MS Memory Limit:32768KB 64bit IO Format:%I64d & %I64u
Description
After 4 years’ waiting, the game “Chinese Paladin 5” finally comes out. Tomato is a crazy fan, and luckily he got the first release. Now he is at home, ready to begin his journey.
But before starting the game, he must first activate the product on the official site. There are too many passionate fans that the activation server cannot deal with all the requests at the same time, so all the players must wait in queue. Each time, the server deals with the request of the first player in the queue, and the result may be one of the following, each has a probability:
1. Activation failed: This happens with the probability of p1. The queue remains unchanged and the server will try to deal with the same request the next time.
2. Connection failed: This happens with the probability of p2. Something just happened and the first player in queue lost his connection with the server. The server will then remove his request from the queue. After that, the player will immediately connect to the server again and starts queuing at the tail of the queue.
3. Activation succeeded: This happens with the probability of p3. Congratulations, the player will leave the queue and enjoy the game himself.
4. Service unavailable: This happens with the probability of p4. Something just happened and the server is down. The website must shutdown the server at once. All the requests that are still in the queue will never be dealt.
Tomato thinks it sucks if the server is down while he is still waiting in the queue and there are no more than K-1 guys before him. And he wants to know the probability that this ugly thing happens.
To make it clear, we say three things may happen to Tomato: he succeeded activating the game; the server is down while he is in the queue and there are no more than K-1 guys before him; the server is down while he is in the queue and there are at least K guys before him.
Now you are to calculate the probability of the second thing.

Input
There are no more than 40 test cases. Each case in one line, contains three integers and four real numbers: N, M (1 <= M <= N <= 2000), K (K >= 1), p1, p2, p3, p4 (0 <= p1, p2, p3, p4 <= 1, p1 + p2 + p3 + p4 = 1), indicating there are N guys in the queue (the positions are numbered from 1 to N), and at the beginning Tomato is at the Mth position, with the probability p1, p2, p3, p4 mentioned above.

Output
A real number in one line for each case, the probability that the ugly thing happens.
The answer should be rounded to 5 digits after the decimal point.

Sample Input
2 2 1 0.1 0.2 0.3 0.4
3 2 1 0.4 0.3 0.2 0.1
4 2 3 0.16 0.16 0.16 0.52

Sample Output
0.30427
0.23280
0.90343
题目大意：有个人想玩仙5，但是需要激活。
在官网上，目前一共有n个人在排队，Tomato排在第N位。给你一个数k。
然后有4种事件可能发生
事件1：激活失败，然后会直接尝试再次激活，概率为p1
事件2：失去连接，这时候，当前激活者将失去排队，然后马上连接上，重新排队，概率为p2
事件3：激活成功，当前激活者出队，开玩，概率为p3。
事件4：官网挂掉，所有排队者退出，概率为p4
问题就是求出，tomato前面至少有k-1个人的时候官网挂掉的概率

分析：
由题建立状态，dp[i][j]表示当前一共i个人在排队，tomato在第j个位置。分析可知，只有当tomato排在第一个的时候，才可能被扔到队尾
那么，且他不能被成功激活。
dp[i][1] = p1*dp[i][1] + p2*dp[i][i] + p4;
当他不在队首的时候，且前面至少有k-1个人，可以由以下转移
dp[i][j] = p1 * dp[i][j] + p2*dp[i][j-1] + p3*dp[i-1][j-1] + p4;
当他前面有大于k-1个人的时候：
dp[i][j] = p1 * dp[i][j] + p2*dp[i][j-1] + p3*dp[i-1][j-1];
(这时候不能当机，所以不能加p4)
整理以上式子
我们设
p22 = p2 / (1 - p1)
p33 = p2 / (1 - p1)
p44 = p4 / (1 - p1)
然后得到
dp[i][1] = p22 * dp[i][i] + p44
当j <= k
dp[i][j] = p22 * dp[i][j-1] + p33 * dp[i-1][j-1] + p44
否则
dp[i][j] = p22 * dp[i][j-1] + p33 * dp[i-1][j-1]
至此如果你就去写代码的话，恭喜你 你会T
在计算机，小数位数很多，概率很接近0的时候，计算的代价是可怕的！
所以我们开一个c数组，来保存上次计算的值
有c[0] = 1;
c[1] = p44;
当j <= k
c[i] = p33 * dp[i-1][j-1] + p44
当j > k
c[i] = p33 * dp[i-1][j-1]
然后原式变成
dp[i][1] = p22 * dp[i][i] + c[1]
……
dp[i][j] = p22 * dp[i][j-1] + c[j]
然后我们发现，这尼玛能解？
当然，可以用高斯消元，但是那玩意儿是n^3的肯定不行
然后我们发现 可以简单的迭代
把第1个式子带入第二个式子
在把第2个式子带入第三个式子
…
一直带入下去，我们发现
最后变成了: (p22简写为p, dp[i][j]简写为a[j])
a[i] = a[i] * p^i + c[1] * p^(i-1) + c[2] * p^(i-2) + … + p * c[i-1] + c[i]
然后化简下就可以求出a[i]
然后算出a[1]
然后根据方程依此算出a[2] ~ a[i-1]
至此如果你交，还是会T
必须特判下p4 接近 0的情况，不然计算代价奇高
最后还的循环数组优化，防止爆空间
能ac真尼玛不容易…
具体见代码

////  Created by Running Photon on 2015-08-27//  Copyright (c) 2015 Running Photon. All rights reserved.////#include <algorithm>#include <cctype>#include <cmath>#include <cstdio>#include <cstdlib>#include <cstring>#include <iomanip>#include <iostream>#include <map>#include <queue>#include <string>#include <sstream>#include <set>#include <vector>#include <stack>#define ALL(x) x.begin(), x.end()#define INS(x) inserter(x, x,begin())#define ll long long#define CLR(x) memset(x, 0, sizeof x)#define MAXN 9999#define MAXSIZE 10#define DLEN 4using namespace std;const int inf = 0x3f3f3f3f;const int MOD = 1e9 + 7;const int maxn = 1e6 + 10;const int maxv = 2e3 + 10;const double eps = 1e-5;inline int read() {    char c = getchar();    int f = 1;    while(!isdigit(c)) {    if(c == '-') f = -1;    c = getchar();    }    int x = 0;    while(isdigit(c)) {        x = x * 10 + c - '0';        c = getchar();    }    return x * f;}//int L[maxv], R[maxv], l[maxv], r[maxv];//vector <int> xx;//void iscretization() {//  s = 0;//  xx.push_back(0);//  sort(ALL(xx));//  xx.resize(unique(ALL(xx)) - xx.begin());//  for(int i = 0; i < n; i++) {//      l[i] = lower_bound(ALL(xx), L[i]) - xx.begin();//      r[i] = lower_bound(ALL(xx), R[i]) - xx.begin();//  }//  t = xx.size();//}int n, m, k;double p1, p2, p3, p4;double c[maxv];double dp[2][maxv];double p[maxv];int main() {#ifdef LOCAL    freopen("in.txt", "r", stdin);//  freopen("out.txt","w",stdout);#endif    while(scanf("%d%d%d%lf%lf%lf%lf", &n, &m, &k, &p1, &p2, &p3, &p4) != EOF) {        double p22 = p2 / (1.0 - p1);        double p33 = p3 / (1.0 - p1);        double p44 = p4 / (1.0 - p1);        if(p4 < eps) {            puts("0.00000");            continue;        }        CLR(dp);        CLR(p);        p[0] = 1.0;        for(int i = 1; i <= n; i++) p[i] = p22 * p[i-1];        dp[1][1] = p44 / (1.0 - p22);        for(int i = 2; i <= n; i++) {            for(int j = 1; j <= i; j++) {                if(j <= k) c[j] = p33 * dp[(i-1) & 1][j-1] + p44;                else c[j] = p33 * dp[(i-1) & 1][j-1];            }            double tmp = 0;            for(int j = 1; j <= i; j++) {                tmp += p[i-j] * c[j];            }            dp[i&1][i] = tmp / (1.0 - p[i]);            dp[i&1][1] = p22 * dp[i&1][i] + p44;            for(int j = 2; j < i; j++) {                dp[i&1][j] = p22 * dp[i&1][j-1] + c[j];            }        }        printf("%.5f\n", dp[n&1][m]);    }    return 0;}