矩阵树定理

bzoj4031 小z的房间

题目大意：给定n*m的网格，有一些障碍点，求生成树的个数。

思路：矩阵树定理直接求解就可以了。但是因为要取模所以要稍微修改一下高斯消元，行列式交换行或列的时候行列式的值要变号。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define p 1000000000
#define LL long long
#define maxm 105
using namespace std;
char map[maxm][maxm];
int dx[4]={1,0,-1,0},dy[4]={0,1,0,-1},bi[maxm][maxm]={0};
LL ai[maxm][maxm]={0};
char in(){
    char ch;
    while(scanf("%c",&ch)==1){
        if (ch=='.'||ch=='*') return ch;
    }
}
int work(int n){
    int i,j;LL x,y,k,ans=1,f=1,t;
    for (i=1;i<=n;++i)
      for (j=1;j<=n;++j) ai[i][j]=(ai[i][j]+p)%p;
    for (i=1;i<=n;++i){
        for (j=i+1;j<=n;++j){
            x=ai[i][i];y=ai[j][i];
            while(y!=0){
                k=x/y;x%=y;swap(x,y);f=-f;
                for (t=i;t<=n;++t) 
                  ai[i][t]=(ai[i][t]-k*ai[j][t]%p+p)%p;
                for (t=i;t<=n;++t) swap(ai[i][t],ai[j][t]);
            }
        }if (!ai[i][i]) return 0;
        ans=ans*ai[i][i]%p;
    }return (int)(ans*f%p+p)%p;
}
int main(){
    int n,m,i,j,k,x,y,tot=0,u,v;scanf("%d%d",&n,&m);
    for (i=1;i<=n;++i){
        for (j=1;j<=m;++j){map[i][j]=in();if (map[i][j]=='.') bi[i][j]=++tot;}
    }for (i=1;i<=n;++i)
      for (j=1;j<=m;++j){
          if (map[i][j]!='.') continue;
          for (k=0;k<4;++k){
              x=i+dx[k];y=j+dy[k];
              if (x<1||x>n||y<1||y>m||map[x][y]!='.') continue;
              u=bi[i][j];v=bi[x][y];
              ++ai[u][u];--ai[u][v];
          }
      }printf("%d
",work(tot-1));
}

bzoj3534 重建（！！！）

题目大意：给定n个城市和他们之间边通行概率，求通行道路恰是生成树的概率。

思路：设矩阵树的数组是G，G[i][j]=Pe/(1-Pe),G[i][i]=-sigma G[i][j],ci=∏(1-Pe)，所以G的n-1阶主子式*tmp就是答案了。

其实原来矩阵树定理中邻接矩阵的1的意义可以看做是概率，这里的是Pe，但并不能简单的改为Pe，而是Pe/(1-Pe)，因为可以发现这样G的主子式是sigma(T)∏(e∈T)Pe/(1-Pe)，乘上ci之后，不是树边的相当于乘了(1-Pe)，是树边的相当于乘了Pe正好是我们想要（一棵生成树的概率就是树边概率*(1-非树边概率)）的答案。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define N 55
#define eps 1e-9
#define LD double
using namespace std;
LD ai[N][N]={0.},ci;int n;
int cmp(LD x,LD y){
    if (x-y>eps) return 1;
    if (y-x>eps) return -1;
    return 0;}
LD gauss(){
    int i,j,k;LD x;--n;
    for (i=1;i<=n;++i){
        if (cmp(ai[i][i],0.)==0)
            for (j=i+1;j<=n;++j)
                if (cmp(ai[j][i],0.)!=0){
                    for (k=1;k<=n;++k) swap(ai[i][k],ai[j][k]);
                    break;}
        for (j=i+1;j<=n;++j){
            if (cmp(ai[j][i],0.)==0) continue;
            x=ai[j][i]/ai[i][i];
            for (k=i;k<=n;++k) ai[j][k]-=ai[i][k]*x;
        }
    }for (x=1.,i=1;i<=n;++i) x*=ai[i][i];
    return fabs(x);}
int main(){
    int i,j;scanf("%d",&n);
    for (ci=1.,i=1;i<=n;++i)
      for (j=1;j<=n;++j){
          scanf("%lf",&ai[i][j]);
          if (i==j) continue;
          if (ai[i][j]>1.-eps) ai[i][j]-=eps;
          if (i<j) ci*=1.-ai[i][j];
          ai[i][j]/=1.-ai[i][j];
      }for (i=1;i<=n;++i)
        for (j=1;j<=n;++j)
          if (i!=j) ai[i][i]-=ai[i][j];
    printf("%.10f
",gauss()*ci);
}