概率DP

HDU 4089 Activation(推导公式)

Activation

Time Limit: 20000/10000 MS (Java/Others) Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 5214 Accepted Submission(s): 1732

Problem 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

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

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

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

可以循环 $i=1-\&gt;n$i=1−>n递推求解 $dp\left[i\right]$dp[i].在求解 $dp\left[i\right]$dp[i]的时候 $dp[i-1]$dp[i−1]就相当于常数了。
在求解 $dp\left[i\right]\left[1i\right]$dp[i][1 i]时等到下列i个方程 //利用消元法，手动求解方程组
$j=1:dp\left[i\right]\left[1\right]=p\ast dp\left[i\right]\left[i\right]+a\left[1\right];$j=1:dp[i][1]=p∗dp[i][i]+a[1];
$2\&lt;=j\&lt;=k:dp\left[i\right]\left[j\right]=p\ast dp\left[i\right][j-1]+a\left[j\right];$2<=j<=k:dp[i][j]=p∗dp[i][j−1]+a[j];
$k\&lt;j=i:dp\left[i\right]\left[j\right]=p\ast dp\left[i\right][j-1]+a\left[j\right];$k<j=i:dp[i][j]=p∗dp[i][j−1]+a[j];
其中 $a\left[j\right]$a[j]都是常数了。上述方程可以解出 $dp\left[i\right]$dp[i]了。
首先是迭代得到 $dp\left[i\right]\left[i\right]$dp[i][i].然后再代入就可以得到所有的 $dp\left[i\right]$dp[i]了。

注意特判一种情况。就是 $p4\&lt;eps$p4<eps时候，就不会崩溃了，应该直接输出 $0$0 //很神奇

#include "bits/stdc++.h"

using namespace std;

typedef long long ll;

const int maxn = 2000+10;
double dp[maxn][maxn], pp[maxn], a[maxn];
double p1, p2, p3, p4, p21, p31, p41;
int N, M, K;

int main() {
    while(cin>>N>>M>>K>>p1>>p2>>p3>>p4) {
        if(p4<1e-6) {
            printf("0.00000\n");
            continue;
        } //这里还不清楚为什么非要加上去。
        p21=p2/(1-p1), p31=p3/(1-p1), p41=p4/(1-p1); //这几个变量是为了表达方便
        dp[1][1]=p41/(1-p21); //这个初始化为突破口
        a[1]=p41; pp[0]=1; //初始化
        for(int i=1; i<=N; ++i) pp[i]=p21*pp[i-1]; //把p21的N次幂都给预处理好，后面使用比较频繁
        for(int i=2; i<=N; ++i) {
            for(int j=2; j<=min(i,K); ++j) a[j]=p31*dp[i-1][j-1]+p41;
            for(int j=K+1; j<=i; ++j) a[j]=p31*dp[i-1][j-1]; // 这两个循环先求出常数项
            double tmp=0; //方程组消元时的常数项积累
            for(int j=1; j<=i; ++j) tmp+=pp[i-j]*a[j];
            dp[i][i]=tmp/(1-pp[i]); //dp[i][i]的积累
            dp[i][1]=p21*dp[i][i]+p41;
            for(int j=2; j<i; ++j) dp[i][j]=p21*dp[i][j-1]+a[j];
        }
        printf("%.5f\n", dp[N][M]);
    }
}

HDU 2262 Where is the canteen(BFS+高斯消元法)

题意： $n\ast m$n∗m的地图，有一个起点，有多个出口，上下左右走，有的格子不能走，求从起点走到一个出口的期望步数是多少。

思路如下：

地图 $n\ast m$n∗m看成一个个格子，编号 $0,1,2,\mathrm{3...}n\ast m-1，$0,1,2,3...n∗m−1， 那么坐标 $i,j$i,j的格子编号为 $i\ast m+j$i∗m+j
$EK$EK表示编号为K的格子到一个出口的期望步数
那么 $EK=0,K$EK=0,K这里为出口的编号
我们要求的则是 $EK,K$EK,K这里为起点的编号
由一个状态可以走向其他状态，假设当前K可以向三个方向走，
那么 $EK=\left(EK\right(a)+EK(b)+EK(c\left)\right)/3+1$EK=(EK(a)+EK(b)+EK(c))/3+1
一般的 $EK=(Enext1+Enext2+Enext3+......Ecnt)/cnt+1$EK=(Enext1+Enext2+Enext3+......Ecnt)/cnt+1
整理得 $Enext1+Enext2+Enext3+...Ecnt-EK\ast cnt=-cnt$Enext1+Enext2+Enext3+...Ecnt−EK∗cnt=−cnt
对每个 $EK(K=0,1,2,\mathrm{3...}n\ast m-1)，$EK(K=0,1,2,3...n∗m−1)，都建立这样的一个等式，那么可以列出 $n\ast m$n∗m个方程，然后采用高斯消元，求出 $E$E（起点）
其中有 $n\ast m$n∗m个方程， $n\ast m$n∗m个变量
$a\left[i\right]\left[j\right]$a[i][j]代表第 $i$i个式子（从 $0$0开始）,第 $j$j个未知数的系数（从 $0$0开始）,其中 $EK$EK为未知数
$a\left[0\right]\left[0\right]\ast E0+a\left[0\right]\left[1\right]\ast E1+...a\left[0\right][n\ast m-1]\ast E(n\ast m-1)=a\left[0\right][n\ast m]$a[0][0]∗E0+a[0][1]∗E1+...a[0][n∗m−1]∗E(n∗m−1)=a[0][n∗m]
$a\left[1\right]\left[0\right]\ast E0+a\left[1\right]\left[1\right]\ast E1+...=a\left[1\right][n\ast m]$a[1][0]∗E0+a[1][1]∗E1+...=a[1][n∗m]
$...$...
$...$...
$a[n\ast m-1]\left[0\right]\ast E0+a[n\ast m-1]\left[1\right]\ast E1+...=a[n\ast m-1][n\ast m]$a[n∗m−1][0]∗E0+a[n∗m−1][1]∗E1+...=a[n∗m−1][n∗m]
所以关键要求 $a\left[i\right]\left[j\right]$a[i][j]，如果 $a$a数组全部求出来了，然后直接带入高斯消元模板就可以了。
预处理：从每个出口进行 $bfs$bfs，把能到达的位置用 $can\left[i\right]\left[j\right]=1$can[i][j]=1标记。
遍历每个格子，对每个格子建立一个方程，求出每个方程每个未知数的系数
用高斯消元求解。
如果起点可以访问到( $flag\left[i\right]\left[j\right]==1$flag[i][j]==1且高斯消元有解），那么输出唯一解，即 $E（sx\ast m+sy）$E（sx∗m+sy）其中， $sx,sy$sx,sy为起点的坐标
否则输出 $-1.$−1.

Where is the canteen

Time Limit: 10000/5000 MS (Java/Others) Memory Limit: 65536/32768 K (Java/Others)
Total Submission(s): 2501 Accepted Submission(s): 810

Problem Description
After a long drastic struggle with himself, LL decide to go for some snack at last. But when steping out of the dormitory, he found a serious problem : he can’t remember where is the canteen… Even worse is the campus is very dark at night. So, each time he move, he check front, back, left and right to see which of those four adjacent squares are free, and randomly walk to one of the free squares until landing on a canteen.

Input
Each case begin with two integers n and m ( n<=15,m<=15 ), which indicate the size of the campus. Then n line follow, each contain m characters to describe the map. There are 4 different type of area in the map:

‘@’ is the start location. There is exactly one in each case.
‘#’ is an impassible square.
‘$’ is a canteen. There may be more than one in the campus.
‘.’ is a free square.

Output
Output the expected number of moves required to reach a canteen, which accurate to 6 fractional digits. If it is impossible , output -1.

Sample Input
1 2
@$
2 2
@.
.$
1 3
@#$

Sample Output
1.000000
4.000000
-1

#include "bits/stdc++.h"

using namespace std;

typedef long long ll;
typedef pair<int,int> pr;
const double eps = 1e-9;

int dx[]={1,-1,0,0};
int dy[]={0,0,1,-1};
double a[250][250], x[250];
int n, m, sx, sy;
bool can[20][20];
char mp[20][20];
queue<pr> q;

void bfs() {
    while(q.size()) {
        int x=q.front().first;
        int y=q.front().second;
        q.pop();
        for(int i=0; i<4; ++i) {
            int xx=x+dx[i];
            int yy=y+dy[i];
            if(xx>=0&&xx<n&&yy>=0&&yy<m&&mp[xx][yy]!='#'&&can[xx][yy]==0) {
                can[xx][yy]=1;
                q.push(pr(xx,yy));
            }
        }
    }
}

bool Gauss(int equ, int var) {
    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][col];
                for(j=col+1; j<var; ++j) a[i][j]-=a[k][j]*a[i][col];
                a[i][col]=0;
            }
        }
    }
    return 1;
}

int main() {
    ios::sync_with_stdio(false);
    while(cin>>n>>m) {
        memset(a,0,sizeof(a));
        memset(can,0,sizeof(can));
        memset(x,0,sizeof(x));
        while(q.size()) q.pop();
        for(int i=0; i<n; ++i) cin>>mp[i];
        for(int i=0; i<n; ++i) {
            for(int j=0; j<m; ++j) {
                if(mp[i][j]=='@') sx=i, sy=j;
                else if(mp[i][j]=='$') {
                    q.push(pr(i,j));
                    can[i][j]=1;
                }
            }
        }
        bfs();
        for(int i=0; i<n; ++i) {
            for(int j=0; j<m; ++j) {
                if(mp[i][j]=='$'||!can[i][j]) a[i*m+j][i*m+j]=1;
                else if(can[i][j]) {
                    int cnt=0;
                    for(int t=0; t<4; ++t) {
                        int x=i+dx[t];
                        int y=j+dy[t];
                        if(x>=0&&y>=0&&x<n&&y<m&&can[x][y]) {
                            cnt++;
                            a[i*m+j][x*m+y]=1;
                        }
                    }
                    a[i*m+j][i*m+j]=-cnt;
                    x[i*m+j]=-cnt;
                }
            }
        }
        if(can[sx][sy]&&Gauss(n*m,n*m))
            cout<<setiosflags(ios::fixed)<<setprecision(6)<<x[sx*m+sy]<<endl;
        else cout<<-1<<endl;
    }
}```