UVA 10828 Back to Kernighan-Ritchie【期望+高斯消元】

题意：

给出类似程序控制流程图，从每个结点出发到每个后继结点的概率均相等。当执行到没有后继结点的结点时，程序终止。程序总是从1号结点开始，给定若干查询，求每个结点的期望执行次数。

分析：

大白书上的例题。
设$dp[i]$表示期望执行次数，$out[i]$表示该点的出度，那么一个子节点的期望次数就是父节点的期望数次乘上从父节点走到该结点的概率即出度的倒数。
那么我们就可以得到各个相连结点之间的线性关系，不能直接递推计算，我们可以将他们转化成矩阵的形式，然后进行消元即可。
注意：
1. 程序总是从结点$1$开始，即我们有$1$的概率要进入结点$1$，所以假设一个结点$0$，连向结点$1$，出度即为$1$，期望执行次数为$1$。
2. 如何处理无穷大？为了方便处理，避免无穷大对回代过程中的影响，我们采用高斯-约当消元法，直接化成对角矩阵，避免回代。
3. 采用高斯-约当方法还可以避免多余方程。遇到$a[i][i] = 0$的情况，如果$a[i][n] = 0$，即为期望为0的点，如果$a[i][n] != 0$即为期望为无穷的点。最后处理一下所有和无穷大的点有关系的点，也为无穷大， 即在一个循环中。
这两篇介绍高斯消元和高斯约当消元，很清晰
矩阵分析-线性系统-2 高斯消元法、高斯－若尔当消元法
选主元的高斯-约当（Gauss-Jordan）消元法解线性方程组/求逆矩阵

代码：

/*************************************************************************
    > File Name: 10828.cpp
    > Author: jiangyuzhu
    > Mail: 834138558@qq.com 
    > Created Time: 2016/7/27 20:27:00
 ************************************************************************/

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<set>
#include<map>
#include<vector>
#include<algorithm>
using namespace std;
typedef pair<int, int>p;
#define sa(a) scanf("%d", &(a))
const int maxn = 1e2 + 5;
vector<int>v[maxn];
double a[maxn][maxn];
int out[maxn];
bool ans[maxn];
double eps = 1e-8;
void gauss_jordan(int n)
{
    int i, j, k;
    int maxr; // 当前这列绝对值最大的行.
    for(i = 0; i < n; i++){
        maxr = i;
        for(j = i + 1; j < n; j++){
            if(fabs(a[j][i]) > fabs(a[maxr][i])) maxr = j;
        }
        if(fabs(a[maxr][i]) < eps) continue;//放弃这一行，直接处理下一行
        if(maxr != i) for(j = 0; j <= n; j++) swap(a[maxr][j], a[i][j]);
    //与除了第i行外的其他行进行消元
        for(k = 0; k < n; k++) if(k != i) 
            for(j = n; j >= i; j--) a[k][j] -= a[k][i] / a[i][i] * a[i][j];
    }
}
int main (void)
{
    int n;
    int kas = 0;
    while(~scanf("%d", &n) && n){
        memset(out, 0, sizeof(out));
        for(int i = 0; i < n; i++) v[i].clear();
        int x, y;
        while(scanf("%d%d", &x, &y) && (x + y)){
            v[y - 1].push_back(x - 1);
            out[x - 1]++;
        }
        memset(a, 0, sizeof(a));
        a[0][n] = 1.0;
        for(int i = 0; i < n; i++){
            a[i][i] = 1.0;
            for(int j = 0; j < v[i].size(); j++){
                int u = v[i][j];
                a[i][u] -= 1.0 / (double)out[u];
            }
        }
        gauss_jordan(n);
        memset(ans, false, sizeof(ans));
        for(int i = n - 1; i >= 0; i--){
            if(fabs(a[i][i]) < eps && fabs(a[i][n]) > eps) ans[i] = true;
            for(int j = i + 1; j < n; j++){
                if(fabs(a[i][j]) > eps && ans[j]) ans[i] = true;
            }
        }
        int q;scanf("%d", &q);
        printf("Case #%d:\n" ,++kas);
        for(int i = 0; i < q; i++){
            scanf("%d", &x);x--;
            if(ans[x]) puts("infinity");
            else printf("%.3f\n", fabs(a[x][x]) < eps?0.000: a[x][n] / a[x][x]);
        }
    }
    return 0;
}