用高斯消元法(gauss)求解一类期望问题

http://acm.zjut.edu.cn/ShowProblem.aspx?ShowID=1423

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<iostream>
#include<string>
#include<cmath>
#include<algorithm>
using namespace std;
#define MAXN 800

double mat[MAXN][MAXN];
char map[40][40];
int num[40][40];

struct node
{
    int x;
    int y;
};
node Star,end;
int n,m,len;
int di[4][2]={1,0,0,1,-1,0,0,-1};
int cut;

int isok(int x,int y)
{
    return x >= 0 && x < n && y >= 0 && y < m && map[x][y]!='X';
}

node queue[2*MAXN];
int head,tail;

void fillnum(int x,int y) //对地图进行重新标号，标号从0开始
{
    num[x][y]=++cut;
    tail=head=0;
    node e,q;
    e.x=x;e.y=y;
    queue[tail++]=e;
    while(head<tail)
    {
        e=queue[head++];
        for(int i=0;i<4;i++)
        {
            q.x=e.x+di[i][0];
            q.y=e.y+di[i][1];
            if(isok(q.x,q.y) && num[q.x][q.y]==-1)
            {
                num[q.x][q.y]=++cut;
                queue[tail++]=q;
            }
        }
    }
}

bool gauss(int n) 
{ 
    int i, j, row, idx; 
    double buf, maxx; 
    for(row = 0; row < n; row ++) 
    { 
        for(maxx = 0, i = row; i < n; i ++)
        { 
            if(maxx < fabs(mat[i][row])) 
            { 
                maxx = fabs(mat[i][row]); 
                idx = i; 
            } 
        } 
        if(maxx == 0)
            continue;
        if(idx != row) 
        { 
            for(i = row; i <= n; i ++) 
                swap(mat[row][i], mat[idx][i]); 
        } 
        for(i = row + 1; i < n; i ++)
        { 
            if(fabs(mat[i][row])<1e-8)
                continue;
            buf = mat[i][row] / mat[row][row]; 
            for(j = row; j <= n; j ++) 
                mat[i][j] -= buf * mat[row][j]; 
        } 
    } 
    for(i = n - 1; i >= 0; i --) 
    { 
        for(j = i + 1; j < n; j ++) 
            mat[i][n] -= mat[i][j] * mat[j][j]; 
        mat[i][i] = mat[i][n] / mat[i][i]; 
    } 
    return true; 
} 

int main()
{
    int i,j,k,x,y;
    freopen("D:\\in.txt","r",stdin);
    while(scanf("%d%d",&n,&m)!=EOF)
    {
        len=0;
        for(i=0;i<n;i++)
        {
            scanf("%s",map[i]);
            for(j=0;j<m;j++)
            {
                if(map[i][j]=='D')
                {
                    Star.x=i;
                    Star.y=j;
                }
                else if(map[i][j]=='E')
                {
                    end.x=i;
                    end.y=j;
                }
            }
        }
        memset(num,255,sizeof(num));
        cut=-1;
        fillnum(Star.x,Star.y);
        if(num[end.x][end.y]==-1)
        {
            printf("tragedy!\n");
            continue;
        }
        memset(mat,0,sizeof(mat));
        for(i=0;i<n;i++)  //建立N元一次方程组
        {
            for(j=0;j<m;j++)
            {
                if(num[i][j]!=-1)
                {
                    int now=num[i][j];
                    int count=0;
                    for(k=0;k<4;k++)
                    {
                        x=i+di[k][0];y=j+di[k][1];
                        if(isok(x,y))
                        {
                            mat[now][num[x][y]]=-1;
                            count++;
                        }
                        mat[now][now]=count;
                        mat[now][cut+1]=count;
                    }
                }
            }
        }
        x=num[end.x][end.y];
        memset(mat[x],0,sizeof(mat[x]));
        mat[x][x]=1;
        if(gauss(cut+1))
        {
            if(mat[num[Star.x][Star.y]][num[Star.x][Star.y]]<=1000000)
                printf("%0.2lf\n",mat[num[Star.x][Star.y]][num[Star.x][Star.y]]);
            else
                printf("tragedy!\n");
        }
        else
        {
            printf("tragedy!\n");
        }
    }
    return 0;
}

http://acm.zjut.edu.cn/ShowProblem.aspx?ShowID=1317

Description:

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

Input:

每行有两个整数m (2<m<=100)，n (0 < n < m)。

这题我们以两个飞盘的距离为状态进行转移。

那么E[n]=E[n+2]/4+E[n-2]/4+E[n]/2+1，化简成:2E[n]-E[n+2]-E[n-2]=4。

首先对于两个飞盘给定的起始距离n，我们可以先搜索一下可否到达状态0，如果不行，则直接输出INF。在搜索的过程中，顺便把重新标号也进行了。为什么这题也要重新标号呢？

因为该题中，假设m是偶数，那么对于任意的n，n+1和n-1都是不可达的状态。请一定记得，如果方程组中有不可达的点的话，就会导致无解的情况！

接下来对每一个可达的状态建立一个如上的方程，最后用高斯消元法解除E[No[n]]即可！

#include<iostream>
#include<string>
#include<cmath>
#include<stdio.h>
#include<memory.h>
#include<algorithm>
using namespace std;
#define maxn 200

double mat[maxn][maxn];
int num[maxn];
int n,m;
int cut;

int getnum(int x) //获取x的状态，对于一个环来说，x和m-x是一样的，我们取小的
{
    if(x<0)
        x=-x;
    if(x>n/2)
        x=n-x;
    return x;
}

void dfs(int x) //重标号
{
    num[x]=cut++;
    int y=getnum(x+2);
    if(num[y]==-1)
        dfs(y);
    y=getnum(x-2);
    if(num[y]==-1)
        dfs(y);
}

bool gauss(int n) 
{ 
    int i, j, row, idx; 
    double buf, maxx; 
    for(row = 0; row < n; row ++) 
    { 
        for(maxx = 0, i = row; i < n; i ++)
        { 
            if(maxx < fabs(mat[i][row])) 
            { 
                maxx = fabs(mat[i][row]); 
                idx = i; 
            } 
        } 
        if(maxx == 0) return false; 
        if(idx != row) 
        { 
            for(i = row; i <= n; i ++) 
                swap(mat[row][i], mat[idx][i]); 
        } 
        for(i = row + 1; i < n; i ++)
        { 
            buf = mat[i][row] / mat[row][row]; 
            for(j = row; j <= n; j ++) 
                mat[i][j] -= buf * mat[row][j]; 
        } 
    } 
    for(i = n - 1; i >= 0; i --) 
    { 
        for(j = i + 1; j < n; j ++) 
            mat[i][n] -= mat[i][j] * mat[j][j]; 
        mat[i][i] = mat[i][n] / mat[i][i]; 
    } 
    return true; 
} 

int main()
{
    int i,len;
    freopen("D:\\in.txt","r",stdin);
    while(scanf("%d%d",&n,&m)==2)
    {
        memset(num,255,sizeof(num));
        if(m>n/2)
            m=n-m;
        cut=0;
        dfs(m);
//         for(i=0;i<=n/2;i++)
//             cout<<num[i]<<" ";
//         cout<<endl;
        if(num[0]==-1)
        {
            printf("INF\n");
            continue;
        }
        memset(mat,0,sizeof(mat));
        len=n/2;
        int now,j;
        for(i=0;i<=len;i++)
        {
            if(num[i]!=-1)
            {
                now=num[i];
                mat[now][now]=2;
                mat[now][cut]=4;
                j=getnum(i-2);
                mat[now][num[j]]-=1;
                j=getnum(i+2);
                mat[now][num[j]]-=1;
            }
        }
        j=num[0];
        memset(mat[j],0,sizeof(mat[j]));
        mat[j][j]=1;
//         for(i=0;i<cut;i++)
//         {
//             for(j=0;j<=cut;j++)
//                 cout<<mat[i][j]<<" ";
//             cout<<endl;
//         }
        if(gauss(cut))
        {
            printf("%0.2lf\n",mat[num[m]][num[m]]);
        }
        else
        {
            printf("INF\n");
        }
    }
    return 0;
}