【HDU4089】Activation-概率DP好题

测试地址：Activation
题目大意： Tomato要在服务器上激活一个游戏，一开始服务器序列中有N个人，他排在第M位，每次服务器会对序列中第一位的玩家进行激活，有四种结果：
1.有p1的概率会激活失败，这时候序列的状态是不变的。
2.有p2的概率第一位的玩家会连接错误，这时候序列中第一位的玩家会成为最后一位，其他玩家相对位置不变。
3.有p3的概率第一位的玩家激活成功，这时候第一位的玩家会离开序列。
4.有p4的概率服务器崩溃，发生这件事之后所有玩家都不能激活了。
求Tomato遇到服务器崩溃并且在服务器崩溃时处于前K位的概率。
做法：这一道题是概率DP，而且是一道非常好的题目。
题意比较复杂，但是看到概率，我们就应该想到用DP来求概率。设f(i,j)为从序列中有i人，Tomato在第j位的状态开始到达目标状态（遇到服务器崩溃并且在服务器崩溃时处于前K位）的概率，有以下状态转移方程（不懂怎么推出来的可以结合题意理解一下）：
当j=1，f(i,j)=f(i,j)×p1+f(i,i)×p2+p4；
当1<j≤k，f(i,j)=f(i,j)×p1+f(i,j−1)×p2+f(i−1,j−1)×p3+p4；
当k<j≤i，f(i,j)=f(i,j)×p1+f(i,j−1)×p2+f(i−1,j−1)×p3；
注意到有些状态之间是会循环引用的，那是不是就不能用DP来解了呢？不是。事实上，在这些状态中进行概率转移的过程称为马尔可夫过程，而马尔可夫过程平衡的条件是所有状态转移方程被同时满足，我们要求的就是平衡状态下的f(n,m)。
于是我们把上列方程化简，方程两侧同减f(i,j)×p1，然后同除1−p1，就得到新的方程组。根据上列方程我们知道f(1,1)=p4/(1−p1−p2)，因此可以在f(2,x)的方程组中作为常数计算。那么怎么解f(2,x)这样循环引用的方程组呢？其实很简单，我们已经有了f(2,1)关于f(2,2)的表达式，那么我们只要迭代推出f(2,2)关于f(2,1)的表达式，就可以把这两个数解出来。对于f(i,x)的方程组也是同理，因为f(i−1,x)已经求出，可以在这个方程组中作为常数使用，那么我们就可以迭代求出f(i,i)关于f(i,1)的表达式，然后解出f(i,1)（或f(i,i)），再通过DP推出所有的f(i,x)即可。这样一步一步解下去，最后就可以求出f(n,m)的值了。
这样我们就解决了这一题，时间复杂度为O(N2)。需要注意一点，如果不做任何特判，上述计算过程可能会溢出，事实证明，溢出就代表着最终概率无限接近0，而最终概率无限接近0就代表p4<eps，所以把这种情况判掉之后，其余情况都可以利用上述计算过程得出。
以下是本人代码：

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#define eps 1e-10
using namespace std;
int n,m,k;
double p1,p2,p3,p4,f[2010][2010];

int main()
{
    while(scanf("%d%d%d%lf%lf%lf%lf",&n,&m,&k,&p1,&p2,&p3,&p4)!=EOF)
    {
        if (p4<eps) {printf("0.00000
");continue;}
        f[1][1]=p4/(1-p1-p2);
        for(int i=2;i<=n;i++)
        {
            double s1=1,s2=0;
            for(int j=2;j<=i;j++)
            {
                if (j<=k)
                {
                    s1=s1*p2/(1-p1);
                    s2=(s2*p2+f[i-1][j-1]*p3+p4)/(1-p1);
                }
                else
                {
                    s1=s1*p2/(1-p1);
                    s2=(s2*p2+f[i-1][j-1]*p3)/(1-p1);
                }
            }
            s1=1-s1*p2/(1-p1);
            s2=(p4+s2*p2)/(1-p1);
            f[i][1]=s2/s1;
            for(int j=2;j<=i;j++)
            {
                if (j<=k) f[i][j]=(f[i][j-1]*p2+f[i-1][j-1]*p3+p4)/(1-p1);
                else  f[i][j]=(f[i][j-1]*p2+f[i-1][j-1]*p3)/(1-p1);
            }
        }
        printf("%.5lf
",f[n][m]);
    }

    return 0;
}