# museum(Gauss-Jordan/Gauss对比-4倍常数)

【题一】博物馆(museum .cpp/c/pas)

厦门一中  刘定峰

输入文件名：museum.in 输出文件名：museum.out源程序文件名：museum.c/cpp/pas

时间限制：4s 内存限制：128MB

             ［题目描述］

　　有一天Petya和他的朋友Vasya在进行他们众多旅行中的一次旅行，他们决定去参观一座城堡博物馆。这座博物馆有着特别的样式。它包含由m条走廊连接的n间房间，并且满足可以从任何一间房间到任何一间别的房间。
　　两个人在博物馆里逛了一会儿后两人决定分头行动，去看各自感兴趣的艺术品。他们约定在下午六点到一间房间会合。然而他们忘记了一件重要的事：他们并没有选好在哪儿碰面。等时间到六点，他们开始在博物馆里到处乱跑来找到对方（他们没法给对方打电话因为电话漫游费是很贵的）
　　不过，尽管他们到处乱跑，但他们还没有看完足够的艺术品，因此他们每个人采取如下的行动方法：每一分钟做决定往哪里走，有的概率在这分钟内不去其他地方（即呆在房间不动），有的概率他会在相邻的房间中等可能的选择一间并沿着走廊过去。这里的i指的是当期所在房间的序号。在古代建造是一件花费非常大的事，因此每条走廊会连接两个不同的房间，并且任意两个房间至多被一条走廊连接。
　　两个男孩同时行动。由于走廊很暗，两人不可能在走廊碰面，不过他们可以从走廊的两个方向通行。（此外，两个男孩可以同时地穿过同一条走廊却不会相遇）两个男孩按照上述方法行动直到他们碰面为止。更进一步地说，当两个人在某个时刻选择前往同一间房间，那么他们就会在那个房间相遇。
　　两个男孩现在分别处在a，b两个房间，求两人在每间房间相遇的概率。

             ［输入格式］

        第一行包含四个整数，n表示房间的个数;m表示走廊的数目;a,b (1 ≤ a, b ≤ n),表示两个男孩的初始位置。
　　之后m行每行包含两个整数，表示走廊所连接的两个房间。
　　之后n行每行一个至多精确到小数点后四位的实数表示待在每间房间的概率。
　　题目保证每个房间都可以由其他任何房间通过走廊走到。

             ［输出格式］

       输出一行包含n个由空格分隔的数字，第i个数字代表两个人在第i间房间碰面的概率（输出保留6位小数）

             ［样例输入］

2 1 1 2
1 2
0.5
0.5

             ［样例输出］

0.5000000.500000

             ［数据规模］

对于30%的数据有 n <= 5

对于60%的数据有 n <= 15

对于100%的数据有 n <= 20，n-1 <= m <= n(n-1)/2

             ［题目来源］

CodeForces113D

一眼Gauss消元……

对没错这题确实是消元。

经典的路径概率问题是设出每点到终点的概率

但是这题没给终点……所以要设

我们假设终点t

那么Pt,t=1

Pi,i=0 (i≠t)

Pi,j=∑Px,y*G(i,x)*G(j,y)    //G(i,j)表示1步从x->i的概率,G(i,i)=pi……自己算

于是枚举t,得到O(n*(n^2)^3) Gauss_Jordan做法

呵呵了……

那么看看朴素的Gauss

喂……这不科学……

好吧不得不承认Gauss确实比+Jordan优化快(Jordan编程复杂度低……)

下面这个程序用了2种方法，要看的自己换函数.

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<iostream>
#include<cmath>
using namespace std;
#define For(i,n) for(int i=1;i<=n;i++)
#define Fork(i,k,n) for(int i=k;i<=n;i++)
#define Rep(i,n) for(int i=0;i<n;i++)
#define ord(i,j) ord[i][j]
#define MAXN (20+3)
#define eps (1e-10)
enum {MAXM=MAXN*MAXN};
int n,m,outdegree[MAXN]={0};
int ord[MAXN][MAXN];
void calc_ord()
{
  For(i,n) For(j,n) ord[i][j]=(i-1)*n+j;
}
bool map[MAXN][MAXN]={0};
double p[MAXN];
double mat[MAXM][MAXM*2],g[MAXN][MAXN];
double f[MAXM][MAXM*2],tmp[MAXM*2];

void print()
{
  For(i,n*n)
  {
    For(j,n*n+1) printf("%.2lf ",f[i][j]);
    cout<<endl;
  }
  cout<<endl;
}

void gauss(int n)
{
  memcpy(f,mat,sizeof(f));
  For(i,n)
  {
  //  print();
    int p=i;
    
    //Fork(j,i+1,n*n) if (fabs(f[j][i])>fabs(f[p][i])) p=j;
    
    while (fabs(f[p][i])<eps&&p<n) p++; 
    if (fabs(f[p][i])<eps) continue;
    if (p^i)
    {
      copy(f[p]+1,f[p]+1+n+1,tmp+1);
      copy(f[i]+1,f[i]+1+n+1,f[p]+1);
      copy(tmp+1,tmp+1+n+1,f[i]+1);
    }
  //  print(); 
    For(j,n)
    {
      if (i==j) continue;
      double p2=f[j][i]/f[i][i];//cout<<"p="<<p<<endl;
      For(k,n+1) f[j][k]-=f[i][k]*p2;
    }
  }
//  For(i,n*n) f[i][n*n+1]/=f[i][i],f[i][i]=1;
//  print();
}
void gauss_no_yuedan(int n)
{
  memcpy(f,mat,sizeof(f));
  for(int i=1;i<=n;i++)
    {
      int pos;
      for(pos = i;pos < n && abs(f[pos][i]) <= eps;pos++);
      if (pos != i)
        {
    copy(f[i]+1,f[i] + n + 2,tmp+1);
    copy(f[pos],f[pos] + n + 2,f[i]+1);
    copy(tmp+1,tmp + n + 2,f[pos]+1);
        }
      for(int j=i+1;j<=n;j++)
        {
    double mul = (f[j][i] / f[i][i]);
    for(int k=i;k<=n+1;k++)
      f[j][k] -= f[i][k] * mul;
        }
    }
//  print();
  for(int i=n;i;i--)
  {
    Fork(j,i+1,n) f[i][n+1]-=f[i][j]*f[j][n+1],f[i][j]=0;
    if (f[i][i]!=0) f[i][n+1]/=f[i][i],f[i][i]=1;
  }
//   print();
}
int main()
{
  freopen("museum.in","r",stdin);
  freopen("museum.out","w",stdout);
  int x,y;
  cin>>n>>m>>x>>y;calc_ord();
  For(i,m)
  {
    int u,v;
    scanf("%d%d",&u,&v);
    map[u][v]=map[v][u]=1;
    outdegree[v]++;outdegree[u]++;
  }
  For(i,n) scanf("%lf",&p[i]),map[i][i]=1,g[i][i]=p[i];
  For(i,n) 
    For(j,n) 
      if (i^j&&map[i][j]) g[i][j]=(1-p[i])/(double)outdegree[i];
  For(i,n)
    For(j,n)
    if (i^j)
    {
      For(x,n)
        For(y,n)
          if (map[i][x]&&map[j][y]) mat[ord(i,j)][ord(x,y)]=g[i][x]*g[j][y];
          else mat[ord(i,j)][ord(x,y)]=0;
      mat[ord(i,j)][ord(i,j)]-=1;
    }
    else mat[ord(i,i)][ord(i,i)]=1; 
  For(i,n*n) mat[i][n*n+1]=0;
//  For(i,n*n) mat[i][n*n+i]=1; 
  
  
  For(i,n)
  {
    mat[ord(i,i)][n*n+1]=1;
  //  print();
    gauss_no_yuedan(n*n);
    mat[ord(i,i)][n*n+1]=0;
    printf("%.6lf ",f[ord(x,y)][n*n+1]/f[ord(x,y)][ord(x,y)]);
  //  if (i^n) printf(" ");else puts("");
  }  
  puts("");
  return 0;
}

PS:正解O(n^6)，方法是把等式左边不确定的那几个方程改成e1*x1+e2*x2+...+en*xn 

ex表示x的情况……ex是坐边系数矩阵第x项