BZOJ 3143 HNOI2013 游走高斯消元期望

这道题是我第一次使用高斯消元解决期望类的问题，首发A了，感觉爽爽的....

不过笔者在做完后发现了一些问题，在原文的后面进行了说明。

中文题目，就不翻大意了，直接给原题：

　　一个无向连通图，顶点从1编号到N，边从1编号到M。

　　小Z在该图上进行随机游走，初始时小Z在1号顶点，每一步小Z以相等的概率随机选择当前顶点的某条边，沿着这条边走到下一个顶点，获得等于这条边的编号的分数。当小Z 到达N号顶点时游走结束，总分为所有获得的分数之和。

　　现在，请你对这M条边进行编号，使得小Z获得的总分的期望值最小。

　　输出最小的总分期望值。

Solution:

　　这题贪心很明显，哪条边走过次数的期望最大，它就应该获得最小的编号。

　　所以假设我们已经求出了每条边走过的期望，我们就可以给它们并编上号了。

　　怎么算出每条边走过的期望呢？

　　每条边连接着两个点u,v，很明显的，当我们经过这条边，一定是从两个点中的某一个进入。

　　所以走过边l的期望=走过u点的期望次数*从u点走到l上的概率+走过v点的期望次数*从v点走到l上的概率（其中从i点走到它连接边的概率为1/d[i]，d[i]为i的度数）

　　即：E[l]=e[u]/d[u]+e[v]/d[v]

　　可是我们只知道e[n]=0。但我们还知道这些点之间哪些是连通的，从而可以得出它们之间的关系：

　　我们就可以利用这些点之间的关系建立起方程组，从而使用高斯消元求解。

　　别忘了，点求解完还要带回到每条边上去哦....

　　附Bzoj上的AC代码（codevs上过不了...我也不知道为什么...）

  1 /*
  2   Problem : Bzoj 3143 概率 & 高斯消元 
  3   Author : Robert Yuan
  4   Memory : 15604 kb
  5   Time : 628 MS    
  6   Result : Accept
  7 */
  8 #include <cmath> 
  9 #include <cstdio>
 10 #include <cstring>
 11 #include <cstdlib>
 12 #include <algorithm>
 13 
 14 using namespace std;
 15 
 16 #define maxn 520
 17 
 18 struct Node{
 19     int data,next;
 20 }node[maxn*maxn<<1];
 21 
 22 struct Edge{
 23     int u,v;
 24     double w;
 25 }edge[maxn*maxn<<1];
 26 
 27 #define now node[point].data
 28 #define then node[point].next
 29 
 30 int n,m,cnt;
 31 int head[maxn],deg[maxn];
 32 const double eps=1e-6;
 33 double w[maxn][maxn],rec_x[maxn],ans;
 34 
 35 bool cmp(const Edge A,const Edge B){
 36     return A.w>B.w;
 37 }
 38 
 39 inline int in(){
 40     int x=0;char ch=getchar();
 41     while(ch>'9' || ch<'0') ch=getchar();
 42     while(ch>='0' && ch<='9') x=x*10+ch-'0',ch=getchar();
 43     return x;
 44 }
 45 
 46 void add(int u,int v){
 47     node[++cnt].data=v;node[cnt].next=head[u];deg[u]++;head[u]=cnt;
 48     node[++cnt].data=u;node[cnt].next=head[v];deg[v]++;head[v]=cnt;
 49 }
 50 
 51 void prework(){
 52     n=in();m=in();
 53     int u,v;
 54     for(int i=1;i<=n;i++) head[i]=-1;
 55     for(int i=1;i<=m;i++)
 56         u=in(),v=in(),edge[i].u=u,edge[i].v=v,add(u,v); 
 57     int point;
 58     for(int i=1;i<=n;i++){
 59         w[i][i]=1;
 60         point=head[i];
 61         while(point!=-1){
 62             w[i][now]=-(double)1/deg[now];
 63             point=then;
 64         }
 65     }
 66     w[1][n+1]=1;
 67 }
 68 
 69 void Swap(int i,int j,int x){
 70     double t;
 71     for(int k=x+1;k<=n+1;k++)
 72         t=w[i][k],w[i][k]=w[j][k],w[j][k]=t;
 73 }
 74 
 75 void gauss(){
 76     int i,j;
 77     for(i=1,j=1;i<=n && j<=n;i++,j++){
 78         int max_r=i;
 79         for(int k=i+1;k<=n;k++)
 80             if(fabs(w[max_r][j])+eps<fabs(w[k][j]))
 81                 max_r=k;
 82         if(fabs(w[max_r][j])<eps){i--;continue;}
 83         if(max_r!=i) Swap(i,max_r,j);
 84         for(int k=i+1;k<=n;k++){
 85             double rate=w[k][j]/w[i][j];
 86             w[k][j]=0;
 87             for(int l=j+1;l<=n+1;l++)
 88                 w[k][l]-=w[i][l]*rate;
 89         }
 90     }
 91     
 92     for(int i=n;i>=1;i--)
 93         if(fabs(w[i][i])>eps){
 94             double ans_c=w[i][n+1];
 95             for(int k=i+1;k<=n;k++)
 96                 ans_c-=w[i][k]*rec_x[k];
 97             rec_x[i]=ans_c/w[i][i];
 98         }
 99 }
100 
101 void mainwork(){
102     gauss();
103     for(int i=1;i<=m;i++){
104         edge[i].w=rec_x[edge[i].u]/deg[edge[i].u]+rec_x[edge[i].v]/deg[edge[i].v];
105     }    
106     sort(edge+1,edge+m+1,cmp);
107     for(int i=1;i<=m;i++)
108         ans+=edge[i].w*i;
109     printf("%.3lf",ans);
110 }
111 
112 int main(){
113 #ifndef ONLINE_JUDGE
114     freopen("x.in","r",stdin);
115 #endif
116     prework();
117     mainwork();
118     return 0;
119 }

View Code

[以下是笔者后来发现的问题]

　　首先感谢某位不愿意透露姓名的人堆同学复制了我的代码，然后代入了样例，结果：

　　然后第三行是无解？可是答案却能跑出来...于是我傻了...开始胡乱吹逼[毕竟这么久没打了...]

　　好吧，然后各种逼都被打败了...(●—●)

　　只能认真看看到底出了什么问题，于是发现这个式子有奥妙：

　　左边的e[i]表示走到i点的期望，右边的e[j]表示走出j点的期望。

　　“走到”和“走出”却并不是一样的！

　　我们设e[i]表示走到i点的概率,e'[i]表示走出i点的概率。

　　如果说i!=n那么走到就能走出,e[i]=e'[i]；

　　如果i==n那么就有e'[n]=0,e[n]=1他们俩不同...尽管都已知，而我们列式子的时候，将e'[n]作为未知数带进别的点的式子里，但在n自己的式子中却用的是e[n]，导致两个变量混淆。

　　所以鉴于e'[n]=0，就将建立方程部分修改了一下：

　　于是现在的式子就发生了变化，最后化出来的矩阵也变成了正常的样子：

　　这样解出来的就是e[n]是到达n点的期望=1，当然在给边设定权值的时候，我们用的都是e'[i]，所以我们手动修改一下e'[n]=0就好了...

　　下面是修改过后的代码：

#include <cmath> 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>using namespace std;#define maxn 520struct Node{int data,next;
}node[maxn*maxn<<1];struct Edge{int u,v;double w;
}edge[maxn*maxn<<1];#define now node[point].data
#define then node[point].nextint n,m,cnt;
int head[maxn],deg[maxn];
const double eps=1e-6;
double w[maxn][maxn],rec_x[maxn],ans;bool cmp(const Edge A,const Edge B){return A.w>B.w;
}inline int in(){int x=0;char ch=getchar();while(ch>'9' || ch<'0') ch=getchar();while(ch>='0' && ch<='9') x=x*10+ch-'0',ch=getchar();return x;
}void add(int u,int v){node[++cnt].data=v;node[cnt].next=head[u];deg[u]++;head[u]=cnt;node[++cnt].data=u;node[cnt].next=head[v];deg[v]++;head[v]=cnt;
}void prework(){n=in();m=in();int u,v;for(int i=1;i<=n;i++) head[i]=-1;for(int i=1;i<=m;i++)u=in(),v=in(),edge[i].u=u,edge[i].v=v,add(u,v); int point;for(int i=1;i<=n;i++){w[i][i]=1;point=head[i];while(point!=-1){if(now!=n) w[i][now]=-(double)1/deg[now];else w[i][now]=0;point=then;}}w[1][n+1]=1;
}void Swap(int i,int j,int x){double t;for(int k=x+1;k<=n+1;k++)t=w[i][k],w[i][k]=w[j][k],w[j][k]=t;
}void gauss(){int i,j;for(i=1,j=1;i<=n && j<=n;i++,j++){int max_r=i;for(int k=i+1;k<=n;k++)if(fabs(w[max_r][j])+eps<fabs(w[k][j]))max_r=k;if(fabs(w[max_r][j])<eps){i--;continue;}if(max_r!=i) Swap(i,max_r,j);for(int k=i+1;k<=n;k++){double rate=w[k][j]/w[i][j];w[k][j]=0;for(int l=j+1;l<=n+1;l++)w[k][l]-=w[i][l]*rate;}}for(int i=n;i>=1;i--)if(fabs(w[i][i])>eps){double ans_c=w[i][n+1];for(int k=i+1;k<=n;k++)ans_c-=w[i][k]*rec_x[k];rec_x[i]=ans_c/w[i][i];}rec_x[n]=0;
}void mainwork(){gauss();for(int i=1;i<=m;i++){edge[i].w=rec_x[edge[i].u]/deg[edge[i].u]+rec_x[edge[i].v]/deg[edge[i].v];}    sort(edge+1,edge+m+1,cmp);for(int i=1;i<=m;i++)ans+=edge[i].w*i;printf("%.3lf",ans);
}int main(){
#ifndef ONLINE_JUDGEfreopen("x.in","r",stdin);
#endifprework();mainwork();return 0;
}