【原创】概率DP总结 by kuangbin

概率DP主要用于求解期望、概率等题目。

转移方程有时候比较灵活。

一般求概率是正推，求期望是逆推。通过题目可以体会到这点。

首先先推荐几篇参考的论文：

《信息学竞赛中概率问题求解初探》

《浅析竞赛中一类数学期望问题的解决方法》

《有关概率和期望问题的研究 》

1、POJ 3744 

Scout YYF I

此题是一个用矩阵优化的求概率的题目。

主要思想是分段，根据转移方程用矩阵求解。

题解见 here

/*
POJ 3744

C++  0ms 184K
*/
#include<stdio.h>
#include<string.h>
#include<algorithm>
#include<iostream>
#include<math.h>
using namespace std;

struct Matrix
{
    double mat[2][2];
};
Matrix mul(Matrix a,Matrix b)
{
    Matrix ret;
    for(int i=0;i<2;i++)
      for(int j=0;j<2;j++)
      {
          ret.mat[i][j]=0;
          for(int k=0;k<2;k++)
            ret.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
      }
    return ret;
}
Matrix pow_M(Matrix a,int n)
{
    Matrix ret;
    memset(ret.mat,0,sizeof(ret.mat));
    for(int i=0;i<2;i++)ret.mat[i][i]=1;
    Matrix temp=a;
    while(n)
    {
        if(n&1)ret=mul(ret,temp);
        temp=mul(temp,temp);
        n>>=1;
    }
    return ret;
}

int x[30];
int main()
{
    int n;
    double p;
    while(scanf("%d%lf",&n,&p)!=EOF)//POJ上G++要改为cin输入
    {
        for(int i=0;i<n;i++)
          scanf("%d",&x[i]);
        sort(x,x+n);
        double ans=1;
        Matrix tt;
        tt.mat[0][0]=p;
        tt.mat[0][1]=1-p;
        tt.mat[1][0]=1;
        tt.mat[1][1]=0;
        Matrix temp;

        temp=pow_M(tt,x[0]-1);
        ans*=(1-temp.mat[0][0]);

        for(int i=1;i<n;i++)
        {
            if(x[i]==x[i-1])continue;
            temp=pow_M(tt,x[i]-x[i-1]-1);
            ans*=(1-temp.mat[0][0]);
        }
        printf("%.7lf\n",ans);//POJ上G++要改为%.7f
    }
    return 0;
}

 2、POJ 2096 Collecting Bugs

dp求期望，概率dp入门题。很简单，题解见here

/*
POJ 2096
概率DP
writed by kuangbin

dp求期望
逆着递推求解
题意：（题意看题目确实比较难道，n和s都要找半天才能找到）
   一个软件有s个子系统，会产生n种bug
   某人一天发现一个bug,这个bug属于一个子系统，属于一个分类
   每个bug属于某个子系统的概率是1/s,属于某种分类的概率是1/n
   问发现n种bug,每个子系统都发现bug的天数的期望。

求解：
         dp[i][j]表示已经找到i种bug,j个系统的bug，达到目标状态的天数的期望
         dp[n][s]=0;要求的答案是dp[0][0];
         dp[i][j]可以转化成以下四种状态:
              dp[i][j],发现一个bug属于已经有的i个分类和j个系统。概率为(i/n)*(j/s);
              dp[i][j+1],发现一个bug属于已有的分类，不属于已有的系统.概率为 (i/n)*(1-j/s);
              dp[i+1][j],发现一个bug属于已有的系统，不属于已有的分类,概率为 (1-i/n)*(j/s);
              dp[i+1][j+1],发现一个bug不属于已有的系统，不属于已有的分类,概率为 (1-i/n)*(1-j/s);
        整理便得到转移方程
*/
#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<string.h>
using namespace std;
const int MAXN=1010;
double dp[MAXN][MAXN];

int main()
{
    int n,s;
    while(scanf("%d%d",&n,&s)!=EOF)
    {
        dp[n][s]=0;
        for(int i=n;i>=0;i--)
          for(int j=s;j>=0;j--)
          {
              if(i==n&&j==s)continue;
              dp[i][j]=(i*(s-j)*dp[i][j+1]+(n-i)*j*dp[i+1][j]+(n-i)*(s-j)*dp[i+1][j+1]+n*s)/(n*s-i*j);
          }
        printf("%.4lf\n",dp[0][0]);//POJ上G++要改成%.4f
    }
    return 0;
}

 3、ZOJ 3329  One Person Game

此题的递推方程稍微复杂点，需要转化后求解系数。

题意：有三个骰子，分别有k1,k2,k3个面。
每次掷骰子，如果三个面分别为a,b,c则分数置0，否则加上三个骰子的分数之和。
当分数大于n时结束。求游戏的期望步数。初始分数为0

题解见here

/*
ZOJ 3329
题意：有三个骰子，分别有k1,k2,k3个面。
每次掷骰子，如果三个面分别为a,b,c则分数置0，否则加上三个骰子的分数之和。
当分数大于n时结束。求游戏的期望步数。初始分数为0

设dp[i]表示达到i分时到达目标状态的期望，pk为投掷k分的概率，p0为回到0的概率
则dp[i]=∑(pk*dp[i+k])+dp[0]*p0+1;
都和dp[0]有关系，而且dp[0]就是我们所求，为常数
设dp[i]=A[i]*dp[0]+B[i];
代入上述方程右边得到：
dp[i]=∑(pk*A[i+k]*dp[0]+pk*B[i+k])+dp[0]*p0+1
     =(∑(pk*A[i+k])+p0)dp[0]+∑(pk*B[i+k])+1;
     明显A[i]=(∑(pk*A[i+k])+p0)
     B[i]=∑(pk*B[i+k])+1
     先递推求得A[0]和B[0].
     那么  dp[0]=B[0]/(1-A[0]);
*/
#include<stdio.h>
#include<string.h>
#include<iostream>
#include<algorithm>
using namespace std;

double A[600],B[600];
double p[100];
int main()
{
    int T;
    int k1,k2,k3,a,b,c;
    int n;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%d%d%d%d%d%d",&n,&k1,&k2,&k3,&a,&b,&c);
        double p0=1.0/k1/k2/k3;
        memset(p,0,sizeof(p));
        for(int i=1;i<=k1;i++)
          for(int j=1;j<=k2;j++)
            for(int k=1;k<=k3;k++)
              if(i!=a||j!=b||k!=c)
                p[i+j+k]+=p0;
        memset(A,0,sizeof(A));
        memset(B,0,sizeof(B));
        for(int i=n;i>=0;i--)
        {
            A[i]=p0;B[i]=1;
            for(int j=1;j<=k1+k2+k3;j++)
            {
                A[i]+=A[i+j]*p[j];
                B[i]+=B[i+j]*p[j];
            }
        }
        printf("%.16lf\n",B[0]/(1-A[0]));
    }
    return 0;
}

4、HDU  4405 Aeroplane chess

这题是2012年网络赛的题目。是很简单的概率dp.转移方程很好想到。

求期望。按照公式从后望前递推就可以得到答案了。

解题报告here

/*

概率DP求期望。
形成一个有向无环图。按照公式递推就可以了。
dp[i]表示i点跳到目标状态的期望步数

*/


#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<string.h>
#include<vector>
using namespace std;

const int MAXN=100010;
double dp[MAXN];
vector<int>vec[MAXN];
bool used[MAXN];
int main()
{
    int n,m;
    int u,v;
    while(scanf("%d%d",&n,&m))
    {
        if(n==0&&m==0)break;
        for(int i=0;i<=n;i++)vec[i].clear();
        memset(dp,0,sizeof(dp));
        while(m--)
        {
            scanf("%d%d",&u,&v);
            vec[v].push_back(u);
        }
        memset(used,false,sizeof(used));
        for(int i=0;i<vec[n].size();i++)
        {
            v=vec[n][i];
            dp[v]=0;
            used[v]=true;
        }
        for(int i=n-1;i>=0;i--)
        {
            if(used[i]==false)
            {
                for(int j=i+1;j<=i+6;j++)dp[i]+=dp[j]/6;
                dp[i]+=1;
                used[i]=true;
            }

            for(int j=0;j<vec[i].size();j++)
            {
                v=vec[i][j];
                dp[v]=dp[i];
                used[v]=true;
            }
        }
        printf("%.4lf\n",dp[0]);
    }
    return 0;
}

 5、HDU 4089 Activation

这题是求概率，但是也有种求期望的感觉，都是要列出公式来，化简，递推出答案。

2011年北京现场赛的题目。再比赛时做出来确实不容易，需要对概率DP很熟悉才能做出来。

解题报告见here

/*
HDU 4098
题意：有n个人排队等着在官网上激活游戏。Tomato排在第m个。
对于队列中的第一个人。有一下情况：
1、激活失败，留在队列中等待下一次激活（概率为p1)
2、失去连接，出队列，然后排在队列的最后（概率为p2）
3、激活成功，离开队列（概率为p3）
4、服务器瘫痪，服务器停止激活，所有人都无法激活了。
求服务器瘫痪时Tomato在队列中的位置<=k的概率

解析：
概率DP；
设dp[i][j]表示i个人排队,Tomato排在第j个位置，达到目标状态的概率(j<=i)
dp[n][m]就是所求
j==1:    dp[i][1]=p1*dp[i][1]+p2*dp[i][i]+p4;
2<=j<=k: dp[i][j]=p1*dp[i][j]+p2*dp[i][j-1]+p3*dp[i-1][j-1]+p4;
k<j<=i:  dp[i][j]=p1*dp[i][j]+p2*dp[i][j-1]+p3*dp[i-1][j-1];
化简：
j==1:    dp[i][1]=p*dp[i][i]+p41;
2<=j<=k: dp[i][j]=p*dp[i][j-1]+p31*dp[i-1][j-1]+p41;
k<j<=i:  dp[i][j]=p*dp[i][j-1]+p31*dp[i-1][j-1];

其中:
p=p2/(1-p1);
p31=p3/(1-p1)
p41=p4/(1-p1)

可以循环i=1->n 递推求解dp[i].在求解dp[i]的时候dp[i-1]就相当于常数了。
在求解dp[i][1~i]时等到下列i个方程
j==1:   dp[i][1]=p*dp[i][i]+c[1];
2<=j<=k:dp[i][j]=p*dp[i][j-1]+c[j];
k<j=i:  dp[i][j]=p*dp[i][j]+c[j];
其中c[j]都是常数了。上述方程可以解出dp[i]了。
首先是迭代得到 dp[i][i].然后再代入就可以得到所有的dp[i]了。

注意特判一种情况。就是p4<eps时候，就不会崩溃了，应该直接输出0
*/
#include<stdio.h>
#include<iostream>
#include<math.h>
#include<algorithm>
#include<string.h>
using namespace std;

const int MAXN=2020;
const double eps=1e-5;
double c[MAXN];
double pp[MAXN];
double dp[MAXN][MAXN];
int main()
{
    int n,m,k;
    double p1,p2,p3,p4;
    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;
        }
        double p=p2/(1-p1);
        double p41=p4/(1-p1);
        double p31=p3/(1-p1);
        pp[0]=1.0;//pp[i]=p^1;
        for(int i=1;i<=n;i++) pp[i]=p*pp[i-1];

        dp[1][1]=p41/(1-p);
        c[1]=p41;
        for(int i=2;i<=n;i++)
        {
            for(int j=2;j<=k;j++)c[j]=p31*dp[i-1][j-1]+p41;
            for(int j=k+1;j<=i;j++) c[j]=p31*dp[i-1][j-1];
            double tmp=c[1]*pp[i-1];
            for(int j=2;j<=k;j++)tmp+=c[j]*pp[i-j];
            for(int j=k+1;j<=i;j++)tmp+=c[j]*pp[i-j];
            dp[i][i]=tmp/(1-pp[i]);
            dp[i][1]=p*dp[i][i]+c[1];
            for(int j=2;j<i;j++)dp[i][j]=p*dp[i][j-1]+c[j];
        }
        printf("%.5lf\n",dp[n][m]);
    }
    return 0;
}

6、HDU 4035 Maze

经典的的概率DP的题目。做了可以体会到dp 求期望的一类的方法。

解题报告见here

/*
HDU 4035
kuangbin
http://www.cnblogs.com/kuangbin/


    dp求期望的题。
    题意：
    有n个房间，由n-1条隧道连通起来，实际上就形成了一棵树，
    从结点1出发，开始走，在每个结点i都有3种可能：
        1.被杀死，回到结点1处（概率为ki）
        2.找到出口，走出迷宫 （概率为ei）
        3.和该点相连有m条边，随机走一条
    求：走出迷宫所要走的边数的期望值。

    设 E[i]表示在结点i处，要走出迷宫所要走的边数的期望。E[1]即为所求。

    叶子结点：
    E[i] = ki*E[1] + ei*0 + (1-ki-ei)*(E[father[i]] + 1);
         = ki*E[1] + (1-ki-ei)*E[father[i]] + (1-ki-ei);

    非叶子结点：（m为与结点相连的边数）
    E[i] = ki*E[1] + ei*0 + (1-ki-ei)/m*( E[father[i]]+1 + ∑( E[child[i]]+1 ) );
         = ki*E[1] + (1-ki-ei)/m*E[father[i]] + (1-ki-ei)/m*∑(E[child[i]]) + (1-ki-ei);

    设对每个结点：E[i] = Ai*E[1] + Bi*E[father[i]] + Ci;

    对于非叶子结点i，设j为i的孩子结点，则
    ∑(E[child[i]]) = ∑E[j]
                   = ∑(Aj*E[1] + Bj*E[father[j]] + Cj)
                   = ∑(Aj*E[1] + Bj*E[i] + Cj)
    带入上面的式子得
    (1 - (1-ki-ei)/m*∑Bj)*E[i] = (ki+(1-ki-ei)/m*∑Aj)*E[1] + (1-ki-ei)/m*E[father[i]] + (1-ki-ei) + (1-ki-ei)/m*∑Cj;
    由此可得
    Ai =        (ki+(1-ki-ei)/m*∑Aj)   / (1 - (1-ki-ei)/m*∑Bj);
    Bi =        (1-ki-ei)/m            / (1 - (1-ki-ei)/m*∑Bj);
    Ci = ( (1-ki-ei)+(1-ki-ei)/m*∑Cj ) / (1 - (1-ki-ei)/m*∑Bj);

    对于叶子结点
    Ai = ki;
    Bi = 1 - ki - ei;
    Ci = 1 - ki - ei;

    从叶子结点开始，直到算出 A1,B1,C1;

    E[1] = A1*E[1] + B1*0 + C1;
    所以
    E[1] = C1 / (1 - A1);
    若 A1趋近于1则无解...

*/
#include<stdio.h>
#include<string.h>
#include<algorithm>
#include<iostream>
#include<math.h>
#include<vector>
using namespace std;
const int MAXN=10010;
const double eps=1e-9;//这里1e-8会WA。设为1e-9和1e-10可以
double k[MAXN],e[MAXN];
double A[MAXN],B[MAXN],C[MAXN];

vector<int>vec[MAXN];//存树

bool dfs(int t,int pre)//t的根结点是pre
{
    int m=vec[t].size();//点t的度
    A[t]=k[t];
    B[t]=(1-k[t]-e[t])/m;
    C[t]=1-k[t]-e[t];
    double tmp=0;
    for(int i=0;i<m;i++)
    {
        int v=vec[t][i];
        if(v==pre)continue;
        if(!dfs(v,t))return false;
        A[t]+=(1-k[t]-e[t])/m*A[v];
        C[t]+=(1-k[t]-e[t])/m*C[v];
        tmp+=(1-k[t]-e[t])/m*B[v];
    }
    if(fabs(tmp-1)<eps)return false;
    A[t]/=(1-tmp);
    B[t]/=(1-tmp);
    C[t]/=(1-tmp);
    return true;
}
int main()
{
   // freopen("in.txt","r",stdin);
   // freopen("out.txt","w",stdout);
    int T;
    int n;
    int u,v;
    int iCase=0;
    scanf("%d",&T);
    while(T--)
    {
        iCase++;
        scanf("%d",&n);
        for(int i=1;i<=n;i++)vec[i].clear();
        for(int i=1;i<n;i++)
        {
            scanf("%d%d",&u,&v);
            vec[u].push_back(v);
            vec[v].push_back(u);
        }
        for(int i=1;i<=n;i++)
        {
            scanf("%lf%lf",&k[i],&e[i]);
            k[i]/=100;
            e[i]/=100;
        }
        printf("Case %d: ",iCase);
        if(dfs(1,-1)&&fabs(1-A[1])>eps)
        {
            printf("%.6lf\n",C[1]/(1-A[1]));
        }
        else printf("impossible\n");
    }
}

 7、HDU 3853 LOOPS

 比较简单的概率DP了，入门基础题。

注意一个小陷进。

解题报告here

/*
HDU 3853

解析：
设dp[i][j]表示(i,j)到(R,C)需要消耗的能量
则：
dp[i][j]=p1[i][j]*dp[i][j]+p2[i][j]*dp[i][j+1]+p3[i][j]*dp[i+1][j]+2;
化简得到：
dp[i][j]=p2[i][j]*dp[i][j+1]/(1-p1[i][j])+p3[i][j]*dp[i+1][j]/(1-p1[i][j])+2/(1-p1[i][j]);
注意一种情况就是p1[i][j]==1的情况。
题目只是保证答案小于1000000.但是有的点可能永远都不可能到达的。
所以这样的点出现p1[i][j]是允许的。
否则就会WA了。
*/
#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<string.h>
#include<math.h>
using namespace std;
const int MAXN=1010;
const double eps=1e-5;
double dp[MAXN][MAXN];
double p1[MAXN][MAXN];
double p2[MAXN][MAXN];
double p3[MAXN][MAXN];

int main()
{
    int R,C;
    while(scanf("%d%d",&R,&C)!=EOF)
    {
        for(int i=1;i<=R;i++)
          for(int j=1;j<=C;j++)
            scanf("%lf%lf%lf",&p1[i][j],&p2[i][j],&p3[i][j]);
        dp[R][C]=0;
        for(int i=R;i>=1;i--)
          for(int j=C;j>=1;j--)
          {
              if(i==R&&j==C)continue;
              if(fabs(1-p1[i][j])<eps)continue;
              dp[i][j]=p2[i][j]/(1-p1[i][j])*dp[i][j+1]+p3[i][j]/(1-p1[i][j])*dp[i+1][j]+2/(1-p1[i][j]);
          }
        printf("%.3lf\n",dp[1][1]);
    }
    return 0;
}

8、POJ 2151 Check the difficulty of problems

此题还不算是概率DP的题目。就是DP题，求概率。

想到转移方程就不难了。

题解见here

/*
POJ 2151
题意：
ACM比赛中，共M道题，T个队，pij表示第i队解出第j题的概率
问 每队至少解出一题且冠军队至少解出N道题的概率。

解析：DP
设dp[i][j][k]表示第i个队在前j道题中解出k道的概率
则：
dp[i][j][k]=dp[i][j-1][k-1]*p[j][k]+dp[i][j-1][k]*(1-p[j][k]);
先初始化算出dp[i][0][0]和dp[i][j][0];
设s[i][k]表示第i队做出的题小于等于k的概率
则s[i][k]=dp[i][M][0]+dp[i][M][1]+``````+dp[i][M][k];

则每个队至少做出一道题概率为P1=(1-s[1][0])*(1-s[2][0])*```(1-s[T][0]);
每个队做出的题数都在1~N-1的概率为P2=(s[1][N-1]-s[1][0])*(s[2][N-1]-s[2][0])*```(s[T][N-1]-s[T][0]);

最后的答案就是P1-P2
*/
#include<stdio.h>
#include<string.h>
#include<algorithm>
#include<iostream>
#include<math.h>
using namespace std;

double dp[1010][50][50];
double s[1010][50];
double p[1010][50];
int main()
{
    int M,N,T;
    while(scanf("%d%d%d",&M,&T,&N)!=EOF)
    {
        if(M==0&&T==0&&N==0)break;
        for(int i=1;i<=T;i++)
          for(int j=1;j<=M;j++)
           scanf("%lf",&p[i][j]);
        for(int i=1;i<=T;i++)
        {
            dp[i][0][0]=1;
            for(int j=1;j<=M;j++)dp[i][j][0]=dp[i][j-1][0]*(1-p[i][j]);

            for(int j=1;j<=M;j++)
              for(int k=1;k<=j;k++)
                dp[i][j][k]=dp[i][j-1][k-1]*p[i][j]+dp[i][j-1][k]*(1-p[i][j]);

            s[i][0]=dp[i][M][0];
            for(int k=1;k<=M;k++)s[i][k]=s[i][k-1]+dp[i][M][k];
        }
        double P1=1;
        double P2=1;
        for(int i=1;i<=T;i++)
        {
            P1*=(1-s[i][0]);
            P2*=(s[i][N-1]-s[i][0]);
        }
        printf("%.3lf\n",P1-P2);
    }
    return 0;
}

 9、Codeforces  148D Bag of mice

抓老鼠问题。转移方程要思考下就出来了。

解题报告here

/*
CF 148D
题意：
原来袋子里有w只白鼠和b只黑鼠
龙和王妃轮流从袋子里抓老鼠。谁先抓到白色老师谁就赢。
王妃每次抓一只老鼠，龙每次抓完一只老鼠之后会有一只老鼠跑出来。
每次抓老鼠和跑出来的老鼠都是随机的。
如果两个人都没有抓到白色老鼠则龙赢。王妃先抓。
问王妃赢的概率。

解析：
设dp[i][j]表示现在轮到王妃抓时有i只白鼠，j只黑鼠，王妃赢的概率
明显 dp[0][j]=0,0<=j<=b;因为没有白色老鼠了
      dp[i][0]=1,1<=i<=w;因为都是白色老鼠，抓一次肯定赢了。
      dp[i][j]可以转化成下列四种状态：
      1、王妃抓到一只白鼠，则王妃赢了，概率为i/(i+j);
      2、王妃抓到一只黑鼠，龙抓到一只白色，则王妃输了，概率为j/(i+j)*i/(i+j-1).
      3、王妃抓到一只黑鼠，龙抓到一只黑鼠，跑出来一只黑鼠，则转移到dp[i][j-3]。
      概率为j/(i+j)*(j-1)/(i+j-1)*(j-2)/(i+j-2);
      4、王妃抓到一只黑鼠，龙抓到一只黑鼠，跑出来一只白鼠，则转移到dp[i-1][j-2].
      概率为j/(i+j)*(j-1)/(i+j-1)*i/(i+j-2);

      当然后面两种情况要保证合法，即第三种情况要至少3只黑鼠，第四种情况要至少2只白鼠
*/

#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<string.h>
#include<math.h>
using namespace std;
const int MAXN=1010;
double dp[MAXN][MAXN];
int main()
{
    int w,b;
    while(scanf("%d%d",&w,&b)!=EOF)
    {
        memset(dp,0,sizeof(dp));
        for(int i=0;i<=b;i++)dp[0][i]=0;
        for(int i=1;i<=w;i++)dp[i][0]=1;
        for(int i=1;i<=w;i++)
          for(int j=1;j<=b;j++)
          {
              dp[i][j]+=(double)i/(i+j);//直接抓到白的
              if(j>=3)//抓到黑的，另外一个人也是抓到黑的，跑出一个黑的
                dp[i][j]+=(double)j/(i+j)*((double)(j-1)/(i+j-1))*((double)(j-2)/(i+j-2))*dp[i][j-3];
              if(j>=2) //抓到黑的，另外一个人也是抓到黑的，跑出一个白的
                dp[i][j]+=((double)j/(i+j))*((double)(j-1)/(i+j-1))*((double)i/(i+j-2))*dp[i-1][j-2];
          }
        printf("%.9lf\n",dp[w][b]);
    }
    return 0;
}

10、POJ 3071 Football

足球赛的淘汰赛问题。问最后胜利的概率最大的球队。

简单概率DP

题解报告见here

/*
POJ 3071
题意：2^n个队进行足球赛，每个队打败另外一个队都有一个概率。
问最后胜利的概率最大的是哪只球队。

概率公式，dp算一下就可以了。
*/
#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<string.h>
using namespace std;

double dp[8][200];//dp[i][j]表示在第i场比赛中j胜出的概率
double p[200][200];
int main()
{
    int n;
    while(scanf("%d",&n)!=EOF)
    {
        if(n==-1)break;
        memset(dp,0,sizeof(dp));
        for(int i=0;i<(1<<n);i++)
          for(int j=0;j<(1<<n);j++)
            scanf("%lf",&p[i][j]);
            //cin>>p[i][j];
        for(int i=0;i<(1<<n);i++)dp[0][i]=1;
        for(int i=1;i<=n;i++)//2^n个人要进行n场比赛
        {
            for(int j=0;j<(1<<n);j++)
            {
                int t=j/(1<<(i-1));
                t^=1;
                dp[i][j]=0;
                for(int k=t*(1<<(i-1));k<t*(1<<(i-1))+(1<<(i-1));k++)
                  dp[i][j]+=dp[i-1][j]*dp[i-1][k]*p[j][k];
            }
        }
        int ans;
        double temp=0;
        for(int i=0;i<(1<<n);i++)
        {
            if(dp[n][i]>temp)
            {
                ans=i;
                temp=dp[n][i];
            }
        }
        printf("%d\n",ans+1);
    }
    return 0;
}

 11、SGU  495 Kids and Prizes 

简单的概率DP。　　O（1）,推公式就可以出来。

解题报告here

/*
SGU 495
题意：n个盒子里装有礼物，m个人随机选择礼物，选完之后空盒子放回
问选中的礼物数的期望。

m个人是独立的。
对于每个礼物不被人选中的概率为((n-1)/n)^m
那么不被选中的礼物数的期望就是 n*((n-1)/n)^m
所以答案就是  n-n*((n-1)/n)^m;

*/
#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<math.h>
using namespace std;
int main()
{
    int n,m;
    while(scanf("%d%d",&n,&m)!=EOF)
    {
        double p=(double)(n-1)/n;
        double ans=n-n*pow(p,m);
        printf("%.10lf\n",ans);
    }
    return 0;
}

 12、ZOJ  3380 Patchouli's Spell Cards

用JAVA的大树做的概率DP。有m个位置，每个位置填入1~n中的一个数，求至少有L个数一样的概率。

解题报告见here

/*
 * ZOJ 3380
 * 题目意思：有m个位置，每个位置填入一个数，数的范围是1~n,问至少有L个位置的数一样的概率
 * 输出要是最简分数的形式，所以用大数JAVA
 * 至少有L个位置一样，就是L，L+1，L+2····m个位置一样。
 * 我们从反面来考虑，总数是n^m，我们求没有L个位置一样的数的概率
 * 设 dp[i][j]表示用前i个数，填充j个位置的方案数（要符合没有L个位置是一样的数)
 * dp[i][j]=dp[i-1][j]+Sigm( dp[i-1][j-k]*C[m-(j-k)][k]  ) k<=j&&k<L
 * 其实就是看第i个数，可以不填，填一个位置，两个位置······这样累加过来。
 * 那么最后的答案就是 (n^m-dp[1~n][m])/(n^m)
 */
import java.util.*;
import java.io.*;
import java.math.*;
public class Main
{
    static BigInteger[][] dp=new  BigInteger[110][110];
    static BigInteger[][] C=new BigInteger[110][110];//组合数
    public static void main(String arg[])
    {
        Scanner cin=new Scanner(new BufferedInputStream(System.in));
        for(int i=0;i<105;i++)
        {
            C[i][0]=C[i][i]=BigInteger.ONE;
            for(int j=1;j<i;j++)
                C[i][j]=C[i-1][j-1].add(C[i-1][j]);
        }
        int N,M,L;
        while(cin.hasNext())
        {
            M=cin.nextInt();
            N=cin.nextInt();
            L=cin.nextInt();
            BigInteger tol=BigInteger.valueOf(N).pow(M);
            if(L>M)
            {
                System.out.println("mukyu~");
                continue;
            }
            if(L>M/2)//这个时候可以直接用组合数求出来
            {
                BigInteger ans=BigInteger.ZERO;
                for(int i=L;i<=M;i++)
                    ans=ans.add(C[M][i].multiply(BigInteger.valueOf(N-1).pow(M-i)));
                ans=ans.multiply(BigInteger.valueOf(N));
                BigInteger gcd=ans.gcd(tol);
                System.out.println(ans.divide(gcd)+"/"+tol.divide(gcd));
                continue;
            }
            for(int i=0;i<=N;i++)
                 for(int j=0;j<=M;j++)
                 {
                     dp[i][j]=BigInteger.ZERO;
                 }
            dp[0][0]=BigInteger.ONE;
            for(int i=1;i<=N;i++)
                for(int j=1;j<=M;j++)
                {
                    for(int k=0;k<L&&k<=j;k++)
                        dp[i][j]=dp[i][j].add(dp[i-1][j-k].multiply(C[M-(j-k)][k]));
                }
           BigInteger ans=BigInteger.ZERO;
           for(int i=1;i<=N;i++)
               ans=ans.add(dp[i][M]);    
           ans=tol.subtract(ans);
           BigInteger gcd=ans.gcd(tol);
           System.out.println(ans.divide(gcd)+"/"+tol.divide(gcd));
        }
    }
}

 13、ZOJ 3640  Help Me Escape

比较简单的概率DP，记忆化搜索很好理解，也很容易写。

解题报告here

/*
题意：
一只吸血鬼，有n条路给他走，每次他随机走一条路，
每条路有个限制，如果当时这个吸血鬼的攻击力大于
等于某个值，那么就会花费t天逃出去，否则，花费1天
的时间，并且攻击力增加，问他逃出去的期望

记忆化搜索做
*/

#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<string.h>
#include<math.h>
using namespace std;
const int MAXN=200010;

double dp[MAXN];

int c[110];
int n;

double solve(int f)
{
    if(dp[f]>0)return dp[f];

    dp[f]=0;
    for(int i=0;i<n;i++)
    {
        if(f>c[i])
        {
            double temp=(1.0+sqrt(5))/2*c[i]*c[i];
            int t=(int)temp;
            dp[f]+=(double)t/n;
        }
        else
        {
            dp[f]+=(1+solve(f+c[i]))/n;
        }
    }
    return dp[f];
}
int main()
{
    int f;
    while(scanf("%d%d",&n,&f)!=EOF)
    {
        for(int i=0;i<n;i++)scanf("%d",&c[i]);
        memset(dp,0,sizeof(dp));
        printf("%.3lf\n",solve(f));
    }
    return 0;
}

 14、HDU 4336 Card Collector

求期望，可以状态压缩概率DP求解。也可以用容斥原理直接求。解题报告here

/*
HDU 4336
题意：
有N(1<=N<=20)张卡片，每包中含有这些卡片的概率为p1,p2,````pN.
每包至多一张卡片，可能没有卡片。
求需要买多少包才能拿到所以的N张卡片，求次数的期望。


可以用容斥原理做。也可以状态压缩进行概率DP
期望DP
*/
#include<stdio.h>
#include<algorithm>
#include<iostream>
#include<string.h>
using namespace std;
const int MAXN=22;
double p[MAXN];
double dp[1<<MAXN];
int main()
{
    int n;
    while(scanf("%d",&n)!=EOF)
    {
        double tt=0;
        for(int i=0;i<n;i++)
        {
            scanf("%lf",&p[i]);
            tt+=p[i];
        }
        tt=1-tt;//tt就表示没有卡片的概率了
        dp[(1<<n)-1]=0;
        for(int i=(1<<n)-2;i>=0;i--)
        {
            double x=0,sum=1;
            for(int j=0;j<n;j++)
            {
                if((i&(1<<j)))x+=p[j];
                else sum+=p[j]*dp[i|(1<<j)];
            }
            dp[i]=sum/(1-tt-x);
        }
        printf("%.5lf\n",dp[0]);

    }
    return 0;
}

/*
HDU 4336
容斥原理
位元素枚举
*/
#include<stdio.h>
#include<string.h>
#include<iostream>
#include<algorithm>
using namespace std;

double p[22];
int main()
{
    int n;
    while(scanf("%d",&n)==1)
    {
        for(int i=0;i<n;i++)scanf("%lf",&p[i]);
        double ans=0;
        for(int i=1;i<(1<<n);i++)
        {
            int cnt=0;
            double sum=0;
            for(int j=0;j<n;j++)
              if(i&(1<<j))
              {
                  sum+=p[j];
                  cnt++;
              }
            if(cnt&1)ans+=1.0/sum;
            else ans-=1.0/sum;
        }
        printf("%.5lf\n",ans);
    }
    return 0;
}

下面介绍的三题是用高斯消元法求解的概率DP

15、ZJUT 1423 地下迷宫

 一个N*M的迷宫，除了障碍外等概率走，求起点走到终点步数的期望。先在起点进行bfs，找出所以可以到达的点并编号，然后建立方程组求解。

/*

地下迷宫

Description:
由于山体滑坡，DK被困在了地下蜘蛛王国迷宫。为了抢在DH之前来
到TFT，DK必须尽快走出此迷宫。此迷宫仅有一个出口，而由于大BOSS
的力量减弱影响到了DK，使DK的记忆力严重下降，他甚至无法记得他
上一步做了什么。所以他只能每次等概率随机的选取一个方向走。
当然他不会选取周围有障碍的地方走。如DK周围只有两处空地，则每
个都有1/2的概率。现在要求他平均要走多少步可以走出此迷宫。

Input:
先是一行两个整数N, M(1<=N, M<=10)表示迷宫为N*M大小，
然后是N行，每行M个字符，'.'表示是空地，'E’表示出口，'D’表示DK，'X’表示障碍。
Output:
如果DK无法走出或要超过1000000步才能走出，输出tragedy!，
否则输出一个实数表示平均情况下DK要走几步可以走出迷宫，四舍五入到小数点后两位。
Sample Input:
1 2
ED
3 3
D.X
.X.
X.E
Sample Output:
1.00
tragedy!


首先对地图节点重新标号。假设E[i]表示DK从i点开始走出迷宫的期望值。
那么E[i]=(E[a1]+E[a2]+E[a3]+...+E[an])/n+1，其中a1...an是i的相邻节点。
那么对于每一个DK可达的节点来说，都可以为它建立这样的一个方程。现
在假设DK可达的点有N个，那么我们最终将会得到N元一次方程组。最后
利用高斯消元解出E[No[S]]。其中S是DK的起点，No[S]是重标号后的起点
这里要重点注意的是，我们联立方程的时候，一定要注意DK可达这个条件，
不然就会导致无解的情况。

*/
#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<string.h>
#include<queue>
#include<math.h>
using namespace std;

#define eps 1e-9
const int MAXN=200;
double a[MAXN][MAXN],x[MAXN];//方程的左边的矩阵和等式右边的值，求解之后x存的就是结果
int equ,var;//方程数和未知数个数

int Gauss()
{
    int i,j,k,col,max_r;
    for(k=0,col=0;k<equ&&col<var;k++,col++)
    {
        max_r=k;
        for(i=k+1;i<equ;i++)
          if(fabs(a[i][col])>fabs(a[max_r][col]))
            max_r=i;
        if(fabs(a[max_r][col])<eps)return 0;
        if(k!=max_r)
        {
            for(j=col;j<var;j++)
              swap(a[k][j],a[max_r][j]);
            swap(x[k],x[max_r]);
        }
        x[k]/=a[k][col];
        for(j=col+1;j<var;j++)a[k][j]/=a[k][col];
        a[k][col]=1;
        for(i=0;i<equ;i++)
          if(i!=k)
          {
              x[i]-=x[k]*a[i][k];
              for(j=col+1;j<var;j++)a[i][j]-=a[k][j]*a[i][col];
              a[i][col]=0;
          }
    }
    return 1;
}

char map[20][20];
int num[20][20];
struct Node
{
    int x,y;
};
int sx,sy,ex,ey;
int n,m;
int dir[4][2]={{0,1},{0,-1},{1,0},{-1,0}};
int cnt;

void bfs()
{
    memset(num,-1,sizeof(num));
    cnt=0;
    num[sx][sy]=cnt++;
    queue<Node>que;
    Node temp;
    Node tt;
    temp.x=sx;temp.y=sy;
    que.push(temp);
    while(!que.empty())
    {
        temp=que.front();
        que.pop();
        for(int i=0;i<4;i++)
        {
            tt.x=temp.x+dir[i][0];
            tt.y=temp.y+dir[i][1];
            if(tt.x>=0&&tt.x<n&&tt.y>=0&&tt.y<m&&map[tt.x][tt.y]!='X'&&num[tt.x][tt.y]==-1)
            {
                num[tt.x][tt.y]=cnt++;
                que.push(tt);
            }
        }
    }
}
int main()
{
    //freopen("in.txt","r",stdin);
    //freopen("out.txt","w",stdout);
    while(scanf("%d%d",&n,&m)!=EOF)
    {
        for(int i=0;i<n;i++)
        {
            scanf("%s",&map[i]);
            for(int j=0;j<m;j++)
            {
                if(map[i][j]=='D')
                {
                    sx=i;sy=j;
                }
                if(map[i][j]=='E')
                {
                    ex=i;ey=j;
                }
            }
        }
        bfs();
        if(num[ex][ey]==-1)
        {
            printf("tragedy!\n");
            continue;
        }
        memset(a,0,sizeof(a));
        memset(x,0,sizeof(x));
        equ=var=cnt;
        for(int i=0;i<n;i++)
         for(int j=0;j<m;j++)
           if(num[i][j]!=-1)
           {
               int now=num[i][j];
               if(map[i][j]=='E')
               {
                   memset(a[now],0,sizeof(a[now]));
                   x[now]=0;
                   a[now][now]=1;
                   continue;
               }

               int Count=0;
               for(int k=0;k<4;k++)
               {
                   int tx=i+dir[k][0];
                   int ty=j+dir[k][1];
                   if(tx>=0&&tx<n&&ty>=0&&ty<m&&map[tx][ty]!='X'&&num[tx][ty]!=-1)
                   {
                       a[now][num[tx][ty]]=-1;
                       Count++;
                   }
                   a[now][now]=Count;
                   x[now]=Count;
               }
           }
        if(Gauss())
        {
            if(x[num[sx][sy]]<=1000000)printf("%.2lf\n",x[num[sx][sy]]);
            else printf("tragedy!\n");
        }
        else printf("tragedy!\n");
    }
    return 0;
}

16、ZJUT 1317 掷飞盘

在一个环上抛掷两个飞盘 ，每个飞盘等概率往左和右走，问两个飞盘走到同一个地方所需要步数的期望。

按照他们的距离表示状态进行概率DP。dp[i]=dp[i-2]/4+dp[i+2]/4+dp[i]/2+1.整理下就出来方程。注意是循环的，要进行处理。

/*
ZJUT 1317
掷飞盘

Description:
m个人位于正m边形的顶点上，彼此抛掷飞盘。他们共有两个飞盘，
且开始时这两个飞盘位于相距为n的两个人的手中（相邻两个人
相距为1，依此类推）。在每次抛掷时两个飞盘被同时抛出，飞盘都
以1/2的概率被抛到掷飞盘的人左边相邻的人，1/2的概率被抛到
右边相邻的人。此过程一直进行，直到两个飞盘被掷到同一个人
手中，求此抛掷飞盘的游戏平均情况下(期望)会在抛掷几次后结束。

Input:
每行有两个整数m (2<m<=100)，n (0 < n < m)。
Output:
对每组数据m,n，输出平均所需步数(四舍五入，保留两位小数），
如果有限步内不可能结束就输出INF。

Sample Input:
3 1
4 1
Sample Output:
4.00
INF

*/
#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<string.h>
#include<queue>
#include<math.h>
using namespace std;

#define eps 1e-9
const int MAXN=200;
double a[MAXN][MAXN],x[MAXN];//方程的左边的矩阵和等式右边的值，求解之后x存的就是结果
int equ,var;//方程数和未知数个数

int Gauss()
{
    int i,j,k,col,max_r;
    for(k=0,col=0;k<equ&&col<var;k++,col++)
    {
        max_r=k;
        for(i=k+1;i<equ;i++)
          if(fabs(a[i][col])>fabs(a[max_r][col]))
            max_r=i;
        if(fabs(a[max_r][col])<eps)return 0;
        if(k!=max_r)
        {
            for(j=col;j<var;j++)
              swap(a[k][j],a[max_r][j]);
            swap(x[k],x[max_r]);
        }
        x[k]/=a[k][col];
        for(j=col+1;j<var;j++)a[k][j]/=a[k][col];
        a[k][col]=1;
        for(i=0;i<equ;i++)
          if(i!=k)
          {
              x[i]-=x[k]*a[i][k];
              for(j=col+1;j<var;j++)a[i][j]-=a[k][j]*a[i][col];
              a[i][col]=0;
          }
    }
    return 1;
}

int n,m;
int num[MAXN];
int cnt;
int getnum(int x)
{
    x=(x%n+n)%n;
    if(x>n-x)x=n-x;
    return x;
}
void dfs(int x)
{
    x=getnum(x);
    num[x]=cnt++;
    int y=getnum(x+2);
    if(num[y]==-1)dfs(y);
    y=getnum(x-2);
    if(num[y]==-1)dfs(y);
}
int main()
{
    while(scanf("%d%d",&n,&m)!=EOF)
    {
        memset(num,-1,sizeof(num));
        cnt=0;
        m=getnum(m);
        dfs(m);
        if(num[0]==-1)
        {
            printf("INF\n");
            continue;
        }
        memset(a,0,sizeof(a));
        memset(x,0,sizeof(x));
        for(int i=0;i<n;i++)
          if(num[i]!=-1)
          {
              int now=num[i];
              a[now][now]=2;
              x[now]=4;
              int t=getnum(i-2);
              a[now][num[t]]-=1;//这里一定要注意的，要减1，不能直接赋值为-1，
              t=getnum(i+2);//因为i-2和i+2是可能一样的
              a[now][num[t]]-=1;
          }
        int t=num[0];
        memset(a[t],0,sizeof(a[t]));
        a[t][t]=1;
        x[t]=0;
        equ=var=cnt;//这个不要忘记了，经常忘掉！！！
        if(Gauss())printf("%.2lf\n",x[num[m]]);
        else printf("INF\n");
    }
    return 0;
}

17、HDU 4418 Time travel

在坐标轴上用高斯消元法求解。注意N=1的时候要特判一下。解题报告here

/*
HDU 4118
题目：给出一个数轴，有一个起点和一个终点，某个人可以
走1，2，3……m步，每一种情况有一个概率，初始有一个方向，
走到头则返回，问到达终点的期望步数为多少。

比较明显的高斯求期望问题

Sample Input
2
4 2 0 1 0
50 50
4 1 0 2 1
100


Sample Output
8.14
2.00

*/
#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<string.h>
#include<queue>
#include<math.h>
using namespace std;

#define eps 1e-9
const int MAXN=220;
double a[MAXN][MAXN],x[MAXN];//方程的左边的矩阵和等式右边的值，求解之后x存的就是结果
int equ,var;//方程数和未知数个数

int Gauss()
{
    int i,j,k,col,max_r;
    for(k=0,col=0;k<equ&&col<var;k++,col++)
    {
        max_r=k;
        for(i=k+1;i<equ;i++)
          if(fabs(a[i][col])>fabs(a[max_r][col]))
            max_r=i;
        if(fabs(a[max_r][col])<eps)return 0;
        if(k!=max_r)
        {
            for(j=col;j<var;j++)
              swap(a[k][j],a[max_r][j]);
            swap(x[k],x[max_r]);
        }
        x[k]/=a[k][col];
        for(j=col+1;j<var;j++)a[k][j]/=a[k][col];
        a[k][col]=1;
        for(i=0;i<equ;i++)
          if(i!=k)
          {
              x[i]-=x[k]*a[i][k];
              for(j=col+1;j<var;j++)a[i][j]-=a[k][j]*a[i][col];
              a[i][col]=0;
          }
    }
    return 1;
}

int num[MAXN];
double p[MAXN];
int cnt;
int n,N;//n=2*N-2
int M;
void bfs(int s)
{
    memset(num,-1,sizeof(num));
    queue<int>que;
    cnt=0;
    num[s]=cnt++;
    que.push(s);
    while(!que.empty())
    {
        int t=que.front();
        que.pop();
        for(int i=1;i<=M;i++)
        {
            if(fabs(p[i])<eps)continue;//这点很重要，这个想到不能达到的点
            int temp=(t+i)%n;
            if(num[temp]==-1)
            {
                num[temp]=cnt++;
                que.push(temp);
            }
        }
    }
}
int main()
{
    //freopen("in.txt","r",stdin);
    //freopen("out.txt","w",stdout);
    int s,e;
    int D;
    int T;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%d%d%d%d",&N,&M,&e,&s,&D);
        for(int i=1;i<=M;i++){scanf("%lf",&p[i]);p[i]/=100;}

        if(e==s)//这个特判一定需要，否则可能N==1,会被0除，RE
        {
            printf("0.00\n");
            continue;
        }

        n=2*(N-1);
        if(D==1)s=n-s;
        bfs(s);
        if(num[e]==-1&&num[n-e]==-1)
        {
            printf("Impossible !\n");
            continue;
        }
        equ=var=cnt;
        memset(a,0,sizeof(a));
        memset(x,0,sizeof(x));
        for(int i=0;i<n;i++)
          if(num[i]!=-1)
          {
              if(i==e||i==n-e)
              {
                  a[num[i]][num[i]]=1;
                  x[num[i]]=0;
                  continue;
              }
              a[num[i]][num[i]]=1;
              for(int j=1;j<=M;j++)
              {
                  int t=(i+j)%n;
                  if(num[t]!=-1)
                  {
                      a[num[i]][num[t]]-=p[j];
                      x[num[i]]+=j*p[j];
                  }
              }
          }
        if(Gauss())printf("%.2lf\n",x[num[s]]);
        else printf("Impossible !\n");
    }
    return 0;
}

 今年概率DP就做到这吧！2012-10-6