HDU 4089 Activation 概率DP 2011年区域赛北京现场赛I题

题目大意：

现在Tomato在等待进入游戏Chinese Paladin, 初始的时候他站在长度为n的队伍的第m个位置

每一秒都可能发生以下事件中的一个：

p1的概率队伍保持不变, p2的概率队伍的第一个人调到队尾, p3的概率队伍的第一个人进入游戏, p4的概率服务器崩了

现在问在服务器崩掉的时候, Tomato站在队伍的前k个位置的概率

大致思路：

计算过程见代码

Result  :  Accepted     Memory  :  31948 KB     Time  :  2500 ms

/*
 * Author: Gatevin
 * Created Time:  2014/12/1 14:19:23
 * File Name: Asuna.cpp
 */
#include<iostream>
#include<sstream>
#include<fstream>
#include<vector>
#include<list>
#include<deque>
#include<queue>
#include<stack>
#include<map>
#include<set>
#include<bitset>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cctype>
#include<cmath>
#include<ctime>
#include<iomanip>
using namespace std;
const double eps(1e-8);
typedef long long lint;

int n, m, k;
double p1, p2, p3, p4;
double dp[2010][2010];

/*
 * 刷概率dp的题又遇到这题,发现之前的理解还是存在问题,现在理解更深刻了
 */

/*
 * 用dp[i][j]表示当前队伍长度为i, 站在位置j时，到达目标状态的概率
 * 显然有：
 * dp[i][j] = p1*dp[i][j] + p2*dp[i][j - 1] + p3*dp[i - 1][j - 1] + p4 (1 < j <= k)
 * dp[i][j] = p1*dp[i][j] + p2*dp[i][j - 1] + p3*dp[i - 1][j - 1] (1 <= k < j)
 * dp[i][1] = p1*dp[i][1] + p2*dp[i][i] + p4 (j = 1)
 * 首先可以发现方程可以化简为:
 * dp[i][j] = p21*dp[i][j - 1] + p31*dp[i - 1][j - 1] + p41 (1 < j <= k) (1)
 * dp[i][j] = p21*dp[i][j - 1] + p31*dp[i - 1][j - 1] (1 <= k < j)       (2)
 * dp[i][1] = p21*dp[i][i] + p41 (j = 1)                                 (3)
 * 其中 p21 = p2/(1 - p2), p31 = p3/(1 - p2), p41 = p4/(1 - p4);
 * 那么由式子(3), dp[1][1] = p41/(1 - p21) = p4/(1 - p1 - p2)
 * 对于i >= 2的式子，当i <= k时， 由(1)式和(3)式可以得到一个共有dp[i][1~i]， i个未知数的方程(dp[i - 1][1 ~ i - 1]视为已知)
 * 对于i > k 由(3)迭代至第二维的数小于k再由(1)式递推即可
 * 这里给出i <= k的迭代计算过程:
 * dp[i][1] = p21*dp[i][i] + p41 (3)
 *          = p21(dp[i][i - 1] + p31*dp[i - 1][i - 1] + p41) + p41 (2)
 *          = p21*(p21*dp[i][i - 2] + p31*dp[i - 1][i - 2] + p41) + p21*p31*dp[i - 1][i - 1] + p21*p41 + p41
 *          = ... 继续迭代
 *          = p41*(1 + p21 + p21^2 + p21^3 + .. + p21^(i - 1)) + (p21^i)*dp[i][1] + p31*(p21*dp[i - 1][i - 1] + .. + (p21^(i - 1))*dp[i - 1][1])
 *          = (p41*(1 - (p21^i))/(1 - p21) + p31*(p21*dp[i - 1][i - 1] + .. + (p21^(i - 1))*dp[i - 1][1])) + (p21^i)*dp[i][1]
 * 所以 dp[i][1] = (p41*(1 - (p21^i))/(1 - p21) + p31*(p21*dp[i - 1][i - 1] + .. + (p21^(i - 1))*dp[i - 1][1]))/(1 - (p21^i))
 * 同样的方法，当 i > k时
 * dp[i][1] = ((1 - p21 + p21^(i - k + 1) - (p21^i))/(1 - p21) + p31*(p21*dp[i - 1][i - 1] + .. + (p21^(i - 1))*dp[i - 1][1]))/(1 - (p21^i))
 * 然后利用(1)或(2)式递推算出dp[i][2~i]即可
 * 依次递推出dp[n][m]即为结果
 */

int main()
{
    while(~scanf("%d %d %d %lf %lf %lf %lf", &n, &m, &k, &p1, &p2, &p3, &p4))
    {
        if(p4 < eps)
        {
            printf("0.00000\n");//注意可能计算中途出现nan的情况
            continue;
        }
        memset(dp, 0, sizeof(dp));
        dp[1][1] = p4/(1 - p1 - p2);
        double p21 = p2/(1 - p1);
        double p31 = p3/(1 - p1);
        double p41 = p4/(1 - p1);
        for(int i = 2; i <= n; i++)
        {
            for(int j = 1; j < i; j++)
                dp[i][1] += dp[i - 1][j]*pow(p21, i - j);
            dp[i][1] *= p31;
            if(i <= k) dp[i][1] += p41*(1 - pow(p21, i))/(1 - p21);
            else dp[i][1] += p41*(1 - p21 + pow(p21, i - k + 1) - pow(p21, i))/(1 - p21);
            dp[i][1] /= (1 - pow(p21, i));
            for(int j = 2; j <= i; j++)
                dp[i][j] = p21*dp[i][j - 1] + p31*dp[i - 1][j - 1] + (j <= k ? p41 : 0);
        }
        printf("%.5f\n", dp[n][m]);
    }
    return 0;
}