数学期望/高斯消元/马尔可夫过程

　　刘汝佳老师白书上的例题- -b

　　本体不满足拓扑关系，但马尔可夫过程是可以高斯消元解的……

　　用「高斯·约当消元」更方便！

 1 //UVA 10828
 2 #include<cmath>
 3 #include<vector>
 4 #include<cstdio>
 5 #include<cstring>
 6 #include<cstdlib>
 7 #include<iostream>
 8 #include<algorithm>
 9 #define rep(i,n) for(int i=0;i<n;++i)
10 #define F(i,j,n) for(int i=j;i<=n;++i)
11 #define D(i,j,n) for(int i=j;i>=n;--i)
12 #define pb push_back
13 using namespace std;
14 int getint(){
15     int v=0,sign=1; char ch=getchar();
16     while(ch<‘0‘||ch>‘9‘){ if (ch==‘-‘) sign=-1; ch=getchar();}
17     while(ch>=‘0‘&&ch<=‘9‘){ v=v*10+ch-‘0‘; ch=getchar();}
18     return v*=sign;
19 }
20 const int N=110;
21 const double eps=1e-8;
22 typedef double Matrix[N][N];
23 /******************tamplate*********************/
24
25 void gauss_jordan(Matrix A,int n){
26     int i,j,k,r;
27     rep(i,n){
28         r=i;
29         for(j=i+1;j<n;++j)
30             if (fabs(A[j][i]) > fabs(A[r][i])) r=j;
31         if (fabs(A[r][i]) < eps)continue;
32         if (r!=i) F(j,0,n) swap(A[r][j],A[i][j]);
33
34         rep(k,n) if (k!=i)
35             D(j,n,i) A[k][j]-=A[k][i]/A[i][i]*A[i][j];
36     }
37 }
38 Matrix A;
39 int n,d[N];
40 vector<int> pre[N];
41 bool inf[N];
42
43 int main(){
44 #ifndef ONLINE_JUDGE
45     freopen("10828.in","r",stdin);
46     freopen("10828.out","w",stdout);
47 #endif
48     int cs=0;
49     while(scanf("%d",&n)==1 && n){
50         memset(d,0,sizeof d);
51         rep(i,n) pre[i].clear();
52         int a,b;
53         while(scanf("%d%d",&a,&b)==2 && a && b){
54             a--; b--;
55             d[a]++;
56             pre[b].pb(a);
57         }
58         memset(A,0,sizeof (A));
59         rep(i,n){
60             A[i][i]=1;
61             rep(j,pre[i].size())
62                 A[i][pre[i][j]]-=1.0/d[pre[i][j]];
63             if (i==0) A[i][n]=1;
64         }
65
66         gauss_jordan(A,n);
67         memset(inf,0,sizeof inf);
68         for(int i=n-1;i>=0;i--){
69             if ( fabs(A[i][i])<eps && fabs(A[i][n])>eps ) inf[i]=1;
70             for(int j=i+1;j<n;j++)
71                 if (fabs(A[i][j])>eps && inf[j]) inf[i]=1;
72         }
73
74         int q,u;
75         scanf("%d",&q);
76         printf("Case #%d:\n",++cs);
77         while(q--){
78             scanf("%d",&u); u--;
79             if (inf[u]) printf("infinity\n");
80             else printf("%.3lf\n",fabs(A[u][u]) < eps ? 0.0 : A[u][n]/A[u][u]);
81         }
82     }
83     return 0;
84 }