【HDU4089】Activation-概率DP好题

测试地址：Activation
题目大意： Tomato要在服务器上激活一个游戏，一开始服务器序列中有$N$个人，他排在第$M$位，每次服务器会对序列中第一位的玩家进行激活，有四种结果：
1.有$p_1$的概率会激活失败，这时候序列的状态是不变的。
2.有$p_2$的概率第一位的玩家会连接错误，这时候序列中第一位的玩家会成为最后一位，其他玩家相对位置不变。
3.有$p_3$的概率第一位的玩家激活成功，这时候第一位的玩家会离开序列。
4.有$p_4$的概率服务器崩溃，发生这件事之后所有玩家都不能激活了。
求Tomato遇到服务器崩溃并且在服务器崩溃时处于前$K$位的概率。
做法：这一道题是概率DP，而且是一道非常好的题目。
题意比较复杂，但是看到概率，我们就应该想到用DP来求概率。设$f(i,j)$为从序列中有$i$人，Tomato在第$j$位的状态开始到达目标状态（遇到服务器崩溃并且在服务器崩溃时处于前$K$位）的概率，有以下状态转移方程（不懂怎么推出来的可以结合题意理解一下）：
当$j=1$，$f(i,j)=f(i,j)\times p_1+f(i,i)\times p_2+p_4$；
当$1<j\le k$，$f(i,j)=f(i,j)\times p_1+f(i,j-1)\times p_2+f(i-1,j-1)\times p_3+p_4$；
当$k<j\le i$，$f(i,j)=f(i,j)\times p_1+f(i,j-1)\times p_2+f(i-1,j-1)\times p_3$；
注意到有些状态之间是会循环引用的，那是不是就不能用DP来解了呢？不是。事实上，在这些状态中进行概率转移的过程称为马尔可夫过程，而马尔可夫过程平衡的条件是所有状态转移方程被同时满足，我们要求的就是平衡状态下的$f(n,m)$。
于是我们把上列方程化简，方程两侧同减$f(i,j)\times p_1$，然后同除$1-p1$，就得到新的方程组。根据上列方程我们知道$f(1,1)=p_4/(1-p_1-p_2)$，因此可以在$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(N^2)$。需要注意一点，如果不做任何特判，上述计算过程可能会溢出，事实证明，溢出就代表着最终概率无限接近$0$，而最终概率无限接近$0$就代表$p_4<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\n");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\n",f[n][m]);
    }

    return 0;
}