在本题中,每一步是独立的,因此即可以看作从$s$移动到$t$的期望步数(对于每一对$s$和$t$都求出答案)
令$f_{i,j}$表示当$s=i$且$t=j$时的答案,则有$f_{i,j}=egin{cases}sum_{(i,k)in n}w_{(i,k)}f_{k,j}+1&(i e j)\0&(i=j)end{cases}$
事实上,对于每一个$j$,其仅有在$i=j$时的方程不同,我们需要利用这个公共部分
更具体的,令$g_{i}=sum_{(i,j)in n}w_{(i,j)}f_{j,i}+1$,即$s=t=i$且强制第一次不能结束的期望步数
将其简单变形,即有$f_{i,i}=sum_{(i,j)in n}w_{(i,j)}f_{j,i}+1-g_{i,i}=0$(即恰好为0)
构造矩阵$M_{i,j}=w_{(i,j)}$,$F_{i,j}=f_{i,j}$,$G_{i,j}=egin{cases}g_{i}&(i=j)\0&(i e j)end{cases}$以及全1矩阵$J_{i,j}=1$,那么根据前面的递推式,就可以得到$MF+J-G=F$
进而,即$(M-I)F=G-J$,问题即变为求$G$以及$frac{G-J}{M-I}$
$G$的计算
再构造一个$1 imes |V|$的矩阵$pi$,满足$sum_{i=1}^{|V|}pi_{1,i}=1$且$pi M=pi$
结论5:$forall 1le ile n,g_{i}pi_{1,i}=1$
考虑条件$MF+J-G=F$,两边同时左乘$pi$即$pi MF+pi J-pi G=pi F$
注意到$pi M=pi$,化简后即有$pi J=pi G$
由此,即$(pi G)_{i}=g_{i}pi_{1,i}$,而$(pi J)_{i}=sum_{i=1}^{|V|}pi_{1,i}=1$,即得证
对于$pi$的计算,相当于解一个$|V|$元的方程,使用高斯消元,其中$pi M=pi$中的$|V|$个方程中是线性相关的(所有方程之和为0),删去其中任意一个即可
$frac{G-J}{M-I}$的计算
这四个矩阵通过前面,都可以在$o(|V|^{3})$内求出,但$M-I$并不一定满秩,因此不能用前面矩阵求逆的手段来计算
结论6(矩阵树定理):对于一张带权有向图$G$,令$M_{i,j}=egin{cases}-e_{(i,j)}&(i e j)\sum_{k=1}^{n}e_{(i,k)}&(i=j)end{cases}$,则矩阵$M$中去掉第$i$行和第$i$列,得到矩阵$M'$的行列式即所有以$i$为根的外向树边权乘积和
结论7:矩阵$M-I$的秩为$|V|-1$
将矩阵$M-I$的每一列求和,其结果都恰好为0,即$M-I$的秩小于$|V|$
显然$I-M$的秩与$M-I$相同,根据边权和为1的性质,其恰为$G$的矩阵$M$
根据矩阵树定理,去掉$M-I$的第$i$行和第$i$列,此时其行列式为以$i$为根的外向树边权乘积和,根据其强连通以及边权$in R^{+}$,不难证明其大于0
因此,去掉第$i$行第$i$列后其满秩,因此$M-I$的秩大于等于$|V|-1$
综上,即证明了$M-I$的秩为$|V|-1$
令$A=G-J$,$B=M-I$,即求$AF=B$
考虑令$A$的第$i$行乘上一个非0数或减去$A$的第$j$行,此时不难发现对于$B$的影响(在$F$不变下)也恰恰是让$B$的第$i$行乘上一个非0数或减去第$j$行
由此进行高斯消元,可以使得$A$中第$i$行仅有第$i$和$|V|$个元素非0(第$|V|$行全为0)
注意到$F_{i,i}=0$,因此有$B_{i,j}=sum_{k=1}^{|V|}A_{i,k}F_{k,j}$,对$i$分类讨论:
1.若$1le i<|V|$,即$B_{i,j}=F_{i,j}+A_{i,n}F_{n,j}$
2.若$i=|V|$,由于$A$的第$|V|$行为0,因此$B_{i,j}=0$,且与$F$的值无关
通过$i=j$时,可以解出$F_{n,i}=frac{B_{i,i}}{A_{i,n}}$,进而即可解出$F_{i,j}=B_{i,j}-A_{i,n}F_{n,j}$
综上,我们就得到了一个时间复杂度为$o(|V|^{3})$的做法
(实现中可能因为一些精度误差而挂掉了QAQ)
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define N 405 4 #define eps 1e-8 5 int n,m,q,x,y,z,deg[N]; 6 double ans,M[N][N],A[N][N],B[N][N],f[N][N]; 7 void guess1(){ 8 for(int i=1;i<=n;i++){ 9 int k=-1; 10 for(int j=i;j<=n;j++) 11 if (fabs(A[j][i])>=eps){ 12 k=j; 13 break; 14 } 15 if (k!=i){ 16 for(int j=i;j<=n+1;j++)swap(A[i][j],A[k][j]); 17 } 18 double x=A[i][i]; 19 for(int j=i;j<=n+1;j++)A[i][j]/=x; 20 for(int j=i+1;j<=n;j++) 21 if (fabs(A[j][i])>=eps){ 22 double x=A[j][i]; 23 for(int k=i;k<=n+1;k++)A[j][k]-=A[i][k]*x; 24 } 25 } 26 for(int i=n;i;i--) 27 for(int j=i+1;j<=n;j++){ 28 A[i][n+1]-=A[i][j]*A[j][n+1]; 29 A[i][j]=0; 30 } 31 } 32 void guess2(){ 33 for(int i=1;i<n;i++){ 34 int k=-1; 35 for(int j=i;j<=n;j++) 36 if (fabs(A[j][i])>=eps){ 37 k=j; 38 break; 39 } 40 if (k!=i){ 41 for(int j=1;j<=n;j++){ 42 swap(A[i][j],A[k][j]); 43 swap(B[i][j],B[k][j]); 44 } 45 } 46 double x=A[i][i]; 47 for(int j=1;j<=n;j++){ 48 A[i][j]/=x; 49 B[i][j]/=x; 50 } 51 for(int j=i+1;j<=n;j++) 52 if (fabs(A[j][i])>=eps){ 53 double x=A[j][i]; 54 for(int k=1;k<=n;k++){ 55 A[j][k]-=A[i][k]*x; 56 B[j][k]-=B[i][k]*x; 57 } 58 } 59 } 60 for(int i=n-1;i;i--) 61 for(int j=1;j<i;j++) 62 if (fabs(A[j][i])>=eps){ 63 double x=A[j][i]; 64 for(int k=1;k<=n;k++){ 65 A[j][k]-=A[i][k]*x; 66 B[j][k]-=B[i][k]*x; 67 } 68 } 69 } 70 int main(){ 71 scanf("%d%d%d",&n,&m,&q); 72 for(int i=1;i<=m;i++){ 73 scanf("%d%d",&x,&y); 74 x++,y++; 75 deg[x]++; 76 M[x][y]++; 77 } 78 for(int i=1;i<=n;i++) 79 for(int j=1;j<=n;j++)M[i][j]=1.0*M[i][j]/deg[i]; 80 for(int i=1;i<n;i++) 81 for(int j=1;j<=n;j++)A[i][j]=M[j][i]-(i==j); 82 for(int i=1;i<=n;i++)A[n][i]=1; 83 A[n][n+1]=1; 84 guess1(); 85 for(int i=1;i<=n;i++)B[i][i]=1.0/A[i][n+1]; 86 for(int i=1;i<=n;i++) 87 for(int j=1;j<=n;j++)B[i][j]--; 88 for(int i=1;i<=n;i++) 89 for(int j=1;j<=n;j++)A[i][j]=M[i][j]-(i==j); 90 guess2(); 91 for(int i=1;i<=n;i++) 92 if (fabs(A[i][n])>=eps)f[n][i]=B[i][i]/A[i][n]; 93 for(int i=1;i<n;i++) 94 for(int j=1;j<=n;j++)f[i][j]=B[i][j]-A[i][n]*f[n][j]; 95 for(int i=1;i<=q;i++){ 96 scanf("%d%d",&x,&y); 97 y++; 98 ans=0; 99 for(int j=1;j<x;j++){ 100 scanf("%d",&z); 101 z++; 102 ans+=f[y][z]; 103 y=z; 104 } 105 printf("%.9f ",ans); 106 } 107 }